do the psd using different methods (maybe there is a problem in the normalization9
1 visualización (últimos 30 días)
Mostrar comentarios más antiguos
SYML2nd
el 28 de Jun. de 2019
Comentada: Star Strider
el 29 de Jun. de 2019
Hi,
I would like to verify if the result is the same using the first method described herehttps://it.mathworks.com/help/signal/ug/power-spectral-density-estimates-using-fft.html and the autocorrelation. Unfortunately I think I have some trouble in the normaliza
vel=readtable('v_c.txt');
v_c=table2array(vel);
Tinc=0.001;
Fs=1/Tinc;
df=0.01
L=length(v_c);
f = transpose(Fs*(0:(L/2))/L)
rms2=rms(v_c(:,2:2)).^2
%First method
xdft = fft(v_c(:,2:2));
xdft = xdft(1:L/2+1);
psdx = (1/(Fs*L)) * abs(xdft).^2;
psdx(2:end-1) = 2*psdx(2:end-1);
%test 1st method
Cumulatedpsdx=sum(psdx)*df
%Seconod method using the autocorrelation
Rxx =xcorr((v_c(:,2:2)),(v_c(:,2:2)), length((v_c(:,2:2)))-1);
freq = -Fs/2:Fs/length(Rxx):Fs/2-(Fs/length(Rxx));
delta_f=(freq(2)-freq(1))
PSDfromRxx = abs(fftshift(fft(Rxx)))*(1/(Fs*length((v_c(:,2:2)))))
%test 1st method
CumulatedPSDfromRxx=sum(PSDfromRxx)*delta_f
figure (1)
loglog(f,psdx(1:50001,:),freq,PSDfromRxx );
If you look at the plot, you can notice that the graphs have a slightly different magnitude. It seems to me that the second graph (freq,PSDfromRxx) is almost half the graph (f,psd1(1:50001,:))
Can you help me to give an explanation to this?? Can you suggest how to set correctly the normalization. It is really important to me
Thank you
0 comentarios
Respuesta aceptada
Star Strider
el 29 de Jun. de 2019
Editada: Star Strider
el 29 de Jun. de 2019
I cannot run your code, since I do not have ‘xdft’, and I cannot figure out from your code what it should be. However, that may not be relevant.
Your observation here may be all that is necessary:
‘It seems to me that the second graph (freq,PSDfromRxx) is almost half the graph (f,psd1(1:50001,:))’
It is, because you are calculating the fft, and although you appear to be scaling it for the length of the ‘Rxx’ vector, you have not allowed for the energy in the fft result being evenly divided between the two halves of the fft result. The amplitude of those is halves (‘PSDfromRxx’) is therefore half the amplitude of the original signal.
This is implied, although not specifically stated in the fft documentation.
Also, you do not have to use readtable and table2array to import your data. This single dlmread call will do everything you want (unless you also want the headers, however you do not use them in your code so they are likely not necessary):
v_c = dlmread('v_c.txt', '\t', 1, 0);
EDIT — (29 Jun 2019 at 21:21)
After your latest edit, ‘psd1’ is missing, so I still cannot run your code. That is likely not a problem, because my initial conclusion remains valid.
6 comentarios
Más respuestas (0)
Ver también
Categorías
Más información sobre Parametric Spectral Estimation en Help Center y File Exchange.
Productos
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!