# How exactly does one apply a Hanning Window for a spectral estimate?

## 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.

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.

• windows are applied in a symmetric fashion to the observation interval, the peak being at the middle. So if your signal observation is within [-T,T] then your window (continuous Hann window in this case) should be a symmetric window having its peak at "0" and tails at -T and T... – Fat32 May 17 '19 at 0:50

• Part 1: Glad to hear you are excited about the method and good news: you do not have to keep adding poles! I generate $10^6$ or $10^7$ consecutive values (to avoid any possible 'period exhaustion' of my simulation program' pseudo-random number generator). I specify, at the beginning, the number of poles/decade and number of decades. My sim program (ExtendSim) does a maximum FFT of 32k, so I specify 4 decades, from Nyquist down to about the frequency resolution element determined by the FFT length. But I could specify more decades, if I wanted to, and it is all done up front. – Ed V Jun 13 '19 at 23:04
Winston and I believe that the observation window is from 0 to T in two variables, $$t$$ and $$\beta$$. The integral equation Winston cites is after a transformation to $$\tau = \beta-t$$ (see page 147 in "Random Data" by Bendat and Piersol). One of the things which makes us suspicious is that the result of the integral in Sage has terms with $$i$$ in them, which cancelled in the simpler, rectangular window case. Thus we believe we are doing something wrong.