| Re: red noise ? [message #47942 is a reply to message #47941] |
Sat, 11 March 2006 13:41  |
Craig Markwardt
Messages: 1869
Registered: November 1996
|
Senior Member |
|
|
"R.G. Stockwell" <no@email.please> writes:
> Does anyone have a routine that calculates a red noise random time series?
> (where red noise indicates a slope of 1/f^2 in the power spectrum)
>
> I whipped up an autoregressive routine to do this, but the slope flattens
> out
> at the higher frequencies.
>
> Pink noise would be nice as well. In fact, if someone had code to create
> a time series of random variables with any arbitrary spectral slope, that
> would
> be great!
Sure, here is a function that does that, modulo a normalization, which
you will have to diddle yourself. You enter the frequency and desired
PDS, and the output is one realization of such a power spectrum,
assuming random phases. [ Of course there are an infinite number of
realizations with the same PDS. The PDS discards 50% all phase
information so you can't go backward to a unique time series from it
alone. ]
Craig
; PDS2LC
; FF - input frequency
; PDS - input PDS
; TOTRATE - input scale factor for output time series
; SEED - (optional) random number seed
; TIME - output time samples
; RETURNS: time series sampled at TIME bins
;
function pds2lc, ff, pds, totrate, seed=seed, time=tt
df = ff(1) - ff(0)
texp = 1/df
dt = 0.5/max(ff)
npts = n_elements(ff)
cpds = sqrt(pds) * exp(2d*!dpi*dcomplex(0,randomu(seed,npts)))
cpds = [cpds, 0, reverse(conj(cpds(1:*)))]
cpds = cpds * sqrt(texp*totrate/2)
cpds(0) = 0
clc = fft(cpds,-1)
lc = (double(clc) + totrate*dt)>0
tt = n_elements(lc)*dt
return, lc
end
--
------------------------------------------------------------ --------------
Craig B. Markwardt, Ph.D. EMAIL: craigmnet@REMOVEcow.physics.wisc.edu
Astrophysics, IDL, Finance, Derivatives | Remove "net" for better response
------------------------------------------------------------ --------------
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
