2
$\begingroup$

Wiener–Khinchin theorem states that autocorrelation function and power spectral density are a Fourier-transform pair - see Wikipidia (and lots of other resources).

That means autocorrelation should be able to be obtained by inverse Fourier transform the spectrum. The following code (run within Octave, with "pkg load signal") shows the Fourier transform of the autocorrelation DOES look like the spectrum, but the inverse Fourier transform of the spectrum does not look like the autocrrelation. What did I do wrong?

### A signal's autocorrelation and its Engergy Spectral Density are Fourier transform pairs.

###   signal
Fs = 1000;
t = 0:1/Fs:1-1/Fs;
x = cos(2*pi*100*t)+randn(size(t));
#plot(x);

###    autocorrelation
Rxx = xcorr(x);
#figure();
plot(Rxx); title("Rxx");

###    autocorrelation FT
RxxDftAbs = abs(fftshift(fft(Rxx)));
freq = -Fs/2:Fs/length(Rxx):Fs/2-(Fs/length(Rxx));
figure(); 
plot(freq,RxxDftAbs); title("RxxDftAbs");

###    Energy Spectral Density
xdft = abs(fftshift(fft(x)));
x_esd = xdft.^2;   # ESD is the same as autocorrelation FT. Here for visualization purpose, using absolute values.
freq = -Fs/2:Fs/length(x_esd):Fs/2-(Fs/length(x_esd));
figure();
plot(freq,x_esd); title("x esd");

### ?????????????
### is it possible to get autocorrelation from the ESD by inverse Fourier transform?
### ?????????????

###    IFT of ESD
x_esd_idft_abs = abs(ifft(fftshift(x_esd)));
figure();
plot(x_esd_idft_abs); title("x esd ift");

Thanks in advance.

Improve this question
$\endgroup$
2
$\begingroup$

Your code works fine. But for the sake of demonstration clarity, just get rid of all the fftshift functions and change your frequency range too.

The main problem is that you should use FFT sizes when calling fft() functions, which by default uses signal length as FFT size which was the problem you faced on the following line :

###    Energy Spectral Density
xdft = abs(fftshift(fft(x)));

which is corrected as :

###    Energy Spectral Density
xdft = abs(fftshift(fft(x,length(Rxx))));

Then your code works just fine.

Improve this answer
$\endgroup$
3
$\begingroup$

Fat32's answer is correct and shows a common pitfall.

The reason that you must do this is because recall that the output of the autocorrelation of a signal of length $N$ is $2N - 1$.

You were performing the FFT with the original sample size of $N = 1000$, effectively destroying necessary information to retrieve the autocorrelation of size $2N - 1 = 1999$.

Improve this answer
$\endgroup$
  • $\begingroup$ If you don't mind me asking, why does the length of the autocorrelation from a signal of length N equal to 2N - 1 ? From a statistical standpoint ( I realize that this is different from DSP standpoint which is why I'm asking ), if you have 1000 data points, then it's only possible to estimate 999 autocorrelations ( from lag 1 to lag 999 ). Thanks. $\endgroup$ – mark leeds Jul 27 '20 at 20:58
  • $\begingroup$ @markleeds This is the case because it is how discrete autocorrelation is defined. You mentioned lags 1-999 when calculating the correlation, but that is not correct. Correlation includes negative lags as well. $\endgroup$ – Envidia Jul 27 '20 at 21:11
  • $\begingroup$ The autcorrelation function ( in statistics ) is symmetric so there is no need for negative values. But, if you say there are negative lags in the DSP version , then the whole function is just totally different from the one in statistics so don't worry about it. If I ever need the autocorrelation function in DSP, I'll investigate it. Thanks. $\endgroup$ – mark leeds Jul 28 '20 at 13:39
  • $\begingroup$ @markleeds Negative lags are not just in "the DSP version". Yes, autocorrelation functions are symmetric but the autocorrelation integral is defined from negative infinity to infinity. Symmetry is simply a property and helps in cutting down calculations (among other things). Same with the Fourier transform of a real signal: the result is symmetric but the integral is still from negative infinity to infinity by definition. The math does not discriminate. It is us who know what functions are being operated on and what simplifications to make. $\endgroup$ – Envidia Jul 28 '20 at 15:05
  • $\begingroup$ If you calculate an autocorrelation function in R ( or any other decent statistical programming language ), you'll never see the negative lag side of it. The manner in which it's it's calculated in statistics doesn't use the spectral density so there's no reason to show the negative side of it. Essentially, given the way it's calculated, you'd need to do the same exact calculation twice. What I'm attempting to say is that the negative side doesn't naturally come out of the procedure. But I see where you're coming from regarding DSP. It's gonna get calculated due to the algorithm used. $\endgroup$ – mark leeds Jul 29 '20 at 11:42

Your Answer

By clicking “Post Your Answer�?, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.