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

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

Answer 1

[Star Strider]

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.

[SYML2nd]

sorry for the mistake, I fixed it the problem. Now the code should run, so you can see the two graphs, that are similare but they are not overlapped.

In order to be sure about the amplitude, I also tried to see if the parameter Cumulatedpsdx (I added it in the updated code) is equal to rms2, and if is also CumulatedPSDfromRxx equal to rms2. This check is verified, so the amplitude should be right. So I don't know what I am missing.

[Star Strider]

If you multiply ‘PSDfromRxx’ by 2, the amplitudes are essentially the same:

PSDfromRxx = abs(fftshift(fft(Rxx)))*(2/(Fs*length((v_c(:,2:2)))));

so:

psdxMean = mean(psdx)
PSDfromRxxMean = mean(PSDfromRxx)
PSDfromRxxMean2 = mean(2*PSDfromRxx)

produces:

psdxMean =
        0.0234964069340912

PSDfromRxxMean =
        0.0117483209493056

PSDfromRxxMean2 =
        0.0234966418986113

I stand by my previous conclusion.

[SYML2nd]

Ok...if you multiply ‘PSDfromRxx’ by 2, the graphs are coincident and the check about the mean is verified.

Now my last question is:

I need to verify the check

CumulatedPSDfromRxx = rms2

Why multiplying by 2 this check is not satisfied? It is because if I multiply the results by 2, I need to select just half of the graph PSDfromRxx ??? (This check is really important for me)

T

[Star Strider]

Check to see how rms is calculated. You need to calculate it with respect to your original fft result, not the corrected version to match the amplitudes, for the very reason I already stated. The original fft result will be correct with respect to the RMS calculation, because the energy in the fft will be the same, although it will have one-half of the correct amplitude in the plot.

Both of these are entirely correct, however it depends on the result you want. The RMS value must be calculated with the original signal (or its Fourier transform), however because of the energy distribution in the fft, the fft amplitudes must be multiplied by 2 for the plots to match.

It only seems confusing. It all makes perfect sense.

[SYML2nd]

Ok...Thank you for your help and patience

[Star Strider]

As always, my pleasure.