MATLAB函数处理数据梳理，得到(幅频响应和相频响应），以便于启发STM32实现方法,STM32 FFT测量相位(是重点，在后面)

原创已于 2024-11-11 08:58:22 修改 · 1.5k 阅读

15 ·

本内容遵循CC 4.0 BY-SA版权协议

GEO检测

于 2024-11-11 08:48:10 首次发布

matlab里面的函数定义

Y = fft(x) 和 y = ifft(X)分别用于实现正变换和逆变换，公式描述如下：

函数描述

为了详细描述快速傅里叶变换（FFT）及其逆变换（IFFT），包括函数参数和使用举例，我们首先需要引入傅里叶逆变换的函数 ifft。然后，我们结合原始描述对 fft 和 ifft 进行详细说明。

快速傅里叶变换 (FFT)

函数： `Y = fft(X)`

使用快速傅里叶变换 (FFT) 算法计算 X 的离散傅里叶变换 (DFT)。

参数：
- X：输入数组，可以是向量、矩阵或多维数组。
- 向量：返回该向量的傅里叶变换。
- 矩阵：将 X 的各列视为向量，并返回每列的傅里叶变换。
- 多维数组：将尺寸大小不等于 1 的第一个数组维度的值视为向量，并返回每个向量的傅里叶变换。
举例：
- 若 X 是 2x1 矩阵，则对两行分别进行傅里叶变换。
- 若 X 是 1x2x3 数组，则对尺寸为 2 的维度（即列）进行傅里叶变换。
- 若 X 是 5x6x2 数组，则对尺寸为 5 的维度（即行）进行傅里叶变换。

函数： `Y = fft(X, n)`

返回 n 点 DFT。

参数：
- n：指定 DFT 的点数。
- 若 X 是向量且长度小于 n，则补上尾零。
- 若 X 是向量且长度大于 n，则进行截断。
- 矩阵和多维数组的处理与向量类似，按列或指定维度处理。

函数： `Y = fft(X, n, dim)`

返回沿维度 dim 的傅里叶变换。

参数：
- dim：指定进行傅里叶变换的维度。
- 例如，若 X 是矩阵，fft(X, n, 2) 返回每行的 n 点傅里叶变换。

快速傅里叶逆变换 (IFFT)

函数： `X_reconstructed = ifft(Y)`

使用快速傅里叶逆变换 (IFFT) 算法计算 Y 的逆离散傅里叶变换 (IDFT)，以重建原始信号 X。

参数：
- Y：经过 FFT 变换后的频域信号。
- 返回的是与 Y 大小相同的时域信号。
举例：
- 若 Y 是由 fft(X) 得到的，则 ifft(Y) 将返回原始的 X（在数值精度允许范围内）。

函数： `X_reconstructed = ifft(Y, n)`

返回 n 点 IDFT。

参数：
- n：指定 IDFT 的点数。
- 若 Y 是向量且长度小于 n，则补上尾零。
- 若 Y 是向量且长度大于 n，则进行截断。
- 矩阵和多维数组的处理与向量类似，按列或指定维度处理。

函数： X_reconstructed = ifft(Y, n, dim)

返回沿维度 dim 的逆傅里叶变换。

参数：
- dim：指定进行逆傅里叶变换的维度。

使用举例

matlab代码

% 示例 1：向量 FFT 和 IFFT
X = [1, 2, 3, 4];
Y = fft(X);  % 计算 X 的 FFT
X_reconstructed = ifft(Y);  % 计算 Y 的 IFFT，重建 X

% 示例 2：矩阵 FFT 和指定维度的 IFFT
X_matrix = [1, 2; 3, 4];
Y_matrix = fft(X_matrix, [], 2);  % 对矩阵的每列进行 4 点 FFT
X_reconstructed_matrix = ifft(Y_matrix, [], 2);  % 对每列进行 IFFT

% 示例 3：多维数组 FFT 和 IFFT
X_nd = rand(2, 3, 4);  % 创建一个 2x3x4 的随机数组
Y_nd = fft(X_nd, [], 1);  % 对第一个尺寸不为 1 的维度（尺寸为 2）进行 FFT
X_reconstructed_nd = ifft(Y_nd, [], 1);  % 对该维度进行 IFFT

通过这些详细的描述和示例，可以更清晰地理解 FFT 和 IFFT 的函数参数及使用方法。参考：

【STM32F429的DSP教程】第28章 FFT和IFFT的Matlab实现（幅频响应和相频响应）-腾讯云开发者社区-腾讯云 (tencent.com)

幅频响应

傅里叶变换的一个常见用途就是查找埋藏在噪声信号中的实际信号的频率成分。下面我们考虑一个这样的例子：

采样率是1000Hz ，信号由如下三个波形组成。

（1）50Hz的正弦波、振幅0,7。

（2）70Hz正弦波、振幅1。

（3）均值为0的随机噪声。

实际运行代码如下：

Fs = 1000;             %采样率
T = 1/Fs;              %采样时间单位
L = 1000;              %信号长度
t = (0:L-1)*T;         %时间序列

x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);  %原始信号
y = x + 2*randn(size(t));                  %原始信号叠加了噪声后
plot(Fs*t(1:50),y(1:50));                  %绘制波形
title('原始信号+零均值随机噪声 ');
xlabel('时间单位:ms');

通过上面的截图，我们是很难发现波形中的频率成分，下面我们通过FFT变换，从频域观察就很方便了，Matlab运行代码如下：

Fs = 1000;           %采样率
T = 1/Fs;            %采样时间单位
L = 1000;            %信号长度
t = (0:L-1)*T;       %时间序列

x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);  %原始信号
y = x + 2*randn(size(t));                  %原始信号叠加了噪声后

NFFT = 2^nextpow2(L);            %求得最接近采样点的2^n，由于上面是1000点，那么最近的就是1024点。
Y=fft(y,NFFT)/L;                 % 进行FFT变换，除以总的采样点数，方便观察实际值。                      
f = Fs/2*linspace(0,1,NFFT/2+1); %频率轴，这里只显示Fs/2部分，另一半是对称的。
 
plot(f,2*abs(Y(1:NFFT/2+1)))     %绘制波形
title('幅频相应');
xlabel('频率');
ylabel('幅度');

从上面的幅频相应，我们可以看出，两个正弦波的频谱并不是准确的0.7和1，而是比较接近，这个就是我们在上节教程中所示的频谱泄露以及噪声的干扰。

相频响应（重要）

这里我们以采样两个余弦波组成的信号为例进行说明，并求出其幅频和相频响应。

（1）50Hz的余弦波，初始相位60°，振幅1.5。

（2）90Hz的余弦波、初始相位60°，振幅1。

（3）采样率256Hz，采集256个点。

Matlab上面运行的代码如下：

Fs = 256;              % 采样率
N  = 256;             % 采样点数
n  = 0:N-1;            % 采样序列
t  = 0:1/Fs:1-1/Fs;     % 时间序列
f = n * Fs / N;          %真实的频率

x = 1.5*cos(2*pi*50*t+pi/3)+ cos(2*pi*90*t+pi/3) ;  %原始信号 
y = fft(x, N);    %对原始信号做FFT变换
Mag = abs(y);    %求FFT转换结果的模值
subplot(2,1,1);
plot(f, Mag);       %绘制幅频相应曲线
title('幅频相应');
xlabel('频率/Hz');
ylabel('幅度');

subplot(2,1,2);
plot(f,  angle(y)*180/pi);   %绘制相频响应曲线，注意这将弧度转换成了角度
title('相频响应方式一');
xlabel('频率/Hz');
ylabel('相角');

求出的幅频响应没问题，而相频响应杂乱无章，造成这个问题的根本原因很多频段的幅值非常小，他们的相角可以不显示出来，这样就可以方便的查看相频响应了。基于此，修正后的代码如下：

Fs = 256;              % 采样率
N  = 256;             % 采样点数
n  = 0:N-1;            % 采样序列
t  = 0:1/Fs:1-1/Fs;     % 时间序列
f = n * Fs / N;          %真实的频率

x = 1.5*cos(2*pi*50*t+pi/3)+ cos(2*pi*90*t+pi/3) ;  %原始信号 
y = fft(x, N);    %对原始信号做FFT变换
Mag = abs(y);    %求FFT转换结果的模值
subplot(3,1,1);
plot(f, Mag);       %绘制幅频相应曲线
title('幅频相应');
xlabel('频率/Hz');
ylabel('幅度');

subplot(3,1,2);
plot(f,  angle(y)*180/pi);   %绘制相频响应曲线，注意这将弧度转换成了角度
title('相频响应方式一');
xlabel('频率/Hz');
ylabel('相角');

subplot(3,1,3);
plot(f,  angle(y)*180/pi.*(Mag>=100)); %绘制相频响应曲线，注意这将弧度转换成了角度
title('相频响应方式二');
xlabel('频率/Hz');
ylabel('相角');

上面代码中Mag>=100是关键，仅展示FFT后幅值大于100的相角。下面再来看Matlab的效果：