数字信号的 FFT 分析

Question 1

已知信号



这里，N=25，Q= 0.9+j0.3。可以推导出 ，



首先根据这个式子计算X(k) 的理论值，然后计算输入序列x(n)的32个值，再利用基2时间抽选的FFT算法，计算x(n)的离散傅里叶变换(DFT) X(k)，与X(k) 的理论值比较（要求计算结果最少6位有效数字）。

代码

format long;                            %精度达到15位
Q=complex(0.9,0.3);                     %计算Q
n=0:24;
xn=Q.^n;                                %生成25点x(n)
WN=exp(-j*2*pi/25);                     %计算WN
Xk=(1-Q.^25)./(1-Q.*WN.^n);             %用理论公式计算X(k)
xn=[xn,zeros(1,7)];                     %补零生成32点x(n)
Xk1=fft(xn,32);                         %用基2-时间抽选FFT计算X(k)
subplot(3,1,1);                         %绘制定义式计算的X(k)
stem(n,abs(Xk));
title('N=25 Use Definition');
subplot(3,1,2);                         %绘制DIT-FFT计算的X(k)
stem(0:31,abs(Xk1));
title('N=32 Use DIT-FFT');
Xk=[Xk,zeros(1,7)];
subplot(3,1,3);                         %绘制差值图
stem(0:31,abs(Xk1-Xk));
title('Difference');

图像

由图像可知，FFT计算得到的X(k)与理论值基本一致，但也存在较大误差。这可能与计算机运算的舍入有关。

Question 2

假设信号x(n)由下述信号组成：

x(n)=0.001cos(0.45nπ)+sin(0.3nπ)−cos(0.302π−π4)

这个信号有两根主谱线 0.3π 和0.302π靠的非常近，而另一根谱线0.45π的幅度很小，请选择合适的长度 N和窗函数，用 DFT 分析其频谱，得到清楚的三根谱线。

分析

由于

2πω1=2π0.45π=409

2πω2=2π0.3π=203

2πω3=2π0.302π=1000151

综上，N取LowestCommonMultiple40,20,1000=1000.

此时三根谱线分别位于k=225,150,151处.

下面只需要选择合适的窗函数显示谱线即可。

代码

%进行1000点FFT运算
N=1000;
n=[0:N-1];
xn=0.001*cos(0.45*pi*n)+ sin(0.3*pi*n)-cos(0.302*pi*n-pi/4);
Xk=fft(xn,N);
%绘制X(k)谱线
subplot(3,1,1);
stem(0:499,abs(Xk(1:500)));    %由于镜像对称，画出一半即可
xlabel(‘k’);
title('X(k)');
%绘制0.3pi和0.302pi处谱线
subplot(3,1,2);
stem(140:160,abs(Xk(141:161)));
xlabel(‘k’);
title('对0.3pi和0.302pi谱线的放大  对应k=150,151');
%绘制0.45pi处谱线
subplot(3,1,3);
stem(220:230,abs(Xk(221:231)));
xlabel(‘k’);
title('对0.45pi谱线的放大  对应k=225');

图像

由图像可知，由于0.45π（k=225）处谱线幅值太小，在该标度下难以体现，而0.3π（k=150）和0.302π（k=151）处谱线相距太近，看起来像一条，故直接描绘X(k)难以分清，选择合适的窗范围，则可分别观察到清晰的谱线。

注

数字信号处理

Matlab实验二