MATLAB与现代信号处理.docx
《MATLAB与现代信号处理.docx》由会员分享，可在线阅读，更多相关《MATLAB与现代信号处理.docx（83页珍藏版）》请在冰豆网上搜索。
MATLAB与现代信号处理
第8章参数化建模3
8.1时间域建模3
8.1.1线性预报法（AR模型）3
8.1.2Prony法（ARMA模型）4
8.1.3Steiglity-McBride法（ARMA模型）4
8.2频率域建模5
8.2.1模拟滤波器的s域建模6
8.2.2数字滤波器的z域内建模6
第9章随机信号分析8
9.1随机信号的数字特征8
9.1.1均值、均方值、方差8
9.1.2离散随机信号9
9.1.3估计9
9.2相关函数和协方差10
9.2.1相关函数和自协方差10
9.2.2互相关函数和互协方差11
9.2.3计算相关函数和协方差的MATLAB函数11
9.3功率谱估计16
9.3.1功率谱密度16
9.3.2周期图法16
9.3.3多窗口法26
9.3.4最大熵法（Maxmumentropymethod，MEM法）28
9.3.5多信号分类法30
9.4传递函数估计31
9.5相干函数33
9.6运用功率谱提取地球自由振荡信息35
9.6.1资料35
9.6.2方法36
9.6.3结果37
第10章数字信号处理的几个前沿课题42
10.1时谱（倒谱）分析42
10.2地震观测系统的仿真和地面运动的恢复44
10.3小波分析举例57
10.3.1小波变换与突变58
10.3.2资料62
10.3.3本尼奥夫应变资料的小波分析63
11900~2001年中国大陆地震资料的小波分析63
21500~2001年华北地区地震资料的小波分析64
第8章参数化建模
由前几章的讨论可见，无论是模拟滤波器还是数字滤波器，均是根据输入信号和系统的传递函数求得系统的输出信号。
反过来，能否通过系统的输入信号和输出信号，或滤波器的频率响应或脉冲响应求得系统的传递函数呢？
这在一定的假设下是可以得到的。
这种技术涉及本章要讲的参数化建模。
所谓参数化建模就是根据未知系统的某些信息（如脉冲响应、频率响应或输入输出序列）建立该系统的有理传递函数模型。
这项技术用于求一个信号、系统或过程的模型参数，并广泛应用于语言分析、数据压缩、高分辨谱估计、信号处理等领域。
参数化建模分为时间域建模和频率域建模两类，本章将分别介绍。
8.1时间域建模
时间域建模就是给定系统时间域内的信息，如脉冲响应或系统输入和输出时间序列，求数字滤波器有理传递函数分子和分母多项式系数。
MATLAB中提供了三种时间域建模函数，下面分别叙述。
8.1.1线性预报法（AR模型）
如果一个信号的各个采样值不是独立的，其任一时刻k的采样值x（k）是它过去的n个采样值及一个白噪声序列k时刻的值u（k）线性组合而成（线性预报）,即
x（k）=-a
（2）x（k-1）-a（3）x（k-2）…-a（n）x（k-n-1）-a（n+1）x（k-n）+u（k）（8-1）
如果知道参数个数n，根据这个信号序列可以得到这n个系数的值。
数字信号处理中称这个信号x符合n阶自回归（automaticregression,AR）线性模型，它可用一个全极点IIR滤波器脉冲响应来逼近，全极点滤波器的传递函数具有如下形式：
（8-2）
MATLAB信号处理工具箱函数lpc利用自回归（AR）模型的相关法来求出AR模型系数a，也是全极点滤波器模型。
其调用格式为：
[a,g]=lpc（x,n）
其中，x为信号，是一个实时间序列；n为AR模型阶次；a为AR模型系数，a=[1,a
（2）,…,a（n+1）;g为AR模型的增益。
这就是说，要建立起AR模型，必须知道其阶次。
函数lpc的算法可简要地分为几个步骤：
（1）调用函数xcorr求信号x的相关估计。
（2）调用函数levinson实现Levinson-Durbin递推算法，求出系数a。
函数levinson的调用格式为：
a=levinson（r,n）
其中，r为x的自相关序列；n为AR模型阶次。
【例8-1】设信号x是一个带白噪声的4阶IIR滤波器的脉冲响应，IIR滤波器的传递函数无零点，其分母多项式系数为a=[10.10.10.10.1]，用线性预报法建立其AR模型。
%Samp8_1
randn（'state',0）;%设置随机函数的状态
a=[10.10.10.10.1];%滤波器分母多项式系数
x=impz（1,a,10）+randn（10,1）/20;%求得带噪声的滤波器脉冲响应
%采用第一种方法
r=xcorr（x）;%求信号x的自相关r，参考第9章的该函数用法。
r（1:
length（x）-1）=[];%该序列的前面部分（相当于x的长度-1）受到边界的影响，扣除
aa=levinson（r,4）%用Levison-Durbin递推算法采用自相关序列求出系数aa
%采用第二种方法
a1=lpc（x,4）%直接调用lpc（LinearPredictorCoefficients）函数求得结果
该程序的运行结果为：
aa=1.00000.18490.12790.11140.1839
a1=1.00000.18490.12790.11140.1839
程序中两种方法求得的结果完全相同。
若信号不包含随机噪声，用上面的程序求得信号AR模型和所采用全极点滤波器模型完全相同。
即使加入随机噪声，得到结果也与原模型参数a接近。
这说明，如果一个信号x是IIR滤波器的脉冲响应，可以利用函数lpc建立的向量a作为滤波器参数化模型。
8.1.2Prony法（ARMA模型）
如果已知某滤波器（或系统）的脉冲响应序列，可以采用Prony方法求得该滤波器（IIR）传递函数的系数。
MATLAB信号处理工具箱提供了Prony法建模的函数prony,其调用格式为：
[b,a]=prony（h,nb,na）
式中，h为时间域脉冲响应序列；nb,na分别为滤波器分子和分母多项式阶数；b,a分别为滤波器分子和分母多项式系数向量。
这就是说，要想建立起ARMA模型，必须知道传递函数分子和分母多项式的最高阶次。
滤波器传递函数表示为：
（8-3）
Prony法是用一个IIR滤波器脉冲响应逼近信号h。
【例8-2】建立一个4阶Butterworth数字滤波器，归一化截止频率为0.2，求其脉冲响应，再用脉冲响应系数和prony函数求滤波器模型,并与原始模型进行比较。
%Samp8_2
[b,a]=butter（4,0.2）%设计4阶截止频率为0.2（归一化）的Butterworth滤波器
h=impz（b,a,26）;%求26个点组成的滤波器的脉冲响应
[bb,aa]=prony（h,4,4）%运用脉冲响应h建立4阶系统的分子和分母多项式系数
运行结果为：
b=0.00480.01930.02890.01930.0048
a=1.0000-2.36952.3140-1.05470.1874
bb=0.00480.01930.02890.01930.0048
aa=1.0000-2.36952.3140-1.05470.1874
程序运行结果表明，用函数prony建模和原始滤波器模型完全相同，显然函数prony所建立的新滤波器的脉冲响应和信号h是完全吻合的。
因此，若已知一个系统的脉冲响应和该系统分子分母多项式阶数，可以运用该方法求得该系统的传递函数。
8.1.3Steiglity-McBride法（ARMA模型）
已知某滤波器（或系统）的脉冲响应序列建立滤波器传递函数的另一种方法是Steiglity-McBride法。
这种方法利用Steiglity-McBride迭代建模，使系统b（z）/a（z）的脉冲响应x’和输出信号x的均方差最小，即
（8-4）
MATLAB信号处理工具箱函数stmcb实现这种方法建模，可用于滤波器设计和系统辨识。
其调用格式为：
[b,a]=stmcb（x,nb,na[,niter[,ai]]）
式中,x为系统的脉冲响应；nb,na分别为系统传递函数b（z）/a（z）分子和分母多项式的阶次；niter为迭代计算次数；缺省值为5；ai为分母系数的初步估计，可缺省；b,a分别为系统（IIR型滤波器）分子和分母多项式系数向量。
系统（IIR滤波器）传递函数具有与（8-3）式相同的形式。
函数的另一种调用形式适用于系统的输入和输出信号作为函数输入：
[b,a]=stmcb（x,u,nb,na[,niter[,ai]]）
式中，x和u分别为系统（滤波器）输出和输入信号，其余参数同上。
函数的这种用法可以根据滤波器的输入和输出（而不是前面的脉冲响应）求得滤波器传递函数。
【例8-3】建立一个6阶Butterworth滤波器，截止频率为0.2，求其脉冲响应，并用脉冲响应数据和函数stmcb建立新滤波器模型，比较两个滤波器的频率特性。
%Samp8_3
[b,a]=butter（6,0.2）%产生6阶截止频率为0.2的Butterworth滤波器
h=filter（b,a,[1zeros（1,100）]）;%采用输入为脉冲函数的方法求得系统的脉冲响应
[bb,aa]=stmcb（h,6,6）%采用stmcb方法求得系统传递函数多项式系数
运行结果为：
b=0.00030.00200.00510.00680.00510.00200.0003
a=1.0000-3.57945.6587-4.96542.5295-0.70530.0838
bb=0.00030.00200.00510.00680.00510.00200.0003
aa=1.0000-3.57945.6587-4.96542.5295-0.70530.0838
结果表明，函数stmcb自滤波器脉冲响应建立的模型和原滤波器模型完全相同。
为了比较lpc，prony和stmcb这三个函数,下例给出了三种方法建模的误差
【例8-4】采用例8-1中的滤波器，用阶数3模拟比较函数lpc、prony和stmcb三种建模方法的误差。
%Samp8_4
randn（'seed',0）;%设置随机数种子
x=impz（1,[10.10.10.1],10）+randn（10,1）/20;%
a1=lpc（x,3）;%运用lpc函数建立传递函数系数
[b2,a2]=prony（x,3,3）;%运用prony函数建立传递函数系数
[b3,a3]=stmcb（x,3,3）;%运用stmcb函数建立传递函数系数
%计算误差
disp（'函数LPC输出:
'）
err1=x-impz（1,a1,10）;%x与lpc函数确定模型脉冲响应的每个数的误差
sumerr1=sum（err1.^2）%x与lpc函数确定模型脉冲响应的总误差
disp（'函数PRONY输出:
'）
err2=x-impz（b2,a2,10）;%x与prony函数确定模型脉冲响应的每个数误差
sumerr2=sum（err2.^2）%x与prony函数确定模型脉冲响应的总误差
disp（'函数STMCB输出:
'）
err3=x-impz（b3,a3,10）;%x与stmcb函数确定模型脉冲响应的每个数误差
sumerr3=sum（err3.^2）%x与stmcbc函数确定模型脉冲响应的总误差
程序输出结果为：
函数LPC输出:
sumerr1=0.0219;函数PRONY输出:
sumerr2=0.0147;
函数STMCB输出:
sumerr3=0.0039
该例表明：
在相同的阶数条件下，函数stmcb精度最好，prony次之，lpc较差。
8.2频率域建模
前面一节我们根据系统（滤波器）的时间域信号（脉冲响应或输入、输出信号）建立系统传递函数。
但如果已知某系统（滤波器）的频率响应，能否建立起该系统（滤波器）的模型呢？
这就涉及到频率域建模。
频率域建模就是由滤波器（系统）的频率响应建立其模型。
根据频率域类型分为面向模拟滤波器的s域内建模和面向数字滤波器的z域内建模。
8.2.1模拟滤波器的s域建模
模拟滤波器传递函数的一般形式为：
（8-5）
MATLAB信号处理工具箱函数invfreqs用于由给定的频率响应数据求形如上式的模拟滤波器传递函数。
调用格式为：
[b,a]=invfreqs（h,w,nb,na[,wt[niter]]）
式中，h为复频率响应向量；w为对应的频率向量；nb,na分别为模拟滤波器分子和分母多项式阶次；wt为与w同维数的权系数向量，调整对频率的拟合误差；niter为迭代次数。
若已知复频率响应的幅值向量mag和相位向量phase，在MATLAB中要采用h=mag.*exp（j*phase）将其复合为复数形式。
【例8-5】已知滤波器原始传递函数中b=[1232],a=[12321],试求其频率响应，并利用频率响应数据求原始滤波器的传递函数。
%Samp8_5
disp（'滤波器初始模型:
'）;
b=[1232]
a=[12321]
[h,w]=freqs（b,a,64）;%计算64个点的模拟滤波器频率响应
disp（'INVFREQS函数得到的模型:
'）
[bb,aa]=invfreqs（h,w,3,4）%运用invfreqs函数求滤波器系数
程序的运行结果为：
滤波器初始模型:
b=1232
a=12321
INVFREQS函数得到的模型:
bb=1.00002.00003.00002.0000
aa=1.00002.00003.00002.00001.0000
运行结果表明，由函数invfreqs从频率响应数据辨识的滤波器模型和原始模型完全相同。
8.2.2数字滤波器的z域内建模
数字滤波器传递函数的一般形式为:
（8-6）
MATLAB信号处理工具箱函数invfreqz用于由给定的频率响应数据求形如上式的数字滤波器的传递函数。
调用格式为：
[b,a]=invfreqz（h,w,nb,na[,wt[,niter]]）
式中，h为复频率响应向量；w为对应的频率向量；nb,na分别为数字滤波器分子和分母多项式阶次；wt为与w同维数的权系数向量，调整对频率的拟合误差；niter为迭代次数。
函数invfreqs和invfreqz利用等误差方法由给定的频率响应数据来辨识最佳模型，当niter缺省时，函数利用下面的算法：
（8-7）
式中，
和
分别表示多项式a,b在频率w（k）处的Fourier变换，
为权系数向量，n为频率点数（h和w的长度）。
这种算法不能保证辨识模型的稳定性。
当迭代次数niter确定后，函数用输出误差迭代算法进行运算：
（8-8）
这样可保证辨识模型的稳定性，且提高辨识模型的精度。
【例8-6】已知原始滤波器为4阶Butterworth低通数字滤波器，归一化截止频率为0.4。
用其频率响应数据产生一个三阶滤波器。
试比较采用等误差和输出误差迭代算法的结果。
%Samp8_6
[b,a]=butter（4,0.4）;%产生4阶Butterworth滤波器
[h,w]=freqz（b,a,64）;%计算频率响应
[bb,aa]=invfreqz（h,w,3,3）;%采用invfreqz求解滤波器的系数
wt=ones（size（w））;%产生权向量
niter=30;%迭代次数
[bbb,aaa]=invfreqz（h,w,3,3,wt,niter）;%采用权向量和迭代次数求滤波器系数
[h1,w1]=freqz（bb,aa,w）;%计算模型1的频率响应
[h2,w2]=freqz（bbb,aaa,w）;%计算模型2的频率响应
plot（w/pi,abs（h）,'k',w1/pi,abs（h1）,'k:
',w2/pi,abs（h2）,'b.'）;%绘制频率响应
xlabel（'归一化频率'）;ylabel（'振幅'）
legend（'理想模型','未采用权向量和迭代次数','采用权向量和迭代次数'）;
gridon
sumerr1=sum（abs（h-h1）.^2）%输出模型1误差
sumerr2=sum（abs（h-h2）.^2）%输出模型2误差
图8-1例8-6的两种算法估计模型的幅频特性的比较
程序命令窗口的输出为：
sumerr1=0.0200；sumerr2=0.0096。
程序运行的图形显示为图8-1。
程序运行结果表明，输出误差迭代法（预估算法）具有较好的拟合精度。
第八章习题
1.用MATLAB函数产生一个5阶低通ChebyshevI型数字滤波器（归一化截止频率为0.4）的脉冲响应序列，再用Prony法建模并与原滤波器模型进行比较。
2.用MATLAB函数产生一个5阶低通Butterworth数字滤波器（归一化截止频率为0.4）并计算其脉冲响应序列，再用Steiglitz-McBride法建模并与原滤波器模型进行比较。
3.用MATLAB函数产生一个低通ChebyshevI型高通数字滤波器，设计指标为：
通带截止频率Fc=2.5kHz,通带波纹为1dB，阻带衰减大于20dB。
试求
（1）高通数字滤波器的系统函数；
（2）滤波器的频率响应；（3）由频率响应数据辨识滤波器模型并与原滤波器模型进行比较。
第9章随机信号分析
随机信号和确定信号是两类性质完全不同的信号，对它们的描述、分析和处理方法也不相同。
随机信号是一种不能用确定数学关系式来描述的，无法预测未来某时刻精确值的信号，也无法用实验的方法重复再现。
随机信号分为平稳和非平稳两类。
平稳随机信号又分为各态历经和非各态历经。
本章所讨论的随机信号是平稳的且是各态历经的。
在研究无限长信号时，总是取某段有限长信号作分析。
这一有限长信号称为一个样本（或称子集），而无限长信号x（t）称为随机信号总体（或称集）。
各态历经的平稳随机过程中的一个样本的时间均值和集的平均值相等。
因此一个样本的统计特征代表了随机信号总体，这使得研究大大简化。
工程上的随机信号一般均按各态历经平稳随机过程来处理。
仅在离散时间点上给出定义的随机信号称为离散时间随机信号，即随机信号序列。
随机信号序列可以是连续随机信号的采样结果，也可以是自然界里实际存在的物理现象，即它们本身就是离散的。
平稳随机过程在时间上是无始无终的，即其能量是无限的，本身的Fourier变换也是不存在的；但功率是有限的。
通常用功率谱密度来描述随机信号的频域特征，这是一个统计平均的频谱特性。
平稳随机过程统计特征的计算要求信号x（n）无限长，而实际上这是不可能的，只能用一个样本，即有限长序列来计算。
因此得到的计算值不是随机信号真正的统计值，而仅仅是一种估计。
本章首先介绍随机信号的数字特征，旨在使大家熟悉描述随机信号的常用特征量。
然后介绍描述信号之间关系的相关函数和协方差。
这些是数字信号时间域内的描述。
在频率域内，本章介绍功率谱及其估计方法，并给出了功率谱在传递函数估计方面的应用。
最后介绍描述频率域信号之间关系的函数---相干函数。
9.1随机信号的数字特征
9.1.1均值、均方值、方差
若连续随机信号x（t）是各态历经的，则随机信号x（t）均值可表示为：
（9-1）
均值描述了随机信号的静态（直流）分量，它不随时间而变化。
随机信号x（t）的均方值表达式为：
（9-2）
均方值
表示了信号的强度或功率。
随机信号x（t）的均方根值表示为：
（9-3）
也是信号能量或强度的一种描述。
随机信号x（t）的方差表达式为：
（9-4）
是信号幅值相对于均值分散程度的一种表示，也是信号纯波动（交流）分量大小的反映。
随机信号x（t）的均方差（标准差）可表示为：
（9-5）
其意义与
相同。
9.1.2离散随机信号
若x（n）是离散的各态历经的平稳随机信号序列。
类似连续随机信号，其数字特征可由下面式子来表示：
均值：
（9-6）
均方值：
（9-7）
方差：
（9-8）
式中，
为均方根值，
为均方差值。
9.1.3估计
以上计算都是对无限长信号而言，而工程实际中所取得信号是有限长的，计算中均无法取
或
。
对于有限长模拟随机信号，计算均值时将（9-1）式改写为：
（9-9）
式中，均值
仅仅是对
的估计。
当T足够长时，均值估计
能精确逼近真实均值
。
对于周期信号，常取T为一个周期，估计均值
就能完全代表真实均值
。
对于有限长随机信号序列，计算其均值估计可将（9-6）式改写为：
（9-10）
当序列长度足够长时，均值估计
能精确逼近真实均值
。
类似地，可以写出均方值和方差估计表达式。
在MATLAB工具箱中，没有专门函数来计算均值、均方值和方差。
但随机信号的统计数字特征值计算都可以通过MATLAB编程实现。
在数值计算中，常将连续随机信号离散化，当作随机序列来处理。
【例9-1】计算以长度N=100000的正态分布高斯随机信号的均值、均方值、均方根值、方差和均方差。
%Samp9_1
N=100000;%数据个数
randn（'state',0）;%设置产生随机数的状态
y=randn（1,N）;%产生一个随机序列
disp（'平均值:
'）;
yMean=mean（y）%计算随机序列的均值
disp（'均方值:
'）;
y2p=y*y'/N%计算其均方值,这里利用了矩阵相乘的算法
disp（'均方根:
'）;
ysq=sqrt（y2p）%计算其均方根值
disp（'标准差:
'）;
ystd=std（y,1）%相当于ystd=sqrt（sum（（y-yMean）.^2）/（N-1））
disp（'方差:
'）;
yd=ystd.*ystd
该程序的运行结果为：
平均值：
yMean=0.0032;均方值：
y2p=1.0073
均方根：
ysq=1.0037
标准差：
ystd=1.0037
方差：
yd=1.0073
注意，函数s=std（x,flag）计算标准差。
x为向量或矩阵；s为标准差；flag为控制符，用来控制标准算法。
当flag=0（或缺省）时，按下式计算无偏标准差：
（9-11）
当flag=1时，按下式计算有偏标准差：
（9-12）
9.2相关函数和协方差
9.2.1相关函数和自协方差
对于随机信号x（t），自相关函数为：
（9-13）
式中，
为时移。
若去掉x（t）的均值部分，则相应的自相关函数称为自协方差，即：
（9-14）
（9-15）
对于离散随机信号序列，x（n）的自相关函数和自协方差为：
（9-16）
（9-17）
式中，m为延迟。
9.2.2互相关函数和互协方差
对于两个不同随机信号x（t）,y（t），互相关函数为：
（9-18）
式中，
为时移。
互协方差为：
（9-19）
对于互协方差有：
（9-20）
对于离散随机序列x（n）和y（n），互相关函数和互协方差为：
（9-21）
升级会员 