利用MATLAB实现信号的1½维谱估计

众所周知，高阶谱是由高阶累积量经过多维傅里叶变换，但是存在着计算量大且耗时比较长的缺点。为了使高阶谱像功率谱那样能够适用FFT，计算简单，许多学者对此进行了研究，但是效果却不尽如人意。近年来，有的学者通过“降维”的思想，将双谱投影到一维频率空间上，取得了一定的理论成果，其主要思想是通过取对角切片的方式实现的，即选择特殊的延时${\tau }_1={\tau }_2=\cdots ={\tau }_k=\tau$，进而简化双谱的计算。
已知零均值平稳随机过程$x\left(n\right)$的三阶累积量和三阶矩相等，三阶累积量的定义为
$$c_{3x}\left({\tau }_1,{\tau }_2\right)=E\left[x\left(n\right)x\left(n+{\tau }_1\right)x\left(n+{\tau }_2\right)\right]$$
当${\tau }_1={\tau }_2$得到三阶累积量的主对角切片为
$$c\left(\tau \right)=c_{3x}\left(\tau ,\tau \right)=E\left[x\left(n\right)x\left(n+\tau \right)x\left(n+\tau \right)\right]$$
对$c\left(\tau \right)$进行一维傅里叶变换可以得到$x\left(n\right)$的$1\frac{1}{2}$维谱，即
$$SB\left(f\right)=\sum _{\tau =-\infty }^{\infty }c\left(\tau \right)e^{-j2\pi f\tau }$$
$1\frac{1}{2}$维谱实际是双谱在平面$f_1=f_2$的投影。
下面给出一个关于求$1\frac{1}{2}$维谱的例子
例：估计仿真信号$x\left(n\right)$的$1\frac{1}{2}$维谱
$$x\left(n\right)=\sum _{i=1}^6{A}_i\cos \left(2\pi f_{i}n+{\phi }_i\right)$$
其中$A_i$为确定常数$i=1,2,\cdots ,6$，${\phi }_i$是$\left[1,2\pi \right]$，$i=1,2,\cdots ,5$，上均匀分布的随机变量，取$f_1=0.6381Hz$，$f_2=0.8435Hz$，$f_3=2Hz$，$f_4=0.4909Hz$，$f_5=1.7671Hz$，$f_6=f_4+f_5$，${\phi }_6={\phi }_4+{\phi }_5$

代码如下

clear;
close all;
clc;
%样本数
N=128*64;
%采样频率
fs=8;
%频率参数
f1=0.6381;
f2=0.8435;
f3=2;
f4=0.4909;
f5=1.7671;
f6=f4+f5;
%仿真信号
a=2*pi*rand(5);
x=cos(2*pi*f1/fs*(1:N)+a(1))+cos(2*pi*f2/fs*(1:N)+a(2))+...
  cos(2*pi*f3/fs*(1:N)+a(3))+cos(2*pi*f4/fs*(1:N)+a(4))+...
  cos(2*pi*f5/fs*(1:N)+a(5))+cos(2*pi*f6/fs*(1:N)+a(4)+a(5));
%估计三阶累积量
K=64;
M=fix(N/K);
%分段，去均值
for i=0:K-1
    y=x(i*M+1:(i+1)*M);
    y=y-mean(y);
    xk(i+1,:)=y;
end
nfft=256;
%计算累积量
for t=-(nfft/2-1):nfft/2
    c(t+nfft/2)=third_slice_cumulant(xk,t);
end
%FFT变换
figure(1);
subplot(211);
Y=fft(c,nfft);
P=Y.*conj(Y)/nfft;
f=fs*(0:nfft/2-1)/nfft;
plot(f,P(1:nfft/2));
title('x(n)的1又1/2谱');
xlabel('频率');
ylabel('幅值');
hold on;
subplot(212);
Y=fft(x,nfft);
P=Y.*conj(Y)/nfft;
f=fs*(0:nfft/2-1)/nfft;
plot(f,P(1:nfft/2));
title('x(n)的功率谱');
xlabel('频率');
ylabel('幅值');

其中third_slice_cumulant.m为

function y=third_slice_cumulant(xk,t)
[K,M]=size(xk);
s1=max([1,1-t]);
s2=min([M,M-t]) ;
for i=1:K
templ=xk(i,:);
temp=sum(templ(s1:s2).*templ(s1+t:s2+t).*templ(s1+t:s2+t));
cc(i)=temp/M;
end
y=sum(cc)/K;
end

运行结果为

可以看出$1\frac{1}{2}$维谱是频率的一维函数，类似于功率谱，但是它只在$f_4，f_5，f_6$处有谱线存在，这和功率谱是不一样的，这是因为这三个频率恰好对应发生二次相位耦合关系的频率分量，因此检测出二次相位耦合现象。