第27章 FFT的Matlab实现

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

    27.1 FFT函数

    27.2 IFFT函数

    27.3 FFTSHIFT函数

    27.4 总结

27.1   FFT函数
27.1.1  语法

    Y = fft(x)

    Y = fft(X,n)

    Y = fft(X,[],dim)

    Y = fft(X,n,dim)

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

Y = fft(x)

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

    如果输入X是一个矩阵，Y = fft（X）返回该矩阵的每列的傅里叶变换。

    如果输入X是一个多维数组，实现第一个尺寸不为1的维度的FFT变化，注意这里第一个尺寸不为1是指一个矩阵的第一个尺寸不为1的维。

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

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

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

Y = fft(X,n)

    此函数用于返回n点的DFT。fft(n)和fft(X，n)是等同的，其中n是向量X中第一个尺寸不为1的维度。如果X的长度小于n，则X的长度通过填充零达到长度为n。如果X的长度大于n，序列X被截断。当X是一个矩阵，各列的长度都以相同的方式进行调整。

Y = fft(X,[],dim)

Y = fft(X,n,dim)

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

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

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

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

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

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

实际运行代码如下：

1. Fs = 1000;          %采样率
   
2. T = 1/Fs;           %采样时间单位
   
3. L = 1000;           %信号长度
   
4. t = (0:L-1)*T;        %时间序列
   
5. 
   
6. x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);  %原始信号
   
7. y = x + 2*randn(size(t));     %原始信号叠加了噪声后
   
8. plot(Fs*t(1:50),y(1:50));      %绘制波形
   
9. title('原始信号+零均值随机噪声 ');
   
10. xlabel('时间单位:ms');

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

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

1. Fs = 1000;         %采样率
   
2. T = 1/Fs;           %采样时间单位
   
3. L = 1000;          %信号长度
   
4. t = (0:L-1)*T;       %时间序列
   
5. 
   
6. x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);  %原始信号
   
7. y = x + 2*randn(size(t));    %原始信号叠加了噪声后
   
8. 
   
9. NFFT = 2^nextpow2(L);  %求得最接近采样点的2^n，由于上面是1000点，那么最近的就是1024点。
   
10. Y=fft(y,NFFT)/L;          % 进行FFT变换，除以总的采样点数，方便观察实际值。                     
    
11. f = Fs/2*linspace(0,1,NFFT/2+1); %频率轴，这里只显示Fs/2部分，另一半是对称的。
    
12. 
    
13. plot(f,2*abs(Y(1:NFFT/2+1))) %绘制波形
    
14. title('幅频相应');
    
15. xlabel('频率');
    
16. ylabel('幅度');

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

    这里我们以采样率=256Hz采样信号1.5*sin(2*pi*50*t+pi/3)，并求出其幅频和相频响应，Matlab上面运行的代码如下：

1. Fs = 256;              % 采样率
   
2. N  = 256;             % 采样点数
   
3. n  = 0:N-1;            % 采样序列
   
4. t  = 0:1/Fs:1-1/Fs;     % 时间序列
   
5. f = n * Fs / N;          %真实的频率
   
6. 
   
7. x = 1.5*sin(2*pi*50*t+pi/3) ;  %原始信号
   
8. y = fft(x, N);    %对原始信号做FFT变换
   
9. Mag = abs(y);    %求FFT转换结果的模值
   
10. subplot(2,1,1);
    
11. plot(f, Mag);       %绘制幅频相应曲线
    
12. title('幅频相应');
    
13. xlabel('频率/Hz');
    
14. ylabel('幅度');
    
15. 
    
16. subplot(2,1,2);
    
17. plot(f,  angle(y)*180/pi);   %绘制相频响应曲线，注意这将弧度转换成了角度
    
18. title('相频响应');
    
19. xlabel('频率/Hz');
    
20. ylabel('幅度');

运行后求出的幅频相应和相频响应结果如下：

27.2   IFFT函数
27.2.1  语法

y = ifft(X)

y = ifft(X,n)

y = ifft(X,[],dim)

y = ifft(X,n,dim)

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

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

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运算。

27.2.3    IFFT实例

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

1. Fs = 1000;             %采样率
   
2. T = 1/Fs;               % 采样时间
   
3. L = 1024;               % 信号长度
   
4. t = (0:L-1)*T;           % 时间序列
   
5. 
   
6. y = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);  %50Hz正弦波和120Hz的正弦波的叠加
   
7. 
   
8. subplot(2,1,1);
   
9. plot(Fs*t(1:50), y(1:50));   %绘制原始信号
   
10. title('原始信号');
    
11. 
    
12. Y = fft(y, L);               %分别调用正变换和逆变换               
    
13. Z = ifft(Y);
    
14. subplot(2,1,2);
    
15. plot(Fs*t(1:50), Z(1:50));  %绘制逆变换后的波形
    
16. title('FFT和IFFT转换后的信号');

运行后求出的结果如下：

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

27.3    FFTSHIFT函数

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

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

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

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

1. Fs = 256;              % 采样率
   
2. N  = 256;             % 采样点数
   
3. n  = 0:N-1;            % 采样序列
   
4. t  = 0:1/Fs:1-1/Fs;     % 时间序列
   
5. f = (-N/2:N/2-1) * Fs / N;   %真实的频率
   
6. 
   
7. x = 1.5*sin(2*pi*50*t+pi/3) ;  %原始信号
   
8. y = fft(x, N);      %对原始信号做FFT变换
   
9. Mag = abs(y);    %求FFT转换结果的模值
   
10. subplot(2,1,1);
    
11. plot(f, Mag);       %绘制幅频相应曲线
    
12. title('fft幅频相应');
    
13. xlabel('频率/Hz');
    
14. ylabel('幅度');
    
15. 
    
16. z = fftshift(y);    %对FFT转换后的结果做偏移。
    
17. subplot(2,1,2);
    
18. plot(f, z);        %绘制幅频相应曲线
    
19. title('fftshift幅频相应');
    
20. xlabel('频率/Hz');
    
21. ylabel('幅度');

Matlab的运行结果如下：

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

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

相频相应要怎么分析呢

怎么看呀?。。?