Plotting different PSDs on MATLAB

3 ビュー (過去 30 日間)
Gianmarco Broilo
Gianmarco Broilo 2022 年 2 月 28 日
回答済み: Milan Bansal 2023 年 12 月 24 日
So I am trying to plot on the same graph the PSD of the noise from an accelerometer and an IMU. The PSDs of the noise are: Pacc = 10e-14 + 10e-18*f^-2; Pimu = 10e-12; with the frequency going from 10e-5 to 10Hz. This is what I have now but I am not too sure on the result:
Thank you for any help!
N = 1e6;
dt = 0.2;
PSD = 10e-12; %PSD of the IMU
sigma2 = PSD/dt; %PSD of WGN is PSD = sigma^2*dt
x = randn(N,1)*sqrt(sigma2);
xi = zeros(size(x));
xint(1) = x(1);
for n = 2:length(x)
xint(n) = (x(n)+x(n-1))/2*dt+xint(n-1);
end
NFFT = 1e5;
[px,~] = pwelch(x,hanning(NFFT),0.5,NFFT,1/dt,'twosided');
[pxint,f] = pwelch(xint,hanning(NFFT),0.5,NFFT,1/dt,'twosided');
Pacc = 10e-14 + 10e-18*f.^-2;
Pint = (2*pi*f).^(-1).*Pacc;
loglog(f,px)
hold on
loglog(f,Pint)
hold off
legend('noise IMU','noise Acc')

回答 (1 件)

Milan Bansal
Milan Bansal 2023 年 12 月 24 日
Hi Gianmarco Broilo
I understand that you are trying to plot the Power Spectral Density (PSD) of the noise from an accelerator and IMU, but the plots are not as expected.
Please modify the code as shown in the code snippet below to get the correct plots.
N = 1e6;
dt = 0.2;
PSD = 10e-12; %PSD of the IMU
sigma2 = PSD/dt; %PSD of WGN is PSD = sigma^2*dt
x = randn(N,1)*sqrt(sigma2);
xi = zeros(size(x));
xint(1) = x(1);
for n = 2:length(x)
xint(n) = (x(n)+x(n-1))/2*dt+xint(n-1);
end
NFFT = 1e5;
[px,~] = pwelch(x,hanning(NFFT),0.5,NFFT,1/dt,'twosided');
[pxint,f] = pwelch(xint,hanning(NFFT),0.5,NFFT,1/dt,'twosided');
Pacc = 10e-14 + (10e-18)*(f.^-2); % modified
Pint = Pacc./(2*pi*f); %modified
loglog(f,Pacc) % modified
hold on
loglog(f,Pint)
hold off
legend('noise IMU','noise Acc')
Hope it helps!

カテゴリ

Help Center および File ExchangeParametric Spectral Estimation についてさらに検索

製品


リリース

R2021b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by