第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：
工作空间上移
