Greetings,
I am using Matlab R2012b 64-bit on Windows 7 in order to estimate the power
spectrum of a simple signal that is:
cos(10*t) + sin(20*t) defined in the time interval from 0.0 to 10.0
Here is what I did:
1) There was no predefined value of points so I took 1001 points by forming a
time array and sampled the function at these points.
2) I calculated the autocorrelation function of the signal with the built-in
function "xcorr".
3) I performed FFT directly on the autocorrelation function with no additional
options.
4) The resulting data set had imaginary parts(I think this should not be the
case since the autocorrelation function is even).
5) I took the absolute value of the FFT of ACF and plotted it as the power
spectrum from the ACF after having it shifted with fftshift.
6) I took the FFT of the sinusoidal signal after padding it with 1000 (n-1)
zeros.
7) I calculated the power spectrum from the FFT of the signal as the square of
the absolute value of the frequency coefficients.
8) I plotted both of them and compared them by assigning an array to difference
between them.
9) They turned out to be the same.
Everything seems to be working correctly but my concerns are as follows:
- How can this data be interpreted? I expected the sum and difference
frequencies of the individual sinusoidal functions in the signal.
-How can I normalize the output and represent negative frequencies properly
since they have no physical significance?
-Does Matlab use FFTW libraries? My main task is to achieve the same results in
a C program using these libraries.
I have benefitted from the following URL:
http://dsp.stackexchange.com/questions/1919/efficiently-calculating-autocorrelation-using-ffts?rq=1
in organizing the post and achieving these results.
I would be grateful for an explanation of this phenomenon and normalization
procedure.
Thanks in advance
P.S : The code snippet I have used is as follows and it can reproduce my
results:
t=transpose(0:0.01:10);
n=size(t,1);
f=cos(10*t)+sin(20*t);
acf=xcorr(f,f);
FFTacft(acf);
PSacfs(FFTacf);
FFTpt([f;zeros(n-1,1)]);
PSfTp.*conj(FFTp);
figure;
subplot(2,1,1);
plot(fftshift(PSf));
title('PS from FFT');
subplot(2,1,2);
plot(fftshift(PSacf));
title('PS from ACF');
d=PSacf-PSf;
fprintf('Max error = %6.2f \n', max(abs(d)));
