Autocorrelation in MATLAB

So I’ve asked around the lab and posted on the MATLAB forum, and nobody seems to be able to solve this, even though it seems it should be simple. This isn’t homework, it is for my work (research).

My code (shown below) does the following. Take a function H which is the square root of a Lorentzian. Inverse fourier transform it to yield the function h. Then take the autocorrelation of h, which should be the Fourier transform of H^2; this is the Fourier transform of a Lorentzian, which is an exponential. I then plot the natural log of this, which should be a triangle function (linear decreasing in each direction from 0). All the supporting math for these statements is here.

When I run this code in MATLAB, it curves up near 0, rather than being linear everywhere (plot posted here). Oddly enough, if I remove the square root and just set H = lorentzian, the log of the autocorrelation is linear like it’s supposed to be. Is my math wrong or my code? Any ideas? Here is the code:

f = [-50:.01:50];
alpha = 1;

H = sqrt(lorentzian(f, alpha)); %square-root lorentzian filter (has exponential autocorrelation)

option = ‘unbiased’; %can be ‘unbiased’ or ‘none’
maxlags = 50;

H_unshifted = ifftshift(H); %shift and scale spectrum for IFFT
h = ifft(H_unshifted); %come back to space domain

[Rhh, lags] = xcorr(h, maxlags, option); %autocorrelation

plot (lags, log(Rhh))
xlabel (‘Lag’)
ylabel (‘log R_{hh}’)

I also define

function l = lorentzian (x, a)
l = a./(x.^2 + a.^2);

Thanks for any insights.

The problem is with the FFTSHIFT, or the lack of it. h has its peak at zero, which is to say that the only large values are near the ends of the vector. MATLAB’s xcorr function pads the vectors with zeroes when it’s doing its convolutions, which means that a lot of the biggest values in your vectors are getting multiplied by zeroes. You want
h = fftshift(ifft(H_unshifted));
or something like it, to move the peak to the middle of the vector before doing the convolution.

Disclaimer: I checked this on Octave and it seemed to work, but I had to hack together an xcorr and I’m not sure it’s quite the same as MATLAB’s.

This worked, thanks!