0

I am writing an app to work with synthetic time series data from a physics experiment. In our experiments we always have $1/f$ noise in our time series, but I haven't been able to find code/packages to generate $1/f$ noise in synthetic time series data. Can anyone tell me how to do this? I have a feeling it should be pretty simple and that I'm missing something obvious, but I don't know what it is. Any advice?

Also, I looked through the list of questions and don't see anything too on point apart from this: making pink noise (1/f) using list of frequencies

Community
  • 1
sunny
  • 1
  • 1

1 Answers1

1

this is an old resource, but good. in fact, i have a very old contribution (another 3 decade old design) that i will repeat here at SE:

A method that Orfanidis (Introduction to Signal Processing mentions came from an old comp.dsp post of mine. (here's a pdf, check Problem B.9) It's just a simple "pinking" filter to be applied to white noise. since the rolloff is -3 dB/octave, -6 dB/octave (1st order pole) is too steep and 0 dB/octave is too shallow.

An equiripple approximation to the ideal pinking filter can be realized by alternating real poles with real zeros. A simple 3rd order solution that i obtained is:

$$ H(z) = \frac{(z-q_1)(z-q_2)(z-q_3)}{(z-p_1)(z-p_2)(z-p_3)} $$

where

$p_1$ = 0.99572754

$p_2$ = 0.94790649

$p_3$ = 0.53567505

$q_1$ = 0.98443604

$q_2$ = 0.83392334

$q_3$ = 0.07568359

The response follows the ideal -3 dB/octave curve to within $\pm$0.3 dB error over a 10 octave range from 0.0009$\times \pi$ to 0.9$\times \pi$. ("Nyquist" = $\pi$.)

robert bristow-johnson
  • 20,661
  • 4
  • 38
  • 76
  • Do you know this filter? Its poles and zeros are very close to yours. Not sure how J.O. Smith came up with that filter. How did you optimize the pole and zero locations? – Matt L. Apr 09 '15 at 09:02
  • R B-J's filter worked great for me in a quick and dirty implementation:

    gain = 1.0/23.1188232529392; out[t] = in[t]gain + in[t-1](-1.89404297gain) + in[t-2](0.9585641562817477gain) + in[t-3](-0.0621320035261672gain) - out[t-1]-2.47930908 - out[t-2]1.985012853639686 - out[t-3]-0.505600430025288;

     – Olli Niemitalo Apr 09 '15 at 15:09
  • @MattL., it's possible that JOS got it from me. or independently. i have never "collaborated" with JOS, but i have corresponded with him dozens of times. usually he credits accurately, but i just don't know where he got this. if it's third order and good for 10 octaves (20 to 20kHz), then there's about only one set of real and alternating poles and zeros that does it. when i first did this (in the 80's) i took the -3 dB straight slope and prewarped it (regarding bilinear transform). then i fit asymptotes to it with either -6 dB slope or 0 dB slope to get the analog corner frequencies. – robert bristow-johnson Apr 09 '15 at 15:52
  • one thing that Sophicles points out is that, without edge effects, the frequencies of the alternating poles and zeros are equally spaced in log-frequency. so that might be a good first guess. and then just manually tweek some of the side corner frequencies. if i were to do it over again, i would probably do it with 5 poles and 4 zeros. – robert bristow-johnson Apr 09 '15 at 15:58
  • @robertbristow-johnson: OK, thanks. Do you have a reference for the poles and zeros being equally spaced in log-frequency? – Matt L. Apr 09 '15 at 16:45
  • i think the Orfanidis book. but i can't remember. but isn't it apparent, @MattL.? consider a log-gain vs. log-frequency plot of this filter. it will be -3.01 dB/octave slope forever on both the left and the right. so now think about the corners of a continuous piecewise-linear approximation to that where you are constrained to -6.02 dB/octave slope or 0 slope. those corners would have to be equally spaced in log-frequency. – robert bristow-johnson Apr 09 '15 at 17:43
  • @MattL. i found my original reference to this from 2 decades ago. but i originally worked this out about a decade earlier before the days that i saw any email or internet. i dunno if i had the numbers in a note or if it was on a Mac 3.5 inch floppy. – robert bristow-johnson Apr 09 '15 at 18:02
  • @robertbristow-johnson: Thanks for the reference. And you're right, equally spacing on a log frequency scale makes sense after all. If I have time I'd like to try to numerically optimize the pole/zero locations for arbitrary filter orders. – Matt L. Apr 10 '15 at 07:11
  • @MattL., that would be nice. as i mentioned, if i were to do it again, i would make it 5 real poles with 4 real zeros wedged in between. i think you can start out with equal spacing in log-frequency and adjust the corners (whether they be poles or zeros) at the left and right sides to deal with the edge effects. – robert bristow-johnson Apr 10 '15 at 08:50