在MATLAB中用FFT作谱分析

一、目的

（1）进一步加深DFT算法原理和基本性质的理解(因为PFT只是DFT的一种快速算法,所以FFT的运算结果必然满足DFT的性质)
（2）熟悉FFT算法原理及子程序的应用。
（3）掌握用FFT对连续信号和时域离散信号进行频谱分析的基本方法。了解可能出现的分析误差和原因,以便在实际中正确应用FFT。

二、原理和内容

        如果用FF对模拟信号进行谱分析,首先要把模拟信号转换成数字信号,转换时要求知道模拟信号的最高截止频率,以便选择满足采样定理的采样频率。一般选择采样频率是模拟信号中最高频率的3~4倍。另外要选择对模拟信号的观测时间,如果采样频率和观测时间确定,则采样点数也确定了。这里观测时间和对模拟信号进行谱分析的分辨率有关,最小的观测时间和分辨率成倒数关系。要求选择的采样点数和观测时间大于它的最小值。用FFT作谱分析时,要求做FFT的点数服从2的整数幂,这一点在上面选择采样点数时可以考虑满足,即使满足不了,可以通过在序列尾部加0完成。
    如果要进行谱分析的模拟信号是周期信号,最好选择观测时间是信号周期的整数倍。如果不知道信号的周期,要尽量选择观测时间长一些,以减少截断效应的影响。用FF对模拟信号作谱分析是一种近似的谱分析。首先一般模拟信号(除周期信号外)的频谱是连续频谱,而用PPT作谱分析得到的是数字谱,因此应该取FT的点数多一些,用它的包络作为模拟信号的近似谱。另外,如果模拟信号不是严格的带限信号,会因为频谱混叠现象引起谱分析的误差,这种情况下可以预先将模拟信号进行预滤,或者尽量将采样频率取高一些。
    一般频率混叠发生在折叠频率附近,分析时要注意因频率混叠引起的误差。最后要注意一般模拟信号是无限长的,分析时要截断,截断的长度和分辨率有关,但也要尽量取长一些,取得太短因截断引起的误差会很大。举一个极端的例子,一个周期性正弦波,如果所取观察时间太短,例如取小于一个周期,它的波形和正弦波相差太大,肯定误差很大,但如果取得长一些,即使不是周期的整倍数,这种截断效应也会小一些。

三、结果与分析

（1）

figure(1) 
x1n=ones(1,4);   
X1k8=fft(x1n,8);      
X1k16=fft(x1n,16);     
X1k32=fft(x1n,32); 
subplot(2,2,1);mstem(X1k8); title('(1a) 8点DFT[x_1(n)]');xlabel('
w/Π');ylabel('幅度');
axis([0,2,0,1.2*max(abs(X1k8))])
subplot(2,2,2);mstem(X1k16); 
title('(1b)16点DFT[x_1(n)]');xlabel(' w/Π');ylabel('幅度');
axis([0,2,0,1.2*max(abs(X1k16))])
subplot(2,2,3);mstem(X1k32);
title('(1c)32点DFT[x_1(n)]');xlabel(' w/Π');ylabel('幅度');
axis([0,2,0,1.2*max(abs(X1k32))])

（2）

figure(2)
M=8;xa=1:(M/2);  xb=(M/2):-1:1;
x2n=[xa,xb]; 
X2k8=fft(x2n,8);       
X2k16=fft(x2n,16); 
subplot(2,2,1);mstem(X2k8);
 title('(2a) 8点DFT[x_2(n)]');xlabel('w/Π');ylabel('幅度');
axis([0,2,0,1.2*max(abs(X2k8))])
subplot(2,2,2);mstem(X2k16); 
title('(2b)16点DFT[x_2(n)]');xlabel('w/Π');ylabel('幅度');
axis([0,2,0,1.2*max(abs(X2k16))])

（3）

figure(3)
M=8;xa=1:(M/2);  xb=(M/2):-1:1;
x3n=[xb,xa]; 
X3k8=fft(x3n,8);       
X3k16=fft(x3n,16); 
subplot(2,2,1);mstem(X3k8);
title('(3a) 8点DFT[x_3(n)]');xlabel('w/Π');ylabel('幅度');
axis([0,2,0,1.2*max(abs(X3k8))])
subplot(2,2,2);mstem(X3k16); title('(3b)16点DFT[x_3(n)]');xlabel('w/Π');ylabel('幅度');
axis([0,2,0,1.2*max(abs(X3k16))])

（4）

figure(4)
N=4;
n=0:N-1;
x4n=cos(pi/4*n);    
 X4k4=fft(x4n,4); 
 X4k8=fft(x4n,8);        
X4k16=fft(x4n,16);    
subplot(2,2,1);mstem(X4k4); 
title('(4a) 4点DFT[x_4(n)]');xlabel('w/Π');ylabel('幅度 ');
axis([0,2,0,1.2*max(abs(X4k4))])
subplot(2,2,2);mstem(X4k8); 
title('(4b)8点DFT[x_4(n)]');xlabel(' w/Π');ylabel('幅度 ');
axis([0,2,0,1.2*max(abs(X4k8))])
subplot(2,2,3);mstem(X4k16);
title('(4c)16点DFT[x_1(n)]');xlabel(' w/Π');ylabel('幅度 ');
axis([0,2,0,1.2*max(abs(X4k16))])

（5）

figure(5)
fs=64;N=16;
n=0:N-1;t=n/fs;
x=cos(8*pi*t)+cos(16*pi*t)+cos(20*pi*t);
y=fft(x,N);
mag=abs(y);
f=n*fs/N;
subplot(3,2,1),plot(f,mag);
xlabel('频率/Hz');
ylabel('振幅');title('N=16');grid on;
subplot(3,2,2),plot(f(1:N/2),mag(1:N/2));
xlabel('频率/Hz');
ylabel('振幅');title('N=16');grid on;
 
fs=64;N=32;
n=0:N-1;t=n/fs;
x=cos(8*pi*t)+cos(16*pi*t)+cos(20*pi*t);
y=fft(x,N);
mag=abs(y);
f=n*fs/N;
subplot(3,2,3),plot(f,mag);
xlabel('频率/Hz');
ylabel('振幅');title('N=32');grid on;
subplot(3,2,4),plot(f(1:N/2),mag(1:N/2));
xlabel('频率/Hz');
ylabel('振幅');title('N=32');grid on;
 

fs=64;N=64;
n=0:N-1;t=n/fs;
x=cos(8*pi*t)+cos(16*pi*t)+cos(20*pi*t);
y=fft(x,N);
mag=abs(y);
f=n*fs/N;
subplot(3,2,5),plot(f,mag);
xlabel(频率/Hz');
ylabel('振幅');title('N=64');grid on;
subplot(3,2,6),plot(f(1:N/2),mag(1:N/2));
xlabel('频率/Hz');


ylabel('振幅');title('N=64');grid on;

做FFT分析时，幅值大小与FFT选择的点数有关，但不影响分析结果。在IFFT时已经做了处理。要得到真实的振幅值的大小，只要将得到的变换后结果乘以2除以N即可。