# Power Spectral Density (PSD) from autocorrelation

I have a 1D fractional Brownian motion (fBm) signal, $$u(x)$$, of size $$N$$, generated through a random number generator with a gaussian mean, $$\mu=0$$, and standard deviation, $$\sigma=\sqrt{N}$$. I want to find the power spectral density (PSD) of said signal. Signal data, MATLAB-R2018a code, and citation is included in the Github repository (https://github.com/dustinma324/Buruglence) .

From "Numerical Recipes" by Press 2007, $$Corr(g,h)_{j} = G_{k}H_{k}^*$$ where $$G_{k}$$ and $$H_{k}$$ are the discrete Fourier transforms (DFT) of $$g_{j}$$ and $$h_{j}$$, and the asterisk represents the complex conjugate. Normalizing the above equation by $$N^2$$ returns the PSD.

I am having trouble understanding the difference between the order of performing DFT. I have two road maps in mind to get the PSD. The first being,

$$PSD = fft\Big(u(x)\Big) \cdot fft\Big(u(x)\Big)^{*}$$

where the PSD is found through breaking up the auto-correlation function into the product of the Fourier transform and the complex conjugate of the Fourier transform (convolution) of the two signals. The second,

$$R_{11}(r) = u(x)*u(-x)$$ $$PSD = \bigg|fft\Big(R_{11}(r=0) \Big)\bigg|$$

where the PSD is found by finding the auto-correlation function of the signal first, then taking the magnitude of the Fourier transform of the auto-correlation function.

The dark blue line is the second method, and red is the first. By taking the $$FFT$$ of the auto-correlation function, the general shape is captured, but the latter method shows small amounts of energy accumulation at the higher frequencies. This violates the Schwartz inequality, since the inner product of to vectors cannot be larger than the product of their magnitudes.

Is there a specific mathematical reason as to why the latter method appears to carry more energy at the higher frequencies?

Here's my MATLAB code below for the two different methods.

Method 1 (Convolution)

n = length(x); %Nx
pwr = [];
meanU = mean(u(:,tavg_start:end),2);
for i = tavg_start:1:nt
U_fft = fft(u(:,i)-meanU,n)/n;
pwr = [pwr U_fft.*conj(U_fft)];
end
meanpwR = mean(pwr,2);


Method 2 (Auto-correlation function)

pwR = [];
n = length(x);
meanU = mean(u(:,tavg_start:end),2);
for i = tavg_start:step:nt
U11 = double(u(:,i)-meanU);
Ucorr = xcorr(U11); % Output will be 2*length(Ucorr)
N = 2^nextpow2(length(Ucorr));
R_fft = fft(Ucorr,N)/(N^2);
pwR = [pwR abs(R_fft)];
end
meanpwR = mean(pwR,2);


I have not been able to find anything that explains my particular question. Any inputs/suggestions/discussions would be greatly appreciated.

UPDATE Feb 4th, 2020

Auto-correlation function is analogous to a low-pass filter, meaning low frequencies are getting through, and high frequencies are attenuated out (gradual loss) at a certain cutoff frequency. The attenuated frequencies (unresolved energy) is being aliased into lower level frequencies when performing the Fourier transform of the ACF, resulting in the slight increase.

A transition band was described in this link (https://www.ni.com/en-us/innovations/white-papers/18/anti-aliasing-filters-and-their-usage-explained.html). I believe for me this band ranged from $$k = N/2$$, Nyquist frequency, to $$k=2N/3$$, where Basu dealiased the same signal using the 3/2 rule in the Fourier domain. This band is said to cause aliasing in the signal, thus showing up in my graph above.

Technically both methods shows some form of aliasing. So are there different kinds of aliasing? One is observed to aliased over a range of wavenumbers, while the other shows a spike.

Thanks,

Dustin

• Did you make sure to flip and take the conjugate of one of the signals before you do the convolution? If not, you aren't performing autocorrelation. Jan 30, 2020 at 0:33
• @Envidia In the first method yes. The second method is through finding the auto-correlation function of the signal first, then taking the fft. Jan 30, 2020 at 0:36
• The second equation only makes it clear that you are going to perform regular convolution of $u(x)$ with itself. If you are using a convolution function in MATLAB or some other application, you have to reverse one of the signals in time and take the complex conjugate in order to perform autocorrelation. Jan 30, 2020 at 0:42
• @Envidia To reverse the function $u(x)$, since it is still in real space, can I simply flip the vector left/right or up/down depending on it being a row or column major? Jan 30, 2020 at 0:51
• Yes you can. However, if you are using a function that finds the autocorrelation for you, then you don't need to do this. Jan 30, 2020 at 0:52

Is there a specific mathematical reason as to why the latter method appears to carry more energy at the higher frequencies?

No. The two methods are equivalent and will match to the numerical precision you are using.

Auto-correlation function is analogous to a low-pass filter, meaning low frequencies are getting through, and high frequencies are attenuated out (gradual loss) at a certain cutoff frequency.

I do not understand this claim. How is auto-correlation a low-pass filter? Auto-Correlation is auto-correlation.

Technically both methods shows some form of aliasing. So are there different kinds of aliasing? One is observed to aliased over a range of wavenumbers, while the other shows a spike.

Your input signal is already sampled and you are not changing the sampling rate, so I do not understand this either.

If the PSDs as estimated from two methods do not match, then there is a mistake in your code. Below is an example using octave/matlab code:

N = 1000;
x = rand(N,1);
Nconv = 2*N - 1;

% PSD from autocorrelation.
a = conv(x,flipud(conj(x)));
A = abs( fftshift(fft(a)) );

% PSD from FFT.
X = fftshift(fft(x,Nconv));
A2 = X .* conj(X);

% The residual between the two estimates.
err = 10*log10( norm(A - A2) / norm(A) );