自相关的功率谱密度 (PSD)

信息处理 傅里叶变换 卷积 功率谱密度 自相关
2022-02-06 13:47:50

我有一个大小为 ,它是通过具有高斯平均值和标准偏差的随机数生成器生成的。我想找到所述信号的功率谱密度(PSD)。信号数据、MATLAB-R2018a 代码和引用包含在 Github 存储库 ( https://github.com/dustinma324/Burugence ) 中。u(x)Nμ=0σ=N

摘自 Press 2007 的“Numerical Recipes”, 其中,星号代表复共轭。对上述方程进行归一化返回 PSD。

Corr(g,h)j=GkHk
GkHkgjhjN2

我无法理解执行 DFT 的顺序之间的区别。我有两个路线图来获得 PSD。第一个是,

PSD=fft(u(x))fft(u(x))

其中 PSD 是通过将自相关函数分解为傅里叶变换的乘积和两个信号的傅里叶变换(卷积)的复共轭来找到的。第二,

R11(r)=u(x)u(x)
PSD=|fft(R11(r=0))|

其中 PSD 是通过首先找到信号的自相关函数,然后取自相关函数的傅里叶变换的幅度来找到的。

在此处输入图像描述

深蓝色线是第二种方法,红色是第一种方法。通过采用自相关函数的,可以捕获一般形状,但后一种方法在较高频率处显示少量能量积累。这违反了 Schwartz 不等式,因为 to 向量的内积不能大于它们大小的乘积。FFT

为什么后一种方法似乎在较高频率下携带更多能量,是否有特定的数学原因?

下面是我的 MATLAB 代码,用于两种不同的方法。

方法 1(卷积)

n = length(x); %Nx
pwr = [];
meanU = mean(u(:,tavg_start:end),2);
for i = tavg_start:1:nt
   U_fft = fft(u(:,i)-meanU,n)/n;
   pwr = [pwr U_fft.*conj(U_fft)];
end
meanpwR = mean(pwr,2);

方法二(自相关函数)

pwR = [];
n = length(x);
meanU = mean(u(:,tavg_start:end),2);
for i = tavg_start:step:nt
    U11 = double(u(:,i)-meanU);
    Ucorr = xcorr(U11); % Output will be 2*length(Ucorr)
    N = 2^nextpow2(length(Ucorr));
    R_fft = fft(Ucorr,N)/(N^2);
    pwR = [pwR abs(R_fft)];
end
meanpwR = mean(pwR,2);

我找不到任何可以解释我的特定问题的东西。任何输入/建议/讨论将不胜感激。

更新 2020 年 2 月 4 日

自相关函数类似于低通滤波器,意味着低频通过,而高频在某个截止频率处被衰减(逐渐损失)。在执行 ACF 的傅里叶变换时,衰减的频率(未解析的能量)被混叠为较低级别的频率,导致略有增加。

此链接中描述了一个过渡带(https://www.ni.com/en-us/innovations/white-papers/18/anti-aliasing-filters-and-their-usage-explained.html)。我相信对我来说,这个频段的范围从奈奎斯特频率,其中 Basu 使用傅里叶域中的 3/2 规则对相同的信号进行去混叠。据说这个频带会导致信号混叠,因此显示在我的上图中。k=N/2k=2N/3

从技术上讲,这两种方法都显示出某种形式的混叠。那么有不同种类的混叠吗?观察到一个在一系列波数上混叠,而另一个则显示出尖峰。

谢谢,

达斯汀

1个回答

我不确定我是否遵循您的确切问题。

为什么后一种方法似乎在较高频率下携带更多能量,是否有特定的数学原因?

不会。这两种方法是等效的,并且会与您使用的数值精度相匹配。

自相关函数类似于低通滤波器,意味着低频通过,而高频在某个截止频率处被衰减(逐渐损失)。

我不明白这种说法。自相关如何成为低通滤波器?自相关是自相关。

从技术上讲,这两种方法都显示出某种形式的混叠。那么有不同种类的混叠吗?观察到一个在一系列波数上混叠,而另一个则显示出尖峰。

您的输入信号已经被采样并且您没有更改采样率,所以我也不明白这一点。

如果从两种方法估计的 PSD 不匹配,那么您的代码中有错误。下面是一个使用 octave/matlab 代码的例子:

N = 1000;
x = rand(N,1);
Nconv = 2*N - 1;

% PSD from autocorrelation.
a = conv(x,flipud(conj(x)));
A = abs( fftshift(fft(a)) );

% PSD from FFT.
X = fftshift(fft(x,Nconv));
A2 = X .* conj(X);

% The residual between the two estimates.
err = 10*log10( norm(A - A2) / norm(A) );