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

2 ビュー (過去 30 日間)
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

採用された回答

Star Strider
Star Strider 2019 年 6 月 29 日
編集済み: Star Strider 2019 年 6 月 29 日
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 件のコメント
SYML2nd
SYML2nd 2019 年 6 月 29 日
Ok...Thank you for your help and patience
Star Strider
Star Strider 2019 年 6 月 29 日
As always, my pleasure.

サインインしてコメントする。

その他の回答 (0 件)

カテゴリ

Help Center および File ExchangeParametric Spectral Estimation についてさらに検索

製品

Community Treasure Hunt

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

Start Hunting!

Translated by