Wiener-Khinchin theorem and fluctuations in flicker noise

Question

Suppose that we have a flicker noise process, that takes values $d(t)\in\{-1,+1\}$. We discretize time into $t\in[0,dt,2dt,\ldots, (N-1)dt]$, and $dt = T/N$ with $T$ being the total length of the 'cycle'. At every time step we flip the state of $d$ with probability $W dt/2$. Ultimately we end up with a list of outcomes composed of $\pm1$ of length $N$, e.g. $d = [-1,-1,-1,\ldots, 1,1,\ldots,-1]$.

The process defined in this way, has correlation $\langle d(t) d(t')\rangle = \exp(- W |t-t'|)$ or $Q(k dt) = \langle d(k dt) d(0)\rangle = \exp(- W k |dt|)$ for discrete times and $k=0,\ldots, N-1$. According to Wiener-Khinchin theorem, I could compute fast Fourier transform (FFT) of my list of outcomes $d$ or of a list of $Q = [1, e^{-W dt}, e^{-2 W dt}, \ldots, e^{-(N-1)W dt}]$ and should get the same results for $N\to \infty$. Taking FFT of $Q$ is easy, since we have a geometric series $$ \tilde{Q}(k)=\sum_{n=0}^{N-1}e^{-2\pi i \frac{kn}{N}} e^{-W n dt} = \frac{1-e^{-W N dt}}{1-e^{-2\pi i \frac{k}{N}-W dt}}\approx \frac{N}{T} \frac{W}{W^2 + (2\pi \frac{k}{N})^2}, $$ where we Taylor expanded the denominator and kept the first order correction, and numerator was approximated by 1, which is valid for large $N$ (note that $dt = T/N$).

If I calculate absolute value square of FFT of $d$, i.e. $$ \tilde{V}(k) = \frac{1}{N}\left|\sum_{n=0}^{N-1} e^{-2\pi i \frac{kn}{N}}d(n)\right|^2, $$ where $d(n)$ is $n$-th outcome of my random process, I should (in the limit of $N\to \infty$) obtain the same result as for $\tilde{Q}(k)$ (isn't this the operational meaning of Wiener-Khinchin theorem?). However, what I got is strongly fluctuating picture, with general trend behaving like $\tilde{Q}$ see the plot, where $\omega_k = 2\pi \frac{k}{NT}$.

The questions that I've got:

1. shouldn't we observe much smoother behavior of sampled process (sampling procedure explained below), similar to 'analytical' results?
2. If I repeat the sampling procedure $M$ times, and average all repetitions, I get smooth behavior, which is exactly what I'd expect. Why it doesn't work for a single (large $N$) realization, where 'ergodicity' should smoothen the fluctuations?

Sampling code snippet (in python):

import numpy as np #standard numerical tools 
import matplotlib.pyplot as plt #for plotting
from scipy.fft import fft, ifft #fast fourier transform and it's inverse

'''
Used parameters:
N = 10**7
Trep = 10**5
Wmin = 0.1
cut_off = 10**5  - how many frequencies we're depicting 
'''
#drawing random number for flipping rate
W = np.random.uniform(low = Wmin, high = Wmax ) 

#expressing total cycle time related to flipping rate (i.e. related to correlation length) 
T = Trep/W


t_list = np.linspace(0,T, N+1)[:-1] #list from 0, to (N-1) dt

dt = t_list[1] #defining dt = T/N 

d = [-1] #wlog setting initial value to -1 

for idx in range(N-1):

    I = np.random.rand() < W * dt/2 #defining successful flip
    
    if I :
        dn.append(d[idx] * -1) #flipping the new outcome
    else:
        dn.append(d[idx] * 1) #leaving the new outcome as the previous one

omega_k = np.linspace(0, N , N+1)[:-1]* 2 * np.pi/T #defining frequencies

Q_WK = N/T * np.real(( W  ) /(W**2 +    omega_k**2)) #approximated geometric series

FT = (abs(fft(d))**2)/(N) #calculating FFT for the outcome dn


plt.rcParams.update({'font.size':30})
plt.figure(100,figsize=(24,16))
plt.grid(True)
plt.plot(omega_k[:cut_off],  0.5 * FT[:cut_off], label = 'sampled FT')
plt.plot(omega_k[:cut_off],  Q_WK[:cut_off], label = r'$\frac{W}{W^2 + \omega_k^2}$', linewidth = 4.0)
plt.legend()
plt.show() 

Answer 1

I think you have run into the classic problem of trying to estimate the spectrum from a single time series observation. It is well known that the amplitude spectrum, i.e., the modulus of the periodogram (Finite Fourier transform) of the signal itself is not a consistent estimator (see, Davenport & Root, An Introduction to the theory of random signals and noise, chapter 6). In other words, as the observation interval increases to infinity the variance of the modulus (or the square of the modulus) does not decrease to zero. This is in contrast to what you would normally expect from the law of large numbers; of course, the samples in the Fourier Sum are not independent, so the LLN cannot just simply be applied.

This is the reason for real-life spectrum analyzers always use a so-called "video filter". (The term "video" here is a historical relic referring to a variable bandwidth and variable shape baseband filter.)

To estimate the spectrum one makes several independent observations of the same time series length whose reciprocal is roughly equal to the effective frequency resolution. The output of that filter, in your case the Fourier components themselves, are squared and passed through the video filter forming essentially a weighted average of these periodograms. Once properly averaged the estimator is consistent, it may be a biased estimator but it is now consistent so that longer and longer observations will lead to less and less fluctuations but you must align the IF with the video filter; a good commercial spectrum analyzer will have a large set of calibrated IF/video filter pairs built-in and will automatically suggest their use.

(Recall the so-called Gibbs phenomenon where increasing the number of Fourier terms does not help reducing the error fluctuations at discontinuities. If you want to avoid those nasty ripples you must window/taper the signal appropriately.)