数字信号处理实验三

实验三 离散付里叶变换（DFT）

一、实验目的：

1． 掌握离散付里叶级数 2． 掌握DFT变换。 3． 掌握DFT特性。

4． 掌握利用DFT计算线性卷积。 5． 掌握快速付里叶变换(FFT)。

二、实验原理：

1．离散付里叶级数（DFS） ~x(n)为周期序列，其频率为基本频率（2~数（DFS）为：；X(k)N1/N）的倍数（或谐波）。其离散付叶级

n0j~x(n)e2Nkn,k0,1,,

IDFS为：~x(n)1NN1k0~X(k)ej2Nkn

2．离散付里叶变换（DFT）

x(n)为长度N的有限长序列，其DFT为：x(n)1NN1N1k0X(k)ej2Nkn

IDFT为：X(k)n0x(n)ej2Nkn,k0,1,

3.DFT的特性：

（1） 线性性：DFT[ax1(n)bx2(n)]aDFT[x1(n)]bDFT[x2(n)]

x(0)n01nN1（2） 循环折叠（圆周对称）性：x((n))N**（3） 共轭性: DFT[x(n)]X((k))N

x(Nn)

*（4） 实序列的对称性（圆周共轭对称性）：X(k)X((k))N

x(nm)x((nm))N （5） 序列的圆周移位：~ 12

（6） 频域中的圆周移位：DFT[WNlnx(n)]X((kl))|NRN(k)

（7） 时域循环卷积：DFT[x1(n)x2(n)]X1(k)X2(k) （8） 频域循环卷积（乘法性）：DFT[x1(n)x2(n)]N121NX1(k)X2(k)

（9） 帕塞瓦尔（Parseval）定理：Exn0x(n)1NN1k0X(k)2

4.用DFT计算线性卷积:

设x1(n)为N1点序列，x2(n)为N2点序列，x3(n)为x1(n)和x2(n)的线性卷积，其为N1N21点序列，x4(n)为x1(n)和x2(n)的圆卷积，其长度为

N，当NN1N21时，x3(n)x4(n)。实际中，采用分段卷积法，即重

叠保留法和重叠相加法。需要对数据流进行分块处理，这时直接采用DFT计算线性卷积会产生一些问题，而应该将x(n)通过重复前M-1个取样进行分块，这样可得到正确结果。

5.快速付里叶变换（FFT）：

掌握基2-时域抽取FFT(DIT-FFT)和基2-频域抽取FFT(DIF-FFT)。MATLAB提供fft函数来计算x的DFT。fft函数是用机器语言写的，采用混合基法，其调用形式为：Xfft(x,N)。如N为2的幂，则得到高速的基2-FFT算法；若N不是2的乘方，则将N分解成质数，得到较慢的混合基FFT算法；最后，若N为质数，则fft函数采用的是原始的DFT算法。

三、实验步骤：

1.离散付里叶级数（DFS）

（1）自已动手:编写实现离散付里叶级数和逆离散付里叶级数的函数。

x(n){0,1,2,3,0,1,2,3,0,1,2,3}求其离散付里（2）已知周期性序列如下所示：~叶级数。

2.离散付里叶变换（DFT）

(1) 编写实现DFT和IDFT的函数。

(2) 已知x(n)是一个六点序列，如下所示：

1 x(n)00n5else

要求计算该序列的离散时间的付里叶变换和离散付里叶变换，并绘出它们的幅度和相

13

位。

(3) 序列后面增加零可以提高信号频谱的密度，比较高密度频谱与高分辨率频谱之间

的区别。考虑如下序列：x(n)cos(0.47ncos(0.53n 求其有限各样本的频谱，并要求：

（a） 当取 x(n)(0n10)时，确定出的离散付里叶变换。

（b） 将（a）的x(n)补零加长到0n100,确定出离散付里叶变换。 （c） 当取x(n)(0n500)，确定出的离散付里叶变换。

3.DFT性质：

（1） 设x(n)10(0.8)n,0n10，求y(n)x((n6))15 （2） 设x1(n){1,2,2},x2(n){1,2,5,4}试分别计算下列圆卷积。 （a）y(n)x1(n)x2(n),N5 (b) y(n)x1(n)x2(n),N6

4.利用DFT计算线性卷积：

用重复法求解：x(n)n1，0n9，h(n)1,0,1，求其线卷积y(n)。

5.快速付里叶变换（FFT）

（1）用MATLAB的fft()来求信号的DFT的幅值谱。

已知模拟信号为x(t)2sin(4t)5cos(8t)，以t=0.1n(n=0:N-1)进行取样，求N点DFT的幅值谱。

设N分别为（1）N=45；（2）N=50；（3）N=55；（4）N=60

（2在上题的基础上，N=64，并在信号中加入噪声（正态）w(t)(用函数randn(1,N)) x(t)2sin(4t)5cos(8t)0.8w(t)

试比较有无噪声时的信号谱。并分析在信号的检测的意义上，这种噪声会不会影响信号的检测？

四、思考题

1.通过分析计算机中的wav文件信号来进一步讨论数字信号处理中的信号分析方法。由于本题需要用户电脑具有麦克风和声卡，所以由同学们在课外完成。 （1） 先掌握几个要用到的函数：

[x,fs,bits]=waveread(‘filename’)

其功能是读取wav文件的函数；其中x表示一长串的数据，一般是两列（立体声）；fs是该wav文件在采集时用的采样频率；bits是指在进行A/D转化时用的量化位长（一般是8bits或16bits）. [d]=FFT(w,l)

是matlab中快速付里叶变换函数的一种输入输出形式。其中w是一列波形数据；

14

l是指示用多少点的FFT;d是频域输出. Sound(w,fs,bits)

与waveread的参数一样，它将数列的数据通过声卡转化为声音。

（2） 选择一个wav文件作为分析对象，可以用麦克风录一段用户自己的声音，也可

选择计算机中.wav的文件，在此选择Windows系统中都有的ding.wav(在

Windows\\media目录中).

（3） 编程求将此声音播放出并求出要处理的数据量的大小（答案为201911）。 （4） 分别画出单声道和双声道要处理的信号的波形图及信号频域的幅值（对信号求fft）。

（5） 求出此信号的主频率。

（6） 说明为何双声道能有一种从一侧耳朵穿过另一侧耳朵的感觉。

2.考虑一被噪声污染的信号，分析它所包含的频率分量，如一个由50Hz和120Hz正弦信号构成的信号，受零均值随机噪声的干扰，数据采样率为1000Hz.可通过fft函数来分析其信号频率成分。

3.混合信号成分分析。有一信号x由三种不同频率的正弦信号混合而成，通过得到信号的DFT，确定出信号的频率及其强度关系。

x(t)2sin(2*20*t)5cos(2*30*t)sin(2*pi*45*t)

4.信号在传输过程中，由于受信道或环境影响，在接收端得到的是噪声环境下的信号。我们利用FFT函数对这一信号进行傅里叶分析，从而确定信号的频率。 x(t)2sin(2*50*t)1.2randn(size(t));

5．天文学家记录了300年来太阳黑子的活动情况，我们对这组数据进行傅里叶分析，从而得出太阳黑子的活动周期。（其中记录的数据文件是sunspot.dat,在matlab中已有）.

15

实验三实验结果

1.离散付里叶级数（DFS）

（1） 自已动手:编写实现离散付里叶级数和逆离散付里叶级数的函数。 解：% dfs.m

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

WN=exp(-j*2*pi/N); nk=n'*k; WNnk=WN.^nk; Xk=xn*WNnk; % Idfs.m

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

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

WNnk=WN.^(-nk); xn=(Xk*WNnk)/N;

~x(n){0,1,2,3,0,1,2,3,0,1,2,3}求其离散付里（2） 已知周期性序列如下所示：

叶级数。 xn=[0 1 2 3];

N=4;

xk=dfs(xn,N)

xk = 6.0000 -2.0000 + 2.0000i -2.0000 - 0.0000i -2.0000 - 2.0000I

2.离散付里叶变换（DFT）

(4) 编写实现DFT和IDFT的函数。 解：%dft.m

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

k=0:N-1;

WN=exp(-j*2*pi/N);

16

nk=n'*k; WNnk=WN.^nk; Xk=xn*WNnk;

%idft.m

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

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

WNnk=WN.^(-nk); xn=Xk*WNnk;

(5) 已知x(n)是一个六点序列，如下所示：

1 x(n)00n5else

要求计算该序列的离散时间的付里叶变换和离散付里叶变换，并绘出它们的幅度和相位。 %计算DTFT n=0:5; x=ones(1,6); k=-200:200; w=(pi/100)*k;

X=x*(exp(-j*pi/100)).^(n'*k); magX=abs(X);

angX=angle(X)*180/pi; subplot(2,2,1) plot(w/pi,magX); axis([0 2 0 6]); grid

title('DTFT幅度'); subplot(2,2,2); plot(w/pi,angX);

axis([0 2 -200 200]); grid

xlabel('以pi为单位的频率'); title('DTFT的相位'); %计算DFT N=6;

X=dft(x,N); magX=abs(X)

phaX=angle(X)*180/pi n1=1:1024;

x1=[x,zeros(1,1018)];

17

N1=1024; X1=dft(x1,N1); magX1=abs(X1);

phaX1=angle(X1)*180/pi; subplot(2,2,3) stem(n,magX); hold on

plot(6*n1/1024,magX1,'r--'); title('DFT幅度'); grid hold off

subplot(2,2,4); stem(n,phaX); hold on

plot(6*n1/1024,phaX1,'r--'); title('DTFT相位'); grid hold off magX =

6.0000 0.0000 0.0000 0.0000 0.0000 0.0000 phaX =

0 -165.9638 -126.8699 -90.0000 -68.7495 -43.0908

18

(6) 序列后面增加零可以提高信号频谱的密度，比较高密度频谱与高分辨率频谱之间

的区别。考虑如下序列：x(n)cos(0.47ncos(0.53n

求其有限各样本的频谱，并要求：

（a） 当0n10，确定出的离散付里叶变换。(给序列进行补零加长至n=99) （b） 当0n500，确定出的离散付里叶变换。 解：（a）确定10点的DFT n=0:499;

x=cos(0.47*pi*n)+cos(0.53*pi*n); n1=0:9; y1=x(1:10); subplot(2,2,1); stem(n1,y1);

title('信号x(n) 0<=n<=9'); axis([0 10 -2.5 2.5]); text(10.2,-2.3,'n'); hold on;

plot([0 10],[0 0]); hold off; %进行DFT

Y1=dft(y1,10); magY1=abs(Y1(1:6)); k1=0:5;

w1=2*pi/10*k1; subplot(2,2,2); stem(w1/pi,magY1);

title('DTFT幅度的样本'); xlabel('以pi为单位的频率');

%将10点的DFT补零加长，看谱线是否变细。 n2=0:99;

y2=[x(1:10), zeros(1,90)]; subplot(2,2,3); stem(n2,y2);

19

title('信号x(n),0<=n<9+90 zeros'); %进行DFT

Y2=dft(y2,100); magY2=abs(Y2(1:51)); k2=0:50;

w2=2*pi/100*k2; subplot(2,2,4); stem(w2/pi,magY2);

title('DTFT幅度');

xlabel('以pi为单位的频率'); axis([0 1 0 9]);

从图看出此序列在w0.5处有一主频率，原始序列则没有这一频率，它只有主频为

w0.47和w0.53的两个频率，由于其点数不够，所以出现频谱混叠失真。给

序列尾部补零，提供了更加细化的谱线。 （b）确定x(n)前500点的DFT n=0:500-1;

x=cos(0.47*pi*n)+cos(0.53*pi*n); subplot(2,1,1); stem(n,x);

title('信号x(n) 0<=n<500'); %进行DFT

X=dft(x,500); magX=abs(X(1:251)); k=0:250; w=2*pi/100*k; subplot(2,1,2); plot(w/pi,magX);

title('DFT幅度');

xlabel('以pi为单位的频率');

20

axis([0 5 0 180]);

结果表明两个靠得很近的频率，这是高分辨率的频谱。从上面实验看出：

a. 给序列尾部添零,只能是得到更加光滑的频谱，但不能提供更多的信息。即得到高密度的频谱。

b.增加样本点则可以使频谱的信息量更大，即可得到高分辨频谱。

3.DFT性质：

（1） 设x(n)10(0.8)n,0n10，求y(n)x((n6))15

解：首先编循环移位函数：

%cirshftt.m

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);

其中m为序列x的长度,N为序列的循环移位的点数 主程序调用移位函数:

figure(1)

n=0:10;x=10*(0.8).^n;

y=cirshftt (x,6,15) %求循环移位 n=0:14;x=[x,zeros(1,4)];

subplot(2,1,1);stem(n,x);title( 'Original sequence' ) xlabel('n' );ylabel( 'x(n)');axis([-1,15,-1,11]) subplot(2,1,2);stem(n,y);

title('Circularly shifted sequence,N=15' ) xlabel('n' );ylabel( 'x((n-6) mod 15)' ); axis([-1,15,-1,11])

21

结果为：

y = 1.3422 1.0737 0 0 0 0 10.0000 8.0000 6.4000 5.1200 4.0960 3.2768 2.6214 2.0972 1.6777

（2） 设x1(n){1,2,2},x2(n){1,2,5,4}试分别计算下列圆卷积。 （a）y(n)x1(n)x2(n),N5 (b) y(n)x1(n)x2(n),N6 解：先编循环卷积函数： %circonvt.m

function y=circonvt(x1,x2,N) if length(x1)>N

error('N must be >=the length of x1') end

if length(x2)>N

error('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';

求x1(n)1,2,2，x1(n)1,2,5,4 （1）y1(n)x1(n)x2(n)

N=5 N=6

（2）y2(n)x1(n)程序：

x1=[1,2,2];

x2(n)

22

x2=[1,2,5,4];

y1=circonvt(x1,x2,5) y2=circonvt(x1,x2,6)

结果为：

y1=9 4 11 18 18

y2=1 4 11 18 18 8

4.利用DFT计算线性卷积：

用重复法求解：x(n)n1，0n9，h(n)1,0,1，求其线卷积y(n)。 解：先编一个重叠保留法作分段卷积函数； %ovrlpsav.m

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*L+N);

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

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

再编主程序进行计算； n=0:9; x=n+1;

h=[1,0,-1]; N=6;

y=ovrlpsav(x,h,N) 结果：

y = 1 2 2 2 2 2 2 2 2 2 -9 -10

5.快速付里叶变换（FFT）

（1）用MATLAB的fft()来求信号的DFT的幅值谱。

已知模拟信号为x(t)2sin(4t)5cos(8t)，以t=0.01n(n=0:N-1)进行取样，求N点DFT的幅值谱。

设N分别为（1）N=45；（2）N=50；（3）N=55；（4）N=60

23

subplot(2,2,1) N=45; n=0:N-1; t=0.01*n; q=n*2*pi/N;

x=2*sin(4*pi*t)+5*cos(8*pi*t); y=fft(x,N); plot(q,abs(y)) title('FFT N=45') %

subplot(2,2,2); N=50; n=0:N-1; t=0.01*n; q=n*2*pi/N;

x=2*sin(4*pi*t)+5*cos(8*pi*t); y=fft(x,N); plot(q,abs(y)) title('FFT N=50') %

subplot(2,2,3); N=55; n=0:N-1; t=0.01*n; q=n*2*pi/N;

x=2*sin(4*pi*t)+5*cos(8*pi*t); y=fft(x,N); plot(q,abs(y)) title('FFT N=55') %

subplot(2,2,4); N=60; n=0:N-1; t=0.01*n; q=n*2*pi/N;

x=2*sin(4*pi*t)+5*cos(8*pi*t); y=fft(x,N); plot(q,abs(y)) title('FFT N=60')

（2在上题的基础上，N=64，并在信号中加入噪声（正态）w(t)(用函数randn(1,N)) x(t)2sin(4t)5cos(8t)0.8w(t)

试比较有无噪声时的信号谱。并分析在信号的检测的意义上，这种噪声会不会影响信号的检测？

24

subplot(2,1,1) N=64; n=0:N-1; t=0.01*n; q=n*2*pi/N;

x=2*sin(4*pi*t)+5*cos(8*pi*t); y=fft(x,N); plot(q,abs(y)) title('FFT N=64') %

subplot(2,1,2); N=64; n=0:N-1; t=0.01*n; q=n*2*pi/N;

x=2*sin(4*pi*t)+5*cos(8*pi*t)+0.8*randn(1,N); y=fft(x,N); plot(q,abs(y))

title('FFT N=64(with noise)')

从下图中可以看出，在信号检测的意义上，这种噪声不会影响信号的检测。

25

四、思考题

1．通过分析计算机中的wav文件信号来进一步讨论数字信号处理中的信号分析方法。由于本题需要用户电脑具有麦克风和声卡，所以由同学们在课外完成。 （1）先掌握几个要用到的函数：

[x,fs,bits]=waveread(‘filename’)

其功能是读取wav文件的函数；其中x表示一长串的数据，一般是两列（立体声）；fs是该wav文件在采集时用的采样频率；bits是指在进行A/D转化时用的量化位长（一般是8bits或16bits）. [d]=FFT(w,l)

是matlab中快速付里叶变换函数的一种输入输出形式。其中w是一列波形数据；l是指示用多少点的FFT;d是频域输出. Sound(w,fs,bits)

与waveread的参数一样，它将数列的数据通过声卡转化为声音。

（2）选择一个wav文件作为分析对象，可以用麦克风录一段用户自己的声音，也可选择计算机中.wav的文件，在此选择Windows系统中都有的ding.wav(在盘符:\\Windows\\media目录中).

（3） 编程求将此声音播放出并求出要处理的数据量的大小（答案为20191*2）。 （4） 分别画出单声道和双声道要处理的信号的波形图及信号频域的幅值（对信号求

fft）。

（5） 求出此信号的主频率。

（6） 说明为何双声道能有一种从一侧耳朵穿过另一侧耳朵的感觉。

2.考虑一被噪声污染的信号，分析它所包含的频率分量，如一个由50Hz和120Hz正弦信号构成的信号，受零均值随机噪声的干扰，数据采样率为1000Hz.可通过fft函数来分析其信号频率成分。

2.混合信号成分分析。有一信号x由三种不同频率的正弦信号混合而成，通过得到信号的DFT，确定出信号的频率及其强度关系。 3.

26

4．模拟信号的取样和重构

设xa(t)e1000t，求其付里叶变换Xz(jw)

解：Xz(jw)Xa(t)ejwtdt

01000te0ejwtdte01000tejwtdt

0.0021(w1000)2

27

严格就来，在MATLAB中不使用SYMBOLIC工具箱是不能分析模拟信号的，但当以充分小的时间间隔取样xa(t)时，可产生平滑的图形，当包含足够的时间时，可显示出所有的模式，这样就可以近似地进行分析。

在此中设在105精度下，xa(t)为F0取

t51052000的带限信号，因此若

12200025105

这样就可得到xa(t)的逼近序列xG(m) 程序：

%Analog Signal 模拟信号的重构

Dt=0.00005;t=-0.005:Dt:0.005; xa=exp(-1000*abs(t));

%continuous-time Fourier transform 连续时间的付里叶变换 Wmax=2*pi*2000; % 最大的角频率

K=500;k=0:1:K;W=k*Wmax/K; %某点的角频率

Xa=xa*exp(j*t'*W)*Dt;%某点的频率,其中t’为t的转置矩阵 Xa=real(Xa);%幅频特性

W=[-fliplr(W),W(2:501)];%其中flipri按左右方向翻转矩阵元素,使得W以0为中心对称

Xa=[fliplr(Xa),Xa(2:501)];%使得Xa以最大的Xa为中心对称 figure(1)

subplot(2,1,1);plot(t*1000,xa,'-'); xlabel('t in msec.');ylabel('xa(t)')

gtext('Analog Signal') %用鼠标放置文本文件

subplot(2,1,2);plot(W/(2*pi*1000),Xa*1000,'r'); xlabel('Frequency in KHz');ylabel('Xa(jW)*1000') gtext('Continuous-time Fourier Transform')

自已动手：

以本题中的xa(t)说明取样速率对频域特性的影响：

28

（1）取样频率Fs=5000Hz,绘出X1(jw)曲线 （2）取样频率Fs=1000Hz,绘出X2(jw)曲线

29