Andor wrote:
> On 15 Okt., 18:16, Sampo Niskanen <spnis...@cc.hut.fi> wrote:

>> The simulation time scale, however, is quite short, so I don't
>> want a gust of wind to last the whole simulation duration (that
>> would effectively change the average wind speed). 

> Ok, I see. I think I had a similar discussion with somebody here
> some time ago who wanted to use "pink noise" for audio testing
> purposes. He was also interested in generating series with PSD
> proportional to 
> 
> 1/(1+ K f),
> 
> with a finite cutoff, going flat towards DC, and having the 1/f
> property for audible frequencies.

Another situation where a low-frequency cutoff would be desirable is 
when the generating mechanism is intermittency. Some nonlinear 
systems have an operating regime where the state grows very slowly, 
until it crosses a threshold inducing bursty output, which eventually 
returns the state to the slow region. If the duration of quiescent 
intermissions is unbounded then the system can exhibit true low-
frequency divergence.

Now, feeding the intermittent signal to an allpass filter can turn 
bursts into smoother oscillation suitable for modulating musical 
quantities, like pitch. But we don't want that viola synth to sound 
like an '80s Casio every now and then, so the intermission length 
must not much exceed the effective length of the allpass impulse 
response in this application.

I don't know if anyone has published on this in the musical context, 
but here are two articles describing intermittent iteration 
processes, one free and one not:
http://www.istia.univ-angers.fr/~chapeau/papers/psiplrc2.pdf
http://prola.aps.org/abstract/PRA/v31/i3/p1830_1


Martin

-- 
When he had finally achieved a position that
allowed him to say everything he thought, he
only thought of his position anymore.
--Gabriel Laub

Andor <andor.bariska@gmail.com> wrote:
> Ok, I see. I think I had a similar discussion with somebody here some
> time ago who wanted to use "pink noise" for audio testing purposes. He
> was also interested in generating series with PSD proportional to

> 1/(1+ K f),

> with a finite cutoff, going flat towards DC, and having the 1/f
> property for audible frequencies. Perhaps it would be a suitable
> nomenclature to call that kind (with finite cutoff frequency towards
> DC) of series "pink noise", as compared to 1/f^alpha noise.

That's true.  However, as seen on the DSP generation of Pink Noise web 
site[1] "pink noise" typically refers simply to noise with PSD 
proportional to 1/f instead of 1/f^alpha.  Therefore some distinction 
must be made between these two.  (Of course the correct way would be to 
always talk about 1/f^alpha, with alpha=1 being a special case, but this 
is currently uncommon.)

In any case, my code allows the limiting frequency to be made 
arbitrarily low (though of course this slows down the generation a bit, 
as generating a sample is O(n) for the number of poles).


[1] http://www.firstpr.com.au/dsp/pink-noise/


-- 
  __________________________________________________
 /____\   Sampo Niskanen <=> sampo.niskanen@iki.fi  \
       \       http://www.iki.fi/sampo.niskanen/     \
        \     ________________________________________\___
         \___/___________________________________________/

On 15 Okt., 18:16, Sampo Niskanen <spnis...@cc.hut.fi> wrote:
> Andor <andor.bari...@gmail.com> wrote:
> > This sounds a bit contradicting. Either you are intersted in 1/f^{5/3}
> > noise, or you are worried about values in the series deviating from
> > zero for long periods.
> > You can't have both.
>
> The application in my case was simulating wind turbulence. &#4294967295;Wind
> turbulence freqency has a power spectrum density proportional to
> 1/(1+K*f)^(5/3), which for large f is equal to 1/f^(5/3). &#4294967295;The
> simulation time scale, however, is quite short, so I don't want a gust
> of wind to last the whole simulation duration (that would effectively
> change the average wind speed).
>
> By choosing the number of poles suitably one can choose how much low
> frequency components to include. &#4294967295;The original (empirical) formula
> doesn't go to infinity at f=0 either, so one could estimate it with a
> suitable number of poles. &#4294967295;In my case I wanted even less low-frequency
> components, and using 2 poles with a 20Hz sampling rate yields maximum
> wind gust lengths of approximately 3-5 seconds (the spectrum turns flat
> below 0.3Hz). &#4294967295;So in my application, it is desireable that the noise is
> pink at the high end of the spectrum, but flat at the low end.

Ok, I see. I think I had a similar discussion with somebody here some
time ago who wanted to use "pink noise" for audio testing purposes. He
was also interested in generating series with PSD proportional to

1/(1+ K f),

with a finite cutoff, going flat towards DC, and having the 1/f
property for audible frequencies. Perhaps it would be a suitable
nomenclature to call that kind (with finite cutoff frequency towards
DC) of series "pink noise", as compared to 1/f^alpha noise. By
definition, 1/f^alpha noise has a PSD satisfying

lim_{f->0} ( P(f) / [c |f|^{-alpha}] ) = 1

for some constant c, so only the behaviour towards DC is important (as
opposed to "pink noise", where only the behaviour towards infinity is
important).

Regards,
Andor

Andor <andor.bariska@gmail.com> wrote:
> This sounds a bit contradicting. Either you are intersted in 1/f^{5/3}
> noise, or you are worried about values in the series deviating from
> zero for long periods.

> You can't have both.

The application in my case was simulating wind turbulence.  Wind 
turbulence freqency has a power spectrum density proportional to 
1/(1+K*f)^(5/3), which for large f is equal to 1/f^(5/3).  The 
simulation time scale, however, is quite short, so I don't want a gust 
of wind to last the whole simulation duration (that would effectively 
change the average wind speed).

By choosing the number of poles suitably one can choose how much low 
frequency components to include.  The original (empirical) formula 
doesn't go to infinity at f=0 either, so one could estimate it with a 
suitable number of poles.  In my case I wanted even less low-frequency 
components, and using 2 poles with a 20Hz sampling rate yields maximum 
wind gust lengths of approximately 3-5 seconds (the spectrum turns flat 
below 0.3Hz).  So in my application, it is desireable that the noise is 
pink at the high end of the spectrum, but flat at the low end.


-- 
  __________________________________________________
 /____\   Sampo Niskanen <=> sampo.niskanen@iki.fi  \
       \       http://www.iki.fi/sampo.niskanen/     \
        \     ________________________________________\___
         \___/___________________________________________/

> Note that your alpha is > 1, indicating non-stationary long
> range correlation. For that case long deviations from zero
> is the rule, not the exception.

You can consider alpha = 2 (Wiener process). In that case, the
expected time between zero crossings is infinite. 5/3 isn't much
smaller than 6/3.

Sampo Niskanen wrote:

> I haven't compared the AR and MA methods either. �I'm no DSP expert and
> once I had a working filter that produced good enough results I went
> back to my original project. �:) �In fact, I decided to use only a
> two-pole filter in my application since I didn't want the values
> deviating from zero for long periods of time.

Hi Sampo

This sounds a bit contradicting. Either you are intersted in 1/f^{5/3}
noise, or you are worried about values in the series deviating from
zero for long periods.

You can't have both.

This deviaton for longer periods from zero is one of the very
characteristics that makes 1/f^alpha noise "statistically difficult"
to handle. This phenomenon was called the "Joseph Effect" by
Mandelbrot, because of Joseph's dream that Egypt was to experience
seven years of plentiful harvests followed by seven years of famine
(long periods of deviation from the zero mean). Particularly
climatological and hydrological time series can be modeled very well
by 1/f^alpha noise. The Nile river data has an estimated Hurst
exponent of H = 0.91, which translates to an alpha of 0.82. Note that
your alpha is > 1, indicating non-stationary long range correlation.
For that case long deviations from zero is the rule, not the
exception.

Just a thought.

> Anyway, I hope the code is useful to others as well.

I'm sure it is :-).

Regards,
Andor

Thanks for all the input.


> 1) When alpha = 1, how do your results compare to the Voss-McCartney
> algorithm?

I haven't compared this method to the Voss-McCartney method.  In my 
application I required 1/alpha^(5/3) noise, which cannot be generated by 
the Voss-McCartney algorithm (or at least I couldn't figure out how to 
modify it in such a way).

I haven't compared the AR and MA methods either.  I'm no DSP expert and 
once I had a working filter that produced good enough results I went 
back to my original project.  :)  In fact, I decided to use only a 
two-pole filter in my application since I didn't want the values 
deviating from zero for long periods of time.

Anyway, I hope the code is useful to others as well.

-- 
  __________________________________________________
 /____\   Sampo Niskanen <=> sampo.niskanen@iki.fi  \
       \       http://www.iki.fi/sampo.niskanen/     \
        \     ________________________________________\___
         \___/___________________________________________/

Stacy wrote:
> 1) When alpha = 1, how do your results compare to the Voss-McCartney
> algorithm?
> 
> 2) Is this the formula for the rolloff?
>      rolloff = (alpha * -10dB/decade)?

That's how I interpret it.

> 3) How do you convert from db/decade to db/octave?
>    i.e. -10dB/decade is 3.0102999 dB/octave.

Multiply by .30102999 (.3 is close enough).

Jerry
-- 
Engineering is the art of making what you want from things you can get.
�����������������������������������������������������������������������
** Posted from http://www.teranews.com **

1) When alpha = 1, how do your results compare to the Voss-McCartney
algorithm?

2) Is this the formula for the rolloff?
     rolloff = (alpha * -10dB/decade)?

3) How do you convert from db/decade to db/octave?
   i.e. -10dB/decade is 3.0102999 dB/octave.


Thanks for posting your code.

Hi Sampo,
Thanks for the presenting the material, very interesting and useful.
Steve