Plotting different PSDs on MATLAB

1 view (last 30 days)
Gianmarco Broilo
Gianmarco Broilo on 28 Feb 2022
Answered: Milan Bansal on 24 Dec 2023
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')

Answers (1)

Milan Bansal
Milan Bansal on 24 Dec 2023
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!

Products


Release

R2021b

Community Treasure Hunt

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

Start Hunting!