数字信号处理(DSP)实验——离散傅立叶变换及谱分析

一、实验目的

1.掌握离散傅里叶变换的计算机实现方法。
2.掌握计算序列的圆周卷积的方法。
3.学习用DFT对连续信号和时域离散信号进行谱分析的方法，了解可能出现的分析误差，以便在实际中正确应用DFT。
4. 理解用FFT对周期序列进行频谱分析时所面临的问题并掌握其解决方法。
5.掌握用时域窗函数加权处理的技术。
6.理解用FFT对非周期信号进行频谱分析所面临的问题并掌握其解决方法。

二、实验原理与方法

１. 对周期序列进行频谱分析应注意的问题

对时间序列作FFT时，实际上要作周期延拓（如果取长序列的一段进行计算还要先作截断）。周期序列是无限长时间序列，如果截断区间刚好就是该序列周期的整数倍，那么在进行周期延拓后，将还原出原来的周期序列，由此可以较精确地计算出的该周期序列的频谱。反之，如果截断区间并不是该序列周期的整数倍，那么在进行周期延拓后，就不可能还原出原来的周期序列，由此计算出的频谱与该周期序列的频谱存在误差，而且误差的大小与截断区间的选取直接相关，如图2-1所示，其中幅度频谱的量值表示为20log|X(k)|，以dB(分贝)为单位。这种误差是由于周期序列与矩形截断序列相乘在频域产生二者的频谱卷积形成的。矩形窗的频谱是抽样函数序列
如图2-2所示：

除了k = 0处主瓣内集中了大部分的能量外，两旁的较小峰值处的旁瓣也分散了一部分能量，它与周期序列频谱卷积的结果使原来集中的频谱展宽，称为频率泄漏。
如果对已知周期的信号作频谱分析，在进行时域截断时，完全可以选取其周期的整数倍裁取，从而可以避免这种频率泄漏的发生。不过，通常需要进行频谱分析的信号是周期未知的信号，或随机信号，无法判断它的周期值，为了尽量避免频率泄漏对结果的影响，在作时间截断时，就应选取其频谱的旁瓣较小的截断函数，以减轻泄漏问题。

2、时域窗函数的应用
作为截断函数，矩形窗在作时间截断时，对所截取区间内的信号不加以任何影响，而其它的窗函数都将对所截取区间内的信号作加权处理。除三角窗、Hanning窗和Hamming窗外，常用的窗函数还有很多，例如Parzen窗、Kaiser窗、Chebyshev窗、Tukey窗、Poisson窗、Caushy窗、Gaussian窗和Blackman窗等等。本次实验利用窗函数作时域加权截断。

图 2-3 中给出了采用Hanning窗对正弦函数作非整周期的时域加权截断后的波形和频谱，与图2-1(b)比较，泄漏已明显减少。

3、对非周期序列进行频谱分析应注意的问题
(１)混叠

(２)泄漏

(３)栅栏效应
非周期信号x(t)应具有连续的频谱，在对x(t)作抽样后进行DFT，得到的是离散的频谱。如果排除混叠和泄漏等误差的影响，所得的结果也只是x(t)的连续频谱上的(N/2)+1个样值。这就象通过栅栏上的等间距缝隙观看到的另一边的景象，故此称栅栏效应。被栅栏遮住的景象中有可能存在与显现出的频谱差异较大的变化，即显示信号特征的频谱分量。为了使被栅栏遮住的部分能尽可能地显现出来，可以采用增加频域样点密度的方法，即在不增加信号采样点的情况下，用时域补零加宽变换尺度N来实现，称为补零重构。例如原来信号采样得到12个样点，在其后面再加上4个零，使序列的总长度为16个样点。这样处理的结果原来信号的采样间隔和频率都没有改变，设采样频率为fs，经补零重构之后，采样频率仍然为fs，但是原来频域样点间宽度为 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(e^jw)的幅度谱|H(e^jw)|和相位谱。

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)