数字信号处理实验三

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

一、实验目的：

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

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

二、实验原理：

1．离散付里叶级数（DFS ） )(~n x 为周期序列，其频率为基本频率（N /2π

）的倍数（或谐波）。其离散付叶级

数（DFS ）为：；∑

-=-±==

1

2,,1,0,)(~)(~N n kn

N

j

k e

n x k X π

IDFS 为：∑

-==

1

2)(~

1)(~N k kn

N

j

e

k X N

n x π

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

)(n x 为长度N 的有限长序列，其DFT 为：21

1()()N j

kn

N

k x n X k e

N

π--==

∑

IDFT 为：21

()(),0,1,N j

kn

N

n X k x n e

k π-==

=±∑

3.DFT 的特性：

（1） 线性性：)]([)]([)]()([2121n x bDFT n x aDFT n bx n ax DFT +=+

（2） 循环折叠（圆周对称）性：??

?

-≤≤-==-1

1)

(0)

0())((N n n N x n x n x N

（3） 共轭性: N k X n x DFT ))(()]([*

*-=

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

-=

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

-=-

（6） 频域中的圆周移位：)())(()]([|ln

k R l k X n x W DFT N N N

-=-

（7） 时域循环卷积：)()()]()([2121k X k X n x n x DFT =? （8） 频域循环卷积（乘法性）：)()(1)]()([2121k X k X N

n x n x DFT ?=

（9） 帕塞瓦尔（Parseval ）定理：∑

∑

-=-==

=

1

2

1

2

)

(1)(N k N n x k X N

n x E

4.用DFT 计算线性卷积:

设)(1n x 为1N 点序列，)(2n x 为2N 点序列，)(3n x 为)(1n x 和)(2n x 的线性卷积，其为121-+N N 点序列，)(4n x 为)(1n x 和)(2n x 的圆卷积，其长度为

N ，当121-+=N N N 时，)()(43n x n x =。实际中，采用分段卷积法，即重

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

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

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

三、实验步骤：

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

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

（2）已知周期性序列如下所示：}3,2,1,0,3,2,1,0,3,2,1,0{)(~

↑

=n x 求其离散付里叶级数。

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

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

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

??

?≤≤=e ls e

n n x 0

501

)(

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

位。

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

的区别。考虑如下序列：()n n n x ππ53.0cos(47.0cos()(+= 求其有限各样本的频谱，并要求：

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

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

3.DFT 性质：

（1） 设n n x )8.0(10)(=,100≤≤n ，求15))6(()(-=n x n y （2） 设}4,5,2,1{)(},2,2,1{)(21==n x n x 试分别计算下列圆卷积。 （a ）5),()()(21=?=N n x n x n y (b) 6),()()(21=?=N n x n x n y

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

用重复法求解：1)(+=n n x ，90≤≤n ，{}1,0,1)(-=n h ，求其线卷积)(n y 。

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

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

已知模拟信号为)8cos(5)4sin(2)(t t t x ππ+=，以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)) )(8.0)8cos(5)4sin(2)(t w t t t x ++=ππ

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

四、思考题

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

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

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

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

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

π

2

)

*

2

20

*

(t

pi

+

=π

x+

t

t

t

sin(

*

2

*

*

)

)

sin(

45

2

cos(

5

*

*

30

)

4.信号在传输过程中，由于受信道或环境影响，在接收端得到的是噪声环境下的信号。我们利用FFT函数对这一信号进行傅里叶分析，从而确定信号的频率。

size

(t

)

randn

t

x+

2

t

=π

(

2.1

(

));

sin(

)

*

*

2

50

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

实验三实验结果

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;

（2） 已知周期性序列如下所示：

}3,2,1,0,3,2,1,0,3,2,1,0{)(~ ↑

=n x 求其离散付里叶级数。 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);

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) 已知)(n x 是一个六点序列，如下所示：

??

?≤≤=e ls e

n n x 0

501

)(

要求计算该序列的离散时间的付里叶变换和离散付里叶变换，并绘出它们的幅度和相位。 %计算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)];

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

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

的区别。考虑如下序列：()n n n x ππ53.0cos(47.0cos()(+=

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

（a ） 当100≤≤n ，确定出的离散付里叶变换。(给序列进行补零加长至n=99) （b ） 当5000≤≤n ，确定出的离散付里叶变换。 解：（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);

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

从图看出此序列在π5.0

w处有一主频率，原始序列则没有这一频率，它只有主频为=

π

.0

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

=

53

47

.0

=

w和π

序列尾部补零，提供了更加细化的谱线。

（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为单位的频率');

axis([0 5 0 180]);

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

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

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

3.DFT 性质：

（1） 设n n x )8.0(10)(=,100≤≤n ，求15))6(()(-=n x n y

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

%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])

结果为：

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） 设}4,5,2,1{)(},2,2,1{)(21==n x n x 试分别计算下列圆卷积。 （a ）5),()()(21=?=N n x n x n y (b) 6),()()(21=?=N n x n x n y 解：先编循环卷积函数： %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';

求{}2,2,1)(1=n x ，{}4,5,2,1)(1=n x （1）)()()(211n x n x n y ?=

N=5 （2）)()()(212n x n x n y ?=

N=6

程序：

x1=[1,2,2];

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 计算线性卷积：

用重复法求解：1)(+=n n x ，90≤≤n ，{}1,0,1)(-=n h ，求其线卷积)(n y 。 解：先编一个重叠保留法作分段卷积函数； %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）用MA TLAB 的fft()来求信号的DFT 的幅值谱。

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

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

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

)

w

(t

t

t

+

x+

t

=π

)

)

8.0

(

)

2

8

sin(

cos(

4

5

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

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

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

四、思考题

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.

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

设t a e t x 1000)(-=，求其付里叶变换)(jw X z

解：dt e t X jw X jwt a z ?+∞

∞--=)()(

2

10000

1000)

1000

(

1002.0w dt

e

e

dt e

e

jwt

t

jwt

t

+=

+

=

-+∞

∞--??

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

在此中设在510-精度下，)(t x a 为2000

=F 的带限信号，因此若

取

5

5

10

252000

2110

5--?=?<<

?=?t

这样就可得到)(t x a 的逼近序列)(m x G 程序：

%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')

自已动手：

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

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