FFT频谱分析法
频谱分辨率D
FFT能够实现的频率分辨率是2pi/N
要求2pi/N≤D
误差主要来自于用FFT做频谱分析时,得到的是离散谱,但是信号是连续谱,只有当N较大时,离散谱的包络才能逼近离散谱,因此N要大一些。
为了方便读频率值,最好关于pi归一化,以w/pi作为横坐标
例1
x(n)=R4(n)
选择FFT的变换区间N为8和16进行频谱分析
clcclose all;clear all;xn=[ones(1,4)];X8k=fft(xn,8);n=0:7;wk=2*n/8;subplot(1,2,1);stem(wk,abs(X8k),'.','r');title('8 point DFT[x(n)]');xlabel('\omega/\pi');ylabel('amplitude');axis([0,2,0,1.2*max(abs(X8k))]);n=0:15;wk=2*n/16;X16k=fft(xn,16);subplot(1,2,2);stem(wk,abs(X16k),'.','r');title('16 point DFT[x(n)]');xlabel('\omega/\pi');ylabel('amplitude');axis([0,2,0,1.2*max(abs(X16k))]);
x(n)的频谱函数8点和16点采样(8点DFT和16点DFT)
例2
选择FFT的变换区间N为8和16进行频谱分析
clcclose all;clear all;xa=1:4;xb=4:-1:1;xn=[xa,xb];X8k=fft(xn,8);n=0:7;wk=2*n/8;subplot(1,2,1);stem(wk,abs(X8k),'.','r');title('8 point DFT[x(n)]');xlabel('\omega/\pi');ylabel('amplitude');axis([0,2,0,1.2*max(abs(X8k))]);n=0:15;wk=2*n/16;X16k=fft(xn,16);subplot(1,2,2);stem(wk,abs(X16k),'.','r');title('16 point DFT[x(n)]');xlabel('\omega/\pi');ylabel('amplitude');axis([0,2,0,1.2*max(abs(X16k))]);
例3
选择FFT的变换区间N为8和16进行频谱分析
clcclose all;clear all;xa=4:-1:1;xb=1:4;xn=[xa,xb];X8k=fft(xn,8);n=0:7;wk=2*n/8;subplot(1,2,1);stem(wk,abs(X8k),'.','r');title('8 point DFT[x(n)]');xlabel('\omega/\pi');ylabel('amplitude');axis([0,2,0,1.2*max(abs(X8k))]);n=0:15;wk=2*n/16;X16k=fft(xn,16);subplot(1,2,2);stem(wk,abs(X16k),'.','r');title('16 point DFT[x(n)]');xlabel('\omega/\pi');ylabel('amplitude');axis([0,2,0,1.2*max(abs(X16k))]);
x3(n)和x2(n)的8点DFT的模相等,所以图示相同;但是16点不满足循环移位关系,模不同
例4
x(n)=cos(pi*n/4)
选择FFT的变换区间N为8和16进行频谱分析
clcclose all;clear all;n=0:7;xn=cos(pi*n/4);X8k=fft(xn,8);wk=2*n/8;subplot(1,2,1);stem(wk,abs(X8k),'.','r');title('8 point DFT[x(n)]');xlabel('\omega/\pi');ylabel('amplitude');axis([0,2,0,1.2*max(abs(X8k))]);n=0:15;xn=cos(pi*n/4);X16k=fft(xn,16);wk=2*n/16;subplot(1,2,2);stem(wk,abs(X16k),'.','r');title('16 point DFT[x(n)]');xlabel('\omega/\pi');ylabel('amplitude');axis([0,2,0,1.2*max(abs(X16k))]);
N=8和N=16均是周期的整数倍,只在±0.25pi有1根单一谱线
例5
xn=cos(pin/4)+cos(pin/8)
选择FFT的变换区间N为8和16进行频谱分析
clcclose all;clear all;n=0:7;xn=cos(pi*n/4)+cos(pi*n/8);X8k=fft(xn,8);wk=2*n/8;subplot(1,2,1);stem(wk,abs(X8k),'.','r');title('8 point DFT[x(n)]');xlabel('\omega/\pi');ylabel('amplitude');axis([0,2,0,1.2*max(abs(X8k))]);n=0:15;xn=cos(pi*n/4)+cos(pi*n/8);X16k=fft(xn,16);wk=2*n/16;subplot(1,2,2);stem(wk,abs(X16k),'.','r');title('16 point DFT[x(n)]');xlabel('\omega/\pi');ylabel('amplitude');axis([0,2,0,1.2*max(abs(X16k))]);
周期为16,N=8不是周期整数倍,所以频谱不正确;N=16是周期,得到正确频谱,仅在±0.25pi和±0.125pi处有两根谱线
例6
x(t)=cos(8pit)+cos(16pit)+cos(20pit)
采样频率Fs=64Hz,变换区间N=16 32 64进行谱分析
clcclear all;close all;Fs=64;T=1/Fs;subplot(3,1,1)N=16;n=0:N-1;xnT=cos(8*pi*n*T)+cos(16*pi*n*T)+cos(20*pi*n*T);X16k=fft(xnT);X16k=fftshift(X16k);Tp=N*T;F=1/Tp;k=-N/2:N/2-1;fk=k*F; %频率分辨率Fstem(fk,abs(X16k),'.');box ontitle('16 point |DFT[x(nT)]|');xlabel('f(Hz)');ylabel('amplitude');axis([-N*F/2-1,N*F/2-1,0,1.2*max(abs(X16k))]);subplot(3,1,2)N=32;n=0:N-1;xnT=cos(8*pi*n*T)+cos(16*pi*n*T)+cos(20*pi*n*T);X32k=fft(xnT);X32k=fftshift(X32k);Tp=N*T;F=1/Tp;k=-N/2:N/2-1;fk=k*F; %频率分辨率Fstem(fk,abs(X32k),'.');box ontitle('32 point |DFT[x(nT)]|');xlabel('f(Hz)');ylabel('amplitude');axis([-N*F/2-1,N*F/2-1,0,1.2*max(abs(X32k))]);subplot(3,1,3)N=64;n=0:N-1;xnT=cos(8*pi*n*T)+cos(16*pi*n*T)+cos(20*pi*n*T);X64k=fft(xnT);X64k=fftshift(X64k);Tp=N*T;F=1/Tp;k=-N/2:N/2-1;fk=k*F; %频率分辨率Fstem(fk,abs(X64k),'.');box ontitle('64 point |DFT[x(nT)]|');xlabel('f(Hz)');ylabel('amplitude');axis([-N*F/2-1,N*F/2-1,0,1.2*max(abs(X64k))]);
x(t)=cos(8pit)+cos(16pit)+cos(20pit)有3个成分f1=4hz,f2=8hz,f3=10hz,所以周期是0.5s
1.采样频率Fs=64hz,变换区间N=16时,观察时间Tp=16T=0.25s,不是周期整数倍
2.N=32,64时观察时间Tp=0.5s和1s,是周期整数倍,
频谱正确