MATLAB离散傅里叶变换及应用

MATLAB 离散傅里叶变换及应用

一、DFT 与IDFT 、DFS 、DTFT 的联系

1、 序列的傅里叶变换(DFT)和逆变换(IDFT)

在实际中常常使用有限长序列。如果有限长序列信号为x(n)，则该序列的离散傅里叶变换对可以表示为

1N ,0,1,k ,

W x (n)DFT[x (n)]X(k)1

N 0n nk N -===∑-=Λ (12-1)

1

N ,0,1,n ,

W X(k)N 1IDFT[X(k)]x (n)1N 0

k nk N -===∑-=-Λ (12-2)

已知x(n)＝［0，1，2，3，4，5，6，7］，求x(n)的DFT 和IDFT 。要求：

(1)画出序列傅里叶变换对应的|X(k)|和arg ［X(k)］图形。 (2)画出原信号与傅里叶逆变换IDFT ［X(k)］图形进行比较。 程序源代码：

xn=[0,1,2,3,4,5,6,7];

N=length(xn); n=0:(N-1);k=0:(N-1);

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

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

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

title('x(n)');

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

title('IDFT|X(k)|');

subplot(2,2,3),stem(k,abs(Xk));

title('|X(k)|');

subplot(2,2,4),stem(k,angle(Xk));

title('arg|X(k)|'); 运行图如下：

x(n)

IDFT|X

(k)|

2

4

6

8

|X (k)|

2

4

6

8

arg|X (k)|

从得到的结果可见，与周期序列不同的是，有限长序列本身是仅有N 点的离散序列，相当于周期序列的主值部分。因此，其频谱也对应序列的主值部分，是含N 点的离散序列。

2、 序列DFT 与周期序列DFS

已知周期序列的主值x(n)＝［0，1，2，3，4，5，6，7］，求x(n)周期重复次数为4次时的DFS 。要求： (1)画出原主值和信号周期序列信号。

(2)画出序列傅里叶变换对应的(k)X

~

和(k)]X ~

arg[的图形。 程序源代码： xn=[0,1,2,3,4,5,6,7]; N=length(xn); n=0:4*N-1;k=0:4*N-1;

xn1=xn(mod(n,N)+1); Xk=xn1*exp(-j*2*pi/N).^(n'*k); subplot(2,2,1),stem(xn); title('原主值信号x(n)');

subplot(2,2,2),stem(n,xn1); title('周期序列信号');

subplot(2,2,3),stem(k,abs(Xk));

title('|X(k)|');

subplot(2,2,4),stem(k,angle(Xk)); title('arg|X(k)|'); 运行结果如下：

2468原主值信号x(n)

周期序列信

号

10

20

30

40

|X (k)|

10

20

30

40

arg|X (k)|

由这个周期序列的实验我们可以看出，有限长序列x(n)可以

看成是周期序列(n)x ~

的一个周期；反之，周期序列(n)x ~可以看成是有限长序列x(n)以N 为周期的周期延拓。频域上的情况也是相同的。从这个意义上说，周期序列只有有限个序列值有意义。 3、 序列DFT 与离散时间傅里叶变换DTFT 的联系

求x(n)＝［0，1，2，3，4，5，6，7］，0≤n ≤7的DTFT ，将(－2p ，2p)区间分成500份。要求： (1)画出原信号。

(2)画出由离散时间傅里叶变换求得的幅度谱X(ejw)和相位谱arg ［X(ejw)］图形。 程序源代码： xn=[0,1,2,3,4,5,6,7];

N=length(xn);

n=0:N-1;

w=linspace(-2*pi,2*pi,500);

X=xn*exp(-j*n'*w);

subplot(3,1,1),stem(n,xn,'k');

ylabel('x(n)');

subplot(3,1,2),plot(w,abs(X),'k');

axis([-2*pi,2*pi,1.1*min(abs(X)),1.1*max(abs(X))]); ylabel('幅度谱');

subplot(3,1,3),plot(w,angle(X),'k');

axis([-2*pi,2*pi,1.1*min(angle(X)),1.1*max(angle(X))]); ylabel('相位谱');

运行结果如下：

x (n )

-6-4

-2

2

4

6

10

2030

幅度谱

-6

-4

-2

2

4

6

-202

相位谱

由图12-3可以看出，两者有一定的差别。主要原因在于，该例进行DTFT 时，X(ejw)在单位圆上取250个点进行分割；而图12-1进行DFT 时，X(k)是在单位圆上N ＝8的等间距点上取值，X(k)的序列长度与X(ejw)相比不够长。

4仍然用x(n)＝［0，1，2，3，4，5，6，7］，将x(n)的有限长序列后面补足至N ＝100，求其DFT ，并与例3进行比较。 程序源代码： N=100;

xn=[0,1,2,3,4,5,6,7,zeros(1,N-8)]; n=0:(N-1);k=0:(N-1);

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

subplot(2,1,1),stem(k,abs(Xk)); title('|X(k)|');

subplot(2,1,2),stem(k,angle(Xk)); title('arg|X(k)|'); 运行结果如下：

0102030405060708090100

|X (k)|

10

20

30

40

50

60

70

80

90

100

arg|X (k)|

二、序列的移位和周期延拓运算。

已知

)

()8.0()(8n R n x n =，利用MATLAB 生成并图示序列

),

(),(m n x n x -和

)())((8n R n x N ),

())((8n R m n x N -其中

为周期的延拓。以表示8)())((,0,248n x n x N m N <<=

程序清单如下： N=24; M=8;

m=3;

n=0:N-1;

xn=0.8.^n.*(n>=0 & n

subplot(3,1,1);stem(n,xn,'.');grid;

axis([0 length(xn),0 1]);title('序列x(n)');

xc=xn(mod(n,8)+1);

subplot(3,1,2);stem(n,xc,'.');grid;

axis([0 length(xc),0 1]);title('序列x(n)的周期延拓序列');

xm=[xn(m+1:M) xn(1:m)];

xm=[xm zeros(1,N-length(xm))];

subplot(3,1,3);stem(n,xm,'.');grid;

axis([0 length(xm),0 1]);title('圆周移位序列x(n+m)');

运行结果如下：

序列x(n)

序列x(n)的周期延拓序列

05101520

圆周移位序列x(n+m)

三、利用MATLAB 验证N 点DFT 的物理意义。

ω

ω

j j jw

e e n x DFT e X DTFT n R n x ----=

==11)]([)()

(),()(44其离散时间傅立叶变换已知有限长序列

试绘制出)(ωj e X 幅度频谱和相位频谱，并分别计算N=8和N=16时的DFT 。 程序源代码：

N1=8;N2=16; % 设置两种DFT 的长度 n=0:N1-1; k1=n;k2=0:N2-1; w=(0:2047)*2*pi/2048;

Xw=(1-exp(-j*4*w))./(1-exp(-j*w)); xn=[n>=0 & n<4]; Xk1=fft(xn,N1); Xk2=fft(xn,N2); subplot(3,1,1);

plot(w/pi,abs(Xw));

grid;title('序列x(n)的幅频曲线|X(e^{j\omega})|'); subplot(3,1,2);

stem(k1*2/N1,abs(Xk1),'.'); grid;title('序列x(n)的8点DFT'); subplot(3,1,3);stem(k2,abs(Xk2),'.'); grid;title('序列x(n)的16点DFT');

运行结果如下；

00.20.40.60.81 1.2 1.4 1.6 1.82

02

4序列x(n)的幅频曲线|X (e j )|

序列x(n)的8点

DFT

序列x(n)的16点DFT

四、利用DFT 计算线性卷积

已知序列x1(n)=0.9^n,n=[0:11];h(n)=R 9(n) 求x1(n)*h(n)；x1(n)与h(n)的10点圆周卷积。 程序源代码:

n1=[0:9];n2=[0:11];m=0:N1-1;n=[0:N1-1];N=12;N1=10;x1=0.9.^n2;

x11=0.9.^n1;x2=ones(1,9);x3=conv(x1,x2) x5=[x11,zeros(1,N1-length(x11))]; x6=[x2,zeros(1,N1-length(x2))]; H=zeros(N1,N1);

x6=[x6 zeros(1,N1-length(x6))];

for n=1:N1

H(n,:)=x6(mod(n-m-1,N1)+1);

end

x4=x5*H';

subplot(221),stem(x1,'.');title('原序列x1')

axis([-1,14,0.6*min(x1),1.1*max(x1)]);

subplot(222),stem(x2,'.');title('原序列x2')

axis([-1,11,0,1.1*max(x2)]);

subplot(223),stem(x3,'.');title('x1卷积x2')

axis([-1,25,1.1*min(x3),1.1*max(x3)]);

subplot(224),stem(x4,'.');title('x1 10点圆周卷积x2') axis([-1,15,0.975*min(x4),1.01*max(x4)]);

运行结果如下：

,.

510

0.2

0.40.60.8

1原序列x1

5

10

0.20.40.60.8

1原序列x2

2

4

6

x1卷积x2

51015

5.4

5.65.86

x1 10点圆周卷积x2