4.6  曲线拟合

在上一节，已经介绍了数据插值，它要求原始数据是精确的，或具有较小的误差。事实上，由于种种原因，实验或测量中所获得的数据总会有一定的误差。在这种情况下，如果强求构造的函数(曲线)通过各插值节点，显然是不合理的。为此，人们设想构造一个函数(曲线)y=g(x)去拟合f(x)，但它不必通过各插值节点，而只是使该曲线从这些插值节点中穿过，且使它在某种意义下最优。

MATLAB的曲线拟合是用常见的最小二乘原理，所构造的g(x)是一个次数小于拟合节点个数的多项式。

4.6.1  最小二乘原理及其曲线拟合算法

设测得离散的n+1个节点的数据如下：

构造一个如下的m次拟合多项式函数g(x)为 (m≤n)：

所谓曲线拟合的最小二乘原理，就是使上述拟合多项式在各数据点处的偏差的平方之和达最小。

上式中的均为已知值，而式中的系数为个未知数，故可以将其看做是的函数，即。于是我们可以把上述曲线拟合归结成对多元函数的求极值问题。为使取极小值，必须满足以下方程组：

经过简单的推导，可以得到一个m+1阶线性代数方程组Sa=t，其中S为m+1阶系数矩阵，t为右端项，而a为未知数向量，即欲求的m次拟合多项式的m+1个系数。这个方程组也称为正则方程组。至于正则方程组的具体推导，可参阅有关数值计算方法的教材。

4.6.2  曲线拟合的实现

在MATLAB中，可以用polyfit函数来求最小二乘拟合多项式的系数，另外可以用polyval函数按所得的多项式计算指定值。

polyfit函数的调用语法是：

[p,s]=polyfit(x,y,m)

输入参数x，y为测量而得的原始数据，为向量；m为欲拟合的多项式的次数。polyfit (x,y,m)将根据原始数据x、y得到一个m次拟合多项式P(x)的系数，该多项式能在最小二乘意义下最优地近似函数f(x)，即有p(xi)≈f(xi)≈yi。

返回的结果中p为m次拟合多项式的系数，而s中的数据则是一个结构数组，代入polyval函数后可以得到拟合多项式相关的误差估计。s最常用的写法可以是：p=polyfit(x,y,M)。

polyval的函数功能是按多项式的系数计算指定点所对应的函数值。

【例4-43】  曲线拟合示例。

本例首先在多项式的基础上加入随机噪声，产生测试数据，然后对测试数据进行数据曲线拟合：

>> clear

>> rand('state',0)

>> x=1:1:10;

>> y=-0.9*x.^2+10*x+20+rand(1,10).*5; %  产生测试数据

>> plot(x,y,'o')                  %  绘图并标出原始数据点

>> p=polyfit(x,y,2)

>> xi=1:0.5:10;

>> yi=polyval(p,xi);                     % 计算拟合的结果

>> hold on

>> plot(xi,yi);                   %  绘制拟合结果图

>> hold off

运行以上命令，得到的结果如图4-10所示。另外得到的多项式系数为：

p =

-0.8923    9.8067   23.6003

也就是说通过曲线拟合，得到了多项式。通过比较系数和观察图形，可以看出本次曲线拟合结果的精度是比较高的。

图4-10  曲线拟合

4.7  Fourier分析

傅立叶(Fourier)分析在信号处理领域有着广泛的应用，现实生活中大部分的信号都包含有多个不同的频率组件，这些信号组件频率会随着时间或快或慢的变化。傅立叶级数和傅立叶变换是用来分析周期或者非周期信号的频率特性的数学工具。从时间的角度来看，傅立叶分析包括连续时间和离散时间的傅立叶变换，总共有4种不同的傅立叶分析类型：连续时间的傅立叶级数、连续时间的傅立叶变换、离散时间的傅立叶级数、离散时间的傅立叶变换等。

频谱分析是在数据中识别频率组成的处理过程。对于离散数据，频谱分析的计算基础是离散傅立叶变换(DFT)。DFT将time-based或者space-based数据转换为frequency-based数据。

一个长度为n的向量x的DFT，也是一个长度为n的向量：

其中是n阶复数根：

在此表达式中，i表示虚数单位 。

DFT有一种快速算法FFT，称为快速傅立叶变换。FFT并不是与DFT不同的另一种变换，而是为了减少DFT运算次数的一种快速算法。它是对变换式进行一次分解，使其成为若干个小数点的组合，从而减少运算量。常用的FFT是以2为基数的，其长度用N表示，N为2的整数倍。

MATLAB中采用的就是FFT算法。MATLAB提供了函数fft和ifft等来进行傅立叶分析。

1．函数fft和ifft

函数fft和ifft对数据作一维快速傅立叶变换和傅立叶反变换，函数fft的调用语法有如下几种。

(1)Y=fft(X)：如果X是向量，则采用快速傅立叶变换算法作X的离散傅立叶变换；如果X是矩阵，则计算矩阵每一列的傅立叶变换。

(2)Y=fft(X,n)：用参数n限制X的长度，如果X的长度小于n，则用0补足；如果X的长度大于n，则去掉长出的部分。

(3)Y=fft(X,[ ],n)或Y=fft(X,n,dim)：在参数dim指定的维上进行操作。

函数ifft的用法和fft完全相同。

2．fft2和ifft2

函数fft2和ifft2对数据作二维快速傅立叶变换和傅立叶反变换。数据的二维傅立叶变换fft2(X)相当于fft(fft(X)’)’，即先对X的列做一维傅立叶变换，然后对变换结果的行做一维傅立叶变换。函数fft2的调用语法有如下几种。

(1)Y=fft2(X)：二维快速傅立叶变换。

(2)Y=fft2(X,MROWS,NCOLS)：通过截断或用0补足，使X成为MROWS*NCOLS的矩阵。

函数ifft2的用法和fft2完全相同。

3．fftshift和ifftshift

函数fftshift(Y)用于把傅立叶变换结果Y(频域数据)中的直流分量(频率为0处的值)移到中间位置：

(1)如果Y是向量，则交换Y的左右半边；

(2)如果Y是矩阵，则交换其一三象限和二四象限；

(3)如果Y是多维数组，则在数组的每一维交换其“半空间”。

函数ifftshift相当于把fftshift函数的操作逆转，用法相同。

【例4-44】  生成一个正弦衰减曲线，进行快速傅立叶变换，并画出幅值(amplitude)图、相位(phase)图、实部(real)图和虚部(image)图。

>> tp=0:2048;                             %  时域数据点数N

>> yt=sin(0.08*pi*tp).*exp(-tp/80);      %  生成正弦衰减函数

>> plot(tp,yt), axis([0,400,-1,1]),      %  绘正弦衰减曲线

>> t=0:800/2048:800;                       %  频域点数Nf

>> f=0:1.25:1000;

>> yf=fft(yt);                             %  快速傅立叶变换

>> ya=abs(yf(1:801));                      %  幅值

>> yp=angle(yf(1:801))*180/pi;            %  相位

>> yr=real(yf(1:801));                    %  实部

>> yi=imag(yf(1:801));                     %  虚部

>> figure

>> subplot(2,2,1)

>> plot(f,ya),axis([0,200,0,60])        %  绘制幅值曲线

>> title('幅值曲线')

>> subplot(2,2,2)

>> plot(f,yp),axis([0,200,-200,10])     %  绘制相位曲线

>> title('相位曲线')

>> subplot(2,2,3)

>> plot(f,yr),axis([0,200,-40,40])      %  绘制实部曲线

>> title('实部曲线')

>> subplot(2,2,4)

>> plot(f,yi),axis([0,200,-60,10])      %  绘制虚部曲线

>> title('虚部曲线')

本例首先生成正弦衰减函数yt，绘制的正弦衰减曲线如图4-11所示。然后对yt进行了快速傅立叶变换，结果如图4-12所示。

图4-11  正弦衰减曲线图

图4-12  傅立叶变换结果