PREAMBLE
I am a graduate student researching evolutionary ecology. My supervisor and I have been trying to learn how to simulate colored noise, discretely, by the $\frac{1}{f^\alpha}$ power law for $\alpha \in \mathbb{R}$, specifically on the interval [0,2].
We found an article that seems to cover a method in great detail.
Link here: Discrete simulation of colored noise and stochastic processes and $1/f^\alpha$ power law noise generation (Kasdin, 1995)
On page 806, we were able to prove to ourselves that by a symmetric autocorrelation function, the formula for the sampled spectrum of an arbitrary segment of a process is indeed,
$$\hat{S}(\omega)\triangleq E\{\tilde{S}(\omega)\}=\int_{-T}^{0}\left(1+\dfrac{\tau}{T}\right) \Big[ \dfrac{1}{T} \int_{t_o+\tau/2}^{t_o+T+\tau/2} R(t, \tau)\Big] e^{-j\omega\tau}d\tau + \int_{0}^{T}\left(1-\dfrac{\tau}{T}\right) \Big[ \dfrac{1}{T} \int_{t_o+\tau/2}^{t_o+T+\tau/2}R(t, \tau) \Big] e^{-j\omega\tau}d\tau$$
$$=\int_{-T}^{T}\left( 1-\dfrac{|\tau|}{T}\right)R(t,\tau)e^{-j\omega\tau}d\tau$$
Where,
$$R(t,\tau)\triangleq \{x(t+\frac{\tau}{2})x(t-\frac{\tau}{2})$$
Next, the paper computed the spectral estimate for Brownian motion. That estimate in the paper is given as,
$$\hat{S}_B(\omega)=2(1+\dfrac{t_o}{T})\dfrac{1}{\omega^2} - 2(\dfrac{t_o}{T})\dfrac{cos\omega T}{\omega^2} - \dfrac{2}{T\omega^3}sin\omega T$$
Where for $T$ much larger than $t_o$, this reduces to $\frac{2}{\omega^2}$. We were able to reproduce this result.
QUESTION AND PROBLEM
The paper now makes the claim that we have unknowingly been using a rectangular window. We are told we can improve the spectral estimate by using a Hanning window with normalization $8/3$ (After some reading, I've been able to grasp the concept of spectral leakage).
The new spectral estimate for Brownian Motion, with $t_o$ set to 0 is,
$$\hat{S}_B=\dfrac{1}{3}\{ \dfrac{4}{\omega^2} - \dfrac{4sin \omega T}{T \omega^3} + \dfrac{80\pi^4 - 24\pi^2T^2\omega^2 + T^4\omega^4}{T^4(4\pi^2/T^2-\omega^2)^3)} + \dfrac{48\pi^2T\omega sin \omega T - 4T^3\omega^3sin\omega T}{T^4(4\pi^2/T^2-\omega^2)^3)} \}$$
Try as we might, we have not yet been able to reproduce this result. We tried the following,
$$\int_{-T}^{0}\left(1+\dfrac{\tau}{T}\right) \Big[ \dfrac{1}{T} \int_{t_o+\tau/2}^{t_o+T+\tau/2} R(t, \tau)\Big] e^{-j\omega\tau}d\tau \cdot HW(x=t)\cdot HW(x=(t+\tau)) + \int_{0}^{T}\left(1-\dfrac{\tau}{T}\right) \Big[ \dfrac{1}{T} \int_{t_o+\tau/2}^{t_o+T+\tau/2}R(t, \tau) \Big] e^{-j\omega\tau}d\tau \cdot HW(x=t)\cdot HW(x=(t+\tau))$$
Where we defined the Hanning Window (HW) as,
$$1-cos((\pi x)/T)^2$$
Done with the following code in Sage:
Q, T, t0, w, tau, t, x = var('Q T t0 w tau t x')
hanning_window = (sqrt(8/3))*(1 - cos((pi*x)/T)^2)
rectangular_window = 1
assume(T>0)
window_1 = hanning_window
window_2 = hanning_window
fneghann = Q*(t+tau/2)*exp(-I*w*tau)*window_1.substitute(x=(t))*window_2.substitute(x=(t+tau))
fposhann = Q*(t-tau/2)*exp(-I*w*tau)*window_1.substitute(x=(t))*window_2.substitute(x=(t+tau))
result = (1/T)*((fneghann.integrate(t,t0-tau/2,t0+T+tau/2)).integrate(tau,-T,0) + (fposhann.integrate(t,t0+tau/2,t0+T-tau/2)).integrate(tau,0,T))
f = (result).expand().full_simplify()
view(f)
Computing the limit as $T\rightarrow \infty$ and setting $t_o = 0$ DOES produce the spectral estimate of $\frac{1}{\omega^2}$ as indicated in the linked paper, but we got something dissimilar to the initial spectral estimate submitted in the article, leading us to believe the result of our limit was luck. I've chosen to not post the result, $f$, because of it's size.
Where did we make a mistake, if any? Did we apply the Hanning Window incorrectly? Is there a mistake in our code? Your help is greatly appreciated.
