simulation

LastWave

#
# Procedure to generate a gaussian process
#

setproc gaussprocess {{&signali corr} {&int size 4096} {&int flagZero 0}} \
"{{{<auto-correlation> [<size>=4096] [<flagZero>]} {Generates a realization of size <size> of the stationnary gaussian random process (of mean 0) \
defined by the autocorrelation function <auto-correlation>. The closest power of 2 greater than <size> must be smaller than the size of the <auto-correlation> signal UNLESS \
<flagZero> is 1 in which case the autocorrelation signal is padded with 0's. \
The <auto-correlation> signal \
must code only the positive lags (starting from lag 0, i.e., coding the variance). The algorithm is optimized for <size> which are power of 2. \
The returned realization has the same dx as the <auto-correlation>}}}" \
{

  if (size <= 1) {errorf "Parameter <size> must be strictly greater than 1"}

  # This algorithm builds a random signal of size m+1 using a circulant matrix 
  # of size M = 2m : <r(0),...,r(m),r(m-1),...,r(1)
  # where r is the autocorrelation function
  # M must be a power of 2 so m must be a power of 2
  m = size-1
  m = 2^ceil(log2(m))
  M = 2*m
  if (!flagZero && corr.size <= m) {
    errorf "Sorry the size %d of the autocorrelation signal must be strictly greater than %d" corr.size m
  }
  
  if (flagZero && corr.size <= m) {
    corr = <corr,Zero(m+1-corr.size)>
  }

  # Building the circulant vector of size 2M
  thecorr = <corr[0:m],corr[m-1:-1:1]>

  # The fourier transform of the circulant vector
  fftrcorr = [new &signal]
  ffticorr = [new &signal]
  fft thecorr Zero(thecorr.size) fftrcorr ffticorr

  # Sqrt
  i = find(fftrcorr<0 && fftrcorr>-1e-15)
  if (i.size != 0) {
    fftrcorr[i] := 0
  }
  fftrcorr = sqrt(fftrcorr)
  ffticorr[:] := 0
  
  # Build a complex white noise on [0,2pi]
  fftrcorr = sqrt(2*M)*fftrcorr*Grand

  # Perform the inverse fourier transform
  corr1 = <>
  fft fftrcorr ffticorr corr corr1 -i
  
  # Returns the extracted signal
  return corr[0:!size]
}


#
# Procedure to compute the mean/autocovariance of the omega_l gaussian process which is used to build the MRM measure
#
setproc _mrwomcorr {{&int size} {&float dt .1} {&float T 200} {&float lambda2 0.02}} \
"{{{[<dt>=.1] [<T>=200] [<lambda2>=0.02] [<sigma>=1>]} {Returns a listv made up of the mean value and the correlation function of the omega process of the mrw of parameters \
T and lambda2 on <size> points}}}" \
{
  xx = dt*<1:size>
  cc = lambda2*ln(T/xx)*(xx<T)
  cc = <lambda2*(1+ln(T/dt)),cc>
  cc.dx = dt
  cc.x0 = 0
  c = cc  
  m = -c[0]
  return {m c}
}

#
# Procedure to compute the mean/autocovariance of the aggregated gaussian process which is used to build the MRM measure
#
setproc _mrwomcorrcheat {{&int size} {&float dt .1} {&float T 200} {&float lambda2 0.02}} \
"{{{[<dt>=.1] [<T>=200] [<lambda2>=0.02] [<sigma>=1>]} {Returns a listv made up of the mean value and the correlation function of the omega process of the mrw of parameters \
T and lambda2 on <size> points}}}" \
{
  xx = <2:size>
  cc = (xx<T && xx<50)*(ln(1.0*T/xx)+1.5-0.5*(xx-1)*(xx-1)*ln(1-1.0/xx)-0.5*(xx+1)*(xx+1)*ln(1+1.0/xx)) + (xx<T && xx>=50)*(ln(1.0*T/xx)+1.0/(12.0*xx*xx)+1.0/(60.0*xx*xx*xx*xx))
  cc = lambda2*<ln(T)+1.5,ln(T)+1.5-2*ln(2),cc>
  cc.x0 = 0
  c = cc
  m = -c[0]
  return {m c}
}

#
# Procedure to buil a MRW/MRM
#
setproc mrw {{&int size 4096} {&int subsample 4} {&float T 200} {&float lambda2 0.02} {&float sigma 1} {&wordlist .l}}  \
"{{{[<size>=4096] [<subsample>=4] [<T>=200] [<lambda2>=0.02] [<sigma>=1] [{-m | -M}]} {Generates realization of size <size> (optimized for \
size which are power of 2) of a MRW process (or the correponding MRMmeasure if -m is set or both in a listv if -M) of parameter <T>, <lambda2> and <sigma>. The realization is computed using a time step of 1. \
The (going to 0) parameter l is 1/subsample.}}}" \
{
  flagM = 0
  flagm = 0
  if (l.llength != 0) {
    if (l.llength != 1) {errorf "bad options"}
    if (l[*list,0] == '-m') {
      flagm = 1
    } elseif (l[*list,0] == '-M') {
      flagM = 1
    } else {
      errorf "bad option '%s'" l[*list,0]
    }
  }

  dt = 1/subsample
  N = size-1
  N = 2^ceil(log2(N))
  {m corr} = [_mrwomcorr N*subsample dt T lambda2]

  om = [gaussprocess corr size*subsample]+m

  if (flagm || flagM) {
    mrm = prim(sigma*sigma*exp(2*om)*dt)
  }  
  if (!flagm || flagM) {
    mrw = prim(exp(om)*Grand*sigma*sqrt(dt))
  }

  if (flagm) {
    return mrm[:subsample:]
  }

  if (flagM) {
    return {mrm[:subsample:] mrw[:subsample:]}
  }
  
    return mrw[:subsample:]  
}