matlab在DSP中的应用（六）---离散傅里叶变换的性质

一、实验目的
(1)加深对离散傅里叶变换(DFT)基本性质的理解。

(2)了解有限长序列傅里叶变换(DFT)性质的研究方法。

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

二、实验原理

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。

例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)的线性性质
代码：

xn1=[0,1,2,4]; %建立xn1序列xn2=[1,0,1,0,1];%建立xn2序列N1=length(xn1);N2=length(xn2);N=max(N1,N2);%确定Nif N1>N2    xn2=[xn2,zeros(1,N1-N2)];%对长度短的序列补0elseif N2>N1     xn1=[xn1,zeros(1,N2-N1)];endyn=2*xn1+3*xn2;%计算ynn=0:N-1;k=0:N-1;Yk1=yn*(exp(-1j*2*pi/N)).^(n'*k);%求yn的N点DFTXk1=xn1*(exp(-1j*2*pi/N)).^(n'*k);%求xn1的N点DFTXk2=xn2*(exp(-1j*2*pi/N)).^(n'*k);%求xn2的N点DFTYk2=2*Xk1+3*Xk2;%由Xk1、Xk2求Ykfigure(1)subplot(4,1,1);stem(xn1);title('x1(n)')subplot(4,1,2);stem(xn2);title('x2(n)')subplot(4,1,3);stem(yn);title('y(n)')subplot(4,1,4);stem(k,Yk1);axis([0 4 -20 30])title('DFT(y(n))')figure(2)subplot(3,1,1);stem(k,Xk1);title('Xk1')subplot(3,1,2);stem(k,Xk2);title('Xk2')subplot(3,1,3);stem(k,Yk2);axis([0 4 -20 30])title('Yk')

输出：

2.循环移位性质

如果有限长序列为x(n)，长度为N，将x(n)左移m位，则:y(n)＝x((n＋m)N)RN(n)

x(n)左移m位的过程可由以下步骤获得：

(1)将x(n)以N为周期进行周期延拓，得到＝x((n)N)；
(2)将左移m位，得到；

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

有限长序列的移位也称为循环移位，原因是将x(n)左移m位时，移出的m位又依次从右端进入主值区。下面举例说明。

例2：
已知有限长序列x(n)＝［1，2，3，4，5，6］，求x(n)左移2位成为新的向量y(n)，并画出循环移位的中间过程。

xn=1:6;%xn序列Nx=length(xn);nx=0:Nx-1;nx1=-Nx:2*Nx-1;x1=xn(mod(nx1,Nx)+1);%建立周期延拓序列x1ny1=nx1+2;%将x1左移2位,左加右减y1=xn(mod(ny1,Nx)+1);%建立周期延拓序列y1RN=(nx1>=0)&(nx1<Nx);%设置主值区间subplot(4,1,1);stem(nx1,RN.*x1);title('主值序列')subplot(4,1,2);stem(nx1,x1);title('周期序列')subplot(4,1,3);stem(ny1,y1);title('移位周期序列')subplot(4,1,4);stem(ny1,RN.*y1);title('移位主值序列')

输出：

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)定义为：

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

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

例3：
求x(n)＝［1，2，3，4，5，6，7］，循环长度分别取N＝7，N＝10。
(1)画出x(n)的图形；
(2)画出x(－n)的图形

代码：

N1=7;x1=1:N1;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(4,1,1);stem(n1,x1);title('x(n) N=7')subplot(4,1,2);stem(n1,y1);title('x(-n) N=7')subplot(4,1,3);stem(n2,x2);title('x(n) N=10')subplot(4,1,4);stem(n2,y2);title('x(-n) N=10')

输出：

例4：
求x(n)＝［1，2，3，4，5，6，7］，循环长度取N＝7。求证：在时域循环折叠后的函数x(－n)，其对应的DFT在频域也作循环折叠，并取X(k)的共轭。

代码：

N1=7;x1=1:N1;n1=0:N1-1;k=0:N1-1;y1=x1(mod(-n1,N1)+1);%建立x(－n)，N＝7序列Xk=x1*exp(-1j*2*pi/N1).^(n1'*k)%求x(n)的DFTYk=y1*exp(-1j*2*pi/N1).^(n1'*k)%求x(－n)的DFT

输出：

Xk =  28.0000 + 0.0000i  -3.5000 + 7.2678i  -3.5000 + 2.7912i  -3.5000 + 0.7989i  -3.5000 - 0.7989i  -3.5000 - 2.7912i  -3.5000 - 7.2678iYk =  28.0000 + 0.0000i  -3.5000 - 7.2678i  -3.5000 - 2.7912i  -3.5000 - 0.7989i  -3.5000 + 0.7989i  -3.5000 + 2.7912i  -3.5000 + 7.2678i

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

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

1)时域循环卷积

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

若x(n)、h(n)和y(n)的DFT分别为X(k)、H(k)和Y(k)，则

2)频域循环卷积

利用时域和频域的对称性，可以得到频域卷积特性。若

则：

下面重点讨论时域循环卷积。时域循环卷积的方法有多种：

方法1：直接使用时域循环卷积。

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

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

即先分别求x1(n)、x2(n)的DFTX1(k)、X2(k)，再求Y(k)的IDFT获得y(n)。参考下图：

这部分太简单，不写matlab程序验证了。

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

其中

设DFT［x(n)］＝X(k)＝Re［X(k)］＋j*Im［X(k)］，则有:
　

三、实验任务

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

代码：

xn=[4,0,3,0,2,0,1];%xn序列Nx=length(xn);nx=0:Nx-1;nx1=-Nx:2*Nx-1;x1=xn(mod(nx1,Nx)+1);%建立周期延拓序列x1ny1=nx1+2;%将x1左移2位,左加右减y1=xn(mod(ny1,Nx)+1);%建立周期延拓序列y1RN=(nx1>=0)&(nx1<Nx);%设置主值区间subplot(4,1,1);stem(nx1,RN.*x1);title('主值序列')subplot(4,1,2);stem(nx1,x1);title('周期序列')subplot(4,1,3);stem(ny1,y1);title('移位周期序列')subplot(4,1,4);stem(ny1,RN.*y1);title('移位主值序列')

输出：

(2)已知一个有限长信号序列x(n)＝［8，7，6，5，4，3］，循环长度取N＝10。
求证：在时域循环折叠后的函数x(－n)，其对应的DFT在频域也作循环折叠。

代码：

N=10;x=[8,7,6,5,4,3];x=[x zeros(1,N-length(x))];n=0:N-1;k=0:N-1;y=x(mod(-n,N)+1);%建立x(－n)，N＝10序列Xk=x*exp(-1j*2*pi/N).^(n'*k)%求x(n)的DFTYk=y*exp(-1j*2*pi/N).^(n'*k)%求x(－n)的DFT

输出：

Xk =  Columns 1 through 8  33.0000 + 0.0000i   7.7361 -16.9273i   5.5000 - 3.4410i   3.2639 - 3.9960i   5.5000 - 0.8123i   3.0000 - 0.0000i   5.5000 + 0.8123i   3.2639 + 3.9960i  Columns 9 through 10   5.5000 + 3.4410i   7.7361 +16.9273iYk =  Columns 1 through 8  33.0000 + 0.0000i   7.7361 +16.9273i   5.5000 + 3.4410i   3.2639 + 3.9960i   5.5000 + 0.8123i   3.0000 - 0.0000i   5.5000 - 0.8123i   3.2639 - 3.9960i  Columns 9 through 10   5.5000 - 3.4410i   7.7361 -16.9273i

可以看到：在时域循环折叠后的函数x(－n)，其对应的DFT在频域也作循环折叠。

(3)已知两个有限长序列x1＝［5，4，－3，－2］，x2＝［1，2，3，0］，用DFT求时域循环卷积y(n)并用图形表示。

代码：

xn1=[5,4,-3,-2];xn2=[1,2,3,0];N=length(xn1);n=0:N-1;k=0:N-1;Xk1=xn1*exp((-1j*2*pi/N)).^(n'*k);%由x1(n)的DFT求X1(k)Xk2=xn2*exp((-1j*2*pi/N)).^(n'*k);%由x2(n)的DFT求X2(k)Yk=Xk1.*Xk2;yn=Yk*exp((1j*2*pi/N)).^(n'*k)/N;%由Y(k)的IDFT求y(n)yn=sign(real(yn)).*abs(yn);%取模值，消除DFT带来的微小复数影响yn1=conv(xn1,xn2);yn1(1)=yn1(1)+yn1(1+N);%循环卷积yn1(2)=yn1(2)+yn1(2+N);yn1(3)=yn1(3)+yn1(3+N);subplot(2,1,1);stem(n,yn);title('逆变换后得到的循环卷积')subplot(2,1,2);stem(n,yn1(1:4))title('理论上的循环卷积')

输出：

(4)运行如下的DFT程序，分析所得结果。

figure; imshow(phantom)phantom_freq = fft2(phantom);figure; imshow(ifft2(phantom_freq(1:2:size(phantom,1),:)))

输出结果：

代码分析：
第1行代码：绘制原图（Figure 1）
第2行代码：以图像矩阵为输入，进行fft变换
第3行代码：取fft变换的结果的奇数行的数据进行逆变换，绘制逆变换的结果（Figure 2）

运行结果分析：
fft变换的实质是DFT变换。从第3行代码的分析可知，逆变换是对fft变换结果的行进行操作的，那么可将原图视为一个一维的有限长序列，假设序列有N个值，那么问题就转化为：将有N个值的序列x[n]进行N点DFT变换，再进行N/2点的逆变换（对N点DFT变换的结果进行采样，采样方式是对变换的N点频域序列每两个点取一点）。由于N/2 < N，因此逆变换后的结果会发生重叠现象，具体重叠方式是：x[n]的前半个序列（N/2个点）和后半个序列（N/2个点）对应点相叠加。对应到图像中，则是：图像的上半个图像和下半个图像对应的发生重叠，因此会出现Figure 2的图像。

思考题：
①离散傅里叶变换(DFT)有哪些常用的基本性质？

答：线性性质、循环移位性质、循环折叠性质、圆周卷积性质、循环对称性质。

②简述离散傅里叶变换(DFT)时域循环卷积的基本方法，其与DTFT、DFS有何联系与区别？

答：首先按照 DTFT 的形式计算两个有限序列（假设两个序列长度的最大值为 N） 的卷积， 然后结果数
组中下标超过 N（从 1 开始） 的值 加到结果为 下标取 N 的模加 1 的下标的值上。

联系和区别：DFS是对DTFT进行频域的离散采样得到的，其在频域是周期离散的，而DFT对应的时域和频域信号均是对DFS取一个周期N个点的信号。