实验五 MATLAB 实现DFT
MATLAB 为计算数据的离散快速傅时叶变换，提供了一系列丰富的数学函数，主要有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 ，且()⎪⎪⎩⎪
⎪⎨
⎧==-===3
3
211
20
1n 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)&(nxn=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 on
plot(N1/2*w/pi,abs(Xw)) %图形上叠加连续频谱的幅度曲线 subplot(3,2,4);stem(k1,angle(X1k)); axis([0,N1,-pi,pi]);line([0,N1],[0,0]);
