Finding the auto-correlation of the noise signal x explicitly using the FFT() and IFFT() functions

Hello,

    Been trying for hours to calculate the auto-correlation of a noise signal without using the xcorr operator and not sure if my results are right. My tutor wants us to find the auto-correlation  by calculating the auto-spectra, the auto-covariance from the spectra and finaly the auto-correlation from the auto-covariance without using matlab operators. My plots  show a peak at zero, which I think is positive, with the auto-correlation having a peak at zero with the max value of 1 in the vertical axis, which seems correct. Although my x-axis show the 0:N points and I don't know how to display the x-axis correct. From what I understand should be a time-axis. I'm also not sure my results are correct or it was just a coincidence. Would be good to have some feedback and tips. The code is underneath a bit long but I hope it's clear. thanks you so much. Code:

fs=11025; %sampling frequency

N=1024; %number of points

T=N/fs; %time

x=2*rand(N,1)-1; %generate noise signal

%creating a delay of signal x

d=zeros(64,1);

x_delay=cat(1,d,x);

%fft of signal x and delay

X=fft(x);

X_delay=fft(x_delay);

%windowing signal x

h=hanning(2*N);

hX=X.*h(1:N);

HX=fft(hX);

%windowing signal x_delay

H=hanning(2*length(X_delay));

hX_delay=X_delay.*H(1:length(X_delay));

HX_delay=fft(hX_delay);

%finding the auto-spectra

for i=1:N

Sxx(i)= (HX(i).*conj(HX_delay(i))/T);

end

%finding auto-covariance

Sxxt=ifft(Sxx);

%plotting the auto-covariance

subplot(4,1,2)

plot(abs(Sxxt));

title('auto-covariance')

%finding the boundary of the auto-covariance Sxx(0)

Sxx0=mean(Sxx);

%finding the auto-correlation

Rxx=Sxxt/Sxx0;

%plotting the auto-correlation

subplot(4,1,3);

plot(d,abs(Rxx));

title('auto-correlation')

Answer 1

[Wayne King]

Hi, here is an example:

      x = randn(16,1);
      xpsdmag = abs(fft(x,2^nextpow2(2*length(x)-1))).^2;
      acf = ifft(xpsdmag);
      acf = acf(1:length(x));
      acf = acf./acf(1);
      % compare to xcorr()
      [xc,lags] = xcorr(x,'coeff');
      xc = xc(length(x):end);
      isequal(xc,acf)

[catarina]

hey...thank you so much! so much simple. Although I'm a bit confused with the nextpow2 but will research about it hadn't seen/used the operator before. I know it has relation to the way fft algorithm is calculated in powers of 2, but not really sure, haven't gone deep into that. Also in regards to the x-axis when I plot the correlation how to I set it? at the moment it shows the lenght of (x) cause I'm doing:

plot(acf)

not sure how to set the x-axis correctly..any tips on that....

Thanks you so much for your help....

[Wayne King]

Hi Catarina, I did that just so I could demonstrate perfect agreement with xcorr()

To plot the normalized acf I would use stem()

x = randn(16,1);

xpsdmag = abs(fft(x,2^nextpow2(2*length(x)-1))).^2;

acf = ifft(xpsdmag);

acf = acf(1:length(x));

acf = acf./acf(1);

stem(0:length(acf)-1,acf); set(gca,'ylim',[-1.5 1.5]);

xlabel('Lags');

[Tahariet Sharon]

Hi,

I was just wondering, if you don't square (.^2) in line 2 of your answer, would it make a difference if we are interested in computing the ACF? I have seen this formula used for ACF but without squaring the absolute value. Thanks!