matlab在信号与图像处理中的应用第7章new

1、0 习题7第7章 随机信号处理与谱分析信号从广义上可分为两大类：确定性信号和随机信号，上一章内容所针对的信号形式就为确定性信号，本章将学习随机信号的处理方法，其中谱分析是一个非常重要的部分。这里介绍的谱分析包括功率谱估计及高阶谱估计两大部分，各自又有很多种估计方法。由于随机信号的处理方法是基于统计特性的分析，因此本章将首先给出描述统计特性的一些特征量，之后介绍功率谱估计方法，包括经典谱估计法、改进的直接法、AR模型法等，接下来将介绍高阶谱估计方法，包括非参数化方法及参数模型法，最后给出一个综合实例。学习目标r 掌握随机信号的统计描述和统计特性r 掌握经典功率谱估计方法及改进方法r 了解AR模型

2、功率谱估计方法r 理解高阶累积量和高阶谱定义r 了解高阶谱估计方法练习案例r 估计给定随机序列的均值和方差。r 产生一个相关正态随机序列，并计算它的自相关函数。r 产生两个被白噪声污染的随机序列，并分别计算二者的自协方差及互协方差。r 产生一频率为300Hz的谐波信号（分实数和复数两种情况），叠加高斯白噪声，利用周期图法估计该混合信号的功率谱。r 产生一频率为300Hz的谐波信号（复数形式），叠加高斯白噪声，利用Welch法估计该混合信号的功率谱。r 模拟一个非高斯ARMA随机过程，并估计它的3阶累积量，采用函数cum3x。r 分别利用直接法和间接法估计一个已知的非高斯信号的双谱。r 估计一已

3、知的非高斯信号ARMA模型参数，并利用此模型估计该信号的双谱。r 生成一混合信号，包含频率分别为300Hz和320Hz两个谐波信号，并叠加高斯白噪声，分别利用周期图法、Welch法、AR模型法、Burg法、特征向量法以及MUSIC法估计该混合信号的功率谱，并对结果进行分析比较。7.1 概述信号常用一组变量值表示，如，如果序列在每个时刻的取值不是随机的，而是服从某种固定函数的关系，则称之为确定性信号，比如上一章介绍的阶跃信号、方波信号、正弦波信号等。与确定性信号不同，如果序列在每个时刻的取值是随机变量，则称之为随机信号。随机信号也称为随机过程，具有以下特点：r 随机信号在任何时间的取值都不能先验

4、确定，是随机变化的；r 虽然随机信号取值不能先验确定，但这些取值服从某种统计规律，或者说随机信号可以用概率分布特性统计描述。例如，一个正弦波信号如果式中相位是一个随机变量，假设在-，上均匀分布那么该正弦波信号就是一个随机信号。随机信号又分为平稳随机信号和非平稳随机信号两大类，如果随机信号的概率分布不随时间推移而变化，即信号的统计特性与起始时间无关，只与时间间隔有关，这类随机信号就称为平稳随机信号，否则为非平稳随机信号。因而，平稳信号常被称作时不变信号，而非平稳信号常被称作时变信号。提示：时变与时不变是指随机信号的某个或某几个统计量是否随时间变化，而非针对信号的波形而言。随机信号的一个重要性质是

5、遍历性，它所关心的是从随机信号的一次观测值能否估计信号的统计量。如果一个随机信号的均值和方差都具有遍历性，则称该信号具有各态历经性，此时只通过一次观测值就可估计随机过程的多个统计特性，这在实际当中就非常有意义了。为了简化分析和运算，我们常常假设所获得的随机信号具有各态历经性。7.2 随机信号的统计特性随机信号可用其统计特性进行描述，这些统计特性又可进一步划分为一阶、二阶和高阶（三阶及三阶以上）统计特性。一阶特性主要是指均值函数，即信号的数学期望，二阶特性主要包括方差、相关函数、协方差函数和功率谱密度，高阶特性则指代高阶矩。对于平稳信号而言，二阶特性往往是最重要的。接下来将对随机信号的几个主要统

6、计特性逐一介绍，并给出MATLAB工具箱提供的相应函数。7.2.1 均值和方差均值和方差是常用的描述随机信号统计特性的特征量。假设表示一个各态历经的随机信号，那么该信号的均值和方差的定义式分别为：式中表示信号的概率密度函数。在MATLAB工具箱中，分别提供了两个函数来实现随机序列的均值和方差的计算，即均值函数mean和方差函数var。l m = mean(X)该函数计算均值的公式为：格式中若输入参量X是向量，则返回值m为该向量的均值；若输入参量X为矩阵，则返回值m为一个行向量，其元素分别为矩阵X每一列的均值。l y = var(X)若输入参量X是向量，则返回值y为该向量的方差；若输入参量X为矩

7、阵，则返回值y为一个行向量，其元素分别为矩阵X每一列的均值。例7-1估计给定随机序列的均值和方差。% 随机序列的均值和方差估计x = normrnd(0.1,0.8,1,6) %产生正态分布的随机序列xm = mean(x)xv = var(x)这里x是一个长度为6的正态分布随机序列，结果如下：x = -0.2461 -1.2325 0.2003 0.3301 -0.8172 1.0527xm = -0.1188xv = 0.68407.2.2 相关函数相关函数是最为常用的描述平稳随机信号统计特性的二阶统计量之一，可用来描述随机信号不同时刻间的相关程度。我们首先定义单个平稳随机信号的自相关函数

8、，其定义式为：同理，对于两个平稳随机信号和，其互相关函数可定义为：式中是和 的二维混合概率密度函数。MATLAB工具箱提供了函数xcorr来实现相关函数的估计，其调用格式为：l c = xcorr(x,y);返回随机序列x和y的互相关函数，长为2N-1，假设x和y长度都为N。l c = xcorr(x,y,maxlags);输入参量maxlags为x和y之间的最大延迟，返回值c长度为2maxlags+1，如缺省，c长度默认为2N-1。l c = xcorr(x,y,maxlags,option);这里option包括以下几个选择项：r biased有偏估计，自相关函数和互相关函数的估计表达式分

9、别取：r unbiased无偏估计，自相关函数和互相关函数的估计表达式分别取：r coeff m=0的相关函数值归一化为1。r none 不做归一化处理。作为特例，函数xcorr也可用来计算自相关函数，调用格式为：l c = xcorr(x);l c = xcorr(x,maxlags,option);例7-2产生一个相关正态随机序列，并计算它的自相关函数。% 相关函数估计程序a = 0.85;sigma = 1.5;N=200;x = randn(N,1);y(1) = sigma*x(1)/sqrt(1-a2);for i=2:Ny(i) = a*y(i-1)+sigma*x(i);end

10、c = xcorr(y,coeff);subplot(121);plot(y,LineWidth,2);grid onsubplot(122);plot(c,LineWidth,2);grid on程序中的变量y就是产生的相关正态随机序列，其相关函数表达式为程序运行结果如图7-1所示。图7-1 随机序列波形及其自相关函数7.2.3 协方差函数除了相关函数，协方差也是一个重要的描述随机过程统计特性的二阶统计量。同样，我们也首先定义单个平稳随机信号的自协方差函数，其定义式为：自协方差函数用于表示它的两个随机变量和之间的相关程度。同理，定义两个平稳随机信号互协方差函数，其估计表达式为：互协方差函数表

11、示随机变量和之间的相关程度，实际上，也是表示随机信号和之间的相关程度。MATLAB工具箱提供了函数cov计算随机信号的协方差函数，其调用格式为：l C = cov(X);当X为向量时，返回值C是一个标量值，表示方差；当X为矩阵形式，每一行表示一次观测值，每列对应一个变量，则返回值C是一个协方差矩阵。提示：方差函数等价于协方差函数的对角元素，即var(X) = diag(cov(X)。l C = cov(X,Y);计算输入向量X和Y的互协方差矩阵，X和Y是具有相同长度的列向量。例7-3产生两个被白噪声污染的随机序列，并分别计算二者的自协方差及互协方差。% 计算自协方差及互协方差程序fs = 10

12、0;t =0:1/fs:.1;f1 = 15;f2 = 30;x0 = cos(2*pi*f1*t).;x = x0+randn(size(x0);y0 = cos(2*pi*f2*t).;y = y0+randn(size(y0);Cx = cov(x)Cy = cov(y)Cxy = cov(x,y)程序运行结果如下： Cx = 1.5176Cy = 0.7867Cxy = 1.5176 0.5067 0.5067 0.7867同时，我们计算序列x和y的方差，得到： var(x)ans = 1.5176 var(y)ans = 0.7867可以看出，方差函数确实等价于协方差函数的对角元素。

13、7.3 功率谱估计 功率谱估计是利用已观测到的一定数量样本数据估计一个平稳随机信号的功率谱密度，因其能够分析信号的能量随频率变化的分布特性，在许多实际应用中功率谱的分析与估计已变得越来越重要。例如，强噪声背景下弱信号的检测、雷达和声纳回波信号分析、故障诊断等多个领域都涉及到了功率谱的分析与估计。功率谱密度，简称功率谱，可定义为从表达式中可以看出，随机信号的自相关函数与它的功率谱构成一对离散Fourier变换对，二者分别从时域和频域反映了随机信号的二阶统计特性。功率谱估计方法有很多种，大体上可以分为三大类：r 非参数方法r 参数方法r 子空间方法非参数方法是直接基于观测信号本身进行功率谱估计的一

14、类方法，而不涉及一些模型参数等。其中最简单的一种方法是周期图法，在其基础上的改进方法包括Welch方法等。参数方法是一类基于模型参数的方法，该类方法假设观测信号为一个带有加性白噪声的线性系统的输出，在此架设基础上进行功率谱的估计。代表算法包括Yule-Walker自回归（AR）算法和Burg算法。这类算法首先要估计产生观测信号的线性系统的一组参数（系数）。当样本数较少时，参数方法往往优于经典的非参数方法。子空间方法常常被称作高分辨率方法或超分辨率方法，这类方法基于观测信号相关矩阵的特征分析或特征分解来对功率谱进行估计，例如多信号分类（MUSIC）方法和特征向量（EV）法。这类方法非常适用于分析

15、信号中含有线谱分量的情况，尤其在信噪比很低的情况下，子空间方法就显示出了优越性。MATLAB信号处理工具箱针对以上三类功率谱估计方法分别提供了相应的函数，线将这些方法汇总在一个表中，见表7-1。表7-1 功率谱估计方法及MATLAB工具箱函数类别方法函数非参数方法周期图法periodogram,Welch方法pwelch参数方法AR模型法spectrum.yulear, pyulearBurg方法spectrum.burg, pburg子空间法MUSIC方法spectrum.music, pmusic特型向量法spectrum.eigenvector, peig7.3.1 非参数方法7.3.1

16、.1 周期图法非参数方法中最简单的是周期图法，它是把随机序列x(n)的N个观测数据视为一能量有限的序列，直接计算x(n)的离散傅立叶变换，然后再取其幅值的平方，并除以N，作为序列x(n)真实功率谱的估计，即这里周期图法的优点是计算效率高，但缺点是频率分辨率较低。MATLAB信号处理工具箱提供了函数periodogram来实现周期图法的功率谱估计，其调用格式包括：l Pxx,w = periodogram(x,window,nfft,range);l Pxx,f = periodogram(x,window,nfft,fs,range);l periodogram();返回值Pxx为利用周期图法

17、估计的功率谱；w和f分别是单位为弧度和Hz的频率向量，与Pxx相同长度；输入参量x为有限长信号向量；window为所选的窗函数，默认矩形窗；nfft是FFT点数，一般取2的整数次幂，默认值256；range包括两种情况：r twosided在频率区间0, fs)上计算双边PSD，当样本数据x为复数时，这是默认情况；r onesided在指定的频率区间上计算单边PSD，当样本数据x为实数时，这是默认情况。提示：range字段可以放在输入参量window之后的任意位置。例7-4产生一频率为300Hz的谐波信号（分实数和复数两种情况），叠加高斯白噪声，利用周期图法估计该混合信号的功率谱。% 周期图法

18、估计信号功率谱Fs = 1000; t = 0:1/Fs:.3;x = cos(2*pi*t*300)+0.1*randn(size(t); % 产生实数混合信号periodogram(x,twosided,512,Fs) % 使用默认窗函数y = exp(j*2*pi*t*300)+0.1*randn(size(t); % 产生复数混合信号periodogram(y,onesided,512,Fs) 程序运行结果如图7-2所示。可以看出，当信号形式为实数时，功率谱图形会产生关于奈奎斯特频率对称的两个波峰，而当信号形式为复数时，仅保留300Hz位置的波峰。图7-2 周期图法估计功率谱密度图（左

19、图为实数信号，右图为复数信号）7.3.1.2 Welch方法周期图法估计出的功率谱性能并不好，当数据长度过大时，谱的曲线起伏剧烈，而当数据长度过小时，谱的分辨率又不好，因此在其基础上又产生了一系列的改进方法，共同的原则是将周期图法进行平滑，使得估计方差减小，其中以Welch方法为代表。Welch方法是一种最为广泛的经典功率谱估计方法，该法在周期图法的基础上主要进行了两处改进：其一是将采样数据进行分段，并允许各段数据有交叠的部分；其二是每段数据加窗处理时，窗函数不局限于矩形窗，也可使用汉宁窗或海明窗等，从而可以改善谱失真的现象，具体步骤如下： 步骤r 先将N个数据分为L段，每段M个数据；r 选择

20、适当的窗函数w(n)，对每段数据依次加权，然后确定每段的周期图；r 对分段周期图进行平均，得到功率谱估计。根据以上步骤，最终得到的功率谱估计表达式为式中为第段的采样数据，为选取的窗函数，是每段数据的长度，是归一化因子，用于保证所得到的谱是渐进无偏估计，表达式为MATLAB信号处理工具箱提供了函数pwelch来实现Welch改进周期图法的功率谱估计，其调用格式为：l Pxx,w = pwelch (x,window,noverlap,nfft);l Pxx,f = pwelch (x,window,noverlap,nfft,fs);l Pxx,f = pwelch (x,window,nove

21、rlap,f,fs);l pwelch (x, );格式中，noverlap为分段重叠点数，必须小于窗长度，默认值为窗长的一半，其他各参量同periodogram。例7-5产生一频率为300Hz的谐波信号（复数形式），叠加高斯白噪声，利用Welch法估计该混合信号的功率谱。% Welch法估计信号功率谱Fs = 1000; t = 0:1/Fs:.3;x = exp(j*2*pi*t*300)+0.1*randn(size(t); pwelch(x,128,256,Fs) 程序运行结果如图7-3所示，与图7-2相比较可以看出，Welch方法估计得到的功率谱图比周期图法估计结果要平滑很多，曲线起

22、伏较为平缓，优于周期图法。图7-3 Welch法估计信号功率谱图7.3.2 参数方法任何具有有理功率谱密度的随机信号都可以看作一白噪声激励一个线性系统而形成，如果根据已观测到得数据估计出这个线性系统地模型参数，就不必假设N个观测值之外的数据都为零，从而克服经典功率谱估计的缺点。因而，参数方法并不直接从观测信号来估计功率谱，而是首先建立一个线性系统模型，将观测信号视为该系统的输出，而后对线性系统的参数进行估计，具体步骤如下： 步骤r 选择合适模型；r 用已观测到的数据估计模型参数；r 将模型参数代入相应公式，得到功率谱估计值。这类方法的优点之一是，当信号采样长度较短时，参数方法能够获得比非参数方

23、法更高的功率谱分辨率。这里我们介绍两种较为典型的参数方法，即AR模型方法和Burg法。7.3.2.1 AR模型法AR模型（自回归模型）使用最为广泛，这是因其参数计算为线性方程，比较简单，同时适用于表示很窄的频谱。AR模型是一个全极点的模型，该模型当前输出等于当前输入与过去输出的加权和，利用差分方程可表示为：式中是高斯白噪声序列，为AR模型阶数，为AR模型系数，。由上面的差分方程，进一步得到系统函数最终得到AR模型的功率谱估计式：待估计参数包括，共个，这些参数可由Yule-Walker方程求得提示：任何ARMA和AR模型均可用一个无限阶的AR模型表示，即使模型选择不当，只要模型阶数足够高，仍能比

24、较好地逼近建模的随机过程。MATLAB信号处理工具箱提供了函数aryule来估计AR模型参数，其调用格式为：l a,e = aryule(x,p)返回值a包含AR模型系数参量，返回值e为方差参数，x为输入信号向量，p为AR模型阶数。这里该函数采用的是Yule-Walker方法估计模型的参数。例7-6估计例7-5中信号的10阶AR模型参量。% AR模型参数估计Fs = 1000;t = 0:1/Fs:.3;x = cos(2*pi*t*300)+0.1*randn(size(t); a,e = aryule(x,10)程序运行结果如下： a = 1.0000 0.1433 0.4348 -0.3

25、201 -0.1534 0.1377 -0.1429 -0.0693 0.0136 0.0276 -0.0281e = 0.0218同时，MATLAB工具箱提供了函数pyulear利用AR模型估计功率谱，即首先由Yule-Walker方法估计AR模型参数，而后将参数代入AR模型的功率谱估计公式，便得到观测信号的功率谱估计值。函数pyulear的具体调用格式包括：l Pxx = pyulear(x,p,nfft,fs)l Pxx,f = pyulear(x,p,nfft,fs,range)l Pxx,w = pyulear(x,p,nfft,range)l pyulear(.)格式中，输入参量p

26、为AR模型阶数，其他参量同函数periodogram例7-7采用AR模型法估计例7-5中信号的功率谱密度。% AR模型法估计功率谱Fs = 1000;t = 0:1/Fs:.3;x = cos(2*pi*t*300)+0.1*randn(size(t); pyulear(x,10);程序运行结果如图7-4所示，注意到，图形中的横坐标为归一化之后的频率，这与上面几个函数如periodogram、pwelch是不同的。图7-4 AR模型法估计功率谱图7.3.2.2 Burg法Burg法估计功率谱又被称为最大熵谱估计，这里首先定义了谱熵的概念，其表达式为：式中为功率谱密度函数。Burg最大熵谱估计可

27、以表述为，求功率谱，使其在约束条件下之下，能够使谱熵达到最大。对这一具有约束条件的优化问题采用Lagrange乘子法求解，得到式中为Lagrange乘子。MATLAB信号处理工具箱提供了函数pburg实现功率谱估计，其调用格式为：l Pxx = pburg (x,p,nfft,fs)l Pxx,f = pburg (x,p,nfft,fs,range)l Pxx,w = pburg (x,p,nfft,range)l pburg (.)格式中，Pxx为功率谱估计值，x是输入信号向量，p为AR模型阶数，nfft是FFT点数，fs为采样频率，range同函数periodogram。例7-8采用Bu

28、rg法估计例7-5中信号的功率谱密度。% Burg法估计功率谱Fs = 1000;t = 0:1/Fs:.3;x = cos(2*pi*t*300)+0.1*randn(size(t); pburg(x,10);程序运行结果如图7-5所示，注意到，图形中的横坐标也为归一化之后的频率，这与函数pyulear相同。图7-5 Burg法估计功率谱图7.3.3 子空间法子空间法主要包括MUSIC法和特征向量法，这类方法都是以自相关矩阵的特征分析为基础，将系统自相关矩阵分成两个子空间：信号子空间和噪声子空间。其中特征向量法主要适用于混有高斯白噪声的正弦信号的功率谱估计，而MUSIC法适合更为普遍情况下正

29、弦信号的功率谱估计。子空间法功率谱估计步骤一般如下： 步骤r 根据N个观测样本值，估计自相关函数，假设和分别为第个复正弦信号的功率和频率，是白噪声方差；r 对进行特征分解，得到个特征值，对应特征向量；r 估计功率谱7.3.3.1 特征向量法在上面的功率谱估计表达式中，如果令，则所得的功率谱估计就称为特征向量估计，此时表达式写作在MATLAB信号处理工具箱中提供了函数peig来实现特征向量法估计功率谱，其调用格式为：l S,w = peig(x,p);l S,f = peig(x,p,f,fs)l S,f = peig(x,p,nfft,fs,nwin,noverlap);l peig()格式中

30、，返回值S为伪谱估计值；w和f分别是单位为弧度和Hz的返回频率值；x为输入信号向量；p用来指定信号子空间中特征向量的数目；nfft为FFT算法的长度，默认值为256；nwin用于指定矩形窗的宽度，信号向量x被分割成长度为nwin的各段，默认值为2*p；noverlap是每段之间重叠的长度，默认值为nwin-1。例7-9试用特征向量法估计例7-5中信号的功率谱密度。% 特征向量法估计功率谱Fs = 1000;t = 0:1/Fs:.3;x = cos(2*pi*t*300)+cos(2*pi*t*250)+0.1*randn(size(t); p=10;peig(x,p)程序运行结果如图7-6所

31、示。7.3.3.2 MUSIC法在上面的功率谱估计表达式中，如果令，则所得的功率谱估计就称为MUSIC法估计，此时表达式写作在MATLAB信号处理工具箱中提供了函数pmusic来实现MUSIC法功率谱估计，其调用格式为：l S,w = pmusic(x,p);l S,f = pmusic(x,p,f,fs)l S,f = pmusic(x,p,nfft,fs,nwin,noverlap);l pmusic()格式中各参数说明同函数peig。例7-10试用MUSIC法估计例7-5中信号的功率谱密度。% MUSIC法估计功率谱Fs = 1000;t = 0:1/Fs:.3;x = cos(2*pi

32、*t*300)+cos(2*pi*t*250)+0.1*randn(size(t); p=10;pmusic(x,p)程序运行结果如图7-7所示。 图7-6 特征向量法估计功率谱图 图7-7 MUSIC法估计功率谱图7.4 高阶谱估计前面介绍的功率谱估计过程中，通常需要假设观测到的数据是由高斯白噪声激励线性系统产生的，而实际当中信号或噪声往往并不服从高斯分布，并且信号的自相关函数只包含幅度信息，而不包含相位信息，这样就使得二阶统计量（时域的相关函数，频域的功率谱）存在一定的局限性，难以满足实际需要。因而研究三阶及更高阶的统计量是非常必要的，将它们统称为高阶统计量。目前，高阶统计量已成为信号处理

33、的一种有力工具，并且在信号分析与处理领域中扮演着越来越重要的角色。本节将首先给出一些相关概念，如高阶累积量和高阶谱，之后具体介绍高阶谱估计的方法，包括非参数化的估计方法以及参数模型法。7.4.1 高阶累积量和高阶谱定义高阶累积量和高阶谱是两个最为常用的高阶统计量，这里首先给出二者的定义。为了方便表示，假设随机过程的均值为0，那么零均值平稳随机过程的二阶、三阶、四阶累积量定义式分别为：这里假设为实信号。根据以上定义式，我们知道二阶累积量其实就是信号的自协方差函数，当均值为零时，也等价于信号的自相关函数。这里我们仅讨论到第四阶累积量，因为实际应用中三阶和四阶累积量是最为常用的，而再高阶的累积量应用

34、场合较少，这里就不再讨论。对上面二阶、三阶、四阶累积量做Fourier变换，有上面三个表达式分别是功率谱、双谱和三谱的定义式。所谓双谱是具有两个频率变量的函数，而三谱是具有三个频率变量的函数。从而，我们可以了解高阶谱的含义。高阶谱也称为多谱，就是有多个频率变量的谱。习惯上，功率谱用表示，双谱用表示，而三谱用表示。高阶累积量和高阶谱具有以下特点：r 假设，其中和是两个相互独立的随机过程，则 ；r 若是高斯随机信号，则，这里，此时如果是非高斯信号，且与独立，则利用高阶累积量可从混合信号中恢复出来；r 若是线性过程，即，其中是独立同分布随机过程，那么相应的功率谱、双谱和三谱的表达式分别如下，其中。从

35、表达式可以看出，功率谱并不包含任何的相位信息，而高阶谱则包含了相位信息。当是非高斯随机过程时，其相位信息就可从高阶谱中恢复。根据以上特点，我们知道相对于二阶统计量，高阶累积量和高阶谱具有其自身的优点，能够弥补二阶统计量在某些方面的不足。具体来说，高阶统计量适用于以下两方面的情况：1) 加性噪声是高斯的，而信号是非高斯的，此时就可利用高阶统计量将信号从混合信号中恢复出来；2) 线性系统不是最小相位情况，或者说随机过程是非线性的，此时高阶谱可以恢复出信号的相位信息。7.4.2 高阶累积量估计对于高阶累积量和高阶谱的估计，我们这里所应用的是高阶谱分析工具箱（High-order Spectral A

36、nalysis Toolbox）。该工具箱功能比较强大，包括了很多先进的信号处理算法，如谱估计、高阶谱估计、时频分布计算等，具体应用领域又包括参数和非参数的盲系统辨识、时延估计、谐波恢复、波达方位估计、Volterra模型参数估计、自适应线性预测。一般观测到的样本数都为有限长度，如，那么利用观测到的样本值来估计二阶、三阶和四阶累积量的公式分别如下：这里 式中和是选取的样本数目，要求当时，只对求和，选取合适的是为了获取无偏估计。令公式中，就会得到自累积量的估计式。工具箱提供了函数cum2x、cum3x、cum4x分别估计二阶、三阶和四阶累积量，函数cumest用来估计自累积量。调用格式分别如下：

37、l cvec = cum2x(x,y,maxlag,samp_seg,overlap,flag);l cvec = cum3x(x,y,z,maxlag,samp_seg,overlap,flag,k1)l cvec = cum4x(w,x,y,z,maxlag,samp_seg,overlap,flag,k1,k2)l cvec = cumest(y,norder,maxlag,samp_seg,overlap,flag,k1,k2);各参数的说明如下：r y是样本矩阵或向量；r norder指定了累积量阶数，可以取2、3、4三个值，默认值为2；r maxlag指定了累积量允许的最大延迟，默

38、认值为0；r samp_seg规定每个数据段的样本数，默认值是样本序列长度；r overlap指定数据段间重叠的百分比，允许的最大值为99，默认值是0；r flag包括biased和unbiased两种情况，如果为biased，则计算有偏样本估计，否则计算无偏样本估计，默认为biased；r 如果y为矩阵，则每一列单独进行累积量估计运算，此时overlap需设为零，samp_seg为行维数，最终累积量为每列估计值的平均结果；r 参数k1和k2决定计算累积量函数的哪个一维切片，默认值都为零。r 产生一频率为300Hz的谐波信号（复数形式），叠加高斯白噪声，试用特征向量法估计混合信号的功率谱密度。

39、例7-11模拟一个非高斯ARMA随机过程，并估计它的3阶累积量，采用函数cum3x。% 函数cum3x估计3阶累积量u = rpiid(1024,exp);n = 25;y = filter(1,-2,1,-1.5,0.8,u);for k = -n:n cmat(:,k+n+1) = cum3x(y,y,y,n,128,20,biased,k);endfigure;contour(-n:n,-n:n,cmat,8)函数rpiid产生一个独立同分布随机序列。这里样本序列被分割成长度128的数据段，重叠百分比为20%，对每个数据段计算有偏的3阶累积量，再计算平均值。程序运行结果如图7-8(a)所

40、示，等高线图显示了三阶累积量的基本对称性。例7-12采用函数cumest估计例11中随机过程的3阶累积量。% 函数cumest估计3阶累积量u = rpiid(1024,exp);n = 25;y = filter(1,-2,1,-1.5,0.8,u);for k = -n:n cmat(:,k+n+1) = cumest(y,3,n,128,20,biased,k);endfigure;contour(-n:n,-n:n,cmat,8)程序运行结果如图7-8(b)所示，与图(a)相比较可以看出，函数cum3x和函数cumest的估计结果基本一致，二者都体现了三阶累积量的基本对称性。(a) 函

41、数cum3x估计结果 (b) 函数cumest估计结果图7-8 一ARMA(2,1)随机过程的三阶累积量估计7.4.3 高阶谱估计最基本的互双谱估计方法是对三阶累积量做Fourier变换，即这里是随机序列，的Fourier变换。上面的估计方法并非一致性估计，通常需要采取平滑的方法获取一致性。当时，就是双谱估计表达式。所谓平滑处理，一般通过窗函数乘以三阶累积量来完成。假设是一个二维窗函数，其二维Fourier变换有界且非负，并且假设同时，窗函数也要满足三阶累积量的对称性。平滑后的互双谱估计式为：高阶谱分析工具箱提供了不同方法实现高阶谱估计的函数，包括：r 函数bispecd，直接法双谱估计r 函

42、数bispeci，间接法双谱估计r 函数bispect，ARMA模型的理论双谱估计对于这些函数及相应的算法，接下来依次进行介绍。7.4.3.1 函数bispecd 函数bispecd实现直接法双谱估计，算法的实现步骤如下： 步骤r 将样本数据分成段，每段含有个样本值，记作，其中，这里允许两段相邻数据之间有重叠；r 计算离散Fourier变换（DFT）系数，其中，；r 计算DFT系数的三重相关，其中，和应选择满足的值；r 段双谱估计的平均值即为样本数据的最终双谱估计值；，函数bispecd的调用格式为：l bspec,waxis = bispecd(y)l bspec,waxis = bispe

43、cd(y,nfft,wind,samp_seg,overlap)各参数的说明如下：r y是样本矩阵或向量；r nfft为计算FFT的长度，一般取2的整数次幂，常规默认值为128；r wind指定频域平滑窗参数；1) 如果wind是标量，则采用长度为此标量值的Rao-Gabr窗；2) 如果wind等于1，则不采用平滑窗；3) 如果wind为向量，则计算二维窗；4) 如果wind为矩阵，则直接指定二维平滑窗；r samp_seg规定每个数据段的样本数，默认值为8；r overlap指定数据段间重叠的百分比，默认值是50；r 返回值bspec为双谱的估计结果，大小为nfftnfft，其原点在中心位置

44、，轴的方向是向下与向右；r 返回值waxis是与bspec矩阵行和列相关的频率向量；7.4.3.2 函数bispeci 函数bispeci实现间接法双谱估计，算法的实现步骤为： 步骤r 将原个样本数据分成段，每段含有个样本值；r 设为第段数据，计算各段的三阶累积量估计值，其中，；r 取所有段的三阶累积量的平均值作为整个观测数据组的三阶累积量估计；r 计算双谱估计，其中，为二维滞后窗函数；函数bispeci的调用格式为：l bspec,waxis = bispeci(y,nlag)l bspec,waxis = bispeci(y,nlag,samp_seg,overlap,flag,nfft,

45、wind)参数说明如下：r y是样本矩阵或向量；r nlag规定了累积量计算延迟数，假设计算三阶累积量，则要求， ，该参数必须指定，一般可将该值取为样本序列 长度的十分之一。r samp_seg规定每个数据段的样本数，默认值为观测序列长度；r overlap指定数据段间重叠的百分比，取值范围0,99，默认值为0；r flag包括biased和unbiased两种情况，如果为biased，则计算有偏样本估计，否则计算无偏样本估计，默认为biased；r nfft为计算FFT的长度，一般取2的整数次幂，常规默认值为128；r wind指定平滑窗参数，当wind = 0时，采用Parzen窗，否则采

46、用六角形窗，默认Parzen窗。r 返回值bspec为双谱的估计结果，大小为nfftnfft，其原点在中心位置，轴的方向是向下与向右；r 返回值waxis是与bspec矩阵行和列相关的频率向量；例7-13分别利用直接法和间接法估计一个已知的非高斯信号的双谱。% 直接法双谱估计load qpc;figure(1)bspecd = bispecd(zmat,128,5,64,50);% 间接法双谱估计figure(2)bspeci = bispeci(zmat,25,64,50,biased,128,5);程序中采用了MATLAB自带的非高斯检测信号qpc.mat，直接法估计函数bispecd的参

47、数设置中，将FFT长度确定为128，每段数据长64，并且重叠百分比为50%，这些与间接法估计函数bispeci的设置相同，同时间接法采取了有偏估计。程序运行结果如图7-9所示，(a)图为直接法估计结果，(b)图为间接法估计结果，可以看出，两种方法得到的峰值位置相同，结果基本一致。(a)直接法 (b)间接法图7-9 直接法和间接法双谱估计图7.4.3.3 函数bispect 函数bispect采用ARMA模型实现理论双谱的估计。ARMA模型经常用于描述一个随机序列，表达式为式中是独立同分布序列，方差。自回归（AR）多项式和滑动平均（MA）多项式分别定义为，利用ARMA模型估计理论双谱，关键是对模型参数和的估计，之后将ARMA模型参数代入