do the psd using different methods (maybe there is a problem in the normalization9
조회 수: 1 (최근 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
댓글 수: 0
채택된 답변
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
추가 답변 (0개)
참고 항목
카테고리
Help Center 및 File Exchange에서 Parametric Spectral Estimation에 대해 자세히 알아보기
제품
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!