MATLAB高阶累积量函数计算介绍

MATLAB有自带的高阶谱分析工具箱，包括cum3x、cum3est、cum4x、cum4est等，这里对其用法进行简单介绍。
语法

y_cum=cum3x(x,y,z,maxlag,nsamp,overlap,flag,k1);

该函数是用来估计x、y、z的三阶互累积量，当x=y=z时，估计的是x的三阶累积量

参数	意义
x、y、z	分析数据向量或者矩阵，维数必须相同
maxlag	估计的最大延迟量，默认值为0
nsamp	每一段的样本数 ，默认值为数据长度
overlap	数据重叠百分比，默认值为0
flag	'biased’表示有偏估计，'unbiased’为无偏估计
k1	固定的延迟量，C3(m,k1)，k1默认值为0
y_cum	返回估计的互三阶累积量

语法

y_cum=cum3est(y,maxlag,nsamp,overlap,flag,k1);

该函数是用来估计y的三阶累积量

参数	意义
y	分析数据列向量
maxlag	估计的最大延迟量，默认值为0
nsamp	每一段的样本数，默认值为数据长度
overlap	数据重叠百分比，默认值为0
flag	'biased’表示有偏估计，'unbiased’为无偏估计
k1	固定的延迟量，C3(m,k1)，k1默认值为0
y_cum	返回估计的三阶累积量

例子：估计非高斯ARMA（1,2）过程的三阶累积量

clear;
close all;
clc;
%产生非高斯ARMA仿真信号
rand('seed',0);
randn('seed',0);
%非高斯输入
u=rpiid(1024,'exp');
%最大延迟量
maxlag=25;
%ARMA（1,2）过程
sig=filter([1,-2],[1,-1.5,0.8],u);
%显示ARMA过程信号
figure(1);
plot(1:length(sig),sig);
xlabel('n');
ylabel('ARMA(1,2)信号幅值')
%估计三阶累积量
%参数设置
nsamp=128;
overlap=0;
flag='biased';
for k=-maxlag:maxlag
    sig3cum(:,k+maxlag+1)=cum3x(sig,sig,sig,maxlag,nsamp,overlap,flag,k);
end
%显示三阶累积量
figure(2);
%等高线图
subplot(2,2,1);
contour(-maxlag:maxlag,-maxlag:maxlag,sig3cum);
xlabel('延迟量k');
ylabel('延迟量l');
%三维图
subplot(2,2,2);
mesh(-maxlag:maxlag,-maxlag:maxlag,sig3cum);
xlabel('延迟量k');
ylabel('延迟量l');
for k=-maxlag:maxlag
    sig3cum(:,k+maxlag+1)=cum3est(sig,maxlag,nsamp,overlap,flag,k);
end
%显示三阶累积量
%等高线图
subplot(2,2,3   );
contour(-maxlag:maxlag,-maxlag:maxlag,sig3cum);
xlabel('延迟量k');
ylabel('延迟量l');
%三维图
subplot(2,2,4);
mesh(-maxlag:maxlag,-maxlag:maxlag,sig3cum);
xlabel('延迟量k');
ylabel('延迟量l');

结果显示
图中显示的是ARMA(1,2)过程的仿真信号
两种函数计算结果对比

从等高线图和三维图来看，两种不同途径计算的结果是一致的。
关于计算四阶互累积量的函数cum4x和四阶累积量的函数cum4est的各个参数和cum3x、cum3est代表的意义相同，这里不再叙述。
高阶谱分析工具箱下载地址：https://www.mathworks.com/matlabcentral/fileexchange/3013-hosa-higher-order-spectral-analysis-toolbox?s_tid=srchtitle