一、实验目的
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(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 %可以在同一个子图里面叠绘多张图
结果：
(2)已知系统函数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
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) 计算序列的圆周卷积，其中x1(n)=[2,2,2,2,2,2]，x2=[1,1,1,1,1,1,1]，N=8
例2 本例为计算序列的圆周卷积程序
运行之前应在命令窗口输入 x1,x2,N 的值
要求：自己选择2个序列进行计算，将实验结果写出。
% 例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
(5) 补零序列的离散傅立叶变换
序列x(n)=R5(n)，已给出序列的傅立叶变换程序和将原序列补零到10长序列的DFT
要求：
(1)编写将序列补零到20长后计算DFT的程序
(2)给出实验结果
(3)实验结果说明什么（即序列补零后进行DFT,频谱有何变化）
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，频谱的谱线增多，谱间隔变小。
(6) 有一个序列为x(n)=2 cos(0.35πn)+cos(0.5πn)（该序列周期计算可得40）
下面给出有10个有效采样点序列的DFT程序；
请写出将第一问中的10长序列补零到40长，计算其DFT；
采样n=0：39，计算有40个有效采样点的序列的DFT；
要求：
(1) 请编写将有10个有效采样点序列补零到40长后计算DFT的程序
(2) 请编写计算有40个有效采样点的序列的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('信号的频谱');
结果：
(3)实验结果说明高密度谱是信号补零后得到的，虽然谱线相当密，但是因为信号有效长度不变，所以其分辨率也不变，因此还是很难看出信号的频谱成分。高分辨率谱是将信号有效长度加长即增大采样点，因此分辨率提高，可以看出信号的成分。
高分辨率谱与高密度谱差异的比较：
1.高密度谱是在原有的序列后插0。
2.高分辨率谱是增大采样点。
3.高密度谱呈许多谱线型，而且当补充0越多，谱线也越密集。
4.而高分辨率谱则在取样点达到一定程度后，谱线一定了，也没有那种密集度。
(7) 将余弦函数cos(2t)以Ts=1/53 s抽样，对余弦序列做样点数为N=128的FFT，画出频谱曲线，观察并记录频率泄漏现象，然后用Hamming窗和三角窗作加权截断，观察并记录泄漏的衰减。
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
为了尽量避免频率泄漏对结果的影响，在作时间截断时，就应选取其频谱的旁瓣较小的截断函数，以减轻泄漏问题。
为了使被栅栏遮住的部分能尽可能地显现出来，可以采用增加频域样点密度的方法，即在不增加信号采样点的情况下，用时域补零加宽变换尺度N来实现，称为补零重构。
四、实验总结
实验结果说明：
(1)两个长度小于等于N的序列的N点圆周卷积长度仍为N；圆周卷积正是周期卷积取主值序列；
(2)高密度谱是信号补零后得到的，虽然谱线相当密，但是因为信号有效长度不变，所以其分辨率也不变，因此很难看出信号的频谱成分。高分辨率谱是将信号有效长度加长，因此分辨率提高，可以看出信号的成分。
(3)DFT[xec(n)]=Re[X(k)]DFT[xoc(n)]=Im[X(k)]；
(4)序列补零后进行DFT，频谱的谱线增多，谱间隔变小，将会显示出原信号频谱更多的细节；
(5)对非周期序列进行频谱分析时，为了尽量避免频率泄漏对结果的影响，在对非周期函数作时间截断时，除尽量增加截断序列的宽度外，也应选取其频谱的旁瓣较小的截断函数，以减轻泄漏问题。在选取了适当的窗函数后，应当使窗函数的宽度与被处理的序列长度相同，如果作变换前还需要补零，则应将原序列与窗函数相乘后再补零，即补零的样点不用窗函数加权处理。
3321
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
