第3章离散傅里叶变换与MATLAB实现 离散傅里叶变换(Discrete Fourier Transform,DFT)是信号分析的基本方法,傅里叶变换是傅里叶分析的核心,通过它把信号从时域变换到频域,进而研究信号的频谱结构和变化规律。 DFT是傅里叶变换在时域和频域上都呈现离散的形式,将时域信号的采样变换为在离散时间傅里叶变换(DTFT)频域的采样。在形式上,变换两端(时域和频域上)的序列是有限长的,而实际上这两组序列都应当被认为是离散周期信号的主值序列。即使对有限长的离散信号做DFT,也应当将其看作是经过周期延拓成为周期信号再做的变换。 3.1有限长序列处理与MATLAB实现 给定序列的傅里叶变换,以及在同一序列不同长度情况下的DFT分析,程序如下: x = [1,1,1,1]; %信号序列 A = 200; %DTFT近似点数 w1 = 0:2/A:2; w = pi*w1.'; xw = 1+exp(-j*w)+exp(-j*w*2)+exp(-j*w*3)+exp(-j*w*4); %DTFT频谱分析 stem(w1.',abs(xw)) %DTFT频谱绘制 N = 32; y1 = fft(x,N); n = 0:N-1; subplot(3,1,1); %32点DFT频谱分析 stem(n,abs(y1),'ok'); title('N=32'); N = 64; y2 = fft(x,N); n=0:N-1; subplot(3,1,2); %64点DFT频谱分析 stem(n,abs(y2),'ok'); title('N=64'); N = 128; y3 = fft(x,N); n = 0:N-1; subplot(3,1,3); %128点DFT频谱分析 stem(n,abs(y3),'ok'); title('N=128'); 同一序列不同长度的DFT互不相同,如图31所示,图中描述了32点DFT、64点DFT以及128点DFT运行后得到的序列DFT输出图像。 图31序列DFT输出 傅里叶级数的系数计算量包括N次复数相乘以及N-1次复数相加。例如,对于N个点的傅里叶级数计算,DFT的全部处理计算量与N的平方成正比,程序如下: hi = 0.00001; t = -0.005:hi:0.005; %取值台阶 xn = exp(-500*abs(t)); %信号序列 Ww = 2*pi*2000; N = 1000; n = 0:2:N; W = n*Ww/N; %傅里叶变换 Xn = xn*exp(-j*t'*W)*hi; Xa = real(Xn); W = [-fliplr(W),W(2:501)]; Xn = [fliplr(Xn),Xn(2:501)]; subplot(1,2,1); plot(t*1000,xn,'o'); xlabel('t(s)'); ylabel('xn'); title('源信号'); subplot(1,2,2); plot(W/(2*pi*1000),Xn*1000,'o'); xlabel('f(KHz)'); ylabel('jw'); title('对应序列的傅里叶变换'); 程序运行后得到的序列DFT输出图像如图32所示。 图32DFT输出图像 3.2N点离散傅里叶变换与MATLAB实现 选择正弦函数与余弦函数作为基函数是因为它们的正交性,周期函数都可以用正弦函数和余弦函数构成的无穷级数来表示。对于周期为T的函数,n取不同值时的周期信号具有谐波关系(即它们都具有一个共同周期T)。n=0时,对应的这一项称为直流分量,k=1时具有基波频率。N点离散傅里叶变换与MATLAB实现程序如下: N = 32; N1 = 16; n = 0:N-1; k = 0:N1-1; x1n = exp(j*pi*n/8); %x1信号 y1 = fft(x1n,N); %[x1]N点DFT y2 = fft(x1n,N1); %[x1]N1点DFT x2n = cos(pi*n/16); %x2信号 y3 = fft(x2n,N); %[x2]N点DFT y4 = fft(x2n,N1); %[x2]N1点DFT subplot(2,2,1); stem(n,abs(y1),'o'); ylabel('|y1|') title('16点的DFT:x1') subplot(2,2,2); stem(n,abs(y3),'o'); ylabel('|y2|') title('16点的DFT:x2') subplot(2,2,3); stem(k,abs(y2),'>'); ylabel('|y1|') title('8点的DFT:x1') subplot(2,2,4); stem(k,abs(y4),'*'); ylabel('|y2|') title('8点的DFT:x2') x1信号和x2信号的周期延拓输出的为单频序列,其DFT由程序中的N=32与N1=16表示,输出图像如图33所示。 图33周期延拓序列DFT输出图像 有限长序列的频谱分析是基于离散傅里叶变换在时域和频域都有有限长序列,并且是离散的这种情况,程序如下: xn=[0.25 0.5 0.75 1] N=length(xn); n=0:(N-1); k=0:(N-1); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% xk=xn*exp(-j*2*pi/N).^(n'*k) xn0=(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,xk); 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)|'); xk 为序列的DFT,xn0为序列对应的IDFT(离散傅里叶逆变换),最终结论为原始序列xn=[0.25 0.5 0.75 1]与IDFT结果xn0相等。运行输出为: xk = 2.5000-0.5000 + 0.5000i -0.5000 - 0.0000i-0.5000 - 0.5000i xn0 = 0.2500 - 0.0000i0.5000 - 0.0000i 0.7500 - 0.0000i1.0000 + 0.0000i >> 基于mod函数实现圆周位移,可将N点序列循环左移或者右移m个采样周期,例如,对于特定序列,有如下: N = 4; n = 0:N-1; m = 2; xn = [4 2 1 1/2]; nm = mod((n-m),N); xm = xn(nm+1); subplot(1,2,1); stem(xn,'o'); xlabel('n'); ylabel('|xn|'); title('信息序列'); subplot(1,2,2); stem(xm,'o'); xlabel('n'); ylabel('|xm|'); title('位移序列'); xn为原始序列,xm为位移序列,xn与xm运行输出对比如图34所示。 图34基于mod函数实现的圆周位移 3.3序列圆周积分与MATLAB实现 序列圆周积分可以由x1、x2和y的波形分布得到,实现程序如下: n = 0:30; x1 = [0.25 1.5 0.5 1.5 0.75]; x2 = 0.9.^n; N = 30; x1 = [x1 zeros(1, N-length(x1))]; x2 = [x2 zeros(1, N-length(x2))]; for n = 1:N for m = 1:N y(n,m) = x1(m)*x2(mod( (n-m), N)+1 ); end end y = sum(y'); subplot(131); stem( x1 ); title('x1'); subplot(132); stem( x2 ); title('x2'); subplot(133); stem( y ); title('y'); x1与x2为圆周积分序列,y为序列圆周积分结果,程序运行后输出图像如图35所示。 图35序列圆周积分图像 3.4本章小结 DFT的特点之一是隐含的周期性,从表面上看,DFT在时域和频域都是非周期的有限长的序列,但实质上,DFT是从DFS(离散傅里叶级数)引申而来的,它们的本质是一致的,因此DFS的周期性决定了DFT具有隐含的周期性。 在实际工程中经常遇到的模拟信号频谱函数也是连续函数,为了利用DFT对模拟信号进行谱分析,首先对模拟信号进行时域采样得到离散量x(n),再对离散量进行DFT,得到的x(k)是x(n)的傅里叶变换在频率区间[0,单位周期]上的N点等间隔采样,这里x(n)和x(k)都是有限长序列。 然而,时间有限长的信号其频谱是无限宽的,反之,信号的频谱有限,则其持续时间将为无限长。因此,按采样定理采样时,采样序列应为无限长,这不满足DFT的条件。实际工程中,对于频谱很宽的信号,为防止时域采样后产生“频谱混叠”,一般用前置滤波器滤除幅度较小的高频成分,使信号的带宽小于折叠频率; 同样,对于持续时间很长的信号,采样点数太多也会导致存储和计算困难,一般也是截取有限点进行计算。用DFT对模拟信号进行谱分析,只能是近似的,其近似程度取决于信号带宽、采样频率和截取长度。