数字信号处理实验 matlab版 离散傅里叶变换的性质

更新时间:2023-09-05 15:38:01 阅读量: 教育文库 文档下载

说明:文章内容仅供预览,部分内容可能不全。下载后的文档,内容与下面显示的完全一致。下载之前请确认下面内容是否您想要的,是否完整无缺。

数字处理实验 matlab版 山大学生最适用 本人自己写的 因为时间比较久了 不能完全保证出现代码都能运行 但95%还是能保证的 谢谢

实验13 离散傅里叶变换的性质

(完美格式版,本人自己完成,所有语句正确,不排除极个别错误,特别适用于山大,勿用冰点等工具下载,否则下载之后的word格式会让很多部分格式错误,谢谢)

XXXX学号姓名处XXXX

一、实验目的

1 加深对离散傅里叶变换(DFT)基本性质的理解。 2 了解有限长序列傅里叶变换(DFT)性质的研究方法。

3 掌握用MATLAB语言进行离散傅里叶变换性质分析时程序编写的方法。

二、实验内容

1 线性性质。 2 循环移位性质。 3 循环折叠性质。

4 时域和频域循环卷积特性。 5 循环对称性。

三、实验环境

MATLAB7.0

四、实验原理

1 线性性质

如果两个有限长序列分别为x1(n)和x2(n),长度分别为N1和N2,且

y(n)=ax1(n)+bx2(n) (a、b均为常数)

则该y(n)的N点DFT为

Y(k)=DFT[y(n)]=aX1(k)+bX2(k) 0≤k≤N-1

其中:N=max[N1,N2],X1(k)和X2(k)分别为x1(n)和x2(n)的N点DFT。

例13-1 已知x1(n)=[0,1,2,4],x2(n)=[1,0,1,0,1],求: (1)y(n)=2x1(n)+3x2(n),再由y(n)的N点DFT获得Y(k);

(2)由x1(n)、x2(n)求X1(k)、X2(k),再求Y(k)=2X1(k)+3X2(k)。

用图形分别表示以上结果,将两种方法求得的Y(k)进行比较,由此验证有限长序列傅里叶变换(DFT)的线性性质。

解 MATLAB程序如下:

>> xn1=[0,1,2,4]; %建立xn1序列 >> xn2=[1,0,1,0,1]; %建立xn2序列 >> N1=length(xn1);N2=length(xn2); >> N=max(N1,N2); %确定N

>> if N1>N2 xn2=[xn2,zeros(1,N1-N2)]; %对长度短的序列补0 >> elseif N2>N1 xn1=[xn1,zeros(1,N2-N1)]; >> end

>> yn=2*xn1+3*xn2; %计算yn >> n=0:N-1;k=0:N-1;

数字处理实验 matlab版 山大学生最适用 本人自己写的 因为时间比较久了 不能完全保证出现代码都能运行 但95%还是能保证的 谢谢

>> Yk1=yn*(exp(-j*2*pi/N)).^(n'*k); %求yn的N点DFT >> Xk1=xn1*(exp(-j*2*pi/N)).^(n'*k); %求xn1的N点DFT >> Xk2=xn2*(exp(-j*2*pi/N)).^(n'*k); %求xn2的N点DFT >> Yk2=2*Xk1+3*Xk2; %由Xk1、Xk2求Yk >> subplot(4,2,1),stem(n,xn1); >> title('x1(n)');

>> subplot(3,2,2),stem(n,Xk1); >> title('X1(k)');

>> subplot(4,2,3),stem(n,xn2); >> title('x2(n)');

>> subplot(3,2,4),stem(n,Xk2); >> title('X1(k)');

>> subplot(4,2,5),stem(n,yn); >> title('yn');

>> subplot(3,2,6),stem(n,Yk2); >> title('2*Xk1+3*Xk2');

>> subplot(4,2,7),stem(n,Yk1); >> title('DFT[y(n)]'); 求得的Y(k),如下所示: Yk=

23.0000 -7.5902+1.5388i 3.5902-0.3633i 3.5902+0.3633i -7.5902-1.5388i 运行结果如图13-1所示。

图13-1 例13-1有限长序列的傅里叶变换的线性性质

数字处理实验 matlab版 山大学生最适用 本人自己写的 因为时间比较久了 不能完全保证出现代码都能运行 但95%还是能保证的 谢谢

2 循环移位性质

如果有限长序列为x(n),长度为N,将x(n)左移m位,则y(n)=x((n+m)N)RN(n) x(n)左移m位的过程可由以下步骤获得:

(1)将x(n)以N为周期进行周期延拓,得到x(n)=x((n)N);

~~

(2)将x(n)左移m位,得到x(n m);

x(n m)的主值序列,得到x(n)循环移位序列y(n)。 (3)取~

有限长序列的移位也称为循环移位,原因是将x(n)左移m位时,移出的m位又依次从

右端进入主值区。下面举例说明。

例13-2 已知有限长序列x(n)=[1,2,3,4,5,6],求x(n)左移2位成为新的向量y(n),并画出循环移位的中间过程。 解 MATLAB程序如下:

>> xn=[1,2,3,4,5,6]; %建立xn序列 >> Nx=length(xn);nx=0:Nx-1; >> nx1=-Nx:2*Nx-1; %设立周期延拓的范围 >> x1=xn(mod(nx1,Nx)+1); %建立周期延拓序列 >> ny1=nx1-2;y1=x1; %将x1左移2位,得到y1 >> RN=(nx1>=0)&(nx1<Nx); %在x1的位置向量nx1上设置主值窗 >> RN1=(ny1>=0)&(ny1<Nx); %在y1的位置向量ny1上设置主值窗 >> subplot(4,1,1),stem(nx1,RN.*x1); %画出x1的主值部分 >> subplot(4,1,2),stem(nx1,x1); %画出x1 >> subplot(4,1,3),stem(ny1,y1); %画出y1 >> subplot(4,1,4),stem(ny1,RN1.*y1); %画出y1的主值部分 运行结果如图13-2所示。

-6

-8-6

-4

-2

2

4

6

8

10

-4

-2

2

4

6

8

10

12

~

数字处理实验 matlab版 山大学生最适用 本人自己写的 因为时间比较久了 不能完全保证出现代码都能运行 但95%还是能保证的 谢谢

图13-2 例13-2有限长序列的循环移位

3 循环折叠性质

如果要把有限长N点序列x(n)直接进行折叠,则x的下标(-n)将不在0≤n≤N-1区域内。但根据有限长序列傅里叶变换隐含的周期性,可以对变量(-n)进行N求余运算。即在MATLAB中,序列x(n)的折叠可以由y=x(mod(-nx,N)+1)得到。

有限长N点序列x(n)的循环折叠序列y(n)定义为

n=0 x(0)

y(n)=x((-n)N)

x(N-n)1 n N-1

可以想像成,序列x(n)以反时针方向等间隔放置在一个圆周上,则x(-n)是将x(n)沿着

圆周顺时针方向等间隔放置。

循环折叠性质同样适用于频域。经循环折叠后,序列的DFT由下式给出:

Y(k) DFT[x(( n)N)]

X(0)

X (( k)N)

X(N-k)

k=0

1 k N-1

就是说,在时域循环折叠后的函数,其对应的DFT在频域也作循环折叠,并取X(k)的共轭。

例13-3 求x(n)=[1,2,3,4,5,6,7],循环长度分别取N=7,N=10。 (1)画出x(n)的图形; (2)画出x(-n)的图形。 解 MATLAB程序如下:

>> x1=[1,2,3,4,5,6,7]; %建立x(n),N=7序列 >> N1=length(x1);n1=0:N1-1; >> y1=x1(mod(-n1,N1)+1); %建立x(-n),N=7序列 >> N2=10;

>> x2=[x1,zeros(1,N2-N1)]; %建立x(n),N=10序列 >> n2=0:N2-1;

>> y2=x2(mod(-n2,N2)+1); %建立x(-n),N=10序列 >> subplot(2,2,1),stem(n1,x1,'k'); %画x(n),N=7 >> title('x(n),N=7');

>> subplot(2,2,3),stem(n1,y1,'k'); %画x(-n),N =7 >> title('x(-n),N=7');

>> subplot(2,2,2),stem(n2,x2,'k');% 画x(n),N=10 >> title('x(n),N=10');

>> subplot(2,2,4),stem(n2,y2,'k'); %画x(-n),N =10 >> title('x(-n),N=10'); 运行结果如图13-3所示。

数字处理实验 matlab版 山大学生最适用 本人自己写的 因为时间比较久了 不能完全保证出现代码都能运行 但95%还是能保证的 谢谢

x(n),N=7

x(n),N=10

02

x(-n),N=7

46

x(-n),N=10

2

4

6

图13-3 例13-3离散序列的循环折叠

例13-4 如例13-3求x(n)=[1,2,3,4,5,6,7],循环长度取N=7。求证:在时域循环折叠后的函数x(-n),其对应的DFT在频域也作循环折叠,并取X(k)的共轭。 解 MATLAB程序如下:

>> x1=[1,2,3,4,5,6,7]; %建立x(n),N=7序列 >> N=length(x1); >> n=0:N-1;k=0:N-1; >> y1=x1(mod(-n,N)+1); %建立x(-n),N=7序列 >> Xk=x1*exp(-j*2*pi/N).^(n'*k) %求x(n)的DFT >> Yk=y1*exp(-j*2*pi/N).^(n'*k) %求x(-n)的DFT 运行结果: Xk =

Columns 1 through 5

28.0000 -3.5000 + 7.2678i -3.5000 + 2.7912i -3.5000 + 0.7989i -3.5000 - 0.7989i Columns 6 through 7

-3.5000 - 2.7912i -3.5000 - 7.2678i Yk =

Columns 1 through 5

28.0000 -3.5000 - 7.2678i -3.5000 - 2.7912i -3.5000 - 0.7989i -3.5000 + 0.7989i Columns 6 through 7

-3.5000 + 2.7912i -3.5000 + 7.2678i

数字处理实验 matlab版 山大学生最适用 本人自己写的 因为时间比较久了 不能完全保证出现代码都能运行 但95%还是能保证的 谢谢

4 时域和频域循环卷积特性

离散傅里叶变换的循环卷积特性也称为圆周卷积,分为时域卷积和频域卷积两类。

1)时域循环卷积

假定x(n)、h(n)都是N点序列,则时域循环卷积的结果y(n)也是N点序列:

若x(n)、h(n)和y(n)的DFT分别为X(k)、H(k)和Y(k),则Y(k)=X(k) H(k) 2)频域循环卷积

利用时域和频域的对称性,可以得到频域卷积特性。若y(n)=x(n)h(n) 则

y(n)

x(n)h(n)

Y(k)

X(k)H(k)

下面重点讨论时域循环卷积。时域循环卷积的方法有多种: 方法1:直接使用时域循环卷积。

由于有限长序列可以看成是周期序列的主值,因此,时域圆周卷积的结果可以由对应的周期序列卷积和取主值部分获得,请参考例11-4。

方法2:用频域DFT相乘再求逆变换。

即先分别求x1(n)、x2(n)的DFTX1(k)、X2(k),再求Y(k)的IDFT获得y(n)。 方法3:用FFT和IFFT进行循环卷积。

基本思路同方法2,但直接使用了MATLAB提供的fft和ifft子函数来实现。见后面的快速傅里叶变换实验。

例13-5 将例11-4已知的两个时域周期序列分别取主值,得到x1=[1,1,1,0,0,0],x2=[0,1,2,3,0,0],求时域循环卷积y(n)并用图形表示。

解 本例采用方法2。程序如下:

>> xn1=[0,1,2,3,0,0]; %建立x1(n)序列 >> xn2=[1,1,1,0,0,0]; %建立x2(n)序列 >> N=length(xn1); >> n=0:N-1;k=0:N-1;

>> Xk1=xn1*(exp(-j*2*pi/N)).^(n'*k); %由x1(n)的DFT求X1(k) >> Xk2=xn2*(exp(-j*2*pi/N)).^(n'*k); %由x2(n)的DFT求X2(k) >> Yk=Xk1.*Xk2; %Y(k)=X1(k)X2(k)

>> yn=Yk*(exp(j*2*pi/N)).^(n'*k)/N; %由Y(k)的IDFT求y(n) >> yn=abs(yn) %取模值,消除DFT带来的微小复数影响

>> subplot(2,3,1),stem(n,xn1); >> title('x1(n)')

>> subplot(2,3,2),stem(n,xn2); >> title('x2(n)')

>> subplot(2,3,3),stem(n,yn); >> title('y(n)')

>> subplot(2,3,4),stem(n,Xk1); >> title('X1(k)')

>> subplot(2,3,5),stem(n,Xk2); >> title('X2(k)')

>>

subplot(2,3,6),stem(n,Yk);

数字处理实验 matlab版 山大学生最适用 本人自己写的 因为时间比较久了 不能完全保证出现代码都能运行 但95%还是能保证的 谢谢

>> title('Y(k)') 得到: yn =

0.0000 1.0000 3.0000 6.0000 5.0000 3.0000

运行结果如图13-5所示。由y(n)图形可见,与例11-4主值区域的卷积结果相同。

X1(k)

X2(k)

Y(k)

5

5

05

图13-5 例13-5离散序列时域循环卷积的结果

5 循环对称性

由于序列x(n)及其离散傅里叶变换X(k)的定义在主值为0~N-1的区间,因此DFT的循环对称性对时间序列是指关于n=0和n=N/2的对称性,对频谱序列是关于数字频率为0和p的对称性。

本实验重点分析实序列的循环对称性。

实序列x(n)可以分解为循环偶序列xe(n)和循环奇序列xo(n): x(n)=xe(n)+xo(n) 0≤n≤N-1 其中:

11

xe(n) [x(n) x( n)] , xo(n) [x(n) x( n)]

22

设DFT[x(n)]=X(k)=Re[X(k)]+j*Im[X(k)],则有

DFT[xe(n)] Re[X(k)]

DFT[xo(n)] j Im[X(k)]

即实序列中的偶序列xe(n)对应于x(n)的离散傅里叶变换X(k)的实部,而实序列中的奇

数字处理实验 matlab版 山大学生最适用 本人自己写的 因为时间比较久了 不能完全保证出现代码都能运行 但95%还是能保证的 谢谢

序列xo(n)对应于x(n)的离散傅里叶变换X(k)的虚部。

例13-6 已知一个定义在主值区间的实序列x=[ones(1,4),zeros(1,4)],试将其分解成为偶对称序列和奇对称序列,并求它们的DFT,验证离散傅里叶变换的循环对称性。 解 程序如下:

>> x=[ones(1,5),zeros(1,5)] %建立x(n)序列 >> N=length(x);

>> n=0:N-1;k=0:N-1;

>> xr=x(mod(-n,N)+1); %求x(-n) >> xe=0.5*(x+xr) %求x(n)的偶序列 >> xo=0.5*(x-xr) %求x(n)的奇序列

>> X=x*(exp(-j*2*pi/N)).^(n'*k); %由x(n)的DFT求X(k) >> Xe=xe*(exp(-j*2*pi/N)).^(n'*k); %由xe(n)的DFT求 Xe(k) >> Xo=xo*(exp(-j*2*pi/N)).^(n'*k); %由xo(n)的DFT求Xo(k)

>> error1=(max(abs(real(X)-Xe))) %计算X(k)的实部与Xe(k)的差值 >> error2=(max(abs(j*imag(X)-Xo))) %计算X(k)的虚部与Xo(k)的差值 >> subplot(2,4,1),stem(n,x); >> title('x(n)')

>> subplot(2,4,2),stem(n,xr); >> title('x(-n)')

>> subplot(2,4,3),stem(n,xe); >> title('xe(n)')

>> subplot(2,4,4),stem(n,xo); >> title('xo(n)')

>> subplot(2,4,5),stem(n,real(X)); >> title('X(k)的实部')

>> subplot(2,4,6),stem(n,imag(X)); >> title('X(k)的虚部')

>> subplot(2,4,7),stem(n,Xe); >> title('Xe(k)=DFT(xe(n))') >> subplot(2,4,8),stem(n,Xo); >> title('Xo(k)=DFT(xo(n))') 运行结果显示: x =

1 1 1 1 1 0 0 0 0 0 xe =

Columns 1 through 8

1.0000 0.5000 0.5000 0.5000 0.5000 0 0.5000 0.5000 Columns 9 through 10 0.5000 0.5000 xo =

Columns 1 through 8

0 0.5000 0.5000 0.5000 0.5000 0 -0.5000 -0.5000 Columns 9 through 10

数字处理实验 matlab版 山大学生最适用 本人自己写的 因为时间比较久了 不能完全保证出现代码都能运行 但95%还是能保证的 谢谢

-0.5000 -0.5000 error1 =

4.0932e-015 error2 =

3.9475e-015

x(n)

0-15

2

Xe(k)=DFT(xe(n))

0246810

02

x(-n)

46

Xo(k)46

xe(n)

X(k)的虚部

810

0246810

810

图13-6 例13-6验证离散实序列的循环对称性

由以上输出数据和图形可知:

(1)xe(n)具有循环对称性。对称中心在n=0和n=5处。

(2)xo(n)具有循环反对称性。对称中心亦在n=0和n=5处。

(3)从图上看,Xe(k)与X(k)的实部相等,Xo(k)与X(k)的虚部相等;从输出数据也可见,error1和error2的差约为0。即可证明,时域的偶、奇分量的确对应于频域的离散傅里叶变换的实部和虚部。

五、实验过程

2 已知有限长序列x(n)=[4,0,3,0,2,0,1],求x(n)右移2位成为新的向量y(n),并画出循环移位的中间过程。

解 MATLAB程序如下:

>> xn=[4 0 3 0 2 0 1];

>> Nx=length(xn);nx=0:Nx-1; >> nx1=-Nx:2*Nx-1;

>> x1=xn(mod(nx1,Nx)+1); >> ny1=nx1+2;y1=x1;

>> RN=(nx1>=0)&(nx1<Nx);

数字处理实验 matlab版 山大学生最适用 本人自己写的 因为时间比较久了 不能完全保证出现代码都能运行 但95%还是能保证的 谢谢

>> RN1=(ny1>=0)&(ny1<Nx);

>> subplot(4,1,1),stem(nx1,RN.*x1); >> title('主值序列');

>> subplot(4,1,2),stem(nx1,x1); >> title('周期序列');

>> subplot(4,1,3),stem(nx1,y1); >> title('移位周期序列');

>> subplot(4,1,4),stem(nx1,RN.*y1); >> title('移位主值序列'); 运行结果如图13-7所示。

主值序列

周期序列

移位周期序列

图13-7

3 已知一个有限长序列x(n)=[8 7 6 5 4 3],循环长度取N=10。求证:在时域循环折叠后的函数x(-n),其对应的DFT在频域也作循环折叠。

解 MATLAB程序如下:

>> N=10;

>> x1=[8,7,6,5,4,3,zeros(1,N-6)]; >> n=0:N-1;k=0:N-1; >> y1=x1(mod(-n,N)+1);

>> Xk=x1*exp(-j*2*pi/N).^(n'*k) >> Yk=y1*exp(-j*2*pi/N).^(n'*k) 运行结果: Xk =

Columns 1 through 6

33.0000 7.7361 -16.9273i 5.5000 - 3.4410i 3.2639 - 3.9960i 5.5000 - 0.8123i

数字处理实验 matlab版 山大学生最适用 本人自己写的 因为时间比较久了 不能完全保证出现代码都能运行 但95%还是能保证的 谢谢

3.0000 - 0.0000i

Columns 7 through 10

5.5000 + 0.8123i 3.2639 + 3.9960i 5.5000 + 3.4410i 7.7361 +16.9273i Yk =

Columns 1 through 6

33.0000 7.7361 +16.9273i 5.5000 + 3.4410i 3.2639 + 3.9960i 5.5000 + 0.8123i 3.0000 - 0.0000i Columns 7 through 10

5.5000 - 0.8123i 3.2639 - 3.9960i 5.5000 - 3.4410i 7.7361 -16.9273i

4 已知两个有限长序列x1=[5 4 -3 -2],x2=[1 2 3 0],用DFT求时域循环卷积y(n)并用图形表示。

解 MATLAB程序如下:

>> xn1=[5,4,-3,-2];xn2=[1,2,3,0]; >> N=length(xn1); >> n=0:N-1;k=0:N-1;

>> Xk1=xn1*(exp(-j*2*pi/N)).^(n'*k); >> Xk2=xn2*(exp(-j*2*pi/N)).^(n'*k); >> Yk=Xk1.*Xk2;

>> yn=Yk*(exp(j*2*pi/N)).^(n'*k)/N; >> yn=abs(yn)

>> subplot(2,3,1),stem(n,xn1); >> title('x1(n)')

>> subplot(2,3,2),stem(n,xn2); >> title('x2(n)')

>> subplot(2,3,3),stem(n,yn); >> title('y(n)')

>> subplot(2,3,4),stem(n,Xk1);

>> title('X1(k)') >> subplot(2,3,5),stem(n,Xk2);

>> title('X2(k)') >> subplot(2,3,6),stem(n,Yk); >> title('Y(k)')

运行结果:

yn = 024

2

8.0000 8.0000 20.0000 4.0000

X2(k)

Y(k)

运行结果如图13-8所示。

02402

图13-8

4

4

数字处理实验 matlab版 山大学生最适用 本人自己写的 因为时间比较久了 不能完全保证出现代码都能运行 但95%还是能保证的 谢谢

六、实验感想

通过此次实验中练习使用matlab语言进行离散傅里叶级数变换的性质实验,更为熟悉的掌握了matlab的功能,在实验过程中也遇到很多小问题,并通过仔细检查和查阅相关书籍解决此类问题,让我深刻认识到,细节的重要性。在使用help过程中,深切体会到良好的英语基础和充实的课堂知识的重要性。

本文来源:https://www.bwwdw.com/article/6o3i.html

Top