第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)}  因为MATLAB不允许零下标，所以移动了一个下标值。  相应的逆变换为：  为了说明FFT的使用，考虑估计连续信号的富里哀变换的问题。  解析上，该富里哀变换为：  虽然在这种情况下，由于知道了富里哀变换的解析结果，再运用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;  图14.1 富里哀变换两种结果的比较 仅从F中取正频率分量，并且乘以采样间隔计算F(w)。 >>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)的富里哀级数表示形式。 给出富里哀级数的复指数形式为：  式中的富里哀级数的系数是：  且基频为。式中T0满足f(t+ T0)=f(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的行向量，向量中的元素为复指数的富里哀级数的系数。该向量包含以升序排列的富里哀级数系数，即： 为了说明这些函数的使用，考虑寻找锯齿信号的富里哀级数表示（见图14.2）。  图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 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 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  图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<alpha<1; alpha = 0.5 is default 'poi' = Poisson, 0<alpha<inf; alpha = 1 is default 'cau' = Cauchy, 1<alpha<inf; alpha = 1 is default 'gau' = Gaussian, 1<alpha<inf; alpha = 1 is default Reference: F. J. Harris,”On the Use of Windows for Harmonic Analysis with the Discrete Fourier Transform,” IEEE Proc. Vol. 66, no. 1, Jan. 1978, pp.51-83 >>Fnh=Fn.*fswindow(Fn, ‘ han ‘); % apply Hanning window >>f=fseval(Fnh, t, wo); % evaluate windowed FS >>plot(t, f) % plot results  图14.4 用汉宁窗口滤波后的富里哀级数 因为乘以窗口函数，改变了富里哀级数的系数，所以，现在图形（见图14.4）看起来好多了。 接下来，让我们把锯齿波加到传递函数为的线性系统中。 函数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<4, wo=1; end N=(length(un)-1)/2; % highest harmonic jnWo=sqrt(-1)*(-N : N)*wo; % frequencies of harmonics y=(polyval(num, jnWo) ./ polyval(den, jnWo)).*un; % output is Yn=H(jnWo)*Un 运用fsresp，我们得到： >>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。  图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)); elseif nargin==3 % trig from -> complex exp form nd=length(d); if nc~=nd; error(‘ A and B must be the same length ‘) end a=0.5*(c-sqrt(-1)*d); a=[conj(a(nc : -1 : 1)) co(1) a]; b=-nc : nc; else error(‘ Improper number of input arguments. ‘) end >>[An, Bn, Ao]=fsform(Fn) An = Columns 1 through 7 -0.2000 -0.2000 -0.2000 -0.2000 -0.2000 -0.2000 -0.2000 Columns 8 through 14 -0.2000 -0.2000 -0.2000 -0.2000 -0.2000 -0.2000 -0.2000 Columns 15 through 21 -0.2000 -0.2000 -0.2000 -0.2000 -0.2000 -0.2000 -0.2000 Columns 22 through 25 -0.2000 -0.2000 -0.2000 0 Bn = Columns 1 through 7 -3.1789 -1.5832 -1.0484 -0.7789 -0.6155 -0.5051 -0.4250 Columns 8 through 14 -0.3638 -0.3151 -0.2753 -0.2418 -0.2130 -0.1878 -0.1655 Columns 15 through 21 -0.1453 -0.1269 -0.1100 -0.0941 -0.0792 -0.0650 -0.0514 Columns 22 through 25 -0.0382 -0.0253 -0.0126 0 Ao = 4.9000 忽略平均值，这个锯齿应展示了奇对称特性，所有An系数应当为0。显然，它们不为0。发生这种情况的原因是因为我们使用了FFT来近似系数。此外，Ao平均值有一点儿误差，它的值应当是5.0。 14.3 小结 下面的两个表总结了本章所讨论的函数。 表14.3 信号处理函数  conv 卷积  conv2 2维卷积  fft 快速富里哀变换  fft2 2维快速富里哀变换  ifft 逆快速富里哀变换  ifft2 2维逆快速富里哀变换  filter 离散时间滤波器  filter2 2维离散时间滤波器  abs 幅值  angle 四个象限的相角  unwrap 在360°边界清除相角突变  fftshift 把FFT结果平移到负频率上  nextpow2 2的下一个较高幂次   表14.4 精通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）现象极小   第15章 低级文件I/O 对于大多数用户，MATLAB函数load和save为装载和存储数据提供了足够的工具。利用以扩展名为.mat结尾的文件名，load和save假定数据是以与平台无关的二进制格式保存，或者用称之为flat的简单的ASCII文件格式保存。当flat ASCII或.mat这两种格式还不够时，MATLAB提供了基于C语言实现的低级文件I/O函数。用这些低级文件I/O函数，MATLAB可以读写你所知道的任意文件格式。例如，知道电子文件或数据库程序所使用的格式，就可以把这些数据读进MATLAB的矩阵中去。类似地，也可以创建电子文件或数据库文件。 MATLAT中这种基本的低级文件I/O命令如下： 表15.1 MATLAB低级文件I/O函数  fclose： 关闭文件  feof： 测试文件结束  ferror： 查询文件I/O的错误状态  fgetl： 读文件的行，忽略回行符  fgets： 读文件的行，包括回行符  fpoen： 打开文件  fprintf： 把格式化数据写到文件或屏幕上  fread： 从文件中读二进制数据  frewind： 返回到文件开始  fscanf： 从文件中读格式化数据  fseek： 设置文件位置指示符  ftell： 获取文件位置指示符  fwrite： 把二进制数据写到文件里   除了这些函数外，所具有的MATLAB版本可能为一个或多个公用软件包提供读写文件的特定函数M文件。有关这些函数的进一步的信息，请使用在线帮助：>>help iofun 第16章 调试工具 在开发函数M文件过程中，不可避免地出现错误，即故障。MATLAB提供了很多函数和方法，帮助调试函数。 在MATLAB表达式中，有两类错误：语法错误和运行错误。当MATLAB计算一个表达式的值或一个函数被编译到内存时会发现语法错误。一旦发现语法错误，MATLAB立即标志这些错误，并提供有关所遇到的错误类型，以及发生错误处M文件的行数。给定这些反馈信息，就很容易纠正这些错误。 而另一方面，即使MATLAB标志了运行错误，但找出错误一般比较困难。当发现运行错误时，MATLAB把控制权返回给命令窗口和MATLAB的工作空间。失去了对发生错误的函数空间的访问权，因此，用户不能询问函数工作空间中的内容排除问题。 根据作者的经验，当一些操作结果导致空矩阵或NaNs时，最容易发生运行错误。所有有关NaNs的操作都返回NaNs值。因此，如果有可能出现NaNs结果，则当出现NaNs时，最好运用逻辑函数isnan来执行一些缺省操作。因为空矩阵为零维，所以对空矩阵寻址常常导致错误。函数find表示了可产生空矩阵结果的一般情况。如果函数find的空矩阵输出用于索引其它数组，所返回的值也将是空的。这样，空矩阵具有传播性质。例如： >>x=pi*(1 : 4) % example data >>i=find(x>20) % use find function >>y=2*x(i) % propagate the empty matrix 很清楚，当希望y具有有限维数和值时，可能发生运行错误。当执行一个操作或使用可返回空结果的函数时，逻辑函数isempty有利于为空矩阵定义一个缺省值，这样避免运行错误。 有几种调试函数M文件的方法。对于简单的问题，可直接使用下列的方法组合： 1、去掉文件中所选择的行的分号，以便中间结果显示在命令窗口中。 2、在文件中加入显示感兴趣的变量的语句 3、把keyboard命令放在文件中所选择的地方，给键盘暂时控制权。这样，可以查询函数空间并按需要改变其值。 4、在M文件开始，在function语句前加上%，将函数M文件改变为脚本M文件。当MATLAB执行该脚本M文件时，该空间就是MATLAB工作空间。这样，当发生错误时可以询问。 当M文件大，递归调用或者多次嵌套（即调用其它M文件函数，被调用M文件函数又调用其它M文件函数，等等）时，用MATLAB的调试函数会更方便。与上述所列方法相反，这些调试函数不要求将有问题的M文件进行编辑。表16.1所给出的这些函数类似于其它高级编程语言中所提供的函数。有关进一步的信息，以及它们的使用实例，参阅《MATLAB用户指南》。 表16.1 MATLAB调试函数  dbclear： 取消断点  dbcont： 在断点后恢复运行  dbdown： 工作空间下移  dbquit： 退出调试模式  dbstack： 列出谁调用谁  dbstatus： 列出所用的断点  dbstep： 执行一行或多行  dbstop： 设置断点  dbtype： 列出带行号的M文件  dbup： 工作空间上移