使用FFT估计信号互功率谱密度

前言

       关于信号自谱的估计可以参考我的另一篇文章：使用FFT估计信号功率谱、功率谱密度。有人在评论区问互谱怎么估计，故决定再出一篇估计互谱的文章，供大家参考。

互功率谱密度（Cross Power Spectrum Density，CPSD）

       互功率密度谱是在频域内描述两个不同信号之间统计相关程度的一种方法。自谱用自相关函数表达，那互谱自然就要用两个信号的互相关函数来表达。计算互谱与自谱的思路基本一致，本文同样针对频率成分不随时间变化的平稳信号，采用FFT方法计算两个信号的互功率谱密度。

1）仿真生成两个信号x(t)、y(t)

fs = 32000;
N = 1024;
t = 0:1/fs:(N-1)/fs;
freqx = 4000;
freqy = 8000;

xt = cos(2*pi*freqx*t) + randn(size(t))/sqrt(2);
yt = cos(2*pi*freqy*t) + randn(size(t))/sqrt(2);

2）加窗并修正，再傅里叶变换

win = hamming(N);
xw = 1.586*xt.*win.'; 
yw = 1.586*yt.*win.'; 

xdft = fft(xw);
xdft = xdft(1:N/2+1);
ydft = fft(yw);
ydft = ydft(1:N/2+1);

3）计算互功率谱密度

Pxy = (1/(fs*N))*abs(xdft.*conj(ydft));
Pxy(2:end-1) = 2*Pxy(2:end-1);

4）绘图，并与MATLAB自带接口cpsd的结果进行对比

fb = 0:fs/N:fs/2;
figure
subplot(211)
plot(fb/1e3,10*log10(Pxy),'r-')
grid on
grid minor   
xlabel("Frequency (kHz)")
ylabel("Power/Frequency (dB/Hz)")
title("Cross Power Spectrum Density Using FFT")
subplot(212)
cpsd(xt,yt,win,0,N,fs);
grid minor

两者曲线几乎一致，验证通过！

结语

        在对比过程中，要注意输入参数是否一致，如窗函数、窗长度等，不同参数的结果区别还是比较明显的。

总结不易，点个赞吧(●'◡'●)