2013数字信号处理实验讲义 - 实验二
实验二 离散傅立叶变换及谱分析
一、 实验目的
1．掌握离散傅里叶变换的计算机实现方法。
2．检验序列离散傅里叶变换的性质。
3．掌握计算序列的循环卷积的方法。
4．学习用DFT 对连续信号和时域离散信号进行谱分析的方法，了解可能出现的分析误差，以便在实际中正确应用DFT 。
二、 实验内容
1．实现序列的离散傅里叶变换并对结果进行分析。（自己选择序列，要求包括复序列，实序列，实偶序列，实奇序列，虚奇序列）
2．计算序列循环卷积。
3．计算实序列离散傅里叶变换并检验DFT 性质。
4．实现补零序列的离散傅里叶变换。
6．实现高密度谱和高分辨率谱，并比较二者的不同。
三、 实验报告要求
见各程序要求
%以下为4个扩展函数
N -1
% （1）离散傅立叶变换 X (k ) =
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;
k =0∑x (n ) W N nk 采用矩阵相乘的方法
1%（2）逆离散傅立叶变换 x (n ) =N
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) 实序列的分解 N -1k =0∑X (k ) W N -nk
% 实序列可分解为共扼对称分量 xec =(1/2)*[x(n)+x((-n))N ]
% 和共扼反对称分量 xoc =(1/2)*[x(n)-x((-n))N ]
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 本例检验实序列的性质DFT[xec(n)]=Re[X(k)] DFT[xoc(n)]=Im[X(k)]
% 设 x(n)=10*(0.8).^n 0
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');
% 例2 本例为计算序列的循环卷积程序
% 运行之前应在命令窗口输入 x1,x2,N 的值
%实验报告要求：自己选择2个序列进行计算，将实验结果写出
if length(x1)>N
error('N must be >= the length of x1')
end
if length(x2)>N
error('N must be >= 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' %用矩阵相乘的方法得到结果
%例4 本例说明补零序列的离散傅立叶变换
%序列x (n ) =R 5(n ) ，已给出序列的傅立叶变换程序和将原序列补零到10长序列的DFT
%实验报告要求： （1）编写将序列补零到20长后计算DFT 的程序
(2)将实验结果画出
(3)实验结果说明什么（即序列补零后进行DFT, 频谱有何变化）
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(3,2,1);stem(n,x);ylabel('x(n)'); %画出矩形序列
subplot(3,2,2);plot(w/pi,Xe);ylabel('|X(ejw)|'); %画出离散时间傅立叶变换
N=10;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(3,2,3);stem(n,x);ylabel('x(n)'); %画出补零序列
subplot(3,2,4);stem(k,magX); %画出DFT 结果
axis([0,10,0,5]);ylabel('|X(k)|');
%例5 本题说明高密度谱和高分辨率谱之间的区别，高密度谱是信号补零后得到的，虽然谱线相当 %密但是因为信号有效长度不变，所以其分辨率也不变，因此还是很难看出信号的频谱成分。高分辨 %率谱是将信号有效长度加长，因此分辨率提高，可以看出信号的成分。
%有一个序列为x (n ) =2cos(0. 35πn ) +cos(0. 5πn ) （该序列周期计算可得40）
%（1）下面给出有10个有效采样点序列的DFT 程序
%（2）请写出将第一问中的10长序列补零到40长，计算其DFT
%（3）采样n=0：39，计算有40个有效采样点的序列的DFT
%实验报告要求： （1）请编写将有10个有效采样点的序列补零到40长后计算DFT 的程序
(2) 请编写计算有40个有效采样点的序列的DFT 的程序
(3) 将实验结果画出并分析实验结果说明什么
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(' 信号的频谱' );