数字信号处理(DSP)实验——离散傅立叶变换及谱分析
一、实验目的
- 学习如何在计算机上实现离散傅里叶变换的方法。
- 学习如何计算序列的圆周卷积的方法。
- 通过学习使用DFT对连续信号和时域离散信号进行谱分析的方法，并了解可能出现的分析误差来正确应用DFT于实际中。
- 深入理解使用FFT分析周期序列频谱面临的挑战并掌握解决方法。
- 学习如何通过加权处理法进行时域窗函数的应用。
- 深入理解使用FFT分析非周期信号频谱面临的挑战并掌握解决方法。
二、实验原理与方法
１. 对周期序列进行频谱分析应注意的问题
对时间序列进行快速傅里叶变换（FFT）时，默认会对其进行无限周延处理（若选择处理较长的时间段还需先对其进行有限窗口截取）。无限时间序列是一种无始无终的数据流，在特定条件下可以通过有限窗口采样得到其片段特征。当所选采样区间恰好是该信号基本周长的整数倍时，在完成周延处理后，则能够还原出原始的时间域信号及其完整的频率特性信息；反之，则当采样区间并非该周长的基本整数倍时，在完成周延处理后，则无法恢复原始的时间域信号及其完整的频率特性信息；此时所得到的结果中存在频率信息偏差，并且这种偏差程度与所选取的时间窗长度直接相关（如图2-1所示），其中幅度频谱量值表示为20log|X(k)|（单位：dB），这种现象源于时域上的矩形窗函数与原信号进行乘积操作而产生的频域卷积效应；而用于描述这一现象的技术称谓则为矩形窗函数
如图2-2所示：
除了在k=0处主波内集中了绝大部分的能量外，在两侧较低峰值区域出现了一些能量分布.这些余波与周期序列频谱卷积的结果导致原频谱带宽增加.这被称为频率泄露现象.对于已知具有明确周期性的信号进行频谱分析时,可以选择整数个基波周期作为截取长度,从而能够有效消除这一现象的发生.然而,由于缺乏先验信息难以预判其确切周期,因此为了尽量避免频率泄露对结果的影响,就需要选择其频谱旁瓣较小的部分作为截取函数来进行处理.
- 时域窗函数的应用
被用作时间截断工具的矩形窗口，在执行时间截断操作时，并不影响该时间段内的信号；而其他类型的窗口会对处理对象进行加权处理。除了三角形窗口（Triangular window）、汉宁（Hanning）窗口和哈明（Hamming）窗口之外，在工程应用中还有许多其他常用的窗口类型可供选择。例如帕塞瓦尔（Parzen）窗口、凯泽（Kaiser）窗口以及切比雪夫（Chebyshev）等带 Reject 窗口等。本次实验中采用了一系列窗口作为时间加权截断的方法。
图 2-3 展示了在非整周期时域上对正弦函数施加 Hanning 窗进行加权截断后得到的波形及其频谱分析结果。相较于图 2-1(b)，泄漏现象已显著减弱。
3、对非周期序列进行频谱分析应注意的问题
(１)混叠
(２)泄漏
（３）栅栏效应
非周期信号x(t)的频谱应为连续的。在对x(t)进行抽样并应用DFT后得到的是离散化的频谱。若忽略混叠和泄漏等误差的影响，则计算结果仅给出x(t)连续频谱上的(N/2)+1个采样点。这类似于通过均匀间隔的栅栏缝隙观察另一边景象的效果。因此我们称这种现象为栅栏效应。被栅栏遮蔽的部分可能存在与显式显示差异较大的变化特征即信号频谱中的某些特定分量可能被遮蔽而不明显显现出来的情况为了尽可能地显现这些被遮蔽的部分可以通过增加时域内的零填充来实现这一目的具体操作方法是在不改变原始信号时域抽样点数量的前提下向时间序列末尾添加若干个零值从而增大时域序列总长度N这种方法被称为补零重构例如原始信号抽样得到12个数据点在其后面附加4个零值使总数据长度变为16个这样处理不会改变原始信号的时间间隔及其频率特性但会使频域中的相邻采样点之间的间距由原来的fs/12变为fs/16这样一来通过补零重构能够提高频域内的采样密度从而使得原本不易于观察到的一些频率位置上的频谱特征得以显现
二、实验内容
1、预先准备：先给出四个扩展函数
(1) 离散傅立叶变换，采用矩阵相乘的方法
function [Xk]=dft(xn,N)
n=[0:1:N-1];
k=[0:1:N-1];
WN=exp(-j*2*pi/N);
nk=n'*k;
WNnk=WN.^(nk);
Xk=xn*WNnk;(2) 逆离散傅立叶变换
function [xn]=idft(Xk,N)
n=[0:1:N-1];
k=[0:1:N-1];
WN=exp(-j*2*pi/N);
nk=n'*k;
WNnk=WN.^(-nk);
xn=(Xk*WNnk)/N;(3) 实序列的分解
% 实序列可分解为共扼对称分量
%和共扼反对称分量
function [xec,xoc]=circevod(x)
N=length(x);
n=0:(N-1);
xec=0.5*(x+x(mod(-n,N)+1)); %根据上面的公式计算，mod函数为取余
xoc=0.5*(x-x(mod(-n,N)+1));(4) 序列的循环移位y(n)=x((n-m))N
function y=cirshftt(x,m,N)
if length(x)>N
error('N mustbe >= the length of x') %要求移位周期大于信号长度
end
x=[x zeros(1,N-length(x))];
n=[0:1:N-1];
n=mod(n-m,N);
y=x(n+1);三、实验内容及步骤
(1) 绘制H(ejw)的幅度谱|H(ejw)|和相位谱。
clc;%本语句的作用是清除命令执行界面中所有的输出信息
clear ;%清除workspace中所有的变量
clf;%清除所有的绘图内容(如果本次程序执行前已经有绘图窗口存在，则可能将本程序将要绘制的图形绘制到之前的窗口中，可能导致疑惑)
w=0:0.01:2*pi;%将连续w分成极小的间隔
hjw=1+exp(-j*w)+exp(-2*j*w)+exp(-3*j*w)+exp(-4*j*w)+exp(-5*j*w)+exp(-6*j*w;
mag=abs(hjw);%取模/幅度值
ang=angle(hjw);%相位
figure(1);%第一个窗口
subplot(211)%将第一个窗口分成2行1列的子图，分的方法是从左到右，从上到下
plot(w,mag);%以w为横坐标，幅度为纵坐标绘制幅度谱
set(gca,'Xtick',[0:pi/8:2*pi]);
% set函数 将当前图形（gca）的x轴坐标刻度（xtick）标志为：
set(gca,'XtickLabel',{'0','pi/8','pi/4','3*pi/8','pi/2','5*pi/8','3*pi/4','7*pi/8','pi','9*pi/8','10*pi/8','11*pi/8','3*pi/2','13*pi/8','7*pi/4','15*pi/8','2*pi'});%此处对X轴标注
title('幅度谱');%设定第1幅子图的标题
subplot(212);%第2个子图
plot(w,ang);%显示相位
set(gca,'ytick',[-pi:pi/4:pi]);
set(gca,'Xtick',[0:pi/8:2*pi]);
set(gca,'XtickLabel',{'0','pi/8','pi/4','3*pi/8','pi/2','5*pi/8','3*pi/4','7*pi/8','pi','9*pi/8','10*pi/8','11*pi/8','3*pi/2','13*pi/8','7*pi/4','15*pi/8','2*pi'});
title('相位谱');%设定第2幅子图的标题
grid on %可以在同一个子图里面叠绘多张图结果：
给定系统函数 H(z) = 1 - Z^{-N}, 利用 MATLAB 绘制 8 阶系统的零极点分布图、幅频特性曲线以及相位频率特性曲线。
b=[1 0 0 0 0 0 0 0 -1]; %H(z)的分子多项式系数矢量
a=1; %H(z)的分母多项式系数矢量
subplot(1,3,1);,
zplane(b,a); %绘制H(z)的零极点图
[H,w]=freqz(b,a); %计算系统的频率响应
subplot(1,3,2);
plot(w/pi,abs(H)); %绘制幅频响应曲线
axis([0,1,0,2.5]);
xlabel('\omega/pi');
ylabel('|H(e^j^\omega)|');
subplot(1,3,3);
plot(w/pi,angle(H)); %绘制相频响应曲线
xlabel('\omega/pi');
ylabel('\phi(\omega)');结果：
(3) 分析实序列的特性
例1 本例用于分析实序列的特性, 设定x(n)=10×(0.8)^n
其中, n满足条件: 0≤n≤10
将x(n)分为其共轭对称分量及共轭反对称分量
%例1 本例检验实序列的性质DFT[xec(n)]=Re[X(k)]DFT[xoc(n)]=Im[X(k)]
%设 x(n)=10*(0.8).^n 0<=n<=10 将x(n)分解为共扼对称及共扼反对称部分
n=0:10;
x=10*(0.8).^n;
[xec,xoc]=circevod(x);
subplot(2,1,1);stem(n,xec); %画出序列的共扼对称分量
title('Circular -even component')
xlabel('n');ylabel('xec(n)');axis([-0.5,10.5,-1,11])
subplot(2,1,2);stem(n,xoc); %画出序列的共扼反对称分量
title('Circular -odd component')
xlabel('n');ylabel('xoc(n)');axis([-0.5,10.5,-4,4])
figure(2)
X=dft(x,11); %求出序列的DFT
Xec=dft(xec,11); %求序列的共扼对称分量的DFT
Xoc=dft(xoc,11); %求序列的共扼反对称分量的DFT
subplot(2,2,1);stem(n,real(X));axis([-0.5,10.5,-5,50])
title('Real{DFT[x(n)]}');xlabel('k'); %画出序列DFT的实部
subplot(2,2,2);stem(n,imag(X));axis([-0.5,10.5,-20,20])
title('Imag{DFT[x(n)]}');xlabel('k'); %画出序列DFT的虚部
subplot(2,2,3);stem(n,Xec);axis([-0.5,10.5,-5,50])
title('DFT[xec(n)]');xlabel('k');
subplot(2,2,4);stem(n,imag(Xoc));axis([-0.5,10.5,-20,20])
title('DFT[xoc(n)]');xlabel('k');结果：
(4) 完成两个序列的圆周卷积运算，在此过程中使用以下参数：x₁(n) = [2,…]、x₂ = [1,…]、N = 8
示例说明的是用于完成两个信号之间关系分析的算法设计。
在运行程序之前，请在命令窗口中输入x₁,x₂和N的具体数值。
请自行选择两个不同的信号进行计算，并记录实验结果。
% 例2 本例为计算序列的圆周卷积程序
% 运行之前应在命令窗口输入 x1,x2,N 的值
% 要求：自己选择2个序列进行计算，将实验结果写出。
clear;
clc;
clf;
x1=[2,2,2,2,2,2];
x2=[1,1,1,1,1,1,1];
N=8;
if length(x1)>N
error('N mustbe >= the length of x1')
end
if length(x2)>N
error('N mustbe >= the length of x2')
end
x1=[x1 zeros(1,N-length(x1))]; %将x1,x2补0成为N长序列
x2=[x2 zeros(1,N-length(x2))];
m=[0:1:N-1];
x2=x2(mod(-m,N)+1); %该语句的功能是将序列翻褶，延拓，取主值序列
H=zeros(N,N);
for n=1:1:N %该for循环的功能是得到x2序列的循环移位矩阵
H(n,:)=cirshftt(x2,n-1,N); %和我们手工计算圆周卷积得到的表是一致的
end
y=x1*H' %用矩阵相乘的方法得到结果输出结果为：
y= 10 10 10 10 10 12 12 10基于补零技术的离散傅里叶变换分析
clear;
clc;
clf;
n=0:4;
x=[ones(1,5)]; %产生矩形序列
k=0:999;w=(pi/500)*k;
X=x*(exp(-j*pi/500)).^(n'*k); %计算离散时间傅立叶变换
Xe=abs(X); %取模
subplot(2,2,1);stem(n,x);ylabel('x(n)'); %画出矩形序列
subplot(2,2,2);plot(w/pi,Xe);ylabel('|X(ejw)|'); %画出离散时间傅立叶变换
N=20;x=[ones(1,5),zeros(1,N-5)]; %将原序列补零为10长序列
n=0:1:N-1;
X=dft(x,N); %进行DFT
magX=abs(X);
k=(0:length(magX)'-1)*N/length(magX);
subplot(2,2,3);stem(n,x);ylabel('x(n)'); %画出补零序列
subplot(2,2,4);stem(k,magX); %画出DFT结果
axis([0,20,0,5]);ylabel('|X(k)|');结果：
实验结果说明：通过将序列进行延展处理后执行L点快速傅里叶变换（DFT），其结果等价于将原信号频谱在[0, 2π)区间内均匀地采样了L次，并因此提高了采样密度并重新定位了采样位置；这使得能够更清晰地反映原信号频谱的细节特征。另外一种方法是直接对延展后的序列执行DFT运算，在此过程中会显著增加高频分量的数量，并进一步缩小频率分辨率
该序列为x(n)=2 cos(0.35πn)+cos(0.5πn)，经周期计算得出其实为40点序列；
列出包含10个有效采样点的DFT程序；
完成对第一题中长度为10的有效采样点序列进行补零至40点，并计算其DFT；
当n=0时取值为39，请完成对包含4O个有效采样点序列进行DFT运算；
要求：
（1）请完成将长度为1O的有效采样点序列扩展至4O后进行DFT运算的任务；
（2）请完成对包含4O个有效采样点进行直接DFT运算的任务；
（3）将所得实验曲线绘制出来，并分析实验结果所反映的现象
%有10个有效采样点序列的DFT程序如下：
M=10;
n=0:M-1;
x=2*cos(0.35*pi*n)+cos(0.5*pi*n);
subplot(2,1,1);stem(n,x);title('没有足够采样点的信号');
Y=dft(x,M);
k1=0:M-1;w1=2*pi/M*k1;
subplot(2,1,2);stem(w1/pi,abs(Y));title('信号的频谱');
M=10;(1)将有10个有效采样点序列补零到40长后计算DFT
%10个有效采样点序列补零到40后的DFT
clear
clc
clf
M=10;
n=0:M-1;
b=40;
s=0:b-1;
x=[2*cos(0.35*pi*n)+cos(0.5*pi*n), zeros(1,30)];
subplot(2,1,1);stem(s,x);title('补零到40点的信号');
Y=dft(x,b);
k1=0:b-1;w1=2*pi/b*k1;
subplot(2,1,2);stem(w1/pi,abs(Y));title('信号的频谱');结果：
(2)40个有效采样点的序列的DFT
%40个有效采样点的序列的DFT
clear
clc
clf
M=40;
n=0:M-1;
x=2*cos(0.35*pi*n)+cos(0.5*pi*n);
subplot(2,1,1);stem(n,x);title('采样点为40的信号');
Y=dft(x,M);
k1=0:M-1;w1=2*pi/M*k1;
subplot(2,1,2);stem(w1/pi,abs(Y));title('信号的频谱');结果：
实验结果表明，在信号补零之后生成了高密度谱，在这种情况下尽管频谱线非常密集但因为信号的有效长度并未改变所以其分辨能力也未提升因此难以明确识别出信号的频谱成分而相比之下通过增加采样点来构建高分辨率谱则能够显著提升分辨能力从而清晰地识别出信号的组成成分
对余弦函数cos(2t)进行采样操作时，在Ts=1/53秒的时间间隔下完成N=128点的快速傅里叶变换（FFT），随后绘制频谱曲线图以分析信号频谱特性。通过施加哈明窗和三角形窗两种加权处理方法对频域数据进行截断处理，并观察分析频率泄露现象的变化特征及衰减过程。
clear
clc
clf
n=[0:127];
N=128;
Ts=1./53;
t=n.*Ts;
xn=cos(2.*pi.*t);
w1=hamming(N);
w2=bartlett(N);
H=xn.*w1';
B=xn.*w2';
subplot(3,2,1);
stem(n,xn,'.');
title('xn函数')
subplot(3,2,2);
stem(abs(fft(xn,N)),'.');
title('xn频谱曲线')
subplot(3,2,3);
stem(n,H,'.');
title('Hamming窗加权')
subplot(3,2,4);
stem(abs(fft(H,N)),'.');
title('Hamming窗加权频谱曲线')
subplot(3,2,5);
stem(n,B,'.');
title('三角窗加权')
subplot(3,2,6);
stem(abs(fft(B,N)),'.');
title('三角窗加权频谱曲线')结果：
N=128
N=106
N=1000
为了避免频率泄漏对结果产生的显著影响，在进行时间截断时应选择其频谱旁瓣相对较小的截断函数以减小泄漏影响。
为了尽可能使被栅栏遮挡的部分显现出来可以通过不增加信号采样点而增加频域样点密度的方法实现即通过时域补零扩展变换窗口宽度称为补零重构。
四、实验总结
实验结果说明：
(1)两个长度不大于N的序列进行N点圆周卷积后所得序列长度仍为N；圆周卷积即为周期卷积取主值序列；
(2)高密度谱是通过在信号补零后获得的一种频谱表示形式，在这种情况下频谱线较为密集但无法提高信号有效频带宽度；因此难以准确识别信号频谱中的各个成分；而高分辨率谱则是通过延长信号有效时长来实现频率分辨率提升从而能够更清晰地辨识信号中的频谱成分；
(3)对偶序列xec(n)和xoc(n)，其DFT结果分别对应复数X(k)的实部与虚部；
(4)在对信号进行补零操作后再执行DFT变换时频谱线数量会增加且频谱间隔减小这将有助于更详细地观察原始信号频谱中的细节信息；
(5)当对非周期性信号进行傅里叶分析时为了避免频率泄露带来的影响建议在截断非周期函数时选择具有较小旁瓣宽度的时间窗函数以减少泄漏影响；在此基础上应选择与被分析信号长度相同的窗函数宽度并在做变换前需补零的情况下建议将原信号与窗函数相乘后再进行补零操作即不采用加权处理对被补零的部分进行处理