# do the psd using different methods (maybe there is a problem in the normalization9

2 views (last 30 days)
SYML2nd on 28 Jun 2019
Commented: Star Strider on 29 Jun 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
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

Star Strider on 29 Jun 2019
Edited: Star Strider on 29 Jun 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.
Star Strider on 29 Jun 2019
As always, my pleasure.