数字信号处理（一）利用FFT对信号进行频谱分析

数字信号处理（一）利用FFT对信号进行频谱分析

1．实验目的

(1) 进一步加深DFT算法原理和基本性质的理解(因为FFT只是DFT的一种快速算法，所以FFT的运算结果必然满足DFT的基本性质)。
(2) 熟悉FFT算法原理和FFT程序的应用。
(3) 学习利用FFT对离散时间信号进行频谱分析的方法，了解可能出现的误差及其原因，以便在实际中正确应用FFT。

2．实验原理

本实验的原理DFT算法及其相关的基本性质，同时学习DFT的快速算法FFT。对于在数字信号处理之中占有重要地位的有限长序列来说，可以利用傅里叶变换和Z变换来进行处理。除此之外，针对有限长这个特点，可以导出一种更有利的工具—离散傅里叶变换
离散傅里叶变换从理论上解决了傅里叶变换应用于实际的可能性，但若直接按DFT公式来计算，运算量实在太大（与N的平方成比例）。快速傅里叶变换是离散傅里叶变换的快速算法，它大大减少了离散傅里叶变换的运算量，一般可以缩短一二个数量级，且N越大改善效果越明显，从而真正的让DFT的运算在实际工作中真正广泛应用。
对于所给出的有限长序列，逐个进行相应点数的离散傅立叶变换，可以得到该序列的频谱图，并进行相应的理论分析。

3．实验内容

（1） 分别做8点DFT和16点DFT，结果如下

图一 x1的8点DFT结果

如图为X1n的8点DFT频谱分析，图中从上到下分别显示的是x1n的离散时域信号，8点DFT结果以及8点DFT的包络，补了4个0之后做的8点DFT。

图二 x1的16点DFT结果
如图所示是对X1n的16点DFT，由于X1n只有4个数据点，相当于在离散时域信号后面补了12个0之后做的16点DFT。
通过实验所得图像分析可得， 为矩形窗函数，其离散傅立叶变换为sa函数的包络，且图像关于中心对称分布。

（2） 分别做8点和16点DFT，结果如下
图三 x2的8点DFT


图四 x2的16点DFT
通过实验所得图像分析可得， 为三角波，其离散傅立叶变换为sa函数平方的包络，且图像关于中心对称分布。

（3） 的8点DFT和16点DFT对于所给函数进行离散傅立叶变化，画出信号的时域和频域波形如

图五 x3的8点DFT

图六 x3的16点DFT

（4） 的16点DFT和32点DFT

图七 x4的16点DFT

图八 x4的32点DFT
通过实验所得图像分析可得， 为余弦波，其离散傅立叶变换为频率为关于 对称的两个冲激。
（5） 的16和32点DFT，结果如下

图九 x5的16点 32点DFT

（6）
处理方式一：采样率fs=64Hz, 分别做N=16, 32, 64点DFT ，结果如下：

观察我们的信号，不难看出应该是有三个频率分量，4Hz,8Hz,10Hz
由于我们选择相同的采样率，都取128个数据点，因此三种情况下我们抽样所得到的离散时域信号是相同的（反映在图中是左侧三个图像一致的），
A. 当做16点DFT时，由于选择的采样率是64Hz,则频域频谱间隔为4Hz,可以观察到图中二图的第二，三根谱线，对应的就是4Hz和8Hz,而10Hz的频率并不在我们的频率刻度线上，因此实际上是观测不到的，但我们确实看到第4，5，6根谱线确实不为0，是因为他们反映了图形的包络，其实就是上课的时候讲到的假包络
B. 当我们做32点DFT时，由于选择的采样率是64Hz,则对应的频谱间隔为2Hz,这样我们可以在DFT结果之中真实准确的看到在前半个周期内第2，4，5根谱线是有值的，后半个周期是对称的；
C. 当我们做64点DFT时，由于选择的采样率是64Hz,则对应的频谱间隔为1Hz,这样我们可以在DFT结果之中真实准确的看到在前半个周期内第4，8，10根谱线是有值的，后半个周期是对称的；

（7） 令x(n)=x4(n)+jx5(n)，
对x(n)计算16点、32点和 64点 (时域补零)的快速傅里叶变换，即X(k)=FFT［x(n)］

通过实验所得图像分析可得， 为离散复信号，用欧拉公式可化简为 。而可知复信号只有正频率，没有负频率，因此对其进行离散傅立叶变化后得到中心为 的sa函数冲激。且对函数的补零越多，其sa函数包络越明显。

4．思考题

问题一： 在N=8时， x2(n)和x3(n)的幅频特性会相同吗? 为什么? N=16时呢?

当N=8时， 的 幅频特性不完全相同，它们的幅度特性相同，但是相位相差 相角。因为 进行8点周期延拓后的周期函数，平移4点后截短就是 。而时域的平移对应频域的相位线性变化，所以它们会有 的相位差。
当N=16是， 的 进行了时域补零，所对应的周期函数没有什么关系。所以它们的幅频特性不同。

问题二：(2) 对一个周期序列分别取一个周期或多个周期所作的频谱分析结果是否有区别?
有区别，因为对于一个周期和多个周期的频谱分析，其取样时矩形函数的宽度不同，则频域所对应的包络的宽度也不同。当所取的周期越长，其取样时矩形函数的宽度越宽，频域卷积的sa函数的宽度也越窄，越接近理想频谱。

5．个人感想与总结：

1. 通过“利用FFT对信号进行频谱分析”的实验，我掌握了利用Matlab对信号做FFT进行频域分析的方法。通过实验中，不同的函数离散傅立叶变换的结果，并对结果进行理论分析，我对于离散傅立叶变换的性质有了更加深入的理解。对于离散傅立叶变化与连续时间傅立叶变换的关系，也更加清晰了。
2. 同时，在实验过程之中，在编程实际中，通过学习书后的代码，和网上共享的资料，我渐渐地对MATLAB编程熟悉起来，而且编程思维和编程能力都有大幅提高，以后在漫长的学习生活之中，相信对matlab的实用会更加驾轻就熟，使得工程实际的效率进一步的提高。
3. 在实验后，对所得到的实验结果进一步分析，尤其是对不同采样率，不同DFT点数，不同数据点数下DFT的结果的认识更加清楚，对上课的时候老师讲的三个同学也有更深的体会。另外对于周期函数的抽样时间长度对于频域的卷积sa函数的宽度的影响，抽样后的周期函数与理想的周期函数的傅立叶变换之间的联系与不同等等。
4. 这次实验让我独立认识问题，思考问题，解决问题，攻关克难的能力都有大的提升。

6. 附录：程序代码：

x1=[1 1 1 1];
N=8;
xk=fftshift(fft(x1,N));
subplot(3,1,1);%用于分割
stem(0:N-1,[x1 zeros(1,N-4)]);
xlabel('n');
ylabel('x(n)');
title('x1n的离散时域信号');
 subplot(3,1,2);
stem(0:N-1,abs(xk));
xlabel('k');
ylabel('xk');
title('x1的N点DFT的结果');
subplot(3,1,3);
plot(0:N-1,abs(xk));
xlabel('w');
ylabel('X(jw)');
title('xk的包络');

x2=[1 2 3 4 4 3 2 1];
N=8;
xk=fftshift(fft(x2,N));
subplot(3,1,1);%用于分割
stem(0:N-1,[x2 zeros(1,N-8)]);
xlabel('n');
ylabel('x(n)');
subplot(3,1,2);
title('x2n的离散时域信号');
stem(0:N-1,abs(xk));
xlabel('k');
ylabel('xk');   
title('x2的N点DFT的结果');
subplot(3,1,3);
plot(0:N-1,abs(xk));
xlabel('w');
ylabel('X(jw)');
title('xk的包络');

x3=[4 3 2 1 1 2 3 4];
N=8;
xk=fftshift(fft(x3,N));
subplot(3,1,1);%用于分割
stem(0:N-1,[x3,zeros(1,N-8)]);
xlabel('n');
ylabel('x(n)');
title('x3n的离散时域信号');
subplot(3,1,2);
stem(0:N-1,abs(xk));
xlabel('k');
ylabel('xk');
title('x3的N点DFT的结果');
subplot(3,1,3);
plot(0:N-1,abs(xk));
xlabel('w');
ylabel('X(jw)');
title('xk的包络');

n=[0:31];
x4=cos(pi/8*n);
N=16;
xk=fftshift(fft(x4,N));
subplot(3,1,1);%用于分割
stem(0:31,x4);
xlabel('n');
ylabel('x(n)');
title('x4的离散时域信号');
subplot(3,1,2);
stem(0:N-1,abs(xk));
xlabel('k');
ylabel('xk');
title('x4的N点DFT的结果');
subplot(3,1,3);
plot(0:N-1,abs(xk));
xlabel('w');
ylabel('X(jw)');
title('xk的包络');

n=[0:31];
x4=cos(pi/8*n);
N=32;
xk=fftshift(fft(x4,N));
subplot(3,1,1);%用于分割
stem(0:31,x4);
xlabel('n');
ylabel('x(n)');
title('x4的离散时域信号');
subplot(3,1,2);
stem(0:N-1,abs(xk));
xlabel('k');
ylabel('xk');
title('x4的N点DFT的结果');
subplot(3,1,3);
plot(0:N-1,abs(xk));
xlabel('w');
ylabel('X(jw)');
title('xk的包络');

n=[0:31];
x5=sin(pi/8*n);
xk1=fftshift(fft(x5,16));
xk2=fftshift(fft(x5,32));
subplot(2,3,1);
stem(0:31,x5);
xlabel('n');
ylabel('x(n)');
title('x5的离散时域信号');
subplot(2,3,2);
stem(0:15,abs(xk1));
xlabel('k');
ylabel('xk1');
title('x5的16点DFT的结果');
subplot(2,3,3);
plot(0:15,abs(xk1));
xlabel('w');
ylabel('X(jw)');
title('xk的包络');

subplot(2,3,4);
stem(0:31,x5);
xlabel('n');
ylabel('x(n)');
title('x5的离散时域信号');
subplot(2,3,5);
stem(0:31,abs(xk2));
xlabel('k');
ylabel('xk2');
title('x5的32点DFT的结果');
subplot(2,3,6);
plot(0:N-1,abs(xk2));
xlabel('w');

fs=64;
delta_t=1/fs;
N=128;
t=[0:N-1]*delta_t;
x6=cos(8*pi*t)+cos(16*pi*t)+cos(20*pi*t);
xk1=fft(x6,16);
xk2=fft(x6,32);
xk3=fft(x6,64);
subplot(3,3,1);
stem(0:127,x6);
xlabel('n');
ylabel('x(n)');
title('x6的离散时域信号');
subplot(3,3,2);
stem(0:15,abs(xk1));
xlabel('k');
ylabel('xk1');
title('x6的16点DFT的结果');
subplot(3,3,3);
plot(0:15,abs(xk1));
xlabel('w');
ylabel('X(jw)');
title('xk的包络');

subplot(3,3,4);
stem(0:127,x6);
xlabel('n');
ylabel('x(n)');
title('x6的离散时域信号');
subplot(3,3,5);
stem(0:31,abs(xk2));
xlabel('k');
ylabel('xk2');
title('x6的32点DFT的结果');
subplot(3,3,6);
plot(0:31,abs(xk2));
xlabel('w');
ylabel('X(jw)');
title('xk的包络');

subplot(3,3,7);
stem(0:127,x6);
xlabel('n');
ylabel('x(n)');
title('x6的离散时域信号');
subplot(3,3,8);
stem(0:63,abs(xk3));
xlabel('k');
ylabel('xk2');
title('x6的64点DFT的结果');
subplot(3,3,9);
plot(0:63,abs(xk3));
xlabel('w');
ylabel('X(jw)');

n=[0:31];
x7=cos(pi/8*n)+j*sin(pi/8*n);
xk1=fftshift(fft(x7,16));
xk2=fftshift(fft(x7,32));
xk3=fftshift(fft(x7,64));
subplot(3,3,1);%用于分割
stem(0:31,x7);
xlabel('n');
ylabel('x(n)');
title('x7的离散时域信号');
subplot(3,3,2);
stem(0:15,abs(xk1));
xlabel('k');
ylabel('xk');
title('x7的16点DFT的结果');
subplot(3,3,3);
plot(0:15,abs(xk1));
xlabel('w');
ylabel('X(jw)');
title('xk的包络');

subplot(3,3,4);%用于分割
stem(0:31,x7);
xlabel('n');
ylabel('x(n)');
title('x7的离散时域信号');
subplot(3,3,5);
stem(0:31,abs(xk2));
xlabel('k');
ylabel('xk');
title('x7的32点DFT的结果');
subplot(3,3,6);
plot(0:31,abs(xk2));
xlabel('w');
ylabel('X(jw)');
title('xk的包络');

subplot(3,3,7);%用于分割
stem(0:31,x7);
xlabel('n');
ylabel('x(n)');
title('x7的离散时域信号');
subplot(3,3,8);
stem(0:63,abs(xk3));
xlabel('k');
ylabel('xk');
title('x7的64点DFT的结果');
subplot(3,3,9);
plot(0:63,abs(xk3));
xlabel('w');

      x2=[1 2 3 4 4 3 2 1];
x3=[4 3 2 1 1 2 3 4];
N1=8;
N2=16;
xk21=fftshift(fft(x2,N1));
xk31=fftshift(fft(x3,N1));
xk22=fftshift(fft(x2,N2));
xk32=fftshift(fft(x3,N2));

subplot(2,3,1);%用于分割
stem(0:N1-1,[x2,zeros(1,N1-8)]);
xlabel('n');
ylabel('x2(n)');
title('x2n的离散时域信号');

subplot(2,3,2);
stem(0:N1-1,abs(xk21));
xlabel('k');
ylabel('xk2');
title('x2的8点DFT的结果');

subplot(2,3,3);
stem(0:N2-1,abs(xk22));
xlabel('k');
ylabel('xk22');
title('x2的16点DFT的结果');

subplot(2,3,4);
stem(0:N1-1,[x3,zeros(1,N1-8)]);
xlabel('n');
ylabel('x3(n)');
title('x3n的离散时域信号');

subplot(2,3,5);
stem(0:N1-1,abs(xk31));
xlabel('k');
ylabel('xk3');
title('x2的8点DFT的结果');

subplot(2,3,6);
stem(0:N2-1,abs(xk32));
xlabel('k');
ylabel('xk3');