【STM32H7的DSP教程】第28章 FFT和IFFT的Matlab实现（幅频响应和相频响应）

完整版教程下载地址：http://www.armbbs.cn/forum.php?mod=viewthread&tid=94547

第28章       FFT和IFFT的Matlab实现（幅频响应和相频响应）

本章主要讲解fft，ifft和fftshift在matlab上的实现。

28.1 初学者重要提示

28.2 Matlab的FFT函数

28.3 Matlab的IFFT函数

28.4 Matlab的FFTSHIFT函数

28.5 总结

28.1   初学者重要提示

1.   求解FFT相频时的修正比较重要，本章做了一个简易修正方法。重点看28.2.5小节。

28.2  Matlab的FFT函数

28.2.1 函数语法

Y = fft(x)

Y = fft(X,n)

Y = fft(X,n,dim)

28.2.2 函数定义

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

 

28.2.3 函数描述

Y = fft(X)

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

•   如果 X 是向量，则 fft(X) 返回该向量的傅里叶变换。
•   如果 X 是矩阵，则 fft(X) 将 X 的各列视为向量，并返回每列的傅里叶变换。
•   如果 X 是一个多维数组，则 fft(X) 将尺寸大小不等于 1 的第一个数组维度的值视为向量，并返回每个向量的傅里叶变换。

注意这里第一个尺寸不为1是指一个矩阵的第一个尺寸不为1的维。

比如一个矩阵是2*1，那么第一个尺寸不为1的维就是行（尺寸为2）。

X是 1*2*3表示第一个尺寸不为1的维就是列（尺寸为2）。

X为维数5*6*2的话，第一个尺寸不为1的维就是行（尺寸为5）。

Y = fft(X, n)

返回 n 点 DFT。如果未指定任何值，则 Y 的大小与 X 相同。

•   如果 X 是向量且 X 的长度小于 n，则为 X 补上尾零以达到长度 n。
•   如果 X 是向量且 X 的长度大于 n，则对 X 进行截断以达到长度 n。
•   如果 X 是矩阵，则每列的处理与在向量情况下相同。
•   如果 X 为多维数组，则大小不等于 1 的第一个数组维度的处理与在向量情况下相同。

 

Y = fft(X, n, dim)

返回沿维度 dim 的傅里叶变换。例如，如果 X 是矩阵，则 fft(X,n,2) 返回每行的 n 点傅里叶变换。

28.2.4 FFT实例一：幅频响应

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

采样率是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');

运行Matlab后，显示波形如下：

 

通过上面的截图，我们是很难发现波形中的频率成分，下面我们通过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，而是比较接近，这个就是我们在上节教程中所示的频谱泄露以及噪声的干扰。

28.2.5 FFT实例二：相频响应（重要）

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

（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的效果：

 

可以看到已经完全没问题了，求出了频率50Hz的余弦初相为60°左右，频率90Hz的余弦初相也是60°。

28.3 Matlab的IFFT函数

28.3.1 函数语法

y = ifft(X)

y = ifft(X,n)

y = ifft(X,[],dim)

y = ifft(X,n,dim)

y = ifft(..., 'symmetric')

y = ifft(..., 'nonsymmetric')

28.3.2 函数描述

y = ifft(X)

此函数用于返回向量X的离散傅立叶变换（DFT）逆变换结果，计算时使用快速傅里叶算法（Fast Fourier transform (FFT)）。

y = ifft(X,n)

此函数用于返回n点的IDFT。

y = ifft(X,[],dim)

y = ifft(X,n,dim)

上面两个函数用于实现指定维度的IFFT运算。

28.3.3 IFFT实例

下面我们对信号：0.7*sin(2*pi*50*t) + sin(2*pi*120*t)求FFT和IFFT，并绘制原始信号和转换后的信号。Matlab上运行的代码如下：

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

y = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);  %50Hz正弦波和120Hz的正弦波的叠加

subplot(2,1,1); 
plot(Fs*t(1:50), y(1:50));   %绘制原始信号
title('原始信号');
 
Y = fft(y, L);               %分别调用正变换和逆变换                 
Z = ifft(Y);
subplot(2,1,2); 
plot(Fs*t(1:50), Z(1:50));   %绘制逆变换后的波形
title('FFT和IFFT转换后的信号');

运行后求出的结果如下：

 

通过上面的运行结果可以看出，转换后的波形与原始的波形基本是一样的。

28.4 Matlab的FFTSHIFT函数

fftshift的作用正是让正半轴部分和负半轴部分的图像分别关于各自的中心对称。因为直接用fft得出的数据与频率不是对应的，fftshift可以纠正过来

以下是Matlab的帮助文件中对fftshift的说明：

Y = fftshift(X) rearranges the outputs of fft, fft2, and fftn by moving the zero-frequency component to the center of the array. It is useful for visualizing a Fourier transform with the zero-frequency component in the middle of the spectrum. For vectors, fftshift(X) swaps the left and right halves of X。

下面我们在Matlab上面实现一个如下的代码来说明fftshift的使用：

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

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

z = fftshift(y);    %对FFT转换后的结果做偏移。
subplot(2,1,2);
plot(f, z);        %绘制幅频相应曲线
title('fftshift幅频相应');
xlabel('频率/Hz');
ylabel('幅度');

Matlab的运行结果如下：

 

通过上面的运行结果我们可以看到，经过fftshift的调节后，正弦波的中心频率正好对应在了相应的50Hz频率点。使用fftshift还有很多其它的好处，有兴趣的可以查找相关的资料进行了解。

28.5 总结

本章节主要讲解了fft，iff和fftshift的基本用法，如果要深入了解，一定要多练习，多查资料和翻阅相关书籍。