华南理工大学信号与系统实验二

实验二利用DFT 分析离散信号频谱
实验日期：评 分：
一、实验目的
应用离散傅里叶变换(DFT)，分析离散信号x [k ]。

深刻理解利用DFT 分析离散信号频谱的原理，掌握改善分析过程中产生的误差的方法。

二、实验原理
根据信号傅里叶变换建立的时域与频域之间的对应关系，可以得到有限长序列的离散傅里叶变换(DFT)与四种确定信号傅里叶变换的之间的关系，实现由DFT 分析其频谱。

三、实验内容
1.利用FFT 分析信号
的频谱；(1)确定DFT 计算的参数；
(2)进行理论值与计算值比较，讨论信号频谱分析过程中误差原因及改善方法。

【答题】
（1）角频率πω83=对应的周期为k k N 3
162==ωπ，取3=k 得到16=N ，所以fft 取n=16
（2）计算值的相角与理论值不符合，经过分析，我找到了修补的方法。

下面的代码中使用将微小数值置零的方法，消除计算机对浮点数存储误差的影响。

实验结果的图片（图1）就是正确的频谱图。

如果不使用这种方法，由于辐角对模长很小的复数存储所带来的精度误差敏感，将呈现出错误的辐角数值（见图2）。

【代码】
%[1]
k = 0:31;
x = cos(3*pi/8 * k);
X = fft(x,16);
subplot(3,1,1)
stem(k,x)
xlabel('k');ylabel('Signal x[k]')
subplot(3,1,2)
stem(-8:7, abs(fftshift(X)))
xlabel('k');ylabel('Spectrum |X[k]|')
subplot(3,1,3)
a = fftshift(X);
a(abs(a)<1e-5)=0;
a = angle(a);
a(a<1e-5)=0;
stem(-8:7, a)
xlabel('k');ylabel('Spectrum arg(X[k])')
angle(fftshift(X))【结果截图】
图1第一题的正确结果图
图2第二题的错误结果图
【结果分析】
观察图1，得知信号)83cos(][k k x π=只具有三次谐波分量。

基波频率为8
0πω=。

2.利用FFT 分析信号
的频谱；(1)确定DFT 计算的参数；
(2)进行理论值与计算值比较，讨论信号频谱分析过程中误差原因及改善方法。

【答题】
（1）能量应不低于90%，即9.0)21()21(0
0max ≥∑∑∞+==k k k n k ，解得3max ≥n ，不妨取n=10（2）理论上计算所得的频谱应该是连续的，但算出来却是离散的。

图3中红线是理论频谱曲线，棒图示计算所得的频谱。

造成这一现象的原因是n 值取得不够大。

当n 趋向于无穷时，采样所得的信号变为非周期的，那么频谱就变为连续的了，这一现象如图4中所示（取n=100）。

【代码】% [2]
n = 10;
k = 0:n-1;
x = 0.5.^k;
subplot(3,1,1)
stem(k,x)
xlabel('k');ylabel('Signal x[k]')
subplot(3,1,2)
w_ = -pi:0.01:pi*(1-2/n);
a_ = 1./(1-0.5*exp(-1j*w_));
hold on
stem(k-n/2,abs(fftshift(fft(x))))
plot(w_/(2*pi/n), abs(a_),'r')
xlabel('k');ylabel('Spectrum |X[k]|')
legend('FFT','predicted')
hold off
subplot(3,1,3)
hold on
stem(k-n/2,angle(fftshift(fft(x))))
plot(w_/(2*pi/n), angle(a_),'r')
legend('FFT','predicted')
hold off
xlabel('k');ylabel('Spectrum arg(X[k])')
【结果截图】
图3第二题结果图（n=10）
图4第二题结果图（n=100）
【结果分析】
当窗口取得很大，如n=100时，计算所得频谱趋向于连续。

该信号具有较多的低频分量。

3.有限长脉冲序列，利用FFT分析其频谱，并绘出其幅度谱与相位谱。

【代码】
% [3]
k = 0:5;
x = [2 3 3 1 0 5];
subplot(3,1,1)
stem(k,x)
xlabel('k');ylabel('Signal x[k])')
X = fftshift(fft(x));
subplot(3,1,2)
stem(k-3,abs(X))
xlabel('k');ylabel('Spectrum |X[k]|')
subplot(3,1,3)
stem(k-3,angle(X))
xlabel('k');ylabel('Spectrum arg(X[k])')【结果截图】
图5第三题结果图
【结果分析】
上图5就是该信号的频谱图。

其中n=6。

4.某周期序列由3个频率组成：，
利用FFT 分析其频谱。

如何选取FFT 的点数N ？此3个频率分别对应FFT 计算结果X [m ]中的哪些点？若选取的N 不合适，FFT 计算出的频谱X [m ]会出现什么情况？
【答题】
N 的选取分析：因为)168cos()169cos()167cos(
][k k k k x πππ++=，故可以选取基波频率160πω=。

那么周期3216
/2==
ππN 【代码】
% [4]
N = 32;
k = 0:N-1;
x = cos(7*pi/16 * k) + cos(9*pi/16 * k) + cos(pi/2 * k);
X = fftshift(fft(x));
subplot(3,1,1)
stem(k,x)
xlabel('k');ylabel('Signal x[k]')
subplot(3,1,2)
stem(k-N/2,abs(X))
xlabel('k');ylabel('Spectrum |X[k]|')
subplot(3,1,3)
stem(k-N/2,angle(X))
xlabel('k');ylabel('Spectrum arg(X[k])')
【结果截图】
图6第四题频谱结果图（N=32）
图7第四题当N不恰当选取时的错误频谱（N=48）
【结果分析】
此三个频率分别对应k=±7、±8、±9共6个点。

若选取的N不恰当，那么频谱将不是正确的频谱，如图7所示。

N必须为信号周期的整数倍。

5.某离散序列由3个频率组成:
利用FFT 分析其频谱。

(1)对x [k ]做64点FFT ，绘出信号频谱，能分辨出其中的两个频率吗？
(2)对x [k ]补零到256点后计算FFT ，能分辨出其中的两个频率吗？
(3)选用非矩形窗计算FFT ，能够分辨出其中的两个频率吗？
(4)若不能够很好地分辨出其中的两个频谱，应采取哪些措施？
【答题】
（1）有两组频率分量特别突出，是4次和5次谐波。

计算所得对应的频率为846424ππω=⨯=和32
556425ππω=⨯=，分别与理论值的相对偏差为%25.6-和1.90%+。

较能符合实际，但仍有偏差。

（2）发现多个频率都有具有较大分量，不能分辨其中的两个频率。

（3）有两组频率分量特别突出，是4次和5次谐波。

与第一题相似，计算所得对应的频率为846424ππω=⨯=和32
556425ππω=⨯=，分别与理论值的相对偏差为%25.6-和1.90%+。

较能符合实际，但仍有偏差。

（4）如果不能增加采样点数，那么N=64做FFT 就是最优解。

如果可以增加采样点数，那么采样周期的倍数个点（300的正整数倍）就能得到最优解。

【代码】
（1）
% [5.1]
N = 64;
k = 0:N-1;
x = cos(2*pi/15 * k) + 0.75*cos(2.3*pi/15 * k);
X = fftshift(fft(x));
subplot(3,1,1)
stem(k,x)
xlabel('k');ylabel('Signal x[k]')
subplot(3,1,2)
stem(k-N/2,abs(X))
xlabel('k');ylabel('Spectrum |X[k]|')
subplot(3,1,3)
stem(k-N/2,angle(X))
xlabel('k');ylabel('Spectrum arg(X[k])')
（2）
% [5.2]
N = 64;
k = 0:N-1;
x = cos(2*pi/15 * k) + 0.75*cos(2.3*pi/15 * k); X = fftshift(fft(x,256));
subplot(3,1,1)
stem(k,x)
xlabel('k');ylabel('Signal x[k]')
subplot(3,1,2)
stem(-128:127,abs(X))
xlabel('k');ylabel('Spectrum |X[k]|')
subplot(3,1,3)
stem(-128:127,angle(X))
xlabel('k');ylabel('Spectrum arg(X[k])')
（3）
% [5.3]
N = 64;
k = 0:N-1;
x = cos(2*pi/15 * k) + 0.75*cos(2.3*pi/15 * k); % w = 1-abs(k)/(N+1); % bartlett window
w = (1+cos(pi*k/N))/2; % hanning window
x_ = x.*w;
X = fftshift(fft(x_));
subplot(5,1,1)
stem(k, x)
xlabel('k');ylabel('Signal x[k]')
subplot(5,1,2)
stem(k, w)
xlabel('k');ylabel('Window w[k]')
subplot(5,1,3)
stem(k, x_)
xlabel('k');ylabel('Windowed Signal x[k]*w[k]') subplot(5,1,4)
stem(k-N/2,abs(X))
xlabel('k');ylabel('Spectrum |X[k]|')
subplot(5,1,5)
stem(k-N/2,angle(X))
xlabel('k');ylabel('Spectrum arg(X[k])')
【结果截图】
图8第五题第一问结果图（N=64）
图9第五题第二问结果图（N=64，补零到256）
图10第五题第二问结果图（加Hanning窗）
【结果分析】
造成不能正确的区分两个频率的根本原因是信号点的数量不足，采样点数至少为300才能正确区分两个频率。

6.已知序列
利用FFT分析下列信号的幅频特性，频率范围为，N=500点。

(1)
(2)
(3)若将上述x[k]乘以cos(p k/2)，重做(1)和(2)。

【答题】
【代码】
（1和2小问）
% [6.1&6.2]
k = -50:50;
N = 500;
x = exp(-(0.1*k).^2/2);
M = 150;
x = [zeros(1,M) x zeros(1,M)];
for i = -50:50
y(i+51) = x(2*i+M+50+1);
g(i+51) = x(4*i+M+50+1);
end
Y = fftshift(fft(y,N));
G = fftshift(fft(g,N));
subplot(4,1,1)
stem(k, x(M:M+100),'.')
xlabel('k');ylabel('Signal x[k]')
subplot(4,2,3)
stem(k, y,'.')
xlabel('k');ylabel('Signal y[k]')
subplot(4,2,4)
stem(k, g,'.')
xlabel('k');ylabel('Signal g[k]')
subplot(4,2,5)
stem((-N/2:N/2-1)*(2*pi/N), abs(Y),'.') xlabel('w');ylabel('Spectrum |Y[w]|') subplot(4,2,6)
stem((-N/2:N/2-1)*(2*pi/N), abs(G),'.') xlabel('w');ylabel('Spectrum |G[w]|') subplot(4,2,7)
stem((-N/2:N/2-1)*(2*pi/N), angle(Y),'.') xlabel('w');ylabel('Spectrum arg(Y[w])') subplot(4,2,8)
stem((-N/2:N/2-1)*(2*pi/N), angle(G),'.') xlabel('w');ylabel('Spectrum arg(G[w])')
（3）
% [6.3]
k = -50:50;
N = 500;
x_ = exp(-(0.1*k).^2/2);
x = x_ .* cos(pi*k/2);
M = 150;
x = [zeros(1,M) x zeros(1,M)];
for i = -50:50
y(i+51) = x(2*i+M+50+1);
g(i+51) = x(4*i+M+50+1);
end
Y = fftshift(fft(y,N));
G = fftshift(fft(g,N));
subplot(5,1,1)
stem(k, x_,'.')
xlabel('k');ylabel('Signal x[k]')
subplot(5,1,2)
stem(k, x(M:M+100),'.')
xlabel('k');ylabel('Signal x[k]*cos(pi*k/2)') subplot(5,2,5)
stem(k, y,'.')
xlabel('k');ylabel('Signal y[k]')
subplot(5,2,6)
stem(k, g,'.')
xlabel('k');ylabel('Signal g[k]')
subplot(5,2,7)
stem((-N/2:N/2-1)*(2*pi/N), abs(Y),'.') xlabel('w');ylabel('Spectrum |Y[w]|') subplot(5,2,8)
stem((-N/2:N/2-1)*(2*pi/N), abs(G),'.') xlabel('w');ylabel('Spectrum |G[w]|') subplot(5,2,9)
stem((-N/2:N/2-1)*(2*pi/N), angle(Y),'.') xlabel('w');ylabel('Spectrum arg(Y[w])') subplot(5,2,10)
stem((-N/2:N/2-1)*(2*pi/N), angle(G),'.') xlabel('w');ylabel('Spectrum arg(G[w])')
【结果截图】
图11第六题第一二问结果图
图12第六题第三问的结果图
图13相位图的局部放大图案
【结果分析】
信号变化较快的地方对应较高的频率，通常在πω±=处，；信号变化缓慢的地方对应着较低的频率，通常在0=ω处。

且原信号与其频谱存在某种相反的关系，原信号宽度越大，则频谱图宽度缩小；频谱图宽度越大，原信号宽度越小。

观察频谱图预测，形如)exp(2x -的信号的频谱图的模长图具有与原信号相似的形式，相位图则具有锯齿波的样式（见图13）。

四、实验思考题
1.既然可直接由DTFT 定义计算序列DTFT ，为何利用DFT 分析序列的频谱？
答：由DTFT 的定义计算频谱，理论上得到的频谱图是连续的，而计算机只能处理离散信号，所以不能在计算机上实现。

而DFT 是有限区间上的变换，得到的频谱是离散的，在计算机上能轻易实现。

当周期N 取较大的值使得所占能量不低于90%就能很好地近似无穷区间的信号。

2.若序列持续时间无限长，且无解析表达式，如何利用DFT 分析其频谱？
答：用窗函数截取序列，得到有限周期的序列，再对这个序列应用DFT 即可近似分析出序列的频谱。

3.在利用DFT 分析离散信号频谱时，会出现哪些误差？如何克服或改善？
答：当窗函数的窗口宽度2M+1低于序列的周期、或不为周期的整数倍时，频谱不能很好地表征实际频谱。

采样所得序列长度小于序列周期时，也不能很好地表征实际频谱。

当频谱的模长很小时（如低于1e-8），应当把它置零，因为计算机对实数有存储误差，会造成相位角计算不准确。

4.在利用DFT 分析离散信号频谱时，如何选择窗函数？
答：最常用的窗函数是矩形窗，但可能会带来较明显的吉伯斯现象。

当不能很好地辨析出谱线时，可以尝试使用三角形窗或海宁窗，它们使得信号的变化平缓得多，从而实现增加一点失真作为代价来减少频谱的起伏的目的。

5.序列补零和增加序列长度都可以提高频谱分辨率吗？两者有何本质区别？
答：序列补零可以使得频谱线变密，但不能提高分辨率。

以增加采样点数的方式增长序列长度既可以使谱线变密又可以提高频谱的分辨率。

序列补零使得包络线不变，只是增加了谱线的密度。

而增加序列长度可以改变包络线的形状，提高谱线的分辨度。

五、实验收获
经过此次实验，使我认识到如何使用MATLAB 进行信号的离散傅里叶变换分析。

MATLAB 中的fft 函数提供了方便快捷的变换方法，使用十分方便简捷。