第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>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 >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);
if nargin==3 ffeval(fun, t) ffeval(fun, t, P) Fn=fft(f(I: n)) Fn=[o conj (Fn(N: -1: 2)) Fn(1: N)0/n; nwo-2*pi/T(N: N 根据上述描述,必须先建立一个函数以计算周期内锯齿信号值。这样,我们有: function f=sawtooth(t, T) SAWTOOTH Sawtooth Waveform Generation SAWTOOTH(t, T) computes values of a sawtooth having a period T at the values in t Used in Fourier series examples f10*rem(t, D/T 现在求近似富里哀级数的系数。 >>T= 2: desired >>N=25: number of harmonics >>Fn=fsfind( sawtooth T, N, T); compute Fourier series coefficients 因为用FFT求近似系数,所以得到的系数并不精确。不过,如果有足够多的系数,特 别是如果时间函数本身,或者其一阶,或者二阶导数连续(锯齿信号不连续),误差就比较 小。如果需要更高的精度,可调用积分程序来计算每一个富里哀级数系数的积分关系式。后 种方法需要很长的时间,因为FFT一次要近似所有的系数。不过,尽管有误差,对于许 多应用,FTT的近似结果已足够准确了 用函数 fseval可实现富里哀级数的计算 function y=fseval(kn, t, wo) FSEVaL Fourier Series Function Evaluation FSEVAL(Kn, t, Wo)computes values of a real valued function given its complex exponential Fourier series coefficients Kn, at the o points given in t where the fundamental frequency is Wo rad/s K contains the Fourier coefficentsin ascending order Kn= N-N+10 if Wo is not given, Wo=l ia assumed
if nargin==3 f=feval(fun, t); else f=feval(fun, t, P); end Fn=fft(f(1 : n)); Fn=[0 conj(Fn(N : -1 : 2)) Fn(1 : N) 0]/n; nwo=2*pi/T*(-N : N); 根据上述描述，必须先建立一个函数以计算周期内锯齿信号值。这样，我们有： function f=sawtooth(t, T) % SAWTOOTH Sawtooth Waveform Generation. % SAWTOOTH(t, T) computes values of a sawtooth having a % period T at the values in t. % % Used in Fourier series examples. f=10*rem(t, T)/T; 现在求近似富里哀级数的系数。 >>T= .2; % desired period >>N=25; % number of harmonics >>Fn=fsfind(‘ sawtooth ‘, T, N, T); % compute Fourier series coefficients 因为用 FFT 求近似系数，所以得到的系数并不精确。不过，如果有足够多的系数，特 别是如果时间函数本身，或者其一阶，或者二阶导数连续（锯齿信号不连续），误差就比较 小。如果需要更高的精度，可调用积分程序来计算每一个富里哀级数系数的积分关系式。后 一种方法需要很长的时间，因为 FFT 一次要近似所有的系数。不过，尽管有误差，对于许 多应用，FTT 的近似结果已足够准确了。 用函数 fseval 可实现富里哀级数的计算。 function y=fseval(kn, t, wo) % FSEVAL Fourier Series Function Evaluation. % FSEVAL(Kn, t, Wo) computes values of a real valued function given % its complex exponential Fourier series coefficients Kn, at the % points given in t where the fundamental frequency is Wo rad/s. % K contains the Fourier coefficentsin ascending order: % Kn = [k k ... k ... k k] % -N -N+1 0 N-1 N % if Wo is not given, Wo=1 ia assumed
Note: this function creates a matri rows= length(t)and columns=(length(K )-1)/2 Copyrigth(c)1996 by Prentice-Hall, Inc =2. w0=1: end length(kn) if rem(nk-1, 2)I(nk==1) erro( Number of element in K must be odd and greater than 1) n=0.5*(nk-1); highest harmonic nwo=wo*(1: n); harmonic frequencies ko=kn(n+1); average value kn=kn(n+2: nk), positive frequency coefs y=ko+2*(real(exp(*t(: )*nwo)*kn)) 遵照所需的语法,我们得到: >>n=-N: n: harmonic index >>Fna=j*5 /(pi*n);% actual Fourier series coefficients >>Fna(N+1)=5; poke in average value >>t=linspace(0,. 4);% time points to evaluate functions >>w=2*pi/T; fundamental frequency >>f-fseval(Fn, t, wo); evaluate approximated Fourier series >>fa=fseval( Fna, t, wo) evaluate actual fourier series >>plot(t, f, t, fa) plot resul 0.2
% Note: this function creates a matrix of size: % rows = length(t) and columns = (length(K)-1)/2 % Copyrigth (c) 1996 by Prentice-Hall,Inc. if nargin==2, wo=1; end nk=length(kn); if rem(nk-1, 2) | (nk==1) erro(‘ Number of element in K must be odd and greater than 1 ‘) end n=0.5*(nk-1);% highest harmonic nwo=wo*(1 : n); % harmonic frequencies ko=kn(n+1); % average value kn=kn(n+2 : nk) ‘; % positive frequency coefs y=ko+2*(real(exp(j*t( : )*nwo)*kn))' ; 遵照所需的语法，我们得到： >>n=-N : N; % harmonic index >>Fna=j*5 ./ (pi*n); % actual Fourier series coefficients >>Fna(N+1)=5; % poke in average value >>t=linspace(0 , .4); % time points to evaluate functions >>wo=2*pi/T; % fundamental frequency >>f=fseval(Fn, t, wo); % evaluate approximated Fourier series >>fa=fseval(Fna, t, wo); % evaluate actual Fourier series >>plot(t, f, t, fa) % plot results for comparison 0 0.1 0.2 0.3 0.4 -2 0 2 4 6 8 10 12
图14.3富里哀级数的结果比较 注意上述两个富里哀级数给出了相似的结果。从图143中可知,在锯齿信号的不连续 点的周围出现了明显的吉布现象波纹。用函数 window,将一个窗口加于富里哀级数系数, 就能减少这种波纹 >>help fswindow FSWINDOW Generate Window Functions FSWINDOW(N, TYPE)creates a window vector of type TYPE having a length equal to the scale N FSWINDOW(X, TYPE)creates a window vector of type TYPE having a length and orientation the same as the vector X FSWINDOW(X, TYPE, alpha) porvides a parameter alpha as required for some window type FSWINDOW with no input arguments returns a string matrix whose i-th row is the i-th TYPE given below TYPE is a string designating the window type desired rec= Re lar or box 'tri'= Triangular or bartlet han'= Hann or hanning ham'= hamming ' bla= blackman common coefs 'blx'= blackman exact coefs rie=Rieman (sin(x)/x) tuk=Tukey, O>Fnh=Fn. 'fswindow(Fn, han); apply Hanning window >>f-fseval( Fnh, t, wo) evaluate windowed FS >>plot(t, f) plot results
图 14.3 富里哀级数的结果比较 注意上述两个富里哀级数给出了相似的结果。从图 14.3 中可知，在锯齿信号的不连续 点的周围出现了明显的吉布现象波纹。用函数 fswindow，将一个窗口加于富里哀级数系数， 就能减少这种波纹。 >>help fswindow FSWINDOW Generate Window Functions. FSWINDOW(N, TYPE) creates a window vector of type TYPE having a length equal to the scale N. FSWINDOW(X, TYPE) creates a window vector of type TYPE having a length and orientation the same as the vector X. FSWINDOW(X, TYPE, alpha) porvides a parameter alpha as required for some window types. FSWINDOW with no input arguments returns a string matrix whose i-th row is the i-th TYPE given below. TYPE is a string designating the window type desired: 'rec' = Rectangular or Boxcar 'tri' = Triangular or Bartlet 'han' = Hann or hanning 'ham' = Hamming 'bla' = blackman common coefs. 'blx' = Blackman exact coefs. 'rie' = Rieman {sin(x)/x} 'tuk' = Tukey, 0>Fnh=Fn.*fswindow(Fn, ‘ han ‘); % apply Hanning window >>f=fseval(Fnh, t, wo); % evaluate windowed FS >>plot(t, f) % plot results
4 图144用汉宁窗口滤波后的富里哀级数 因为乘以窗口函数,改变了富里哀级数的系数,所以,现在图形(见图144)看起来 好多了 接下来,让我们把锯齿波加到传递函数为H(s)= 的线性系统中 s+60 函数 fares专门解决此类问题。 function y=fsresp(num, den, un, wo) %o FSRESP Fourier Series Linear System Response FSRESP(N, D, Un, Wo)returns the comples exponential FS of the output of a linear system when the input is given by a FS N and d are the numerator and denominator coefficients espectively of the system transfer function Un is the complex exponential Fourier Series of the system input Wo is the fundamental frequency associated with the input Copyright(c)1996 by Prentice-Hall, Inc if nargin>Yn=fsresp(60,[1 601, Fn, wo); find response coefficients >>y=fseval(Yn, t, wo); evaluate output 结果见图145
0 0.1 0.2 0.3 0.4 0.5 0 2 4 6 8 10 图 14.4 用汉宁窗口滤波后的富里哀级数 因为乘以窗口函数，改变了富里哀级数的系数，所以，现在图形（见图 14.4）看起来 好多了。 接下来，让我们把锯齿波加到传递函数为 H s s ( ) = + 60 60 的线性系统中。 函数 fsresp 专门解决此类问题。 function y=fsresp(num, den, un, wo) % FSRESP Fourier Series Linear System Response % FSRESP(N, D, Un, Wo) returns the comples exponential FS of the % output of a linear system when the input is given by a FS % N and D are the numerator and denominator coefficients % respectively of the system transfer function. % Un is the complex exponential Fourier Series of the system input. % Wo is the fundamental frequency associated with the input. % Copyright (c) 1996 by Prentice-Hall, Inc. if nargin>Yn=fsresp(60, [1 60], Fn, wo); % find response coefficients >>y=fseval(Yn, t, wo); % evaluate output >>plot(t, f, f, y) % plot system input and output 结果见图 14.5
0 0.2 04 图145线性系统的输入输出曲线 作为最后一个例子,运用 fsforn把锯齿信号的富里哀级数转换为三角形式 Inction [a, b, ao]=fsform(c, d, co) FSFORM Fourier Series Format Conversion KN=FSFORM(An, Bn, Ao) converts the trigonometric FS with An being the cSoine and Bn being the sine coefficients to the complex exponential FS with coefficients Kn Ao is the DC component and An, Bn and Ao are assumed to be real [Kni]=FSFORM(An, Bn, Ao) returns the index vector i that identifies the harmonic number of each element of Kn [An, Bn, Ao]=FSFORM(Kn) does the reverse format conversion Copyrigth(c)1996 by Prentice-Hall, Inc c=length(c) if nargin==1% complex exp-> trig form if rem(nc-1, 2)I(nc=1) error( Number of elements in K must bo odd and greater than 1) nn=(nc+3)/2; a=2*real(c(nn: nc)) b=-2*imag(c(nn: nc) ao=real(c(nn-D))
0 0.1 0.2 0.3 0.4 -2 0 2 4 6 8 10 12 图 14.5 线性系统的输入输出曲线 作为最后一个例子，运用 fsform 把锯齿信号的富里哀级数转换为三角形式。 function [a, b, ao]=fsform(c, d, co) % FSFORM Fourier Series Format Conversion % KN=FSFORM(An, Bn, Ao) converts the trigonometric FS with % An being the CSOINE and Bn being the SINE coefficients to % the complex exponential FS with coefficients Kn. % Ao is the DC component and An,Bn and Ao are assumed to be real. % % [Kni]=FSFORM(An, Bn, Ao) returns the index vector i that % identifies the harmonic number of each element of Kn. % % [An, Bn, Ao]=FSFORM(Kn) does the reverse format conversion. % Copyrigth (c) 1996 by Prentice-Hall,Inc. nc=length(c); if nargin==1 % complex exp -> trig form if rem(nc-1, 2)|(nc==1) error(‘ Number of elements in K must bo odd and greater than 1 ‘) end nn=(nc+3)/2; a=2*real(c(nn : nc)); b=-2*imag(c(nn : nc)); ao=real(c(nn-1));