实验五 MATLAB 实现DFTMATLAB 为计算数据的离散快速傅时叶变换，提供了一系列丰富的数学函数，主要有fft 、ifft 、fft2、ifft2和czt 等。

当所处理的数据的长度为2的幂次时，采用基-2算法进行计算，计算速度会显著增加。

所以，要尽可能使所要处理的数据长度为2幂次或者用添零的方式来添补数据使之成为2的幂次。

1．fft 和ifft 函数 调用格式是： （1）()X fft Y =如果X 是向量，则采用傅时叶变换来求解X 的离散傅里叶变换；如果X 是矩阵，则计算该矩阵每一列的离散傅里叶变换；如果X 是()D N *维数组，则是对第一个非单元素的维进行离散傅里叶变换。

（2）()N X fft Y ,=N 是进行离散傅里叶变换的X 的数据长度，可以通过对X 进行补零或截取来实现。

（3）[]()dim ,,X fft Y =或()dim ,,N X fft Y =在参数dim 指定的维上进行离散傅里叶变换；当X 为矩阵时，dim 用来指定变换的实施方向：dim=1，表明变换按列进行；dim=2，表明变换按行进行。

函数ifft 的参数应用与函数fft 完全相同。

2．fft2和ifft2函数调用格式是： （1）()X fft Y 2=如果X 是向量，则此傅里叶变换即变成一维傅里叶变换fft ；如果X 是矩阵，则是计算该矩阵的二维快速傅里叶变换；数据二维傅里叶变换fft 2（X ）相当于()()''X fft fft ，即先对X 的列做一维傅里叶变换，然后再对变换结果的行做一维傅里叶变换。

（2）()N M X fft Y ,,2=通过对X 进行补零或截断，使得X 成为()N M *的矩阵。

函数ifft2的参数应用与函数fft2完全相同fftn 、ifftn 是对数据进行多维快速傅立变换，其应用与fft2、ifft2类似；在此，不再叙述。

3．czt 函数调用格式是：()A W M X czt X ,,,=式中X 是待变换的时域信号()n x ，其长度设为N ，M 是变换的长度，W 确定变换的步长，A 确定变换的起点。

若M=N ，A=1，则CZT 变成DFT 。

4．fftsfift Y=fftshift(X)用来重新排列X=fft(x)的输出，当X 为向量时，把X 的左右两半进行交换，从而将零频分量移至频谱的中心；如果X 为二维傅立叶变换的结果，它同时将X 的左右和上下部分进行交换。

5．fftfilt y=fftfilt(b,x)采用重叠相加法FFT 对信号向量x 快速滤波，得到输出序列向量y ，向量b 为FIR 滤波器的单位脉冲响应，h(n)=b(n+1),n=0,1,…,length(b)-1。

y=fftfilt(b,x,N)自动选取FFT 长度NF=2^nextpow2(N),输入数据x 分段长度M=NF-length(b)+1，其中nextpow2(N)函数求得一个整数，满足：2^(nextpow2(N)-1)≤N ≤2^nextpow2(N)N 缺省时，fftfilt 自动选择合适的FFT 长度NF 和对x 的分段长度M 。

一、利用MATLAB 实现的快速傅里的变换例1 已知有限长序列()n x 的长度4=N ，且()⎪⎪⎩⎪⎪⎨⎧==-===33211201n n n n n x ，用FFT 求()k X ，再用IFFT 求()n x 。

解：利用快速傅里叶变换函数求解的MA TLAB 实现程序清单如下clear xn=[1,2,-1,3]; X=fft(xn) x=ifft(X)程序运行结果如下X =5.0000 2.0000 + 1.0000i -5.0000 2.0000 - 1.0000i x =1 2 -1 3例2 设()n x 是由两个正弦信号及白噪声的叠加，试用FFT 文件对其作频谱分析。

解：程序清单如下% 产生两个正弦加白噪声；N=256;f1=.1;f2=.2;fs=1; a1=5;a2=3; w=2*pi/fs;x=a1*sin(w*f1*(0:N-1))+a2*sin(w*f2*(0:N-1))+randn(1,N); % 应用FFT 求频谱； subplot(2,2,1); plot(x(1:N/4)); title('原始信号'); f=-0.5:1/N:0.5-1/N; X=fft(x); y=ifft(X); subplot(2,2,2);plot(f,abs(X)); %也可求出将零频分量移至频谱的中心的频谱幅值plot(f,fftshift(abs(X))) title('频域信号'); subplot(2,2,3);plot(real(x(1:N/4))); %y=x title('时域信号');运行结果如图1所示，该程序同时完成了傅里叶变换与傅里叶反变换。

图1 例2运行结果二、 MATLAB 实现序列的移位和周期延拓运算例3：已知()()n R n x n 88.0=,利用MATLAB 生成并图示(),n x ()m n x -，()()()n R n x N 8，其中N=24，N m <<0，()()8n x 表示()n x 以8为周期的延拓。

图2 序列的移位和周期延拓解：程序清单如下N=24;M=8;m=3; %设移位值为3 n=0:N-1;x1=0.8.^n;x2=[(n>=0)&(n<M)];xn=x1.*x2; %产生x(n) xm=xn; nm=n+m %产生x(n-m)xc=xn(mod(n,8)+1); %产生x(n)的周期延拓,求余后加1是因为MA TLAB 向量下标从开始 xcm=xn(mod(n-m,8)+1); %产生x(n)移位后的周期延拓 subplot(2,2,1);stem(n,xn,'.');axis([0,length(n),0,1]);title('x(n)')subplot(2,2,2);stem(nm,xm,'.');axis([0,length(nm),0,1]);title('x(n-m)')subplot(2,2,3);stem(n,xc,'.');axis([0,length(n),0,1]);title('x((n)的周期延拓') subplot(2,2,4);stem(n,xcm,'.');axis([0,length(n),0,1]);title('x(n)的循环移位')程序运行结果如图2所示。

二、 MATLAB 验证N 点DFT 的物理意义例4：已知()()n R n x 4=，ωωωj j j e e n x FT e X ----==11)]([)(4，绘制相应的幅频和相频曲线，并计算图示N=8和N=16时的DFT 。

解：程序清单如下N1=8;N2=16; % 两种FFT 的变换长度 n=0:N1-1;k1=0:N1-1; k2=0:N2-1; w=2*pi*(0:2047)/2048;Xw=(1-exp(-j*4*w))./( 1-exp(-j*w));%对x(n)的频谱函数采样2048个点可以近似的看作是连续的频谱 xn=[(n>=0)&(n<4)]; %产生x(n)X1k=fft(xn,N1); %计算N1=8点的X 1(k) X2k=fft(xn,N2); %计算N2=16点的X 2(k) subplot(3,2,1);plot(w/pi,abs(Xw));xlabel(‘w/π’)subplot(3,2,2);plot(w/pi,angle(Xw));axis([0,2,-pi,pi]);line([0,2],[0,0]); xlabel(‘w/π’)subplot(3,2,3);stem(k1,abs(X1k),’.’);xlabel(‘k(ω=2πk/N1)’);ylabel(‘|X1(k)|’);hold onplot(N1/2*w/pi,abs(Xw)) %图形上叠加连续频谱的幅度曲线 subplot(3,2,4);stem(k1,angle(X1k)); axis([0,N1,-pi,pi]);line([0,N1],[0,0]);xlabel(‘k(ω=2πk/N1)’) ;ylabel(‘Arg[X1(k)]’);hold onplot(N1/2*w/pi,angle(Xw)) %图形上叠加连续频谱的相位曲线 subplot(3,2,5);stem(k2,abs(X2k),’.’);xlabel(‘k(ω=2πk/N2)’);ylabel(‘|X2(k)|’);hold on plot(N2/2*w/pi,abs(Xw))subplot(3,2,6);stem(k2,angle(X2k),’.’); axis([0,N2,-pi,pi]);line([0,N2],[0,0]);xlabel(‘k(ω=2πk/N2)’) ;ylabel(‘Arg[X2(k)]’);hold on plot(N2/2*w/pi,angle(Xw)) 程序运行结果如图3所示。

图3 离散傅里叶变换和傅里叶变换的采样关系前面的理论已经知道，序列x(n)的N 点DFT 的物理意义是对)(ωj e X 在[0，2π]上进行N 点的等间隔采样。

图3可以直观的看到X(k)与)(ωj e X 之间的采样关系。

实验六：MATLAB 实现DFT 的应用 一．MATLAB 实现快速卷积例1用快速卷积法计算下面两个序列的卷积)()4.0sin()(15n R n n x =，)(9.0)(20n R n h n =快速卷积计算框图如图所示解：程序清单如下M=15;N=20;nx=1:15;nh=1:20; xn=sin(0.4*nx);hn=0.9.^nh;L=pow2(nextpow2(M+N-1));%L 变为2的整数次幂 Xk=fft(xn,L); Hk= fft(hn,L); Yk=Xk.*Hk;yn=ifft(Yk,L);ny=1:Lsubplot(3,1,1);stem(nx,xn,'.');title('x(n)'); subplot(3,1,2);stem(nh,hn,'.');title('h(n)');subplot(3,1,3);stem(ny,real(yn),'.');title('y(n)');程序运行结果如图1所示。