FFT频谱分析法
频谱分辨率D
FFT能够实现的频率分辨率是2pi/N
要求2pi/N≤D
误差主要来自于用FFT做频谱分析时,得到的是离散谱,但是信号是连续谱,只有当N较大时,离散谱的包络才能逼近离散谱,因此N要大一些。
为了方便读频率值,最好关于pi归一化,以w/pi作为横坐标
例1
x(n)=R4(n)
选择FFT的变换区间N为8和16进行频谱分析
clc
close 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进行频谱分析
clc
close 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进行频谱分析
clc
close 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进行频谱分析
clc
close 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进行频谱分析
clc
close 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进行谱分析
clc
clear 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; %频率分辨率F
stem(fk,abs(X16k),'.');
box on
title('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; %频率分辨率F
stem(fk,abs(X32k),'.');
box on
title('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; %频率分辨率F
stem(fk,abs(X64k),'.');
box on
title('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,是周期整数倍,
频谱正确