实验二 DFT 和FFT

实验二 DFT 和FFT

一、实验目的

1）认真复习周期序列 DFS、有限长序列DFT 的概念、旋转因子的定义、以及DFS 和DFT的性质等有关内容；复习基2-FFT 的基本算法，混合基-FFT 的基本算法、Chirp-Z 变换的算法等快速傅立叶变换的方法。

2）掌握有限长序列的循环移位、循环卷积的方法，对序列共轭对称性的含义和相关内容加深理解和掌握，掌握利用DFT 分析序列的频谱特性的基本方法。

3）掌握 FFT 算法的基本原理和方法、Chirp-Z 变换的基本原理和方法，掌握利用FFT 分析序列的频谱特性的方法。

4）熟悉利用 MATLAB 进行序列的DFT、FFT 的分析方法。 二、实验内容 1）设周期序列( ) { xn =……,0,1,2,3,0,1,2,3,0,1,2,3,……}，求该序列的离散傅立叶级数X (k) = DFS[x(n)]，并画出DFS 的幅度特性。

主程序：

clc; N=4; n=0:N-1; k=0:N-1; x=[0 1 2 3];

xk=x*exp(-j*2*pi/N).^(n'*k); stem(k,abs(xk)); xlabel('k');

title('xk的幅度性')

分析：由定义可知，对于周期序列，根据离散傅里叶级数公式即可求出，此实验中显示了一个周期的傅里叶级数。 2）设周期方波序列为

其中 N 为基波周期，L/N 是占空比。 (1) 用L 和N求| X(k) |的表达式；

信息科学与工程学院 数字信号处理实验报告

(2) 当L 和N 分别为：L=5，N=20；L=5，N=40；L=5，N=60 以及L=7，N=60 时画出DFS 的幅度谱；

(3) 对以上结果进行讨论，总结其特点和规律。

主程序：1

clc;

%fudupu(5,20) L=5;N=20; n=0:4*N-1; k=0:4*N-1;

x=(ones(1,l),zroes(1,N-L)); x=[x x x x ];

xk=x*exp(-j*2*pi/N).^(n'*k); stem(k,abs(xk)); xlabel('k');

title('xk的幅度性')

结果：（修改代码中的L和N（x(n)），可以得到其他占空比时DFS的幅度谱）

2

clc;

xk1=fudupu(5,20); figure

xk2=fudupu(5,40); figure;

xk3=fudupu(5,60); figure

xk4=fudupu(7,60) 结果：

信息科学与工程学院 数字信号处理实验报告

求幅度谱的函数

function xk=fudupu(L,N) n=0:4*N-1; k=0:4*N-1;

x=[ones(1,L),zeros(1,N-L)]; x=[x x x x ];

xk=x*exp(-j*2*pi/N).^(n'*k); stem(k,abs(xk)); xlabel('k');

title('xk的幅度性')

分析：由四组图对比可知,N越大，其频域抽样间隔越小，N为频域的重复周期。占空比L/N主要决定第一零点带宽（在一个周期内）。

3）设有限长序列x(n) = {0,1,2,3}，计算DTFT[x(n)]=X(ejω)，并画出它的幅度谱；然后利用kω=2*pi*k/4,k=0,1,2,3对X(ejω)进行采样，并证明它等于实验a 中的X(k)。

主程序：

clc; N=4;

x=[0 1 2 3]; n=0:N-1; k=0:N-1; m=0:2047;

w=2*pi*m/2048;

xw=exp(-j*w)+exp(-j*2*w)+exp(-j*3*w); subplot(3,1,1);

plot(w/pi*2,abs(xw)); axis([0,4, 0, 6]) xk=fft(x,N); subplot(3,1,2); stem(k,abs(xk)); hold on ;

plot(w/pi*2,abs(xw));

信息科学与工程学院 数字信号处理实验报告

N=4; n=0:N-1; k=0:N-1; x=[0 1 2 3];

xk=x*exp(-j*2*pi/N).^(n'*k);

subplot(3,1,3); hold on ;

stem(k,abs(xk));axis([0,4, 0, 6]) xlabel('k');

title('xk的幅度性')

分析：对比第一题的结果可以看出，对离散傅里叶变换的频谱进行抽样，在满足采样定理的条件下，可以无失真的恢复原来的波形。

4）序列x(n)=R4(n)，计算DTFT[x(n)]=X(ejω)，并绘制其幅度和相位谱。 (1) 计算x(n)的4 点DFT，并绘制DFT 的幅度与相位谱；

(2) 将x(n)补零形成8 点序列，计算8 点DFT，并绘制幅度与相位谱，求频率分辨率； (3) 将x(n)补零形成16 点序列，计算16 点DFT，并绘制幅度与相位谱，求频率分辨率；

主程序：

N1=8; N2=16; n=0:N1-5; k1=0:N1-1; k2=0:N2-1;

w=2*pi*(0:2047)/2048;

Xw=(1-exp(-j*4*w))./(1-exp(-j*w));

%对x(n)的频谱函数采样2048个点可以看做连续的频谱 xn=[(n>=0)&(n<4)]; X1k=fft(xn,N1); X2k=fft(xn,N2); subplot(321);

plot(w/pi,abs(Xw)); title('连续函数幅度谱'); xlabel('w/pi');

信息科学与工程学院 数字信号处理实验报告

subplot(322);

plot(w/pi,angle(Xw)); title('连续函数相位谱'); axis([0,2,-pi,pi]); line([0,2],[0,0]); xlabel('w/pi'); subplot(323);

stem(k1,abs(X1k),'.'); title('X1(k)幅度谱'); hold on;

plot(N1/2*w/pi,abs(Xw));%图形上叠加连续谱的幅度曲线 subplot(324)

stem(k1,angle(X1k),'.'); axis([0,N1,-pi,pi]); line([0,N1],[0,0]); title('X1(k)相位谱'); hold on;

plot(N1/2*w/pi,angle(Xw));%图形上叠加连续谱的相位曲线 subplot(325);

stem(k2,abs(X2k),'.'); title('X2(k)幅度谱'); hold on;

plot(N2/2*w/pi,abs(Xw)); subplot(326)

stem(k2,angle(X2k),'.'); axis([0,N2,-pi,pi]); line([0,N2],[0,0]); title('X2(k)相位谱'); hold on;

plot(N2/2*w/pi,angle(Xw));

实验结果：

信息科学与工程学院 数字信号处理实验报告

分析：对序列x(n)的N点DFT的物理意义是对其傅里叶变换的频谱在[0,2*pi]上进行N点等间隔的采样。由实验结果可以看出。 5）序列x(n)=cos(0.48πn)+cos(0.52πn)

(1) 求x(n)的10 点DFT，并画出它幅度与相位谱； (2) 求x(n)的100 点DFT，并画出它幅度与相位谱； 根据实验结果，讨论 DFT 进行谱分析的条件。 主程序：1 clc;

x=cos(0.48*pi*n)+cos(0.52*pi*n) N=10;

xk=fdxw(x,N) 2 clc;

xn=cos(0.48*pi*n)+cos(0.52*pi*n) y=fdxw(xn,100) 幅度和相位函数 ： function xk=fdxw(x,N) k=0:N-1;

w=2*pi*(0:2047)/2048; xk=fft(x,N); subplot(211); stem(k,abs(xk),'.'); title('X(k)幅度谱'); subplot(212);

stem(k,angle(xk),'.'); title('X(k)相位谱');

实验结果：

信息科学与工程学院 数字信号处理实验报告

6）序列x(n)=5(0.9)nR11(n)

(1) 求循环反转序列x((-n))11，并绘制x(n)和x((-n))11 的波形；求两序列的DFT，验证DFT 的循环反转性质。

(2) 把序列x(n)分解成圆周共轭奇分量xoc(n)和圆周共轭偶分量xec(n)，并求出对应的DFT，验证DFT 的圆周共轭对称性质。

(3) 绘制x((n+4))11R11(n)、x((n-3))15R15(n)和x((n-6))15 的波形，验证序列圆周移位性质。 主程序：1 clc; N=11; n=0:N-1; r=ones(1,11); x=5*((0.9).^n).*r; subplot(321); stem(n,x,'.'); title('原始序列');

x1=[x,zeros(1,N-length(x))]; y=x1(mod(-n,N)+1); subplot(322); stem(n,y,'.');

title('反转序列'); figure

y1=fdxw(x,11); figure

x1=[x,zeros(1,N-length(x))]; y2=fdxw(y,11);

信息科学与工程学院 数字信号处理实验报告

2 clc; N=11; n=0:N-1; r=ones(1,11); x=5*((0.9).^n).*r;

x1=[x,zeros(1,N-length(x))]; y=x1(mod(-n,N)+1); xe=0.5*(x+y); xo=0.5*(x-y); Xe=fft(xe,N); Xo=fft(xo,N); subplot(321); X=fft(x);

stem(n,abs(X),'.');

title('原始序列DFT幅度谱'); subplot(322);

stem(n,abs(Xo),'.');

title('奇对称序列DFT幅度谱'); subplot(323);

stem(n,abs(Xe),'.');

title('偶对称序列DFT幅度谱'); subplot(324);

stem(n,angle(X),'.');

title('原始序列DFT相位谱'); subplot(325);

stem(n,angle(Xo),'.');

title('奇对称序列DFT相位谱'); subplot(326);

stem(n,angle(Xe),'.');

title('偶对称序列DFT相位谱');

信息科学与工程学院 数字信号处理实验报告

3 clc; N=11; n=0:N-1; r=ones(1,11); x=5*((0.9).^n).*r; n4=0:10; n5=0:14; n6=0:14;

x4=cirshftt(x,4,11); x5=cirshftt(x,-3,15); x6=cirshftt(x,-6,15); subplot(311); stem(n4,x4,'.');

title('X((n+4))U(11)'); subplot(312); stem(n5,x5,'.');

title('X((n-3))U(15)'); subplot(313); stem(n4,x4,'.');

title('X((n-6))U(15)');

循环移位的函数

function y=cirshftt(x,m,N)%求循环移位X((n+m))N R(N) if length(x)>N

error('Nx=[x,zeros(1,N-length(x))]; n=0:N-1;

y=x(mod(n-m,N)+1);

分析：由实验结果可以知道，实域循环反褶，频域对应共轭。

信息科学与工程学院 数字信号处理实验报告

由图可以看出，xe具有循环对称性，xo具有循环反对称性质。由于xe虚部为0，对应其相位为0.

7）序列x1(n)={1,2,1}，x2(n)={1,2,3,2}，

(1) 编制程序在时域中直接计算4 点循环卷积1 4 2，y(n) = x (n)⊗x (n)，并绘制波形图； (2) 编制程序先计算x1(n)和x2(n)的4 点DFT，再利用IDFT 计算1 2 4 x (n)⊗x (n)，比较结果以上两种计算结果是否一致？

(3) 计算1 2 5x (n)⊗x (n)和1 2 6，x (n)⊗x (n)，分析循环卷积点数N 对循环卷积的影响，并比较哪种循环卷积结果与线性卷积是一致的？总结出用循环卷积计算线性卷积的条件。 主程序：

（1） 计算x1和x2四点循环卷积 1 clc;

x1=[1 2 1]; n1=0:2; x2=[1 2 3 2]; n2=0:3;

y1=circonvt(x1,x2,4); subplot(311) stem(n1,x1,'.'); grid on; title('x1'); subplot(312) stem(n2,x2,'.'); grid on; title('x2'); subplot(313) stem(n2,y1,'.'); grid on; title('y1'); 结果：

信息科学与工程学院 数字信号处理实验报告

2 clc;

x1=[1 2 1 0]; x2=[1 2 3 2]; X1=fft(x1,4); X2=fft(x2,4); Y=X2.*X1; stem(Y); y=ifft(Y) subplot(211) stem(y);

title('乘积法球卷积'); n2=0:3;

y1=circonvt(x1,x2,4); subplot(212) stem(n2,y1,'.') 结果：

信息科学与工程学院 数字信号处理实验报告

两者点数样，但是有所移位 3 clc;

x1=[1 2 1 0 0]; n1=0:4;

x2=[1 2 3 2 0]; n2=0:4;

y1=circonvt(x1,x2,5); subplot(211) stem(n2,y1,'.'); axis([0,5,0,15]); title('N=5'); x1=[1 2 1 0 0 0]; n1=0:5;

x2=[1 2 3 2 0 0]; n2=0:5;

y2=circonvt(x1,x2,6); subplot(212) stem(n2,y2,'.'); title('N=6'); 结果：

信息科学与工程学院 数字信号处理实验报告

可以看出6点循环卷积和线性卷积的结果一致，可以看出当圆周卷积点数L>=n1+n2-1时，此时圆周卷积可以代表线性卷积。

循环卷积函数

function fn=circonvt(x1,x2,N)

%fn=sum(x1(m)*x2((n-m) mod N)) if (length(x1)>N|length(x2)>N)

error('N的长度必须大于输入数据的长度'); end

x1=[x1,zeros(1,N-length(x1))]; x2=[x2,zeros(1,N-length(x2))]; m=0:N-1; x=zeros(N,N); for n=0:N-1

x(:,n+1)=x2(mod((n-m),N)+1)'; end;

fn=x1*x; %循环计算卷积

8） 序列x(n)=(n+1)R10(n)，h(n)={1,0,-1}，利用重叠保留法，编制程序用N=6 点的循环卷积计算线性卷积的值y(n)=x(n)*h(n)，并与直接线性卷积结果进行比较。 主程序： clc; n=0:9;

r=ones(1,10); x1=(n+1).*r; h=[1 0 -1]; f1=conv(x1,h); subplot(211); stem(f1);

title('线形卷积'); N=6;

f2=fdconv(x1,h,N);

信息科学与工程学院 数字信号处理实验报告

subplot(212) stem(f2);

title('重叠法求卷积'); 函数

function [y]=fdconv(x,h,N) Lx=length(x); M=length(h);

M1=M-1;L=N-M1;%各段搭接长度M1，有效程度为L h=[h,zeros(1,N-M)];%将h延长至循环长度N; x=[x,zeros(1,N+M1-1)];%在x前面加上M-1个0 K=floor((Lx+M1-1)/(L))+1;%段数 Y=zeros(K+1,N); for k=0:K-1;

xk=x(k*L+1:k*L+N); Y(k+1,:)=circonvt(xk,h,N); end;

Y=Y(:,M:N)';%去掉前面M-1个样本 y=(Y(:))';%将列向量转换成行向量输出

对比两个结果可以看出，自己编写的函数长度为14，产生误差的是最后4个样本，这是由于自己编写的函数输入的长度不合适，在其后面随便补了几个0，从而产生几个样本的误差。 9）序列x(n)=R6(n)，用快速傅立叶变换FFT 计算6 点DFT[x(n)]和8 点DFT[x(n)]，绘制波形图并比较结果。 主程序： clc;

x=ones(1,6); y1=fft(x,6); y2=fft(x,8); subplot(211) stem(y1,'.','r'); title('6点DFT') subplot(212) stem(y2,'.','r'); title('8点DFT')

信息科学与工程学院 数字信号处理实验报告

实验结果：

10） 设两个序列x(n)=(2n+3)R17(n)，h(n)=(n+1)R4(n)，采用重叠相加法，按分段长度L=7的FFT 计算线性卷积：y(n)=x(n)*h(n)，并与直接线性卷积的结果进行比较。 主程序： clc; n=0:16

%r=ones(1,17); x=2*n+3; r1=ones(1,4); h=[1 2 3 4]; y1=conv(x,h); subplot(211); stem(y1,'.','r'); title('线性卷积') y2=fdconv1(x,h,7); subplot(212) stem(y2,'.','r')

title('重叠相加法') 结果：

分析：fdconv函数源代码与第8题一致，对比两个结果可以看出，自己编写的函数长度为24，

信息科学与工程学院 数字信号处理实验报告

产生误差的是最后4个样本，这是由于自己编写的函数输入的长度不合适，在其后面随便补了几个0，从而产生几个样本的误差。

11）设序列x(n)=(0.8)nR15(n)，计算序列x(n)在单位园上的Chirp-z 变换，并与DFT[x(n)]的结果进行比较。 主程序： clc; n=0:14;

x=(0.8.^n).*(ones(1,15)); y1=fft(x,15); subplot(211); stem(y1,'.'); title('DFT');

w=exp(-1*j*2*pi/15); y2=czt(x,15,w,1); subplot(212); stem(y2,'.'); title('chirp_z') 实验结果：

分析：由实验结果可以看出，两种方法求的结果一致，在单位圆上求Z变换可以直接调用ctz函数。

12）设信号x(t)=sin(f1t)+sin(f2t)+sin(f3t)+sin(f4t)+sin(f5t)，其中f1=6Hz，f2=6.5Hz，f3=8Hz和f4=9Hz，f5=10Hz，对信号x(t)进行频率为40Hz 进行抽样，时域抽样400 点。 (1) 用Chirp-z 变换计算DFT[x(n)]； (2) 直接计算DFT[x(n)]； 主程序： clc;

n=0:400; t=1/40; nt=t*n;

f1=6;f2=6.5;f3=8;f4=9;f5=10;

信息科学与工程学院 数字信号处理实验报告

x=sin(f1*nt)+sin(f1*nt)+sin(f1*nt)+sin(f1*nt)+sin(f1*nt); subplot(311) stem(x); subplot(312) Xk=czt(x); Xk1=fft(x); %Xk=czt(x); 结果：

分析：对时域进行采样，相当于将t用nfs代入，由结果可以看出两种方法求得的结果一致。 三、实验总结

本次实验相对于上次实验难度有了提高，主要理解了求ＤＦＴ的方法，以及在处理循环序列方面常用的几种情况比如求循环卷积，周期延拓等，验证ＤＦＴ的性质。实验中加深了对课本知识的理解，起到相辅相成的作用，实验效果比较理想。实验过程中遇到一些小问题，也及时得到了解决。