Plotting different PSDs on MATLAB

6 visualizaciones (últimos 30 días)
Gianmarco Broilo
Gianmarco Broilo el 28 de Feb. de 2022
Respondida: Milan Bansal el 24 de Dic. de 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')

Respuestas (1)

Milan Bansal
Milan Bansal el 24 de Dic. de 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!

Categorías

Más información sobre Spectral Measurements en Help Center y File Exchange.

Productos


Versión

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by