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