第14章富里哀分析 象富里哀级数,富里哀变换以及它们离散时间相应部分构成了信号处理的基础。为了 便于这类问题的分析, MATLAB提供了函数frff,i2和 fftshift。这类函数集执行 维和二维离散富里哀变换及其逆变换。这些函数允许人们完成很多信号处理任务。除此之外, 还可在可选的信号处理工具箱中得到其他扩展的信号处理工具 因为信号处理包含如此广泛的领域,甚至要说明用 MATLAB中离散富里哀变换函数可 解决的这类小问题,就超出了本书的范围。因此,这里将只介绍用函数ft近似连续时间信 号的富里哀变换的一个例子。此外,还将讨论《精通 MATLAB工具箱》中处理富里哀级数 的函数集 14.1快速富里哀变换 在MA∏LAB中,函数ft计算一个信号的离散富里哀变换。在数据的长度是2的幂次 或质因数的乘积的情况下,就用快速富里哀变换(FFT来计算离散富里哀变换。当数据长度 是2的幂次时,计算速度显著增加,因此,只要可能,选择数据长度为2的幂次或者用零来 填补数据,使得数据长度等于2的幂次显得非常重要。在《 MATLAB参考指南》中可找到 有关该问题的讨论 MATLAB中实现的快速富里哀变换,是按照工科教材中常使用的方法。 F(k=FFT(f(n); F(k)=2f(n)e" J2Tnk/Nk=0, 1,-.N-1 因为 MATLAB不允许零下标,所以移动了一个下标值 F(k)=2f(n)e"J2T(n-I(k-I)/N k=1,2, 相应的逆变换为: f(n)=FFT F(K)) n) NkEF(k)e jx(n-1Xk-ID/N n=1,2,……,N 为了说明FFT的使用,考虑估计连续信号的富里哀变换的问题
第 14 章 富里哀分析 象富里哀级数，富里哀变换以及它们离散时间相应部分构成了信号处理的基础。为了 便于这类问题的分析，MATLAB 提供了函数 fft,ifft,fft2,ifft2 和 fftshift。这类函数集执行一 维和二维离散富里哀变换及其逆变换。这些函数允许人们完成很多信号处理任务。除此之外， 还可在可选的信号处理工具箱中得到其他扩展的信号处理工具。 因为信号处理包含如此广泛的领域，甚至要说明用 MATLAB 中离散富里哀变换函数可 解决的这类小问题，就超出了本书的范围。因此，这里将只介绍用函数 fft 近似连续时间信 号的富里哀变换的一个例子。此外，还将讨论《精通 MATLAB 工具箱》中处理富里哀级数 的函数集。 14.1 快速富里哀变换 在 MATLAB 中，函数 fft 计算一个信号的离散富里哀变换。在数据的长度是 2 的幂次 或质因数的乘积的情况下，就用快速富里哀变换（FFT)来计算离散富里哀变换。当数据长度 是 2 的幂次时，计算速度显著增加，因此，只要可能，选择数据长度为 2 的幂次或者用零来 填补数据，使得数据长度等于 2 的幂次显得非常重要。在《MATLAB 参考指南》中可找到 有关该问题的讨论。 MATLAB 中实现的快速富里哀变换，是按照工科教材中常使用的方法。 F(k)=FFT{f(n)} F k f n e j nk N n N ( ) ( ) / = − = − 2 0 1 k = 0,1,,N -1 因为 MATLAB 不允许零下标，所以移动了一个下标值。 F k f n e j n k N n N ( ) ( ) ( )( )/ = − − − = 2 1 1 1 k = 1,2，,N 相应的逆变换为： f n FFT F K f n N F k e j n k N k N ( ) { ( )} ( ) ( ) ( )( )/ = = − − − − = 1 2 1 1 1 1 n = 1,2,,N 为了说明 FFT 的使用，考虑估计连续信号的富里哀变换的问题
f(t t<0 解析上,该富里哀变换为: 12 F(w) 3+w 虽然在这种情况下,由于知道了富里哀变换的解析结果,再运用FFT没有多大的实用 价值,但这个例子说明了对不常见的信号,特别是那些解析上难以找到富里哀变换的信号, 一个估计富里哀变换的方法。下面的 MATLAB语句用FFT估计F(w),并且用图形把所得 到结果与上面的解析表达式的结果进行比较: >>N=128; % choose a power of 2 for speed >>t=linspace(0, 3, N); time points for function evaluation >>f2*exp(-3t); evaluate the function and minimize aliasing: f(3)-0 >>Ts=t(2)-t(1); the sampling period >>Ws=2*pi/Ts, the sampling frequency in rad/sec >>F=fft(f); compute the fft >>Fp=F(1:N2+1)*Ts 0.7 0.6 0.4 0.3 0 80 1201 Frequency, Rad/s 图14.1富里哀变换两种结果的比较 仅从F中取正频率分量,并且乘以采样间隔计算F(w)
f t e t ( ) = − 12 0 3 t 0 t < 0 解析上，该富里哀变换为： F w jw ( ) = + 12 3 虽然在这种情况下，由于知道了富里哀变换的解析结果，再运用 FFT 没有多大的实用 价值，但这个例子说明了对不常见的信号，特别是那些解析上难以找到富里哀变换的信号， 一个估计富里哀变换的方法。下面的 MATLAB 语句用 FFT 估计 F(w)，并且用图形把所得 到结果与上面的解析表达式的结果进行比较： >>N=128; % choose a power of 2 for speed >>t=linspace(0, 3, N); % time points for function evaluation >>f=2*exp(-3*t); % evaluate the function and minimize aliasing:f(3)~0 >>Ts=t(2)-t(1); % the sampling period >>Ws=2*pi/Ts; % the sampling frequency in rad/sec >>F=fft(f); % compute the fft >>Fp=F(1 : N/2+1)*Ts; 0 20 40 60 80 100 120 140 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Frequency, Rad/s |F(w)| 图 14.1 富里哀变换两种结果的比较 仅从 F 中取正频率分量，并且乘以采样间隔计算 F(w)
>>W=Ws*(0:N2)N 它建立了连续频率轴,该轴起始于0,终止于奈魁斯特( Nyquist)频率Ws2 >>Fa=2 /(3+]*w);% evaluate analytical Fourier transform >>plot(W, abs( Fa),W, abs(Fp), +) generate plot, + mark fft results >>xlabel( Frequency, Rad/s), ylabel('IF(w)I) MATLAB提供了大量的完成一般信号处理任务的函数。它们列于表141: 表141 信号处理函数 conv 卷积 2维卷积 快速富里哀变换 fft2 2维快速富里哀变换 ifft 逆快速富里哀变换 fft2 2维逆快速富里哀变换 离散时间滤波器 filter 2 2维离散时间滤波器 幅值 四个象限的相角 unwrap 在360°边界清除相角突变 fftshift 把FFT结果平移到负频率上 nextpow2 2的下一个较高幂次 142富里哀级数 MATLAB本身没有特别关于富里哀级数分析和处理的函数。不过,通过创建M文件 函数,可容易加上这些函数。在这一节,将介绍《精通 MATLAB工具箱》中富里哀级数函 数。在介绍之前,首先定义实周期信号f(t)的富里哀级数表示形式 给出富里哀级数的复指数形式为: f(t)=∑Fnea 式中的富里哀级数的系数是 Fn=中+f( t)e wo dt
>>W=Ws*(0 : N/2)/N 它建立了连续频率轴，该轴起始于 0，终止于奈魁斯特（Nyquist）频率 Ws/2， >>Fa=2./(3+j*w); % evaluate analytical Fourier transform >>plot(W, abs(Fa), W, abs(Fp), ‘ + ‘ ) % generate plot, ‘ + ‘ mark fft results >>xlabel(‘ Frequency, Rad/s ‘),ylabel(‘ |F(w)| ‘) MATLAB 提供了大量的完成一般信号处理任务的函数。它们列于表 14.1： 表 14.1 信号处理函数 conv 卷积 conv2 2 维卷积 fft 快速富里哀变换 fft2 2 维快速富里哀变换 ifft 逆快速富里哀变换 ifft2 2 维逆快速富里哀变换 filter 离散时间滤波器 filter2 2 维离散时间滤波器 abs 幅值 angle 四个象限的相角 unwrap 在 360°边界清除相角突变 fftshift 把 FFT 结果平移到负频率上 nextpow2 2 的下一个较高幂次 14.2 富里哀级数 MATLAB 本身没有特别关于富里哀级数分析和处理的函数。不过，通过创建 M 文件 函数，可容易加上这些函数。在这一节，将介绍《精通 MATLAB 工具箱》中富里哀级数函 数。在介绍之前，首先定义实周期信号 f(t)的富里哀级数表示形式。 给出富里哀级数的复指数形式为： f t F en jnw t n ( ) = = 0 式中的富里哀级数的系数是： F T f t e dt n jnw t t t T = + − 1 0 0 0 ( )
且基频为oo=2π/T。式中T满足f(tT=f(t)e 给出富里哀级数的三角数形式为: f(t)=Ao+E(An cos(noo t)+ B, sin(noot)i 式中的富里哀级数的系数是 It*lo f(t)dt T An==t+To f(t)cos( noot)dt Bn=止f(t)sn(not)dt 且基频为oo=2π/To。式中T满足f(tT0=ft 表14.2 精通 MATLAB的富里哀级数函数 fsderiv(Kn, Wo) 富里哀级数的微分 fseval(Kn, t, Wo) 计算富里哀级数 fsfind( fname,T, N) 寻找时间函数的富里哀级数的系数 [An, Bn, Ao]=fsform(Kn) 富里哀级数不同形式之间的转换 Kn=fsform(An, Bn, Ao) sharm(Kn, 1) 提取特殊的富里哀级数的谐波 fsms(Kn) 计算信号的均方根值 fsresize(Kn, N) 重新调整富里哀级数的系数向量的大小 fsresp(Num, Den, Un, Wo) 线性系统对输入富里哀级数Un的富里哀级数 fsround (Kn) 设置次要的富里哀级数的系数为0 fswindow(N, type 产生一个窗口的函数,使吉布(Gibb)现象极 上述两种形式中,一般复指数的富里哀级数更容易进行解析上处理。而三角富里哀级 数则提供了更多的直观理解,更容易将正弦和余弦波形可视化。由于这个原因,《精通 MATLAB工具箱》中的富里哀级数函数一般假定为复指数形式。不过,该工具箱提供了这 两种富里哀级数形式之间转换的函数。表142总结了《精通MA∏LAB工具箱》中的富里哀 级数函数 因为有无穷多个谐波或富里哀级数的系数,有必要对富里哀级数截尾,只考虑有限个 谐波。如果考虑N个谐波, MATLAB中富里哀级数表示成一个长度为2N+1的行向量,向
且基频为 0 = 2 0 / T 。式中 T0 满足 f(t+ T0)=f(t)。 给出富里哀级数的三角数形式为： f t A A n t B n t n n n ( ) = + { cos( ) + sin( )} = 0 0 0 1 式中的富里哀级数的系数是： A T f t dt t t T 0 0 1 = + 0 ( ) A T f t n t dt n t t T = + 2 0 0 0 ( ) cos( ) B T f t n t dt n t t T = + 2 0 0 0 ( )sin( ) 且基频为 0 = 2 0 / T 。式中 T0 满足 f(t+ T0)=f(t)。 表 14.2 精通 MATLAB 的富里哀级数函数 fsderiv(Kn,Wo) 富里哀级数的微分 fseval(Kn,t,Wo) 计算富里哀级数 fsfind(‘ fname ‘,T,N) 寻找时间函数的富里哀级数的系数 [An,Bn,Ao]=fsform(Kn) Kn=fsform(An,Bn,Ao) 富里哀级数不同形式之间的转换 fsharm(Kn,i) 提取特殊的富里哀级数的谐波 fsmsv(Kn) 计算信号的均方根值 fsresize(Kn,N) 重新调整富里哀级数的系数向量的大小 fsresp(Num,Den,Un,Wo) 线性系统对输入富里哀级数 Un 的富里哀级数 响应 fsround(Kn) 设置次要的富里哀级数的系数为 0 fswindow(N, ‘ type ‘) fswindow(Kn, ‘ type ‘) 产生一个窗口的函数，使吉布（Gibb）现象极 小 上述两种形式中，一般复指数的富里哀级数更容易进行解析上处理。而三角富里哀级 数则提供了更多的直观理解，更容易将正弦和余弦波形可视化。由于这个原因，《精通 MATLAB 工具箱》中的富里哀级数函数一般假定为复指数形式。不过，该工具箱提供了这 两种富里哀级数形式之间转换的函数。表 14.2 总结了《精通 MATLAB 工具箱》中的富里哀 级数函数。 因为有无穷多个谐波或富里哀级数的系数，有必要对富里哀级数截尾，只考虑有限个 谐波。如果考虑 N 个谐波，MATLAB 中富里哀级数表示成一个长度为 2N+1 的行向量，向
量中的元素为复指数的富里哀级数的系数。该向量包含以升序排列的富里哀级数系数,即 F=「FNF FI FoFL 为了说明这些函数的使用,考虑寻找锯齿信号的富里哀级数表示(见图142)。 0.2 图142锯齿信号 虽然很容易找到这个信号的富里哀级数的系数,但可用函数 fsfind近似求解这些系数 function[Fn, nwo, f, t=fsfind(fun, T,N, P) FSFIND Find Fourier Series Approximation %o Fn=FSFIND(FUN, T, P)computes the Complex Exponential Fourier Series of a signal described by the function'FUN FuN is the character string name of a user created M-file function %o the unction is called as f=fUN(t) where t is a vector over the range o<=t<=l The FFT is used Choose sufficient harmonics to minimize aliasing T is the period of the function. N is the number of harmonics Fn is the vector of fS coefficients [ Fn, nWo, t, p]=FSFIND(FUN, T, P)returns the frequencies associated lo with fn in n wo and return values of the function fun in f evaluated at the points int t over the range O<=t<=T FSFIND(FUN, T, P) passes the data in P to function FUN as fFUN(t, p). This allows parameters to be passed to FUN Copyrigth(c)1996 by Prentice-Hall, Inc n=2*N linspace(0, T, n+1)
量中的元素为复指数的富里哀级数的系数。该向量包含以升序排列的富里哀级数系数，即： F = F−N [ F F F F F F ] -N+1 -1 0 1 N-1 N 为了说明这些函数的使用，考虑寻找锯齿信号的富里哀级数表示（见图 14.2）。 0 0.1 0.2 0.3 0.4 0.5 0 2 4 6 8 10 图 14.2 锯齿信号 虽然很容易找到这个信号的富里哀级数的系数，但可用函数 fsfind 近似求解这些系数。 function [Fn, nwo, f, t]=fsfind(fun, T, N, P) % FSFIND Find Fourier Series Approximation. % Fn=FSFIND(FUN, T, P) computes the Complex Exponential % Fourier Series of a signal described by the function ' FUN '. % FUN is the character string name of a user created M-file function. % the unction is called as f=FUN(t) where t is a vector over % the range 0<=t<=T. % % The FFT is used. Choose sufficient harmonics to minimize aliasing. % % T is the period of the function. N is the number of harmonics. % Fn is the vector of FS coefficients. % % [Fn, nWo, t, p]=FSFIND(FUN, T, P) returns the frequencies associated % with Fn in nWo and return values of the function FUN % in f evaluated at the points int t over the range 0<=t<=T. % % FSFIND(FUN, T, P) passes the data in P to function FUN as % f=FUN(t, p). This allows parameters to be passed to FUN. % Copyrigth (c) 1996 by Prentice-Hall,Inc. n=2*N; t=linspace(0, T, n+1);