数字信号处理实验二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 的分析方法。 二．实验内容

a. 设周期序列x(n)={ …,0,1,2,3,0,1,2,3,0,1,2,3,….}，求该序列的离散傅立叶级数X",(k) = DFS[x(n)]，并画出DFS 的幅度特性。 在matlab中新建函数dfs:

function [Xk]=dfs(xn,N) n=0:1:N-1; k=0:1:N-1;

Wn=exp(-j*2*pi/N); nk=n'*k; Wnk=Wn.^nk; Xk=xn*Wnk;

~~在matlab中输入以下代码：

xn=[0,1,2,3]; k=0:1:3; N=4;

－ 1 －

Xk=dfs(xn,N); y=abs(Xk); stem(k,y);

title('周期序列的离散傅立叶级数');

生成图像如下：

由定义可知，对于周期序列，根据离散傅里叶级数公式即可求出，实验中显示了一个周期的傅里叶级数。

b. 设周期方波序列为

x(n)=??1(mN?n?mN?L-1)(m=0,?1,?2,....)

?0(mN?L?n?(m?1)N-1)其中N 为基波周期，L/N 是占空比。 (1) 用L 和N求| X",(k) |的表达式；

(2) 当L 和N 分别为：L=5，N=20；L=5，N=40；L=5，N=60 以及

－ 2 －

L=7，N=60 时画出DFS 的幅度谱；

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

L=5; N=20;

k=[-N/2:N/2];

xn=[ones(1,L),zeros(1,N-L)]; Xk=dfs(xn,N);

y=abs([Xk(N/2+1:N) Xk(1:N/2+1)]); subplot(4,1,1); stem(k,y);

title('L=5,N=20'); L=5; N=40;

k=[-N/2:N/2];

xn=[ones(1,L),zeros(1,N-L)]; Xk=dfs(xn,N);

y=abs([Xk(N/2+1:N) Xk(1:N/2+1)]); subplot(4,1,2); stem(k,y);

title('L=5,N=40'); L=5; N=60;

k=[-N/2:N/2];

xn=[ones(1,L),zeros(1,N-L)]; Xk=dfs(xn,N);

y=abs([Xk(N/2+1:N) Xk(1:N/2+1)]); subplot(4,1,3); stem(k,y);

title('L=5,N=60'); L=7; N=60;

k=[-N/2:N/2];

xn=[ones(1,L),zeros(1,N-L)]; Xk=dfs(xn,N);

y=abs([Xk(N/2+1:N) Xk(1:N/2+1)]); subplot(4,1,4); stem(k,y);

－ 3 －

title('L=7,N=60');

生成图像如下：

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

c. 设有限长序列x(n) = {0,1,2,3}，计算DTFT[x(n)]=X(ejω)，并画出它

2?的幅度谱；然后利用kw1=4~X它等于实验a 中的(k)。

k,k=0,1,2,3对X(ejω)进行采样，并证明

X(ejω)=n????x(n)e??jwn=e-jw+2e-j2w+3e-j3w

X(ej0)=1+2+3=6;

－ 4 －

X(ej2?/4)=-2+2j; X(ej4?/4)=2; X(ej6?/4)=-2-2j; 实验代码：

N=4;

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

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

Xw=exp(-j*w)+2*exp(-j*2*w)+3*exp(-j*3*w); subplot(211); plot(w/pi,abs(Xw)); title('X(ejw)幅度谱'); Xk=fft(xn,N); subplot(212);

stem(k,abs(Xk),'fill'); hold on;

plot(N/2*w/pi,abs(Xw));

生成图像如下：

－ 5 －

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

d. 序列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，并绘制幅度与相位谱，求频率分辨率； 在matlab中新建函数dft:

－ 6 －

function [Xk]=dft(xn,N) n=0:1:N-1; k=0:1:N-1;

Wn=exp(-j*2*pi/N); nk=n'*k; Wnk=Wn.^nk; Xk=xn*Wnk;

function [ xn ] = idft( Xk,N ) n=0:1:N-1; k=0:1:N-1;

Wn=exp(-j*2*pi/N); nk=n'*k;

Wnk=Wn.^(-nk); Xk=(Xk*Wnk)/N; End

（1）在matlab中输入以下序列：

x=[1 1 1 1]; N=4;

X=dft(x,N); subplot(2,1,1); stem(abs(X));

title('x(n)4点DFT对应的幅度谱'); subplot(2,1,2);

stem(angle(X)*180/pi);

title('x(n)4点DFT对应的相位谱');

生成图像如下：

－ 7 －

（2）在matlab中输入以下序列：

x=[1 1 1 1];

x=[x,zeros(1,4)]; N=8;

X=dft(x,N); subplot(2,1,1); stem(abs(X));

title('x(n)补零成8点DFT对应的幅度谱'); subplot(2,1,2);

stem(angle(X)*180/pi);

title('x(n)补零成8点DFT对应的相位谱');

生成图像如下：

－ 8 －

（3）在matlab中输入以下序列：

x=[1 1 1 1];

x=[x,zeros(1,12)]; N=16;

X=dft(x,N); subplot(2,1,1); stem(abs(X));

title('x(n)补零成16点DFT对应的幅度谱'); subplot(2,1,2);

stem(angle(X)*180/pi);

title('x(n)补零成16点DFT对应的相位谱');

－ 9 －

生成图像如下：

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

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

n=0:1:99;

x=cos(0.48*pi*n)+cos(0.52*pi*n); n1=0:1:9; y1=x(1:1:10);

－ 10 －

y2=dfs(y1,10); subplot(2,1,1) stem(n1,abs(y2));

title('x(n)10点DFT对应的幅度谱'); subplot(2,1,2);

stem(n1,angle(y2));

title('x(n)10点DFT对应的相位谱');

生成图像如下：

（2）实验代码：

n=0:1:99;

x=cos(0.48*pi*n)+cos(0.52*pi*n); y1=x(1:1:100); y2=dfs(y1,100); subplot(2,1,1) stem(n,abs(y2));

title('x(n)100点DFT对应的幅度谱'); subplot(2,1,2); stem(n,angle(y2));

－ 11 －

title('x(n)100点DFT对应的相位谱');

生成图像如下：

f. 序列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 的波形，验证序列圆周移位性质。

在matlab中新建函数circevod和函数cirshftt: 函数circevod：

－ 12 －

function [xec,xoc]=circevod(x) if any(imag(x)~=0)

error('x is not a real sequence') end

N=length(x); n=0:(N-1);

xec=0.5*(x+x(mod(-n,N)+1)); xoc=0.5*(x-x(mod(-n,N)+1));

函数cirshftt:

function y=cirshftt(x,m,N) if length(x)>N

error('N must be >= 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）实验代码：

n=0:10;

x=5*(0.9).^n;

y=x(mod(-n,11)+1); figure(1);

subplot(2,1,1); stem(n,x);

title('x(n)时域波形'); subplot(2,1,2); stem(n,y);

title('循环反转序列x((-n)11时域波形');

X=dft(x,11); Y=dft(y,11); figure(2);

subplot(2,2,1); stem(n,real(X));

title('x(n)的DFT实部'); subplot(2,2,2); stem(n,imag(X));

title('x(n)的DFT虚部'); subplot(2,2,3); stem(n,real(Y));

title('循环反转序列x((-n)11的DFT实部');

－ 13 －

subplot(2,2,4); stem(n,imag(Y));

title('循环反转序列x((-n)11的DFT虚部');

生成图像如下：

－ 14 －

（2）实验代码：

n=0:10;

x=5*(0.9).^n;

[xec,xoc]=circevod(x); figure(1);

subplot(2,1,1); stem(n,xec);

title('圆周共轭奇分量'); subplot(2,1,2); stem(n,xoc);

title('圆周共轭偶分量');

X=dft(x,11);

[xec,xoc]=circevod(x); Xec=dft(xec,11); Xoc=dft(xoc,11); figure(2);

subplot(2,2,1); stem(n,real(X));

－ 15 －

title('x(n)的DFT实部'); subplot(2,2,2); stem(n,imag(X));

title('x(n)的DFT虚部'); subplot(2,2,3); stem(n,real(Xec));

title('圆周共轭奇分量的DFT实部'); subplot(2,2,4);

stem(n,imag(Xoc));

title('圆周共轭偶分量的DFT虚部');

生成图像如下：

－ 16 －

（3）实验代码：

n=0:10; n1=0:14; x=5*(0.9).^n; x1=[x,zeros(1,3)]; y1=cirshftt(x,-4,11); y2=cirshftt(x,3,15); y3=cirshftt(x,6,15); subplot(2,2,1); stem(n,x);

title('x(n)序列'); subplot(2,2,2); stem(n,y1);

title('x((n+4))11R11(n)序列'); subplot(2,2,3); stem(n1,y2);

title('x((n-3))15R15(n)序列'); subplot(2,2,4); stem(n1,y3);

－ 17 －

title('x((n-6))15R15(n)序列');

生成图像如下：

由图可以看出，xe具有循环对称性，xo具有循环反对称性质。 g. 序列x1(n)={1,2,1}，x2(n)={1,2,3,2}，

(1) 编制程序在时域中直接计算4点循环卷积1 4 2y(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 6x (n)?x (n)，分析循环卷积点数N 对循环卷积的影响，并比较哪种循环卷积结果与线性卷积是一致的？总结出用循环卷积计算线性卷积的条件。

－ 18 －

在matlab中新建函数circonvl：

function y=circonvl(x1,x2,N) 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))]; 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

H(n,:)=cirshftt(x2,n-1,N); end

y=x1*H';

（1）实验代码：

x1=[1,2,1]; x2=[1,2,3,2];

y=circonvl(x1,x2,4); stem(y);

title('4点循环卷积波形图');

生成图像如下：

－ 19 －

（2）实验代码：

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

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

y1=circonvl(x1,x2,4); x11=[1 2 1 0]; X1=fft(x11,4); X2=fft(x2,4); Y=X1.*X2; y2=ifft(Y,4); subplot(211); stem(n2,y1,'.'); grid on;

title('通过时域卷积求循环卷积') subplot(212); stem(n2,y2,'.'); grid on;

title('通过反变换方法求循环卷积')

－ 20 －

生成图像如下：

（3）实验代码：

x1=[1,2,1]; x2=[1,2,3,2]; subplot(2,1,1)

y=circonvl(x1,x2,5); stem(y);

subplot(2,1,2)

y=circonvl(x1,x2,6); stem(y);

生成图像如下：

－ 21 －

h. 序列x(n)=(n+1)R10(n)，h(n)={1,0,-1}，利用重叠保留法，编制程序用N=6 点的循环卷积计算线性卷积的值y(n)=x(n)*h(n)，并与直接线性卷积结果进行比较。 在matlab中新建函数ovrlpsav：

function [y]=ovrlpsav(x,h,N) Lenx=length(x); M=length(h); M1=M-1; L=N-M1;

h=[h zeros(1,N-M)];

x=[zeros(1,M1),x,zeros(1,N-1)]; K=floor((Lenx+M1-1)/(L)); Y=zeros(K+1,N); for k=0:K

xk=x(k*L+1:k*K+N);

Y(k+1,:)=circonvl(xk,h,N);

－ 22 －

end

Y=Y(:,M:N)'; y=(Y(:))';

实验代码：

n=0:9; N=6; x=n+1;

h=[1,0,-1];

y=ovrlpsav(x,h,N); stem(y);

生成图像如下：

i. 序列x(n)=R6(n)，用快速傅立叶变换FFT 计算6 点DFT[x(n)]和8 点DFT[x(n)]，绘制波形图并比较结果。 实验代码：

－ 23 －

x1=ones(1,6); x2=[x1 0 0]; y1=fft(x1,6); y2=fft(x2,8); subplot(2,1,1); stem(abs(y1)); subplot(2,1,2); stem(abs(y2));

生成图像如下：

j. 设两个序列x(n)=(2n+3)R17(n)，h(n)=(n+1)R4(n)，采用重叠相加法，按分段长度L=7的FFT 计算线性卷积：y(n)=x(n)*h(n)，并与直接线性卷积的结果进行比较。 实验代码：

n1=0:16; n2=0:3;

－ 24 －

x=2*n1+3; h=n2+1;

y=conv(x,h); stem(y);

生成图像如下：

k. 设序列x(n)=(0.8)nR15(n)，计算序列x(n)在单位园上的Chirp-z 变换，并与DFT[x(n)]的结果进行比较。 实验代码：

n=0:14;

x=((0.8).^n).*stepfun(n,0);

w=exp(-j*2*pi/15); %z = a*(w.^-(0:m-1))，所以w=exp(-j*2*pi/15)

y1=czt(x,15,w,1); %用函数czt(x,15,w,1)求Chirp-z变换，15是x的长度 y2=fft(x,15); %求FFT subplot(2,1,1); stem(n,abs(y1));

axis([0,14,0,6]); %定义横轴和纵轴的范围

－ 25 －

title('Chirp-z变换的图形'); subplot(2,1,2); stem(n,abs(y2));

title('DFT变换的图形'); axis([0,14,0,6]);

生成图像如下：

l. 设信号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)]；

(3) 在5—12Hz 的频段范围求Chirp-z 变换。

－ 26 －

实验代码：

N=400; n=0:399; fs=40;

w=exp(-j*2*pi/400); %z = a*(w.^-(0:m-1))，所以w=exp(-j*2*pi/400) xn1=sin(6*n*fs)+sin(6.5*n*fs)+sin(8*n*fs)+sin(9*n*fs)+sin(10*n*fs); y1=czt(xn1,400,w,1); %用函数czt(x,400,w,1)求Chirp-z变换，400是x的长度 y2=fft(xn1,400); %直接求FFT subplot(2,1,1); stem(n,abs(y1));

title('Chirp-z变换的图形'); subplot(2,1,2); stem(n,abs(y2));

title('DFT变换的图形');

生成图像如下：

三．实验心得

通过这次实验，进一步加强了对于课本是上理论知识的理解，熟悉了

－ 27 －

利用MATLAB 进行序列的DFT、FFT 的分析方法。此外，对周期序列DFS、有限长序列DFT 的概念、旋转因子的定义、以及DFS 和DFT的性质有了更深的理解。实验过程中遇到了一些以前没有接触过的函数，通过查阅资料最终掌握了它们的使用方法。

－ 28 －

本文来源：https://www.bwwdw.com/article/m6kr.html