2
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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... $\endgroup$ – Fat32 May 17 at 0:50
0
$\begingroup$

Given the goal stated in the second sentence of your Preamble, this is what is possible using the Corsini and Saletti algorithm:

Low frequency PSDs

More details, and the literature reference, at https://dsp.stackexchange.com/a/56820/41790. Hope this helps!

$\endgroup$
  • 1
    $\begingroup$ Thanks for the reply and the information! Can the algorithm generate noise 'on-the-fly'? Or is the method contingent on a pre-existing sequence of noise? $\endgroup$ – Winston Campeau Jun 10 at 22:33
  • $\begingroup$ The algorithm works on the fly, is fast and is actually quite easy to implement. Thanks for referencing the Kasdin paper: I had been meaning to get it, for the past month or so, and I am glad to finally have it. Of course, Kasdin does not really like ARMA-based models, for the three reasons he gives on page 817 of his paper, but the Corsini and Saletti algorithm has worked well for me since 1990. So horses for courses, as they say. $\endgroup$ – Ed V Jun 10 at 22:47
  • 1
    $\begingroup$ I went over the paper with my supervisor and we're pretty excited about this method! However, we then addressed the passage you mentioned in Kasdin's paper and we would like to clarify something with you, if you can. We will probably be iterating over 10^6 to 10^9 times -- Kasdin says that with longer sequences you must continually add poles in a Gauss-Markov model resulting in an "unwiedly" model. Do you add an additional 3 poles for every decade? And is this feasible over the number of iterations in our model? Thanks! $\endgroup$ – Winston Campeau Jun 13 at 22:33
  • $\begingroup$ 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. $\endgroup$ – Ed V Jun 13 at 23:04
  • 1
    $\begingroup$ That's wonderful! I would very, very, much appreciate a PDF of the code. My profile name without spaces or capitalized letters at gmail dot com might be one avenue. $\endgroup$ – Winston Campeau Jun 20 at 15:49
0
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.