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.