FFT频谱分析
- 格式：doc
- 大小：9.52 MB
- 文档页数：9
实验二 FFT频谱分析
一、实验目的
1.理解FFT算法的编程思想。
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 。
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),对信号进行谱分析。
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是 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实现任意两个序列的快速卷积。
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。
%程序
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;
