Generating 1/f (flicker/pink) noise

I wrote:

> I thought that part in the paper was a bit odd.

I also like this quote from the paper:

"... for large enough values of f_min (no more than 2 or 3 decades
from 0) ... "

:-)

>On 25 Mar, 14:51, dbd <d...@ieee.org> wrote:
>> On Mar 24, 4:51 pm, maxlittle2...@googlemail.com wrote:
>>
>>
>>
>> > On 23 Mar, 18:33, dbd <d...@ieee.org> wrote:
>>
>> > > On Mar 22, 4:03 pm, maxlittle2...@googlemail.com wrote:
>>
>> > > > ...
>> > > > Ah, well I did find those resources, but no simple "one pager"
of
>> > > > Matlab code for generating perfect 1/f^alpha noise.
>> > > > ...
>> > > > Max
>>
>> > > There can be no page that is 'simple' and 'perfect' for all
>> > > applications.
>>
>> > > One characteristic of your implementation choice is that the power
>> > > spectrum does not just have the desired expectation value in any
>> > > frequency region, but it equals the expectation value. While this
is a
>> > > possible power spectum value, it is not a likely one, and is
unlikely
>> > > to repeat. This can be useful in certain testing situations. It
might
>> > > be inappropriate in others. An alternative is to use independent
>> > > Gaussian samples for real and imaginary frequency components and
>> > > weight them with the desired frequency weighting. An example
>> > > implementation of this, modified from your implementation, follows
>> > > below. This could provide more realistic statistical behavior in
some
>> > > applications.
>>
>> > > Dale B. Dalrymplehttp://dbdimages.com
>>
>> > Dale
>>
>> > Thanks, these are excellent points. I've added an option to the
>> > routine so that a choice can be made between deterministic and
>> > stochastic power.  I've also added an option that normalises the
>> > signal amplitude so that it lies in the range -1 to +1. This may be
>> > useful in certain circumstances. Grab the code from:
>>
>> >http://www.eng.ox.ac.uk/samp/powernoise_soft.html
>>
>> > Max
>>
>> Max
>>
>> If you really think a normalization and a peak limiting capability
>> should be included inside the function, I suggest you do so with a
>> gain and a clipping capability instead of putting the mean and scale
>> of the output under the control of the extremal behavior of the noise
>> generated.
>>
>> There are test purposes where the original function or my addition
>> could be called multiple times to generate longer sequences of test
>> data. Even these special cases no longer work when gain and offset
>> change from block to block.
>>
>> Dale B. Dalrymple
>
>Thanks for the thoughtful suggestions Dale.
>
>I think the normalization to the range [-1, 1] is required for some
>applications (particularly nonlinear time series analysis): at least
>one member of this group has found that some gain is desireable, and I
>reason that it's better to have a signal that uses a standard number
>range. Otherwise, as you imply, there probably isn't much point in
>putting this scaling in the function. If by clipping I understand you
>to mean some kind of nonlinear "clamping" of values outside a certain
>range, I don't agree with this because it will distort the spectrum.
>
>This also leads on to your second paragraph. Unfortunately this
>method, as it stands, can't be called multiple times to generate
>longer sequences: there will be several problems. One of them is that
>the signal is not guaranteed to be continuous across multiple calls.
>The resulting discontinuities will introduce spectral distortions.
>
>Of course, one could imagine a solution to fix this problem that
>matches the endpoints of each successive call, but then you'd also
>need to match as many derivatives of the signal as are required to
>ensure the appropriate differentiability for the chosen "roughness" of
>the signal (which, informally, is inversely proportional to alpha). I
>don't think such an approach would be elegant though, and I'm sure
>there's a better one, although I can't think of it right now.
>
>Max
>

Hi Max & Dale,
I have been trying to use your code to create a background noise signal
with 1/f characteristics (between [-1 1]). That is fine and moves very
efficiently. Now, I need to create a signal (1/f again) which is not
normalized between [-1 1] but rather has an average gain defined by factor
K. How do I do that?

Simply adding or multiplying the background noise by the factor of K does
not provide the same 1/f characteristics. What is the right solution?
Best
Nim

>>On 25 Mar, 14:51, dbd <d...@ieee.org> wrote:
>>> On Mar 24, 4:51 pm, maxlittle2...@googlemail.com wrote:
>>>
>>>
>>>
>>> > On 23 Mar, 18:33, dbd <d...@ieee.org> wrote:
>>>
>>> > > On Mar 22, 4:03 pm, maxlittle2...@googlemail.com wrote:
>>>
>>> > > > ...
>>> > > > Ah, well I did find those resources, but no simple "one pager"
>of
>>> > > > Matlab code for generating perfect 1/f^alpha noise.
>>> > > > ...
>>> > > > Max
>>>
>>> > > There can be no page that is 'simple' and 'perfect' for all
>>> > > applications.
>>>
>>> > > One characteristic of your implementation choice is that the
power
>>> > > spectrum does not just have the desired expectation value in any
>>> > > frequency region, but it equals the expectation value. While this
>is a
>>> > > possible power spectum value, it is not a likely one, and is
>unlikely
>>> > > to repeat. This can be useful in certain testing situations. It
>might
>>> > > be inappropriate in others. An alternative is to use independent
>>> > > Gaussian samples for real and imaginary frequency components and
>>> > > weight them with the desired frequency weighting. An example
>>> > > implementation of this, modified from your implementation,
follows
>>> > > below. This could provide more realistic statistical behavior in
>some
>>> > > applications.
>>>
>>> > > Dale B. Dalrymplehttp://dbdimages.com
>>>
>>> > Dale
>>>
>>> > Thanks, these are excellent points. I've added an option to the
>>> > routine so that a choice can be made between deterministic and
>>> > stochastic power.  I've also added an option that normalises the
>>> > signal amplitude so that it lies in the range -1 to +1. This may be
>>> > useful in certain circumstances. Grab the code from:
>>>
>>> >http://www.eng.ox.ac.uk/samp/powernoise_soft.html
>>>
>>> > Max
>>>
>>> Max
>>>
>>> If you really think a normalization and a peak limiting capability
>>> should be included inside the function, I suggest you do so with a
>>> gain and a clipping capability instead of putting the mean and scale
>>> of the output under the control of the extremal behavior of the noise
>>> generated.
>>>
>>> There are test purposes where the original function or my addition
>>> could be called multiple times to generate longer sequences of test
>>> data. Even these special cases no longer work when gain and offset
>>> change from block to block.
>>>
>>> Dale B. Dalrymple
>>
>>Thanks for the thoughtful suggestions Dale.
>>
>>I think the normalization to the range [-1, 1] is required for some
>>applications (particularly nonlinear time series analysis): at least
>>one member of this group has found that some gain is desireable, and I
>>reason that it's better to have a signal that uses a standard number
>>range. Otherwise, as you imply, there probably isn't much point in
>>putting this scaling in the function. If by clipping I understand you
>>to mean some kind of nonlinear "clamping" of values outside a certain
>>range, I don't agree with this because it will distort the spectrum.
>>
>>This also leads on to your second paragraph. Unfortunately this
>>method, as it stands, can't be called multiple times to generate
>>longer sequences: there will be several problems. One of them is that
>>the signal is not guaranteed to be continuous across multiple calls.
>>The resulting discontinuities will introduce spectral distortions.
>>
>>Of course, one could imagine a solution to fix this problem that
>>matches the endpoints of each successive call, but then you'd also
>>need to match as many derivatives of the signal as are required to
>>ensure the appropriate differentiability for the chosen "roughness" of
>>the signal (which, informally, is inversely proportional to alpha). I
>>don't think such an approach would be elegant though, and I'm sure
>>there's a better one, although I can't think of it right now.
>>
>>Max
>>
>
>Hi Max & Dale,
>I have been trying to use your code to create a background noise signal
>with 1/f characteristics (between [-1 1]). That is fine and moves very
>efficiently. Now, I need to create a signal (1/f again) which is not
>normalized between [-1 1] but rather has an average gain defined by
factor
>K. How do I do that?
>
>Simply adding or multiplying the background noise by the factor of K
does
>not provide the same 1/f characteristics. What is the right solution?
>Best
>Nim
>
>
>
>

Actually never mind,
I changed the code to adapt that:

function x = powernoise(alpha,N, varargin) 
% function x = powernoise(alpha,N, 'key','value') 
% Generate samples of power law noise. The power spectrum
% of the signal scales as f^(-alpha).
%
% Useage:
%  x = powernoise(alpha, N)
%  x = powernoise(alpha, N, 
% 'randpower','off','normalize','on','normfactor',n)
%
% Inputs:
%  alpha - power law scaling exponent
%  N     - number of samples to generate
%
% Output:
%  x     - N x 1 vector of power law samples
%
% With no option strings specified, the power spectrum is
% deterministic, and the phases are uniformly distributed in the range
% -pi to +pi. The power law extends all the way down to 0Hz (DC)
% component. By specifying the 'randpower' option string however, the
% power spectrum will be stochastic with Chi-square distribution. The
% 'normalize' option string forces scaling of the output to the range
% [-1, 1], consequently the power law will not necessarily extend
% right down to 0Hz.
% 
% if you decide to normalize, you can also add normfactor option:
% 'normfactor',k will normalize the output between [-k,k]
%
% (cc) Max Little, 2008. This software is licensed under the
% Attribution-Share Alike 2.5 Generic Creative Commons license:
% http://creativecommons.org/licenses/by-sa/2.5/
% If you use this work, please cite:
% Little MA et al. (2007), "Exploiting nonlinear recurrence and fractal
% scaling properties for voice disorder detection", Biomed Eng Online,
6:23
%
% As of 20080323 markup
% If you use this work, consider saying hi on comp.dsp
% Dale B. Dalrymple 
% 
% Modified by Nima Dehghani : added normfactor to scale between
%[-normfactor normfactor]
% Sep 3rd 2009

callbacks = {'normalize','randpower','normfactor'};

  try
    options = varargin;
    for index = 1:length(options)
        if iscell(options{index}) & ~iscell(options{index}{1}),
options{index} = { options{index} }; end;
    end;
    if ~isempty( varargin ), g=struct(options{:}); 
    else g= []; end;
  catch
    fprintf('%s: calling convention {''key'', value, ... }
error',mfilename);
    return;
  end;

try, g.normalize;               catch; g.normalize = 'off'; end;
try, g.randpower;               catch; g.randpower = 'off'; end;
try, g.normfactor;              catch; g.normfactor = 1;    end;


if strcmpi(g.normalize, 'on'),   opt_normal = true;   else optnormal =
false;   end;
if strcmpi(g.normalize, 'on'),   opt_randpow = true;  else opt_randpow =
false; end;



N2 = floor(N/2)-1;
f = (2:(N2+1))';
A2 = 1./(f.^(alpha/2));

if (~opt_randpow)
    p2 = (rand(N2,1)-0.5)*2*pi;
    d2 = A2.*exp(i*p2);
else
    % 20080323
    p2 = randn(N2,1) + i * randn(N2,1);
    d2 = A2.*p2;
end

d = [1; d2; 1/((N2+2)^alpha); flipud(conj(d2))];
x = real(ifft(d));

if (opt_normal) 
    x = ((x - min(x))/(max(x) - min(x)) - 0.5) * g.normfactor*2;
end