自相关和谱密度

信息处理 matlab 信号分析 功率谱密度 八度
2022-02-04 08:03:24

Wiener-Khinchin 定理指出自相关函数和功率谱密度是傅里叶变换对 - 请参阅Wikipidia(以及许多其他资源)。

这意味着自相关应该能够通过傅里叶逆变换光谱来获得。下面的代码(在 Octave 中运行,带有“pkg load signal”)显示了自相关的傅里​​叶变换确实看起来像频谱,但频谱的傅里叶逆变换看起来不像自相关。我做错什么了?

### A signal's autocorrelation and its Engergy Spectral Density are Fourier transform pairs.

###   signal
Fs = 1000;
t = 0:1/Fs:1-1/Fs;
x = cos(2*pi*100*t)+randn(size(t));
#plot(x);

###    autocorrelation
Rxx = xcorr(x);
#figure();
plot(Rxx); title("Rxx");

###    autocorrelation FT
RxxDftAbs = abs(fftshift(fft(Rxx)));
freq = -Fs/2:Fs/length(Rxx):Fs/2-(Fs/length(Rxx));
figure(); 
plot(freq,RxxDftAbs); title("RxxDftAbs");

###    Energy Spectral Density
xdft = abs(fftshift(fft(x)));
x_esd = xdft.^2;   # ESD is the same as autocorrelation FT. Here for visualization purpose, using absolute values.
freq = -Fs/2:Fs/length(x_esd):Fs/2-(Fs/length(x_esd));
figure();
plot(freq,x_esd); title("x esd");

### ?????????????
### is it possible to get autocorrelation from the ESD by inverse Fourier transform?
### ?????????????

###    IFT of ESD
x_esd_idft_abs = abs(ifft(fftshift(x_esd)));
figure();
plot(x_esd_idft_abs); title("x esd ift");

提前致谢。

2个回答

Fat32 的答案是正确的,并且显示了一个常见的陷阱。

您必须这样做的原因是,回想一下长度的信号的自相关的输出是N2N1

您正在使用原始样本大小的自相关的必要信息N=10002N1=1999

您的代码工作正常。但为了演示清楚起见,只需摆脱所有fftshift功能并更改您的频率范围。

主要问题是在调用 fft() 函数时应该使用 FFT 大小,默认情况下使用信号长度作为 FFT 大小,这是您在以下行中遇到的问题:

###    Energy Spectral Density
xdft = abs(fftshift(fft(x))); 

更正为:

###    Energy Spectral Density
xdft = abs(fftshift(fft(x,length(Rxx)))); 

然后你的代码工作得很好。