Power spectral density: why this 2 method are equal

the Power spectral density can be calculated in two way:

1. by doing the fourier transform of the autocorrelation


2. by doing (abs(X(f)).^2 where X(f)=fft(x(t))

Can you explain me passage by passage why these are equal?

I'll try to explain better my problem, because I tried to simplify it, but maybe (since I am not an expert) I have worsen the situation. I had to do Power spectral density of a random signal (the velocity of a cfd simulation). Firstly I did the fft of the autocorrelation:

   v is my signal
   r =xcorr(v,v)
   PSD = abs(fftshift(fft(v)))

After that I realised that in order to be verified the fact that the integral of the PSD is equal to the rms^2 of the signal, my result has to be scaled by the factor 1/(fsN), as it is shown here http://aaronscher.com/Course_materials/Communication_Systems/documents/PSD_Autocorrelation_Noise.pdf Unfortunately I did not understand why this factor was necessary and I did not find any explanation about this online or on my books. So I looked at this video that gives an explanation using the Parseval Theorem https://www.youtube.com/watch?v=D67ZgH8FEAI&t=1s . So I thought if I demonstrate that abs(fftshift(fft(v))) is equal to abs(fft(v)).^2 than I can demonstrate why I need that scaling factor 1/(fsN) in order to satisfy my check.

I hope you can forgive me if I did not explained everything from the beginning but I tried to make it simple (I obviously failed). I hope you can help me, I am struggling with this demonstration from one week and none can help me

Answer

Regarding what we discussed in the comments, here is an explanation of the different concepts I was referring to:

The concrete usage of terms may vary, but my interpretation is this: For a stochastic process $x(t)$, the autocorrelation function is defined as $$r(\tau) = \mathbb{E}\{x(t) \cdot x(t+\tau)\}.$$ Note that the process needs to be stationary, otherwise the autocorrelation function depends on $t$ as well. Its power spectral density is given by $$R(f) = \mathbb{E}\{|X(f)|^2\}.$$ These two quantities are linked to each other via the Wiener-Khinchin-Theorem which shows that $$R(f) = \int_{-\infty}^{\infty} r(\tau) {\rm e}^{-\jmath 2\pi f \tau} {\rm d} \tau,$$ i.e., they are Fourier transform pairs.

Since the expectation is hard to implement in practice (we typically do not have an infinite ensemble available), we need to resort to other methods. When the process is ergodic, ensemble average can be replaced by temporal average and we may, e.g., estimate the autocorrelation function from a window $T$ via $$r_T(\tau) = \frac{1}{T} \int_T x(t) \cdot x(t+\tau){\rm d}t.$$ Note that $r(\tau) = \lim_{T\rightarrow \infty} r_T(\tau)$. On the other hand, considering using $R_T(f) = |X_T(f)|^2$ as an estimate for the power spectral density (referred to as the periodigram), where $X_T(f)$ is the spectum of a $T$-windowed signal leads to an estimate that does not converge since $S(f) \neq \lim_{T\rightarrow \infty} S_T(f)$.

Now, I'm not sure what exactly is the question. Is it the explanation of the Wiener-Khinchin-Theorem for the stochastic signals? Or is it the deterministic equivalent as, e.g., discussed in this answer?

*edit: Regarding your clarification that you posted as a reply to your question(?): My feeling is this question has little to do with ACF and PSD being Fourier transform pairs, instead you seem to be worried how to normalize the result of a Discrete Fourier transform (carried out via an FFT) such that its units make sense physically. It's probably worth asking this in a separate question, but the short version is here. There are two things to be considered:

• There are different conventions how to normalize DFT and IDFT. Let $\mathbf{F}$ be a DFT matrix with elements ${\rm e}^{-\jmath \frac{mn}{N}}$, then you can either (i) define DFT via $\mathbf{D} = \mathbf{F} \cdot \mathbf{d}$ in which case you need to define the IDFT via $\mathbf{d} = {\color{red}{ \frac{1}{N} }} \mathbf{F}^H \cdot \mathbf{D}$ or (ii) define the DFT via $\mathbf{D} = {\color{red}{ \frac{1}{N} }}\mathbf{F} \cdot \mathbf{d}$ which leads to the IDFT $\mathbf{d} = \mathbf{F}^H \cdot \mathbf{D}$ or alternatively (iii) you let the DFT be $\mathbf{D} = {\color{red}{ \frac{1}{\sqrt{N}} }}\mathbf{F} \cdot \mathbf{d}$ which leads to the IDFT $\mathbf{D} = {\color{red}{ \frac{1}{\sqrt{N}} }} \mathbf{F}^H \cdot \mathbf{d}$. Different conventions are in place, Matlab uses the first by default, whereas Maple will prefer the third.


• The Fourier transform changes the units: the ACF has a unit of power (say, Watt), the Fourier transform is an integral over time, hence it adds seconds, so that the unit is Watt seconds or better Watt/Hz. It's a density anyways (how do Watts distribute over Hertz?). If you want to approximate a Fourier transform via a Discrete Fourier transform, you have to be super careful, since the DFT/IDFT does not change the units. If your input to the DFT is a sequence with unit, say, Watt, its DFT coefficients will have the same unit. If you want to mirror what happens in a Fourier transform, you need to multiply the result with the sampling period of the samples in time domain. In other words: If you have a function function $u(t)$ with spectrum $U(f)$ and you want to approximate $U(f)$ with a DFT, you need to consider a finite window (say, length $T$) and a finite sampling interval (say, $t_0 = T/N$). Then, ${\color{red}{t_0}}\cdot D[\mu]$ is an approximation for $U(\mu f_0)$ where $f_0 = 1/T$. It's an approximation since sampling as well as truncation can cause errors (one of them will for sure).

I think the factor you were looking for can very likely be a combination of these two effects, your $1/(f_s N)$ is the same as my $t_0/N$.