matlab 拟合

多项式拟合
p=polyfit(x,y,n)
其中x,y是输入数据，对应自变量和应变量；n表示多项式拟合的阶数。
函数的结果p是多项式的系数；
Y=polyval(p,x1)
根据拟合求出的多项式系数p,求解x1变量对应的y值；

频谱分析

信号的分类

信号分类的方法有很多，这里只给出四种后面将会涉及的分类方式。

按数学关系：确定性信号和随即信号；

按能量特性：能量信号和功率信号；

按取值特征：连续信号和离散信号；

按重复特性，周期性信号和非周期性信号。

随机信号：不能用确定的数学关系来描述的，其值的变动服从统计规律，常用一些数学统计量来描述该类信号的特征。

周期信号：每隔一定时间T，周而复始且无始无终的信号。

以下将重点阐述能量信号和功率信号的概念。我们选择一个具体的物理系统来引出这两种信号的概念。已知阻值为R的电阻上的电压和电流分别为v(t)和i(t)，则通过电阻的瞬时功率为：p(t)=v²(t)/R=i²(t)R。方便起见，假设电阻R为1Ω。故在1Ω电阻上消耗的能量为：

                                                                      

其中f(t)为电压或者电流。

注意积分区间是从，如果积分值存在，则该信号为能量信号；如果该积分不存在，则考虑求该信号的平均功率，如下式：

                                                                  

若该值存在，则该信号即为功率信号。

几点备注：

ü  能量信号的平均功率为0，即E/T（E有限，T无限）为0；

ü  功率信号的能量为无穷，即PT(P有限，T无限)无穷；

ü  信号f(t)可以既非功率信号也非能量信号（如单位斜坡信号f(t)=t）E和P都不存在；

ü  信号f(t)不可能同时既是功率信号也是能量信号（容易理解）；

ü  一般来说，周期信号和随机信号都是功率信号；而非周期的确定信号是能量信号。

ü  除了具有无限能量和功率的信号外，非周期信号要么是能量信号要么是功率信号；

时域和频域

时域和频域是信号通常的两个表示方式。

所谓时域就是自变量是时间，而纵轴是信号的变化，描述了信号在不同时刻取值的函数。

频域的自变量是频率，纵轴是该频率的振幅或相位或功率或能量，描述了一个信号的频谱结构。

如何将信号从时域转换到频域？

一般通过傅里叶变换将信号从时间域转换到频率域。注意，并不是所有的信号都能满足傅里叶变换的条件，通常情况下周期信号满足傅里叶变换的条件。对于测量中经常碰到的随机信号，并不满足傅里叶变换的条件。但是通常会截取其中一段（即样本，用样本的统计值来估计整体信号的频谱特征），而截取的这一段可以进行傅里叶变换，进而可进行功率谱分析。

将信号从时间域转换到频率域后，信号即可表示成很多的谐波之和。

振幅谱：用频域做横坐标，每个谐波的振幅做纵坐标，绘出的图形称为振幅谱；

相位谱：用频域做横坐标，每个谐波的初始相位做纵坐标，绘出的图形称为相位谱。

功率谱和能量谱：

将(1)和(2)式进行傅里叶变换得到信号能量和功率在频域的表示方式，这里不介绍推导过程，只讲述这两个概念的意思。

能量谱密度，也称能量谱，表示信号能量在各个频率的分布情况，信号的总能量等于各个频率分量单独贡献出的能量的连续和。

功率谱密度：也称功率谱，反映了信号功率在各个频率的分布情况，信号的平均功率等于各个频率分量单独贡献出的功率之连续和。

随机信号的特征

随机信号是一类持续时间无限长，具有无限大能量的功率信号，它不满足傅里叶变换条件，而且也不存在解析表达式，因此就不能应用确定信号的频谱计算方法去分析随机信号的频谱。然而，虽然随机信号的频谱不存在，但其相关函数确实是确定的。如果随机信号时平稳的，那么其相关函数的傅里叶变换就是它的功率谱密度函数。

由于在测量中经常用到的是随机信号，这里简单给出描述随机信号的数学指标。这些指标通常包含期望、方差、相关函数、自协方差、互相关函数、互协方差等。

功率谱密度估计

本文以下将详细介绍功率谱估计的方法及其matlab实现。功率谱估计的方法可以分为经典谱估计法和现代谱估计法。

首先，给出一个结论：随机信号自相关函数的傅里叶变换为信号的功率谱密度，互相关函数和互功率谱密度也是一傅里叶变换对。

经典谱估计法概述：

经典谱估计法最常见的有BT(Blackman-Tukey)法和周期图法，它们都是通过傅里叶变换来实现的。BT法是间接法，它先由采样数据估计不同延迟的自相关函数，然后用不同的窗函数对自相关估值开窗加权，再对加权后的自相关估值做傅里叶变换，得到功率谱估计。周期图法是直接法，它是直接将采样数据（也可经窗函数加权）做傅里叶变换，再取其幅度平方而得到功率谱估计。两种方法的结果是一致的。

然而，以上两种方法求出的功率谱分辨率低，故出现了一些改进的方法。大体改进思路是：对数据进行分段，每段求取功率谱密度，继而求平均功率谱密度；数据分段，每段用窗函数进行加权，继而求平均。此外，分割的每段数据还可以有数据重叠。

以上就是经典谱估计方法的发展历程，以下依次用简单的matlab例子进行说明。经典谱估计法的常用matlab函数有periodogram和pwelch函数。

1、 直接法：

直接法又称周期图法，它是把随机序列x(n)的N个观测数据视为一能量有限的序列，直接计算x(n)的离散傅立叶变换，得X(k)，然后再取其幅值的平方，并除以N，作为序列x(n)真实功率谱的估计。
clear;
Fs=1000; %采样频率
n=0:1/Fs:1;

%产生含有噪声的序列
xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));
window=boxcar(length(xn)); %矩形窗
nfft=1024;
[Pxx,f]=periodogram(xn,window,nfft,Fs); %直接法
plot(f,10*log10(Pxx));  %10*log10(Pxx)  %将结果转化为分贝（dB）

为了进一步说明直接法的原理，这里也给出了不用periodogram函数的直接法部分程序：

pxx1=abs(fft(xn(1:256),Nsec).^2)/Nsec;   %根据psd的定义

pxx2=abs(fft(xn(256:512),Nsec).^2)/Nsec;

pxx3=abs(fft(xn(512:768),Nsec).^2)/Nsec;

pxx4=abs(fft(xn(768:1024),Nsec).^2)/Nsec;

pxx=10*log¹⁰((pxx1+pxx2+pxx3+pxx4)/4)  % 转化为分贝

 2、间接法：
间接法先由序列x(n)估计出自相关函数R(n)，然后对R(n)进行傅立叶变换，便得到x(n)的功率谱估计。
clear;
Fs=1000; %采样频率
n=0:1/Fs:1;
%产生含有噪声的序列
xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));
nfft=1024;
cxn=xcorr(xn,'unbiased'); %计算序列的自相关函数
CXk=fft(cxn,nfft);
Pxx=abs(CXk);
index=0:round(nfft/2-1);
k=index*Fs/nfft;
plot_Pxx=10*log10(Pxx(index+1));
plot(k,plot_Pxx);
3、改进的直接法：
对于直接法的功率谱估计，当数据长度N太大时，谱曲线起伏加剧，若N太小，谱的分辨率又不好，因此需要改进。

3.1、Bartlett法
Bartlett平均周期图的方法是将N点的有限长序列x(n)分段求周期图再平均。

clear；
Fs=1000;
n=0:1/Fs:1;
xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));
nfft=1024;
window=boxcar(length(n)); %矩形窗
noverlap=0; %数据无重叠
p=0.9; %置信概率
[Pxx,Pxxc]=pwelch(xn,nfft,Fs,window,noverlap,p);  %该方法采用的主要是矩形窗
index=0:round(nfft/2-1);
k=index*Fs/nfft;
plot_Pxx=10*log10(Pxx(index+1));
plot_Pxxc=10*log10(Pxxc(index+1));
figure(1)
plot(k,plot_Pxx);
pause;
figure(2)
plot(k,[plot_Pxx plot_Pxx-plot_Pxxc plot_Pxx+plot_Pxxc]);
3.2、Welch法
Welch法对Bartlett法进行了两方面的修正，一是选择适当的窗函数w(n)，并再周期图计算前直接加进去，加窗的优点是无论什么样的窗函数均可使谱估计非负。二是在分段时，可使各段之间有重叠，这样会使方差减小。

clear;
Fs=1000;
n=0:1/Fs:1;
xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));
nfft=1024;
window=boxcar(100); %矩形窗
window1=hamming(100); %海明窗
window2=blackman(100); %blackman窗
noverlap=20; %数据无重叠
range='half'; %频率间隔为[0 Fs/2]，只计算一半的频率
[Pxx,f]=pwelch(xn,window,noverlap,nfft,Fs,range);
[Pxx1,f]=pwelch(xn,window1,noverlap,nfft,Fs,range);
[Pxx2,f]=pwelch(xn,window2,noverlap,nfft,Fs,range);
plot_Pxx=10*log10(Pxx);
plot_Pxx1=10*log10(Pxx1);
plot_Pxx2=10*log10(Pxx2);
figure(1)
plot(f,plot_Pxx);
pause;
figure(2)
plot(f,plot_Pxx1);
pause;
figure(3)
plot(f,plot_Pxx2);

分辨率与频谱泄露

在周期图法中，对有限长0-(N-1)的随机数据序列可以看做无限长的随机数据序列经矩形窗开窗截断的结果，我们知道两个时间序列相乘的傅里叶变换等于两个时间序列傅里叶变换的褶积。如果信号的真正功率集中在一个窄的频带内，则该褶积运算将把这个窄带的功率扩展到邻近的频段，这种现象称之为“泄露现象”。弱信号的主瓣很容易被强信号泄露到邻近的副瓣所淹没或畸变。对于短观测数据序列，功率谱轨迹的分辨率是不高的。

结论

在经典谱估计中，总体存在着频率分辨率低、方差性能不好的问题，原因是谱估计时需要对数据进行加窗处理，用有限个数据或其自相关函数来估计无限个数据的功率谱。

GMT系列十 GMT格网化方法

在GMT的很多命令中都会涉及到格网文件，在上一节中已经较为详细的介绍了格网文件的相关内容。本章将详细介绍如何生成格网文件。

在我们所获取的数据里，通常都是(x,y,z)形式，而这些数据有的已经是按规则的格网排列好的，而大多数都是不规则分布的数据形式。对于按规则格网排列好的(x,y,z)数据，只需用xyz2grd命令即可生成格网文件；对于不规则分布或者随意分布的数据，首先需要将这些数据进行格网化，生成按规则格网排列的(x,y,z)数据，然后再转化为格网数据。这里主要讨论随意分布的数据如何格网化并生成格网数据。

数据格网化，常用有两种算法：nearest neighbor gridding和gridding with splines in tenson

1.       nearest neighbor gridding

该方法对应的GMT命令为nearneighbor，该命令的详细用法这里不做介绍，只讲述这种格网化算法。

nearneighbor命令将会指定数据范围和格网间隔，从而即可确定格网结点。此外，该命令还需要给出搜索半径，-Sradius选项。有了结点，有了搜索半径，即可搜索出该结点附近分布的点值了，根据这些点值，利用距离做权重，计算出结点处的值，即为格网点的值了。该算法要求在每个结点的搜索范围内只要有一个点，如果没有点，则该结点值被设为NaN值。

例：nearneighbor –R245/255/20/30 –I5m –S40k –Gship.nc –V ship.xyz

注：当数据较密时，适合用该方法。

2.       gridding with spline in tension

利用全局的数据进行格网化。具体说来就是，将所有的数据投影到一个surface上，surface上对应的格网点值就是格网值。因此，相比这两种方法，nearneighbor是利用局部范围的数据进行格网化，而该方法是用全局的数据进行格网化。

在利用全局的数据进行格网化之前，为了消除aliasing效应，需要用blockmean或blockmedian或blockmode命令对数据进行预处理。blockmean命令适用于较为平滑的数据；blockmedian适用于起伏较大的数据。预处理后的数据，即可用surface命令进行格网化并生成格网数据了。

例：blockmean –R245/255/20/30 –I5m –V ship.xyz > ship_5m.xyz

surface ship_5m.xyz –R245/255/20/30 –I5m –Gship.nc –V

注：相同的数据，surface和nearneighbor命令格网化的结果差别较大。

此外，GMT还提供了两个格网数据提取和信息查询的命令，分别是grdraster和grdinfo，这里就不在详细介绍了。

GMT系列九 GMT文件格式介绍

GMT中会用到几种不同格式的数据文件，这里将详细介绍这几类文件。

1.      table文件

1)       ascii table文件

ascii table文件是最常用的一种输入文件，该文件就是m行n列的形式，文件的第1,2列是x和y的坐标值。ascii table文件具体又可分为单段和多段文件。

ü   单段文件（只有一个文件头）

最为常见的文件形式。如果有文件头，需用-H选项标示，默认文件头只有一行，是以”>”字符开始的行（如果想更改这个默认字符，可用-m选项，如-mA表示以字母A开头的行为文件头），如果文件头的行数大于1，则需要用-Hnrecs设定行数，如-H5表明文件头为5行。

ü   多段文件（有多个文件头）

用-m选项标识输入的将是一个多段文件，并可用该选项更改头记录的标记。

2)       二进制table文件

头记录等同于所有域值为NaN的记录，也可分为单段和多段文件。

ü   单段文件

可用-H选项标示，并指定头记录的行数。

ü   多段文件

用-m选项标示，但后面不许加flags，因为所有的头记录都是NaN值。

2.      格网文件

GMT里存储2D格网值的文件通常为netCDF格式的文件。默认情况下，netCDF文件的前2列变量是Z变量，而x和y的坐标将会根据z变量的维数确定。格网的存储顺序是从上到下（ymax(north) to ymin(south)）、从左到右（xmin(west) to xmax(east)）。有两种方式来确定结点，一种结点是格网线的交叉点，另一种是格网的中心点。这两种确定格网的方法分别为gridline registration和pixel registration。

gridline registration

结点是格网线的交叉点，每个结点值代表以结点为中心，xinc和yinc矩形范围内的值；结点数为

nx=(x_max-x_min)/x_inc+1

ny=(y_max-y_min)/y_inc+1

pixel registration

结点位于每个格网的中心，其值代表每个格网的均值，结点数为：

nx=(x_max-x_min)/x_inc

ny=(y_max-y_min)/y_inc

该方法确定的结点数在x和y方向分别比上一种方法少一个。

3.      sun raster files

该类文件通常可用后缀.ras来标示，可用于定制一些特定的符号或logo等。在GMT参考手册的第七章中会见到该文件的例子。

4.      cpt文件

这里只详细介绍通用的cpt文件，对于用于绘制分类图的cpt文件，这里不做介绍。

对于用RGB方式生成的cpt文件，其格式为：

z₀       R_min  G_min B_min  z₁      R_max G_max B_max [A]   [;label]

…

z_n-2    R_min  G_min B_min  z_n-1   R_max G_max B_max [A]   [;label]

B        R_back G_back B_back

F        R_fore G_fore B_fore

N       R_nan  G_nan  B_nan

z是表示z变量的范围，即从z₀到z₁之间的值被设定为对应的颜色；R, G, B用来设定颜色；可选项A是用来表示当在绘制scale时，annotation放置的位置，A选项的值可以是Low, Upper, Both. 以;开始的label是在绘制scale时设置的label. B设置背景色；F设置前景色；N是那些值为NaN对应的颜色。