FFT频谱分析

一、实验目的

1.理解FFT算法的编程思想。

2.熟练掌握利用FFT对信号作频谱分析。

包括正确地进行参数选择、作频谱图以及读频谱图。 3.了解FFT的应用。 二、实验环境

安装MATLAB6.5以上版本 三、实验原理

1. 谱分析参数选择

1)设信号x(t)最高频率为fc,对其进行取样得x(n),根据取样定理，取样频率fs必须满足： fs>=2fc 。

2)设谱分辨率为F,则最小记录时间tpmin= 1/F ；取样点数N≥ 2fc/F ;为使用快速傅里叶变换(FFT)进行谱分析，N还须满足：N= 2M 。 2.用FFT计算信号x(n)的频谱。[设x(n)为实信号] 1)对信号x(n)作N点FFT,得频谱X(k)(k=0~N-1) Matlab语句：Y=fft(x,N) 其中：x----x(n) Y----X(k) 2)幅频谱:|X(k)|；

Matlab语句：abs(Y) 3)功率谱:PSD(k)=|X(k)|2/N=X(k)X*(k)/N Matlab语句：PSD=Y.*conj(Y)/N 其中：conj(Y)-- X*(k)[X(k)的共轭] 3、用FFT实现线性卷积运算

用FFT 实现y(n)=x(n)*h(n)的步骤为： 1) 设x(n)及h(n)的长度分别为N1和N2。为使循环卷积等于线性卷积，用补0的方法使x(n),h(n)

长度均为N，则N须满足N≥ N1+N2-1 ；为用FFT计算DFT, 则N还须满足N= 2M 。

2) 用FFT计算X(k),H(k); (N点) 3)Y(k)= X(k).* H(k) 。 4) y(n)= IFFT[Y(k)] 。 四、实验内容

理解DIT-FFT算法程序的原理，见附录一。 1.FFT谱分析

设信号为x(t)=sin(2πf1t)+sin(2πf2t)+随机噪声，f1=50Hz, f2=120Hz,以取样频率fs=1kHz对x(t)进行取样,样本长度tp=0.25s,得x(n),对x(n)作256点FFT,得频谱X(k),画原信号x(n),幅频谱|X(k)|以及功率谱PSD(k),对信号进行谱分析。 %程序 pufenxi.m clear clc

fs=1000;

t=0:1/fs:0.25; N=256; f1=50;f2=120; s=sin(2*pi*f1*t)+sin(2*pi*f2*t);

x=s+randn(1,length(t)); %信号+噪声x(n) Y=fft(x,N); PSD=Y.*conj(Y)/N;

f=fs/N*(0:N/2-1); %为什么画图时，f的最大值是fs/2?此题作为实验的思考题。大家可以将这里的自变量的取值改为（0:N-1），并把Y的自变量也改为Y(1:N)，画图看看结果。

subplot(311);plot(x); subplot(312);plot(f,abs(Y(1:N/2))); subplot(313);plot(f,PSD(1:N/2));

1)画出图形窗口显示的图形，并注名每个图形的含义。 2)回答下列问题：

i)观察幅频谱图，可以发现，信号x(n)含有的两个频率分量分别是 50 Hz和 120 Hz。 ii)在该程序中的“f=fs/N*(0:N/2-1);”下添加“k=0:N/2-1;”，另外增加一幅图“plot(k,abs(Y(1:N/2)));”观察f和k的变量下Y的差别；

重新运行该程序并观察幅频谱图，图中两峰值对应的下标分别是

14 和 30 。它们的含义为 所测频率的k值，即k=N*f/ fs 。 iii)再将该程序中的N改为512, 重新运行该程序并观察幅频谱图，这时图中两峰值对应的下标分别是 25 和 60 。结果是否和上面的相同？ 不同 为什么？ N不同，N=512 。

iV)本例的频谱分辨率F是 3.90625 Hz,改变f2=60Hz,问：在幅频谱中，能否分

辨f1和f2对应的频率分量？ 能 。为什么？ |f2-f1|>=F 。

再改变f2=52Hz,问：在幅频谱中，能否分辨f1和f2对应的频率分量？ 不能 。 为什么？ |f2-f1|再改变f2=600Hz，在幅频谱中，f2对应的频率分量出现在 400 Hz；

问：在fs=1000Hz的情况下，能否正确检测f2对应的频率分量？ 不能 。 为什么？ 因为fs<2max(f1,f2) 产生了混叠现象 。

为了正确检测f2对应的频率分量，则fs至少取多少Hz? 1205 Hz。在该程序中改变

fs，验证你的结论。

2. 用FFT计算X(k)=FFT[R4(n)],分别取N=4,8,16和64, 绘出幅频谱|X(k)|。

clear clc

N=input('N= '); % x=input('x= '); % x(length(x)+1:N)=zeros(1,N-length(x)); %补0 x(1:N) l=log2(N); x1=zeros(1,N);

for j1=1:N %倒序 x1(j1)=x(bin2dec(fliplr(dec2bin(j1-1,l)))+1); end

%%FFT(DIT)%% M=2;

while(M<=N)

W=exp(-2*j*pi/M); %旋转因子W V=1;

for k=0:1:M/2-1 %k为每级蝶形运算旋转因子的个数 for i=0:M:N-1 %i为各群的首序号 p=k+i; q=p+M/2; A=x1(p+1); B=x1(q+1)*V;

x1(p+1)=A+B; %本级蝶形运算，x1最终存放X(k) x1(q+1)=A-B; end

V=V*W; %旋转因子W的变化 end

M=2*M; %第M级 end

%%%%%%

subplot(211);stem(x,'.');grid on; % title('x(n)'); % subplot(212);stem(abs(x1),'.');grid on; % title('|X(k)|'); %

3 FFT实现任意两个序列的快速卷积。 %程序 fftjuanji.m clear clc

x1=input('x1=');x2=input('x2='); % N1=length(x1);N2=length(x2); %序列x1(n),x2(n)的长度 E=ceil(log2(N1+N2-1)); %ceil---向+∞方向取整 N=2^E; % x1=[x1,zeros(1,N-N1)]; % x2=[x2,zeros(1,N-N2)];

X1=fft(x1,N); % X2=fft(x2,N);

Y=X1.*X2; % y=ifft(Y,N) %

结果分析：

1)回到MATLAB窗口，键入： x1=[1 1 1], x2=[1 2],回车。

结果:y= 1 3 3 2 2)问：可用Matlab中的什么函数验算上述卷积结果？ Conv 并编写程序验证直接卷积和快速卷积的结果

4.研究高密度频谱（计算分辨力）与高分辨率频谱（物理分辨力） 实验内容：设有连续信号

xa(t)=cos(2*pi*6.5*1000*t)+ cos(2*pi*7*1000*t)+ cos(2*pi*9*1000*t)以采样频率fs=32kHZ对信号xa(t)采样，分析下列几种情况的幅频特性。

1) 采集数据长度N=16点，做N=16的FFT;采集数据长度N=16点，补零到256点，做256

点的FFT；

2) 采集数据长度N=64点，做N=64的FFT;采集数据长度N=64点，补零到256点，做256

点的FFT；

3) 采集数据长度N=256点，做N=256点的FFT。 %程序Untitled31.m clear clc

N=input('N= '); N1=input('N1= '); n=0:N-1; n1=0:N1-1; t=1/32*n; x=input('x= ');

x1=[x(1:1:N),zeros(1,(256-N))];

subplot(311);stem(n,x,'.');title('x(n)'); grid on; X=fft(x,N); X1=fft(x1,N1);

magx=abs(X(1:1:N/2+1));k=[0:1:N/2];w=2*pi/N*k;

magx1=abs(X1(1:1:N1/2+1));k1=[0:1:N1/2];w1=2*pi/N1*k1;

subplot(312);plot(w/pi,magx);title('FFT N=16');grid on; xlabel('频率（单位：pi)'); ylabel('|X|');

subplot(313);plot(w1/pi,magx1);title('FFT N=256');grid on; xlabel('频率（单位：pi)'); ylabel('|X1|');