Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

How to simulate pink noise in R

I know that white noise can be achieved by treating the output of rnorm() as a timeseries. Any suggestions on how to simulate pink noise?

like image 920
Mike Lawrence Avatar asked Jan 02 '12 05:01

Mike Lawrence


1 Answers

As said by @mbq, you can just use p@left to get the vector, instead of saving and reading the wav file. On the other hand, you could directly use the function generating the time serie in tuneR:

TK95 <- function(N, alpha = 1){ 
    f <- seq(from=0, to=pi, length.out=(N/2+1))[-c(1,(N/2+1))] # Fourier frequencies
    f_ <- 1 / f^alpha # Power law
    RW <- sqrt(0.5*f_) * rnorm(N/2-1) # for the real part
    IW <- sqrt(0.5*f_) * rnorm(N/2-1) # for the imaginary part
    fR <- complex(real = c(rnorm(1), RW, rnorm(1), RW[(N/2-1):1]), 
                  imaginary = c(0, IW, 0, -IW[(N/2-1):1]), length.out=N)
     # Those complex numbers that are to be back transformed for Fourier Frequencies 0, 2pi/N, 2*2pi/N, ..., pi, ..., 2pi-1/N 
     # Choose in a way that frequencies are complex-conjugated and symmetric around pi 
     # 0 and pi do not need an imaginary part
    reihe <- fft(fR, inverse=TRUE) # go back into time domain
    return(Re(reihe)) # imaginary part is 0
}

and this works perfectly :

par(mfrow=c(3,1))
replicate(3,plot(TK95(1000,1),type="l",ylab="",xlab="time"))

enter image description here

like image 170
Simon C. Avatar answered Oct 26 '22 00:10

Simon C.