信息处理 - 来自 MATLAB 中自相关的 PSD - 吾爱随笔录

来自 MATLAB 中自相关的 PSD

信息处理信号分析噪音功率谱密度自相关随机过程

2022-02-15 14:05:03

我正在尝试模拟由方程式定义的简单随机过程：

\frac{1}{v} \frac{d b}{d t} + Γ_{0} b = \sqrt{σ} R (t),

$\begin{equation} \frac{1}{v}\frac{db}{dt} +\Gamma_0 b= \sqrt{\sigma}R(t), \end{equation}$ 在哪里

R (t)

$R(t)$ 是零均值白噪声（时间相关的增量），并且

v

$v$ ,

Γ_{0}

$\Gamma_0$ 和

σ

$\sigma$ 都是真正的正参数。我正在使用 Euler-mayurama 方案对此进行模拟，即

b (t + Δ t) = b (t) - v Γ_{0} b (t) Δ t + v \sqrt{σ} \sqrt{Δ t} N (0, 1),

$\begin{equation} b(t+\Delta t) = b(t) -v \Gamma_0 b(t) \Delta t + v\sqrt{\sigma} \sqrt{\Delta t} N(0,1), \end{equation}$ 和

N (0, 1)

$N(0,1)$ 是标准正态随机数，对于有限区间

t \in [0, T]

$t\in [0,T]$ . 我正在尝试使用 Wiener-Khintchine 定理计算功率谱密度

S [f] = \int_{- \infty}^{\infty} e^{- 2 π f τ} ⟨ b (t) b^{*} (t + τ) ⟩ d τ,

$\begin{equation} S[f] = \int_{-\infty}^{\infty} e^{-2\pi f\tau} \left\langle b(t)b^*(t+\tau)\right\rangle d\tau, \end{equation}$ 从自相关函数

⟨ b (t) b^{*} (t + τ) ⟩

$\left\langle b(t)b^*(t+\tau)\right\rangle$ . 从分析上讲，我可以计算出它，它给了我

S [f] = \frac{v σ}{2 Γ_{0}} \frac{2 v Γ_{0}}{(v Γ_{0})^{2} + (2 π f)^{2}},

$\begin{equation} S[f] = \frac{v\sigma}{2\Gamma_0}\frac{2v\Gamma_0}{(v\Gamma_0)^2 + (2\pi f)^2}, \end{equation}$ 但是当我尝试计算相同的数值时，我的函数看起来完全不同。我是否正确执行此操作？（见下面的 Matlab 代码）

%% Autocorrelation test
close all; clear all; clc;

% Inputs
dt = 1e-3;
T = 2;
tv = 0:dt:T;
Lt = length(tv);

v = 3.5;
Gamma0 = 20.3;
sigma = 0.75;
a = v*Gamma0;

% Frequency spectrum properties
Fs = 1/dt; % (Hz) sampling frequency, based on Ben's paper

sims = 1; % number of simulations
U_Walks = zeros(Lt,sims);
for nn = 1:sims
    b = zeros(Lt,1);
    for ii = 1:Lt-1
        b(ii+1) = b(ii) - v*Gamma0*b(ii)*dt + ...
            v*sqrt(sigma)*sqrt(dt).*normrnd(0,1);
    end
    Rxx(:,nn) = xcorr(b);
    U_Walks(:,nn) = b;
    disp(nn);
end

%% Stats
LR = length(Rxx);
N = 2^nextpow2(LR);
Y = fft(Rxx,N); % Taking only one side of the y-axis symmetric FFT
Y = Y(1:N/2+1);
mY = (Y).^2/N;
F = Fs*(0:(N/2))/N; % calculating the frequency range for the x-axis

S_f = 10*log10(mY);

S_f_analytic = 10*log10((v*sigma/2/Gamma0.*(2*a./(a^2 + (2*pi*F).^2))));

%% PLOTS
close all;

figure('units','normalized','outerposition',[0 0 1 1])
set(gcf,'Color','w');
plot(F,S_f,'r',F,S_f_analytic,'b','LineWidth',3);
legend('Numeric','Analytic');
title(['PSD']);
xlabel('f (Hz)'); ylabel('S[f] (dB)');
set(gca,'FontSize',20);

1个回答

您的代码/方法中有两个错误。第一个是术语 $\sqrt{\Delta t}$ 在你的第二个公式中；它应该被替换为 $\Delta t$ . 第二个是根据估计的自相关计算功率谱。您所做的是将 FFT 的结果平方Y以获得mY，但这是不正确的。首先，Y是复值，其次，自相关的傅里叶变换是功率谱，所以不需要对 FFT 的结果求平方，只需要计算其幅度即可。

总之，更换行

v*sqrt(sigma)*sqrt(dt).*normrnd(0,1);

经过

v*sqrt(sigma)*dt.*normrnd(0,1);

和

mY = (Y).^2/N;

经过

mY = abs(Y)/N;

这应该会给你想要的结果，如下图所示：

其它你可能感兴趣的问题

上一篇如何从幅度谱图重建声音？下一篇模拟到离散到模拟系统采样问题