#########################################################################
#       $Log: Util.S,v $
#
#               (c) Copyright  1997
#                          by                                   
#      Author: Rene Carmona,Bruno Torresani,Wen-Liang Hwang,Andra Wang   
#                  Princeton University
#                  All right reserved                           
#########################################################################





smoothwt <- function(modulus, subrate, flag = FALSE)
#########################################################################
#      smoothwt:   
#      --------
#       smooth the wavelet (or Gabor) transform in the b direction
#
#       input:
#       ------
#	 modulus: modulus of the wavelet (or Gabor) transform
#	 subrate: size of the smoothing window
#	 flag: if set to TRUE, subsamples the smoothed transform by
#		subrate.	
#
#       output:
#       -------
#        smoothed transform (2D array)
#
#########################################################################
{
  sigsize <- dim(modulus)[1]
  nscale <- dim(modulus)[2]

  if(flag) {
    smodulus <- matrix(0,sigsize,nscale)
    dim(smodulus) <- c(sigsize * nscale, 1)
  }
  else {
    modsize <- as.integer(sigsize/subrate) + 1
    smodulus <- matrix(0,modsize,nscale)
    dim(smodulus) <- c(modsize * nscale, 1)
  }
	
  dim(modulus) <- c(sigsize * nscale, 1)

  z <- .C("Ssmoothwt",
        sm = as.double(smodulus),
        as.double(modulus),
        as.integer(sigsize),
        as.integer(nscale),
        as.integer(subrate),
	as.integer(flag))

  if(flag)
   dim(z$sm) <- c(sigsize, nscale)
  else 	
   dim(z$sm) <- c(modsize, nscale)
  z$sm
}




smoothts <- function(ts, windowsize)
#########################################################################
# smooth a time series by averaging window
#
# Arguments: 
#	<ts> time series
#       <windowsize> window size
#########################################################################
{
   sigsize <- length(ts)
   sts <- numeric(sigsize)
   nscale <- 1   
   stssize <- sigsize
   dim(sts) <- c(sigsize, 1)
   dim(ts) <- c(sigsize, 1)
   subrate <- 1
   
   z <- .C("Smodulus_smoothing",
        as.double(ts),
        sts = as.double(sts),
	as.integer(sigsize),
        stssize = as.integer(stssize),            
        as.integer(nscale),
        as.integer(subrate))

  sts <- z$sts
  sts
}

adjust.length <- function(inputdata)
#########################################################################
# adjust.length adds zeros to the end of the data if necessary so that 
# its length is a power of 2.  It returns the data with zeros added if
# nessary and the length of the adjusted data.
#
# Arguments: 
#	<inputdata> is either a text file or an S object containing
# 		    data.
#########################################################################
{
  if (is.character(inputdata))
    s <- scan(inputdata)
  else
    s <- inputdata
  np <- length(s)

  pow <- 1
  while ( 2*pow < np )
    pow <- 2*pow
  new.np <- 2*pow

  if ( np == new.np )
    list( signal=s, length=np )
  else
  {
    new.s <- 1:new.np
    new.s[1:new.np] <- 0
    new.s[1:np] <- s

    list( signal=new.s, length=new.np )
  }
}



#########################################################################
#      npl:
#      ---
#       plots with n rows.
#
#       input:
#       ------
#	 nrow: number of rows
#
#########################################################################
npl <- function(nbrow)
{
        par(mfrow = c(nbrow, 1))
        cat("")
}



#########################################################################
#      SpecGen:
#      --------
#       Estimate power spectrum
#
#       input:
#       ------
#	 input: input signal.
#        spans: length of Daniell's smoother.
#        taper: relative size of tapering window (cosine taper).
#
#########################################################################
SpecGen <- function(input,spans = 9,taper = 0.2)
{
   tmp <- spec.pgram(input,spans,taper)
   tmp1 <- tmp$spec
   tmp <- 10^(tmp1/10)
   spec <- tmp
   spec
}



#########################################################################
#      SampleGen:
#      ----------
#       generate a real sample with prescribed spectrum      
#
#       input:
#       ------
#	 nspec: spectral density
#
#########################################################################
SampleGen <- function(spec)
{
   sqspec <- sqrt(spec)
   d <- length(sqspec)
   size <- d 
   tmp <- rnorm(2 * size)
    real <- numeric(size)
   imag <- numeric(size)
   real <- tmp[1:size]
   a <- size+1
   b <- 2 * size
   imag <- tmp[a:b]
   real1 <- real * sqspec[1:size]
   imag1 <- imag * sqspec[1:size]
   real2 <- numeric(2 * size)
   imag2 <- numeric(2 * size)
   real2[1:size] <- real1
   imag2[1:size] <- imag1
   for(i in 2:size) {
      real2[b-i+2] <- real1[i]
      imag2[b-i+2] <- -imag1[i]
   }	
   imag2[1] <- 0
   i <- sqrt(as.complex(-1))
   tmp1 <-  fft(real2 + imag2 * i, inverse=T)
   sample <- Re(tmp1)/sqrt((2*size))
   sample
}



hurst.est <- function(wspec,range,nvoice,plot=T)
#########################################################################
#      hurst.est:   
#      ----------
#       estimate Hurst exponent from (part of) wavelet spectrum
#
#       input:
#       ------
#	 wspec: wavelet spectrum.
#	 nvoice: number of voices.
#	 plot: if set to T, displays regression line on current plot.
#
#       output:
#       -------
#        intercept and slope of regression line. Slope is Hurst exponent.
#
#########################################################################
{
  loctmp <- lsfit(range,log(wspec[range],base=2^(2/nvoice)))
  if(plot) abline(loctmp)
  loctmp$coef
}




wspec.pl <- function(wspec,nvoice)
#########################################################################
#      wspec.pl
#      ----------
#       Displays normalized log of wavelet spectrum
#
#       input:
#       ------
#	 wspec: wavelet spectrum.
#	 nvoice: number of voices.
#
#########################################################################
{
  plot.ts(log(wspec,base=2^(2/nvoice)))
  title("log(wavelet spectrum)", xlab="log(scale)", ylab="V(a)")
  cat(" ")
}








#########################################################################
#      $Log: Crazy_Climbers.S,v $
#########################################################################
#
#               (c) Copyright  1997
#                          by                                   
#      Author: Rene Carmona, Bruno Torresani, Wen-Liang Hwang   
#                  Princeton University
#                  All right reserved                           
#########################################################################






crc <- function(tfrep, tfspec = numeric(dim(tfrep)[2]), bstep = 3,
	iteration = 10000, rate = .001, seed = -7, nbclimb=10,
	flag.int = TRUE, chain = TRUE, flag.temp = FALSE)
#########################################################################
#  crc:   Time-frequency multiple ridge estimation (crazy climbers)
#  ----
# 	use the crazy climber algorithm to evaluate ridges of
#    	   continuous Time-Frequency transform
#
#      input:
#      ------
# 	tfrep: wavelet or Gabor transform
#	tfspec: additional potential (coming from learning the noise)
#       iteration: number of iterations
#       rate: initial value of the temperature
#       seed: initialization for random numbers generator
#       nbclimb: number of crazy climbers
#       bstep: step size for the climber walk in the time direction
#	flag.int: if set to TRUE, computes the integral on the ridge.
#	chain: if set to TRUE, chains the ridges.
#	flag.temp: if set to TRUE, keeps a constant temperature.
#
#      output:
#      -------
#       beemap: 2D array containing the (weighted or unweighted)
#               occupation measure (integrated with respect to time)
#
#########################################################################
{

   tfspectrum <- tfspec 

   d <- dim(tfrep)
   sigsize <- d[1]
   nscale <- d[2]
   beemap <- matrix(0,sigsize,nscale)
   sqmodulus <- Re(tfrep*Conj(tfrep))
   for (k in 1:nscale)
     sqmodulus[,k] <- sqmodulus[,k] - tfspectrum[k]
   dim(beemap) <- c(nscale * sigsize, 1)
   dim(sqmodulus) <- c(nscale * sigsize, 1)
	
   
   z <- .C("Sbee_annealing",
        as.double(sqmodulus),
        beemap= as.double(beemap),
        as.single(rate),
	as.integer(sigsize),
        as.integer(nscale),
        as.integer(iteration),
	as.integer(seed),
        as.integer(bstep),
        as.integer(nbclimb),
	as.integer(flag.int),
	as.integer(chain),
	as.integer(flag.temp))

   beemap <- z$beemap
   dim(beemap) <- c(sigsize,nscale)   
   if(dev.set() != 1) image(beemap)
   beemap
}





cfamily <- function(ccridge, bstep = 1, nbchain = 100, ptile = 0.05)
#########################################################################
#     cfamily:
#     --------
#     chain the ridges obtained by crazy climber.
#
#      input:
#      ------
#       ccridge: unchained ridge (output of Bee_Annealing)
#       bstep: maximal length for a gap in a ridge
#       nbchain: maximum number of chains
#       ptile: relative threshold for the ridges 
#
#      output:
#      ------
#       ordered: ordered map: image containing the ridges
#                (displayed with different colors)
#       chain: 2D array containing the chained ridges, according
#		to the chain data structure:
#		chain[,1]: first point of the ridge
#		chain[,2]: length of the chain
#		chain[,3:(chain[,2]+2)]: values of the ridge
#	nbchain: number of chains.
#
#########################################################################
{
   d <- dim(ccridge)
   sigsize <- d[1]
   nscale <- d[2]
   threshold <- range(ccridge)[2] * ptile
   sz <- sigsize + 2   	
   chain <- matrix(-1,nbchain,sz)
   orderedmap <- matrix(0,sigsize,nscale)	

   dim(chain) <- c(nbchain * sz, 1)
   dim(ccridge) <- c(nscale * sigsize, 1)
   dim(orderedmap) <- c(nscale * sigsize, 1)
	

   
   z <- .C("Scrazy_family",
        as.double(ccridge),
        orderedmap = as.single(orderedmap),
	chain = as.integer(chain),
	chainnb = as.integer(nbchain),
	as.integer(sigsize),
        as.integer(nscale),
        as.integer(bstep),
	as.single(threshold))
	
   orderedmap <- z$orderedmap
   chain <- z$chain

   dim(orderedmap) <- c(sigsize,nscale)   
   dim(chain) <- c(nbchain,sz)   
   nbchain <- z$chainnb
   chain <- chain +1
   chain[,2] <- chain[,2]-1

   list(ordered = orderedmap, chain = chain, nbchain=nbchain)
}   


crfview <- function(beemap,twod=T)
#########################################################################
#     crfview:
#     --------
#     displays a family of chained ridges
#
#########################################################################
{
  if (twod)
    image(beemap$ordered > 0)
  else {
    Dim <- dim(beemap$ordered)
    matmp <- matrix(0,Dim[1],Dim[2])
    matmp[1,1] <- 1
    image(matmp)
    for (k in 1:beemap$nbchain){
      start <- beemap$chain[k,1]
      end <- start + beemap$chain[k,2]  
      lines(start:(end-1),beemap$chain[k,3:(end-start+2)],type="l")
    }
  }
  title(" Chained Ridges")
  cat("")
}

#########################################################################
#      $Log: Crc_Rec.S,v $
# Revision 1.2  1995/04/05  18:56:55  bruno
# *** empty log message ***
#
# Revision 1.1  1995/04/02  01:04:16  bruno
# Initial revision
#
#
#               (c) Copyright  1997                             
#                          by                                   
#      Author: Rene Carmona, Bruno Torresani, Wen-Liang Hwang   
#                  Princeton University
#                  All right reserved                           
#########################################################################


#########################################################################
#
#	Functions to reconstruct a signal from the output of
#	 the crazy climber algorithm.
#
#########################################################################




crcirrec <- function(siginput, inputwt, beemap, noct, nvoice, compr,
	minnbnodes = 2, w0 = 2*pi, bstep = 5,ptile = .01, prob = 0.8,
	epsilon = .5, fast = F, para = 0, real = F, plot=1)
#########################################################################
#     crcirrec:
#     -------
#      Reconstruction of a real valued signal from ridges found by 
#      crazy climbers on a wavelet transform (sampled irregularly).
#
#      input:
#      ------
#      siginput: input signal
#      inputwt: continuous wavelet transform (output of cwt)
#      beemap: output of crazy climber algorithm
#      noct: number of octaves (powers of 2)
#      nvoice: number of different scales per octave
#      compr: subsampling rate for the ridge
#      bstep: used for the chaining
#      ptile: 
#      epsilon: coeff of the Q2 part of reconstruction kernel
#      fast: if set to TRUE, computes the Q2 kernel by Riemann sums
#            if not, uses a Romberg adaptive step quadrature.
#      para:
#      plot: plot the original and reconstructed signal in display
#      prob: the percentile to keep 
#
#      output:
#      -------
#      rec: reconstructed signal
#      ordered: image of the ridges (with different colors)
#      comp: 2D array containing the signals reconstructed from ridges
#
#########################################################################
{

   tmp <- cfamily(beemap,bstep,ptile=ptile)
   image(tmp$ordered)
   chain <- tmp$chain
   nbchain <- tmp$nbchain
   ordered <- tmp$ordered
   sigsize <- length(siginput)
   rec <- numeric(sigsize)
   plnb <- 0

   par(mfrow=c(1,1))
   plot.ts(siginput)
   title("Original signal")

   tmp <- matrix(0,nbchain,length(siginput))

   totnbnodes <- 0
   idx <- numeric(nbchain)	
   p <- 0

   for (j in 1:nbchain){
      phi.x.min <- 2 * 2^(chain[j,3]/nvoice)
      if (chain[j,2] > (para*phi.x.min) ){ 

         cat("Chain number",j)

         phi.x.max <- 2 * 2^(chain[j,(2+chain[j,2])]/nvoice)
         x.min <- chain[j,1]
         x.max <- chain[j,1] + chain[j,2] - 1
         x.min <- x.min - round(para * phi.x.min)
         x.max <- x.max + round(para * phi.x.max)

         tmp2 <- irregrec(siginput[chain[j,1]:(chain[j,1]+chain[j,2]-1)],
            inputwt[chain[j,1]:(chain[j,1]+chain[j,2]-1),],
            chain[j,3:(chain[j,2]+2)], compr,noct,nvoice,
            epsilon, w0 = w0 , prob = prob, fast = fast, para = para,
            minnbnodes = minnbnodes, real = real);


         totnbnodes <- totnbnodes + tmp2$nbnodes

         np <- length(tmp2$sol)
         start <- max(1,x.min)
         end <- min(sigsize,x.min+np-1)
         start1 <- max(1,2-x.min)
         end1 <- min(np, sigsize +1 - x.min)
         end <- end1 - start1 + start
         rec[start:end] <- rec[start:end]+tmp2$sol[start1:end1]

         plnb <- plnb + 1
         p <- p + 1
         idx[p] <- j
  
       }
   }

   if(plot == 1){
      par(mfrow=c(2,1))
      par(cex=1.1)
      plot.ts(siginput)
      title("Original signal")
      title("Reconstructed signal")
   }
   else if (plot == 2){
      par(mfrow=c(plnb+2,1))
      par(mar=c(2,4,4,4))
      par(cex=1.1)
      par(err=-1)
      plot.ts(siginput)
      title("Original signal")

      for (j in 1:p)
         plot.ts(tmp[idx[j],]);

      plot.ts(Re(rec))	
      title("Reconstructed signal")


   }

   cat("Total number of ridge samples used: ",totnbnodes,"\n")

   par(mfrow=c(1,1))

   list(rec=rec, ordered=ordered,chain = chain, comp=tmp)
}



crcirgrec <- function(siginput, inputgt, beemap, nvoice, freqstep, scale,
	compr, prob = 0.8, bstep = 5, ptile = .01, epsilon = .5, fast = T, para = 0,
	minnbnodes=3, hflag = F, plot = 2, real = F)
#########################################################################
#     crcirgrec:
#     ----------
#      Reconstruction of a real valued signal from ridges found by 
#      crazy climbers on the gabor transform (irregularly sampled).
#
#      input:
#      ------
#      siginput: input signal
#      inputgt: continuous gabor transform (output of cgt)
#      beemap: output of crazy climber algorithm
#      nvoice: number of different frequencies
#      freqstep: difference between two consecutive frequencies
#      scale: scale of the window
#      compr: subsampling rate for the ridge
#      bstep: used for the chaining
#      ptile: 
#      epsilon: coeff of the Q2 part of reconstruction kernel
#      fast: if set to TRUE, computes the Q2 kernel by Riemann sums
#            if not, uses a Romberg adaptive step quadrature.
#      para:
#      prob: percentile to keep 
#      minnbnodes: minimal number of nodes for a sampled ridge.
#      hflag: if set to FALSE, uses the identity as first term
#             in the reconstruction kernel. If not, uses Q1 instead. 
#      plot: if set to 1, displays the signal, the components, and
#            signal and reconstruction one after another. If set to
#            2, displays the signal, the components, and the
#            reconstruction on the same page. Else, no plot.
#      real: if set to true, uses only real constraints.
#
#      output:
#      -------
#      rec: reconstructed signal
#      ordered: image of the ridges (with different colors)
#      comp: 2D array containing the signals reconstructed from ridges
#
#########################################################################
{
   tmp <- cfamily(beemap,bstep,ptile=ptile)
   image(tmp$ordered)
   chain <- tmp$chain
   nbchain <- tmp$nbchain
   ordered <- tmp$ordered
   sigsize <- length(siginput)
   rec <- numeric(sigsize)
   plnb <- 0

   par(mfrow=c(1,1))
   plot.ts(siginput)
   title("Original signal")

   totnbnodes <- 0
   tmp <- matrix(0,nbchain,length(siginput))
   for (j in 1:nbchain){
      if (chain[j,2] > scale){
         nbnodes <- round(chain[j,2]/compr);
         if(nbnodes < minnbnodes)
            nbnodes <- minnbnodes

         totnbnodes <- totnbnodes + nbnodes

         cat("Chain number",j)

         phi.x.min <- scale
         phi.x.max <- scale
         x.min <- chain[j,1]
         x.max <- chain[j,1] + chain[j,2] - 1
         x.min <- x.min - round(para * phi.x.min)
         x.max <- x.max + round(para * phi.x.max)
         x.inc <- 1
         np <- as.integer((x.max-x.min)/x.inc) + 1

         tmp2 <- girregrec(siginput[chain[j,1]:(chain[j,1]+chain[j,2]-1)],
            inputgt[chain[j,1]:(chain[j,1]+chain[j,2]-1),],
            chain[j,3:(chain[j,2]+2)], nbnodes,nvoice,freqstep,scale,
            epsilon,fast,prob=prob,para=para, hflag = hflag, real = real);


         start <- max(1,x.min)
         end <- min(sigsize,x.min+np-1)
         start1 <- max(1,2-x.min)
         end1 <- min(np, sigsize +1 - x.min)
         end <- end1 - start1 + start
         rec[start:end] <- rec[start:end]+tmp2$sol[start1:end1]
         tmp[j,start:end] <- tmp2$sol[start1:end1]

         plnb <- plnb + 1
         plot.ts(tmp[j,]);
      }

   }

   if(plot == 1){
      par(mfrow=c(2,1))
      par(cex=1.1)
      plot.ts(siginput)
      title("Original signal")
      plot.ts(rec)
      title("Reconstructed signal")
   }
   else if (plot == 2){
      par(mfrow=c(plnb+2,1))
      par(mar=c(2,4,4,4))
      par(cex=1.1)
      par(err=-1)
      plot.ts(siginput)
      title("Original signal")

      for (j in 1:nbchain)
         plot.ts(tmp[j,]);
        
      plot.ts(rec)
      title("Reconstructed signal")
   }

   cat("Total number of ridge samples used: ",totnbnodes,"\n")

   par(mfrow=c(1,1))
   list(rec=rec, ordered=ordered,chain = chain, comp=tmp)
}















#########################################################################
#      $Log: Crc_Rec.S,v $
#########################################################################
#
#               (c) Copyright  1997                             
#                          by                                   
#      Author: Rene Carmona, Bruno Torresani, Wen-Liang Hwang   
#                  Princeton University
#                  All right reserved                           
#########################################################################


#########################################################################
#
#	Functions to reconstruct a signal from the output of
#	 the crazy climber algorithm.
#
#########################################################################




crcrec <- function(siginput, inputwt, beemap, noct, nvoice, compr,
	minnbnodes = 2, w0 = 2*pi, bstep = 5,ptile = .01,
	epsilon = 0, fast = F, para = 5, real = F, plot=2)
#########################################################################
#     crcrec:
#     -------
#      Reconstruction of a real valued signal from ridges found by 
#      crazy climbers on a wavelet transform.
#
#      input:
#      ------
#      siginput: input signal
#      inputwt: continuous wavelet transform (output of cwt)
#      beemap: output of crazy climber algorithm
#      noct: number of octaves (powers of 2)
#      nvoice: number of different scales per octave
#      compr: subsampling rate for the ridge
#      bstep: used for the chaining
#      ptile: 
#      epsilon: coeff of the Q2 part of reconstruction kernel
#      fast: if set to TRUE, computes the Q2 kernel by Riemann sums
#            if not, uses a Romberg adaptive step quadrature.
#      para:
#      plot: plot the original and reconstructed signal in display
#
#      output:
#      -------
#      rec: reconstructed signal
#      ordered: image of the ridges (with different colors)
#      comp: 2D array containing the signals reconstructed from ridges
#
#########################################################################
{

   tmp <- cfamily(beemap,bstep,ptile=ptile)
   chain <- tmp$chain
   nbchain <- tmp$nbchain
   ordered <- tmp$ordered
   sigsize <- length(siginput)
   rec <- numeric(sigsize)
   plnb <- 0

   if(plot != F){
     par(mfrow=c(2,1))
     plot.ts(siginput)
     title("Original signal")
     image(tmp$ordered)
     title("Chained Ridges")
   }

   tmp <- matrix(0,nbchain,length(siginput))

   totnbnodes <- 0
   idx <- numeric(nbchain)	
   p <- 0

   for (j in 1:nbchain){
      phi.x.min <- 2 * 2^(chain[j,3]/nvoice)
      if (chain[j,2] > (para*phi.x.min) ){ 

         cat("Chain number",j)

         phi.x.max <- 2 * 2^(chain[j,(2+chain[j,2])]/nvoice)
         x.min <- chain[j,1]
         x.max <- chain[j,1] + chain[j,2] - 1
         x.min <- x.min - round(para * phi.x.min)
         x.max <- x.max + round(para * phi.x.max)

         tmp2 <- regrec(siginput[chain[j,1]:(chain[j,1]+chain[j,2]-1)],
            inputwt[chain[j,1]:(chain[j,1]+chain[j,2]-1),],
            chain[j,3:(chain[j,2]+2)], compr,noct,nvoice,
            epsilon, w0 = w0 ,fast = fast, para = para,
            minnbnodes = minnbnodes, real = real);


         if(is.list(tmp2)==T) {
         totnbnodes <- totnbnodes + tmp2$nbnodes

         np <- length(tmp2$sol)
         start <- max(1,x.min)
         end <- min(sigsize,x.min+np-1)
         start1 <- max(1,2-x.min)
         end1 <- min(np, sigsize +1 - x.min)
         end <- end1 - start1 + start
         rec[start:end] <- rec[start:end]+tmp2$sol[start1:end1]

         tmp[j,start:end] <- Re(tmp2$sol[start1:end1])
         }
         plnb <- plnb + 1
         p <- p + 1
         idx[p] <- j
  
       }
   }

   if(plot == 1){
      par(mfrow=c(2,1))
      par(cex=1.1)
      plot.ts(siginput)
      title("Original signal")
      plot.ts(Re(rec))
	
      title("Reconstructed signal")
   }
   else if (plot == 2){
      par(mfrow=c(plnb+2,1))
      par(mar=c(2,4,4,4))
      par(cex=1.1)
      par(err=-1)
      plot.ts(siginput)
      title("Original signal")

      for (j in 1:p)
         plot.ts(tmp[idx[j],]);
      plot.ts(Re(rec))	
      title("Reconstructed signal")


   }

   cat("Total number of ridge samples used: ",totnbnodes,"\n")

   par(mfrow=c(1,1))

   list(rec=rec, ordered=ordered,chain = chain, comp=tmp)
}



gcrcrec <- function(siginput, inputgt, beemap, nvoice, freqstep, scale,
	compr, bstep = 5, ptile = .01, epsilon = 0, fast = T, para = 5,
	minnbnodes=3, hflag = F, real= F, plot = 2)
#########################################################################
#     gcrcrec:
#     -------
#      Reconstruction of a real valued signal from ridges found by 
#      crazy climbers on the gabor transform.
#
#      input:
#      ------
#      siginput: input signal
#      inputgt: continuous gabor transform (output of cgt)
#      beemap: output of crazy climber algorithm
#      nvoice: number of different frequencies
#      freqstep: difference between two consecutive frequencies
#      scale: scale of the window
#      compr: subsampling rate for the ridge
#      bstep: used for the chaining
#      ptile: 
#      epsilon: coeff of the Q2 part of reconstruction kernel
#      fast: if set to TRUE, computes the Q2 kernel by Riemann sums
#            if not, uses a Romberg adaptive step quadrature.
#      para:
#      minnbnodes: minimal number of nodes for a sampled ridge.
#      hflag: if set to FALSE, uses the identity as first term
#             in the reconstruction kernel. If not, uses Q1 instead. 
#      plot: if set to 1, displays the signal, the components, and
#            signal and reconstruction one after another. If set to
#            2, displays the signal, the components, and the
#            reconstruction on the same page. Else, no plot.
#      real: if set to true, uses only real constraints.
#
#      output:
#      -------
#      rec: reconstructed signal
#      ordered: image of the ridges (with different colors)
#      comp: 2D array containing the signals reconstructed from ridges
#
#########################################################################
{
   tmp <- cfamily(beemap,bstep,ptile=ptile)
   chain <- tmp$chain
   nbchain <- tmp$nbchain
   ordered <- tmp$ordered
   sigsize <- length(siginput)
   rec <- numeric(sigsize)
   plnb <- 0

   if(plot != F){
     par(mfrow=c(2,1))
     plot.ts(siginput)
     title("Original signal")
     image(tmp$ordered)
     title("Chained Ridges")
   }

   totnbnodes <- 0
   tmp <- matrix(0,nbchain,length(siginput))
   for (j in 1:nbchain){
      if (chain[j,2] > scale){
         nbnodes <- round(chain[j,2]/compr);
         if(nbnodes < minnbnodes)
            nbnodes <- minnbnodes

         totnbnodes <- totnbnodes + nbnodes

         cat("Chain number",j)

         phi.x.min <- scale
         phi.x.max <- scale
         x.min <- chain[j,1]
         x.max <- chain[j,1] + chain[j,2] - 1
         x.min <- x.min - round(para * phi.x.min)
         x.max <- x.max + round(para * phi.x.max)
         x.inc <- 1
         np <- as.integer((x.max-x.min)/x.inc) + 1

         tmp2 <- gregrec(siginput[chain[j,1]:(chain[j,1]+chain[j,2]-1)],
            inputgt[chain[j,1]:(chain[j,1]+chain[j,2]-1),],
            chain[j,3:(chain[j,2]+2)], nbnodes,nvoice,freqstep,scale,
            epsilon,fast,para=para, hflag = hflag, real = real);


         start <- max(1,x.min)
         end <- min(sigsize,x.min+np-1)
         start1 <- max(1,2-x.min)
         end1 <- min(np, sigsize +1 - x.min)
         end <- end1 - start1 + start
         rec[start:end] <- rec[start:end]+tmp2$sol[start1:end1]
         tmp[j,start:end] <- tmp2$sol[start1:end1]

         plnb <- plnb + 1
         plot.ts(tmp[j,]);
      }

   }

   if(plot == 1){
      par(mfrow=c(2,1))
      par(cex=1.1)
      plot.ts(siginput)
      title("Original signal")
      plot.ts(rec)
      title("Reconstructed signal")
   }
   else if (plot == 2){
      par(mfrow=c(plnb+2,1))
      par(mar=c(2,4,4,4))
      par(cex=1.1)
      par(err=-1)
      plot.ts(siginput)
      title("Original signal")

      for (j in 1:nbchain)
         plot.ts(tmp[j,]);
        
      plot.ts(rec)
      title("Reconstructed signal")
   }

   cat("Total number of ridge samples used: ",totnbnodes,"\n")

   par(mfrow=c(1,1))
   list(rec=rec, ordered=ordered,chain = chain, comp=tmp)
}




scrcrec <- function(siginput, tfinput, beemap, bstep = 5,
	ptile = .01, plot = 2)
#########################################################################
#     scrcrec:
#     --------
#      Simple reconstruction of a real valued signal from ridges found by 
#      crazy climbers.
#
#      input:
#      ------
#      siginput: input signal
#      tfinput: continuous time-frequency transform (output of cwt or cgt)
#      beemap: output of crazy climber algorithm
#      bstep: used for the chaining
#      ptile: 
#      plot: if set to 1, displays the signal, the components, and
#            signal and reconstruction one after another. If set to
#            2, displays the signal, the components, and the
#            reconstruction on the same page. Else, no plot.
#
#      output:
#      -------
#      rec: reconstructed signal
#      ordered: image of the ridges (with different colors)
#      comp: 2D array containing the signals reconstructed from ridges
#
#########################################################################
{
   tmp <- cfamily(beemap,bstep,ptile=ptile)
   chain <- tmp$chain
   nbchain <- tmp$nbchain
   rec <- numeric(length(siginput))

   if(plot != F){
     par(mfrow=c(2,1))
     plot.ts(siginput)
     title("Original signal")
     image(tmp$ordered)
     title("Chained Ridges")
   }

   npl(1)

   tmp3 <- matrix(0,nbchain,length(siginput))

   rdg <- numeric(dim(tfinput)[1])
   
   for (j in 1:nbchain){
       bnext <- chain[j,1]
     for(k in 1:chain[j,2]) {
       rdg[bnext] <- chain[j,2+k]
       bnext <- bnext + 1
     }
     tmp3[j,] <- sridrec(tfinput,rdg)
     rec <- rec + tmp3[j,]
     rdg[] <- 0
     plot.ts(tmp3[j,]);
   }

   if(plot == 1){
      par(mfrow=c(2,1))
      par(cex=1.1)
      plot.ts(siginput)
      title("Original signal")
      plot.ts(rec)
      title("Reconstructed signal")
   }
   else if (plot == 2){
      par(mfrow=c(nbchain+2,1))
      par(mar=c(2,4,4,4))
      par(cex=1.1)
      par(err=-1)
      plot.ts(siginput)
      title("Original signal")

      for (j in 1:nbchain)
         plot.ts(tmp3[j,]);
        
      plot.ts(rec)
      title("Reconstructed signal")
   }
   list(rec=rec, ordered=tmp$ordered,chain = tmp$chain, comp=tmp3)
}











#########################################################################
#       $Log: Cwt_DOG.S,v $
#########################################################################
#
#               (c) Copyright  1997
#                          by                                   
#      Author: Rene Carmona, Bruno Torresani, Wen-Liang Hwang   
#                  Princeton University
#                  All right reserved                           
#########################################################################





DOG <- function(input, noctave, nvoice = 1, moments,
	twoD = TRUE, plot = TRUE)
#########################################################################
#       DOG: Derivative of Gaussian  
#       ----
# 	 continuous wavelet transform function:
#         compute the continuous wavelet transform with (complex-valued)
#	  derivative of gaussians wavelet
#
#       input:
#       ------
# 	 input: input signal (possibly complex-valued)
#	 noctave: number of powers of 2 for the scale variable
#	 nvoice: number of scales between 2 consecutive powers of 2
#        moments: number of vanishing moments for the wavelet
#	 twoD: if set to TRUE, organizes the output as a 2D array 
#			(signal_size X nb_scales)
#		      if not: 3D array (signal_size X noctave X nvoice)
#	 plot: if set to TRUE, displays the modulus of cwt on the graphic
#		device.
#
#       output:
#       -------
#        tmp: continuous (complex) wavelet transform
#
#########################################################################
{
   oldinput <- input
   isize <- length(oldinput)

   tmp <- adjust.length(oldinput)
   input <- tmp$signal
   newsize <- length(input)

   pp <- noctave * nvoice
   Routput <- matrix(0,newsize,pp)
   Ioutput <- matrix(0,newsize,pp)
   output <- matrix(0,newsize,pp)
   dim(Routput) <- c(pp * newsize,1)
   dim(Ioutput) <- c(pp * newsize,1)
   dim(input) <- c(newsize,1)

   z <- .C("Scwt_DOG",
            as.single(Re(input)),
            as.single(Im(input)),
            Rtmp = as.double(Routput),
            Itmp = as.double(Ioutput),
            as.integer(noctave),
            as.integer(nvoice),
            as.integer(newsize),
            as.integer(moments))

   Routput <- z$Rtmp
   Ioutput <- z$Itmp
   dim(Routput) <- c(newsize,pp)
   dim(Ioutput) <- c(newsize,pp)
   i <- sqrt(as.complex(-1))
   if(twoD) {
     output <- Routput[1:isize,] + Ioutput[1:isize,] * i
     if(plot) image(Mod(output))
     output
   } 
   else {
     Rtmp <- array(0,c(isize,noctave,nvoice))
     Itmp <- array(0,c(isize,noctave,nvoice))
     for(i in 1:noctave)
       for(j in 1:nvoice) {
         Rtmp[,i,j] <- Routput[1:isize,(i-1)*nvoice+j]
         Itmp[,i,j] <- Ioutput[1:isize,(i-1)*nvoice+j]
      }
     Rtmp + Itmp * i
   }
}



vDOG <- function(input, scale, moments)
#########################################################################
#       vDOG:   
#       -----
#        continuous wavelet transform on one scale:
# 	 compute the continuous wavelet transform with (complex-valued)
#	  derivative of gaussians wavelet
#
#       input:
#       ------
# 	 input: input signal (possibly complex-valued)
#        scale: value of the scale at which the transform is computed
#        moments: number of vanishing moments
#
#       output:
#       -------
#        Routput + i Ioutput: voice wavelet transform (complex 1D array)
#
#########################################################################
{
   oldinput <- input
   isize <- length(oldinput)

   tmp <- adjust.length(oldinput)
   input <- tmp$signal
   newsize <- length(input)

   Routput <- numeric(newsize)
   Ioutput <- numeric(newsize)
   dim(input) <- c(newsize,1)

   z <- .C("Svwt_DOG",
            as.single(Re(input)),
            as.single(Im(input)),
            Rtmp = as.double(Routput),
            Itmp = as.double(Ioutput),
            as.single(scale),
            as.integer(newsize),
            as.integer(moments))

   Routput <- z$Rtmp
   Ioutput <- z$Itmp
   i <- sqrt(as.complex(-1))

   Routput[1:isize] + Ioutput[1:isize] * i
}






#########################################################################
#      $Log: Cwt_Morlet.S,v $
#########################################################################
#
#               (c) Copyright  1997                             
#                          by                                   
#      Author: Rene Carmona, Bruno Torresani, Wen-Liang Hwang   
#                  Princeton University
#                  All right reserved                           
#########################################################################

cwt <- function(input, noctave, nvoice = 1, w0 = 2*pi,
	twoD = TRUE, plot = TRUE)
#########################################################################
#      cwt:
#      ---
#       continuous wavelet transform function
# 	compute the continuous wavelet transform with (complex-valued)
#			Morlet wavelet
#
#       Input:
#       ------
# 	 input: input signal (possibly complex-valued)
#	 noctave: number of powers of 2 for the scale variable
#	 nvoice: number of scales between 2 consecutive powers of 2
#        w0: central frequency of Morlet wavelet
#	 twoD: if set to TRUE, organizes the output as a 2D array 
#			(signal_size X nb_scales)
#		      if not: 3D array (signal_size X noctave X nvoice)
#	 plot: if set to TRUE, displays the modulus of cwt on the graphic
#		device.
#
#       Output:
#       -------
#        tmp: continuous (complex) wavelet transform
#
#########################################################################
{
   oldinput <- input
   isize <- length(oldinput)

   tmp <- adjust.length(oldinput)
   input <- tmp$signal
   newsize <- length(input)

   pp <- noctave * nvoice
   Routput <- matrix(0,newsize,pp)
   Ioutput <- matrix(0,newsize,pp)
   output <- matrix(0,newsize,pp)
   dim(Routput) <- c(pp * newsize,1)
   dim(Ioutput) <- c(pp * newsize,1)
   dim(input) <- c(newsize,1)

   z <- .C("Scwt_morlet",
            as.single(Re(input)),
            as.single(Im(input)),
            Rtmp = as.double(Routput),
            Itmp = as.double(Ioutput),
            as.integer(noctave),
            as.integer(nvoice),
            as.integer(newsize),
            as.single(w0))

   Routput <- z$Rtmp
   Ioutput <- z$Itmp
   dim(Routput) <- c(newsize,pp)
   dim(Ioutput) <- c(newsize,pp)
   i <- sqrt(as.complex(-1))
   if(twoD) {
     output <- Routput[1:isize,] + Ioutput[1:isize,] * i
     if(plot) {
	image(Mod(output), xlab="Time", ylab="log(scale)")
	title("Wavelet Transform Modulus")
     }
     output
   } 
   else {
     Rtmp <- array(0,c(isize,noctave,nvoice))
     Itmp <- array(0,c(isize,noctave,nvoice))
     for(i in 1:noctave)
       for(j in 1:nvoice) {
         Rtmp[,i,j] <- Routput[1:isize,(i-1)*nvoice+j]
         Itmp[,i,j] <- Ioutput[1:isize,(i-1)*nvoice+j]
      }
     Rtmp + Itmp * i
   }
}




cwtpolar <- function(cwt,threshold=0.0)
#########################################################################
#       cwtpolar:   
#       --------
# 	continuous wavelet transform conversion:
#        converts one of the possible outputs of cwt to modulus and phase
#
#       input:
#       ------
# 	 cwt: 3D array containing a continuous wavelet transform (output
#		of cwt, with twoD=FALSE)
#	 threshold: the phase is forced to -pi if the modulus is 
#		less than threshold.
#
#       output:
#       -------
#        output1: modulus
#        output2: phase      
#
#########################################################################
{
   tmp1 <- cwt
   sigsize <- dim(tmp1)[1] # sig size
   noctave <- dim(tmp1)[2]
   nvoice <- dim(tmp1)[3]

   output1 <- array(0,c(sigsize,noctave,nvoice))
   output2 <- array(0,c(sigsize,noctave,nvoice))
   for(i in 1:noctave)
     for(j in 1:nvoice) {
       output1[,i,j] <- sqrt(Re(tmp1[,i,j])*Re(tmp1[,i,j])+
          Im(tmp1[,i,j])*Im(tmp1[,i,j]))
       output2[,i,j] <- atan(Im(tmp1[,i,j]),Re(tmp1[,i,j]))
    }

   ma <- max(output1)
   rel <- threshold * ma
   dim(output1) <- c(sigsize * noctave * nvoice,1)
   dim(output2) <- c(sigsize * noctave * nvoice,1)
   output2[abs(output1) < rel] <- -3.14159
   dim(output1) <- c(sigsize,noctave,nvoice)
   dim(output2) <- c(sigsize,noctave,nvoice)

   list(modulus=output1,argument=output2)
}
   


cwtimage <- function(input)
#########################################################################
#       cwtimage:   
#       ---------
#        continuous wavelet transform display
# 	 converts the output (modulus or argument) of cwtpolar to a 
#		2D array and displays on the graphic device
#
#       input:
#       ------
# 	 input: 3D array containing a continuous wavelet transform (output
#		of cwtpolar
#       output:
#       -------
#        output: 2D array 
#
#########################################################################
{
  sigsize <- dim(input)[1]
  noctave <- dim(input)[2]
  nvoice <- dim(input)[3]
  
  output <- matrix(0,sigsize, noctave * nvoice)
  for(i in 1:noctave) {
     k <- (i-1) * nvoice
   for(j in 1:nvoice) {
     output[,k+j] <- input[,i,j]
   }
  }
  image(output)
  output
}
   


vwt <- function(input, scale, w0 = 2*pi)
#########################################################################
#      vwt:   voice Morlet wavelet transform   
#      ----
#       continuous wavelet transform on one scale
# 	compute the continuous wavelet transform with (complex-valued)
#			Morlet wavelet
#
#       input:
#       ------
# 	 input: input signal (possibly complex-valued)
#        scale: value of the scale at which the transform is computed
#        w0: central frequency of Morlet wavelet
#
#       output:
#       -------
#        Routput + i Ioutput: voice wavelet transform (complex 1D array)
#
#########################################################################
{
   oldinput <- input
   isize <- length(oldinput)

   tmp <- adjust.length(oldinput)
   input <- tmp$signal
   newsize <- length(input)

   Routput <- numeric(newsize)
   Ioutput <- numeric(newsize)
   dim(input) <- c(newsize,1)

   z <- .C("Svwt_morlet",
            as.single(Re(input)),
            as.single(Im(input)),
            Rtmp = as.double(Routput),
            Itmp = as.double(Ioutput),
            as.single(scale),
            as.integer(newsize),
            as.single(w0))

   Routput <- z$Rtmp
   Ioutput <- z$Itmp
   i <- sqrt(as.complex(-1))

   Routput[1:isize] + Ioutput[1:isize] * i
}




morlet <- function(sigsize, location, scale, w0 = 2*pi)
#########################################################################
#       morlet: Morlet's wavelet  
#       -------
#        Generates a Morlet wavelet for given location and scale
#
#       input:
#       ------
#        sigsize: signal size (dimension of the array)
#        location: location of the wavelet
#        scale: value of the scale at which the transform is computed
#        w0: central frequency of Morlet wavelet
#
#       output:
#       -------
#        z$wavelet.r + z$wavelet.i * i: wavelet (complex 1D array
#          of size sigsize)
#
#########################################################################
{
   wavelet.r <- numeric(sigsize)
   wavelet.i <- numeric(sigsize)


   z <- .C("morlet_time",
            as.single(w0),
            as.single(scale),
            as.integer(location),
            wavelet.r = as.double(wavelet.r),
            wavelet.i = as.double(wavelet.i),
            as.integer(sigsize))

   i <- sqrt(as.complex(-1))

   z$wavelet.r + z$wavelet.i * i
}


vecmorlet <- function(sigsize, nbnodes, bridge, aridge, w0 = 2*pi)
#########################################################################
#       vecmorlet:   vector of Morlet's wavelets
#       ---------
#        Generates morlet wavelets located on the ridge
#
#       input:
#       ------
#        sigsize: signal size (dimension of the array)
#        nbnodes: number of ridge samples
#        bridge: b coordinates of the ridge samples
#        aridge: acoordinates of the ridge samples
#        w0: central frequency of Morlet wavelet
#
#       output:
#       -------
#        z$morlet.r + z$morlet.i * i: 2D array containing the morlet
#           wavelets located on the ridge
#
#########################################################################
{

   morlet.r <- numeric(nbnodes * sigsize)
   morlet.i <- numeric(nbnodes * sigsize)


   z <- .C("vmorlet_time",
            as.single(w0),
            as.single(aridge),
            as.integer(bridge),
            morlet.r = as.double(morlet.r),
            morlet.i = as.double(morlet.i),
            as.integer(sigsize),
            as.integer(nbnodes))

   i <- sqrt(as.complex(-1))

   z$morlet.r + z$morlet.i * i
}












#########################################################################
#      $Log: Cwt_Squeezing.S,v $
#########################################################################
#
#               (c) Copyright  1997
#                          by                                   
#      Author: Rene Carmona, Bruno Torresani, Wen-Liang Hwang   
#                  Princeton University
#                  All right reserved                           
#########################################################################






cwtsquiz <- function(input, noctave, nvoice = 1, w0 = 2*pi,
	twoD = TRUE, plot = TRUE)
#########################################################################
#      cwtsquiz:
#      --------
#       squeezed wavelet transform function
#
#       Input:
#       ------
# 	 input: input signal (possibly complex-valued)
#	 noctave: number of powers of 2 for the scale variable
#	 nvoice: number of scales between 2 consecutive powers of 2
#        w0: central frequency of Morlet wavelet
#	 twoD: if set to TRUE, organizes the output as a 2D array 
#			(signal_size X nb_scales)
#		      if not: 3D array (signal_size X noctave X nvoice)
#	 plot: if set to TRUE, displays the modulus of cwt on the graphic
#		device.
#
#       output:
#       -------
#        tmp: continuous (complex) squeezed wavelet transform
#
#########################################################################
{
   oldinput <- input
   isize <- length(oldinput)

   tmp <- adjust.length(oldinput)
   input <- tmp$signal
   newsize <- length(input)

   pp <- noctave * nvoice
   Routput <- matrix(0,newsize,pp)
   Ioutput <- matrix(0,newsize,pp)
   output <- matrix(0,newsize,pp)
   dim(Routput) <- c(pp * newsize,1)
   dim(Ioutput) <- c(pp * newsize,1)
   dim(input) <- c(newsize,1)

   z <- .C("Scwt_squeezed",
            as.single(input),
            Rtmp = as.double(Routput),
            Itmp = as.double(Ioutput),
            as.integer(noctave),
            as.integer(nvoice),
            as.integer(newsize),
            as.single(w0))

   Routput <- z$Rtmp
   Ioutput <- z$Itmp
   dim(Routput) <- c(newsize,pp)
   dim(Ioutput) <- c(newsize,pp)
   i <- sqrt(as.complex(-1))
   if(twoD) {
     output <- Routput[1:isize,] + Ioutput[1:isize,] * i
     if(plot){
        image(Mod(output),xlab="Time", ylab="log(scale)")
	title("Squeezed Wavelet Transform Modulus")
     }
     output
   } 
   else {
     Rtmp <- array(0,c(isize,noctave,nvoice))
     Itmp <- array(0,c(isize,noctave,nvoice))
     for(i in 1:noctave)
       for(j in 1:nvoice) {
         Rtmp[,i,j] <- Routput[1:isize,(i-1)*nvoice+j]
         Itmp[,i,j] <- Ioutput[1:isize,(i-1)*nvoice+j]
      }
     Rtmp + Itmp * i
   }
}






#########################################################################
#       $Log: Cwt_Thierry.S,v $
#########################################################################
#
#               (c) Copyright  1997                             
#                          by                                   
#      Author: Rene Carmona, Bruno Torresani, Wen-Liang Hwang   
#                  Princeton University
#                  All right reserved                           
#########################################################################





cwtTh <- function(input, noctave, nvoice = 1, moments,
	twoD = TRUE, plot = TRUE)
#########################################################################
#       cwtTh:   Cauchy's wavelet transform
#       -----
# 	 continuous wavelet transform function:
#         compute the continuous wavelet transform with (complex-valued)
#	  Cauchy's wavelet
#
#       input:
#       ------
# 	 input: input signal (possibly complex-valued)
#	 noctave: number of powers of 2 for the scale variable
#	 nvoice: number of scales between 2 consecutive powers of 2
#        moments: number of vanishing moments for the wavelet
#	 twoD: if set to TRUE, organizes the output as a 2D array 
#			(signal_size X nb_scales)
#		      if not: 3D array (signal_size X noctave X nvoice)
#	 plot: if set to TRUE, displays the modulus of cwt on the graphic
#		device.
#
#       output:
#       -------
#        tmp: continuous (complex) wavelet transform
#
#########################################################################
{
   oldinput <- input
   isize <- length(oldinput)

   tmp <- adjust.length(oldinput)
   input <- tmp$signal
   newsize <- length(input)

   pp <- noctave * nvoice
   Routput <- matrix(0,newsize,pp)
   Ioutput <- matrix(0,newsize,pp)
   output <- matrix(0,newsize,pp)
   dim(Routput) <- c(pp * newsize,1)
   dim(Ioutput) <- c(pp * newsize,1)
   dim(input) <- c(newsize,1)

   z <- .C("Scwt_thierry",
            as.single(Re(input)),
            as.single(Im(input)),
            Rtmp = as.double(Routput),
            Itmp = as.double(Ioutput),
            as.integer(noctave),
            as.integer(nvoice),
            as.integer(newsize),
            as.integer(moments))

   Routput <- z$Rtmp
   Ioutput <- z$Itmp
   dim(Routput) <- c(newsize,pp)
   dim(Ioutput) <- c(newsize,pp)
   i <- sqrt(as.complex(-1))
   if(twoD) {
     output <- Routput[1:isize,] + Ioutput[1:isize,] * i
     if(plot) image(Mod(output),xlab="Time",ylab="log(scale)")
     title("Wavelet Transform Modulus")
     output
   } 
   else {
     Rtmp <- array(0,c(isize,noctave,nvoice))
     Itmp <- array(0,c(isize,noctave,nvoice))
     for(i in 1:noctave)
       for(j in 1:nvoice) {
         Rtmp[,i,j] <- Routput[1:isize,(i-1)*nvoice+j]
         Itmp[,i,j] <- Ioutput[1:isize,(i-1)*nvoice+j]
      }
     Rtmp + Itmp * i
   }
}



vwtTh <- function(input, scale, moments)
#########################################################################
#       vwtTh:   
#       -----
#        continuous wavelet transform on one scale:
# 	 compute the continuous wavelet transform with (complex-valued)
#	  Cauchy's wavelet
#
#       input:
#       ------
# 	 input: input signal (possibly complex-valued)
#        scale: value of the scale at which the transform is computed
#        moments: number of vanishing moments
#
#       output:
#       -------
#        Routput + i Ioutput: voice wavelet transform (complex 1D array)
#
#########################################################################
{
   oldinput <- input
   isize <- length(oldinput)

   tmp <- adjust.length(oldinput)
   input <- tmp$signal
   newsize <- length(input)

   Routput <- numeric(newsize)
   Ioutput <- numeric(newsize)
   dim(input) <- c(newsize,1)

   z <- .C("Svwt_Thierry",
            as.single(Re(input)),
            as.single(Im(input)),
            Rtmp = as.double(Routput),
            Itmp = as.double(Ioutput),
            as.single(scale),
            as.integer(newsize),
            as.integer(moments))

   Routput <- z$Rtmp
   Ioutput <- z$Itmp
   i <- sqrt(as.complex(-1))

   Routput[1:isize] + Ioutput[1:isize] * i
}



#########################################################################
#      $Log: Cwt_phase.S,v $
#########################################################################
#
#               (c) Copyright  1997
#                          by                                   
#      Author: Rene Carmona, Bruno Torresani, Wen-Liang Hwang   
#                  Princeton University
#                  All right reserved                           
#########################################################################





cwtp <- function(input, noctave, nvoice = 1, w0 = 2*pi,
	twoD = TRUE, plot = TRUE)
#########################################################################
#       cwtp:   
#       ----
#        continuous wavelet transform function:
# 	  compute the continuous wavelet transform with Morlet wavelet,
#    	  and its phase derivative (minus the wavelet frequency)
#
#       input:
#       ------
# 	 input: input signal (possibly complex-valued)
#	 noctave: number of powers of 2 for the scale variable
#	 nvoice: number of scales between 2 consecutive powers of 2
#        w0: central frequency of Morlet wavelet
#	 twoD: if set to TRUE, organizes the output as a 2D array 
#			(signal_size X nb_scales)
#		      if not: 3D array (signal_size X noctave X nvoice)
#	 plot: if set to TRUE, displays the modulus of cwt on the graphic
#		device.
#
#       output:
#       -------
#        wt: continuous (complex) wavelet transform
#        f: derivative of wavelet transform phase
#
#########################################################################
{
   oldinput <- input
   isize <- length(oldinput)

   tmp <- adjust.length(oldinput)
   input <- tmp$signal
   newsize <- length(input)

   pp <- noctave * nvoice
   Routput <- matrix(0,newsize,pp)
   Ioutput <- matrix(0,newsize,pp)
   Foutput <- matrix(0,newsize,pp)
   output <- matrix(0,newsize,pp)
   dim(Routput) <- c(pp * newsize,1)
   dim(Ioutput) <- c(pp * newsize,1)
   dim(Foutput) <- c(pp * newsize,1)
   dim(input) <- c(newsize,1)

   dyn.load("/usr/people/bruno/Swave/wave.a");
   dyn.load("/usr/people/whwang/Swave/wave.a");

   z <- .C("Scwt_phase",
            as.single(input),
            Rtmp = as.double(Routput),
            Itmp = as.double(Ioutput),
            Ftmp = as.double(Foutput),
            as.integer(noctave),
            as.integer(nvoice),
            as.integer(newsize),
            as.single(w0))

   Routput <- z$Rtmp
   Ioutput <- z$Itmp
   Foutput <- z$Ftmp
   dim(Routput) <- c(newsize,pp)
   dim(Foutput) <- c(newsize,pp)
   dim(Ioutput) <- c(newsize,pp)

   if(twoD) {
     Ftmp <- Foutput[1:isize,]
     output <- Routput[1:isize,] + Ioutput[1:isize,] * i
     if(plot) image(Mod(output))
     list(wt=output,f=Ftmp)
   } 
   else {
     Rtmp <- array(0,c(isize,noctave,nvoice))
     Itmp <- array(0,c(isize,noctave,nvoice))
     Ftmp <- array(0,c(isize,noctave,nvoice))
     for(i in 1:noctave)
       for(j in 1:nvoice) {
         Rtmp[,i,j] <- Routput[1:isize,(i-1)*nvoice+j]
         Itmp[,i,j] <- Ioutput[1:isize,(i-1)*nvoice+j]
         Ftmp[,i,j] <- Foutput[1:isize,(i-1)*nvoice+j]
         }       
     i <- sqrt(as.complex(-1))
     list(wt=Rtmp + Itmp * i,f=Ftmp)
   }
}





#########################################################################
#      $Log: Gabor.S,v $
#########################################################################
#
#               (c) Copyright  1997
#                          by                                   
#      Author: Rene Carmona, Bruno Torresani, Wen-Liang Hwang   
#                  Princeton University 
#                  All right reserved                           
#########################################################################





cgt <- function(input, nvoice, freqstep = (1/nvoice),
	scale = 1, plot = TRUE)
#########################################################################
#       cgt:   
#       ---
#        continuous Gabor transform function:
# 	  compute the continuous Gabor transform with gaussian window.
#
#       input:
#       ------
# 	 input: input signal (possibly complex-valued)
#	 nvoice: number of frequency bands
#        freqstep: sampling rate for the frequency axis
#        scale: size of the window
#	 plot: if set to TRUE, displays the modulus of cwt on the graphic
#		device.
#
#       output:
#       -------
#        output: continuous (complex) gabor transform
#
#########################################################################
{
   oldinput <- input
   isize <- length(oldinput)

   tmp <- adjust.length(oldinput)
   input <- tmp$signal
   newsize <- length(input)

   pp <- nvoice
   Routput <- matrix(0,newsize,pp)
   Ioutput <- matrix(0,newsize,pp)
   output <- matrix(0,newsize,pp)
   dim(Routput) <- c(pp * newsize,1)
   dim(Ioutput) <- c(pp * newsize,1)
   dim(input) <- c(newsize,1)


   z <- .C("Sgabor",
            as.single(input),
            Rtmp = as.double(Routput),
            Itmp = as.double(Ioutput),
            as.integer(nvoice),
	    as.single(freqstep),
            as.integer(newsize),
            as.single(scale))

   Routput <- z$Rtmp
   Ioutput <- z$Itmp
   dim(Routput) <- c(newsize,pp)
   dim(Ioutput) <- c(newsize,pp)

       
   i <- sqrt(as.complex(-1))
   output <- Routput[1:isize,] + Ioutput[1:isize,] * i
   if(plot) {
      image(Mod(output),xlab="Time",ylab="Frequency")
      title("Gabor Transform Modulus")
   }
   output
}



vgt <- function(input, frequency, scale, plot = F)
#########################################################################
#      vgt:   
#      ---
#       continuous Gabor transform on one frequency:
# 	 compute the continuous Gabor transform with (complex-valued)
#	 gaussian window
#
#       input:
#       ------
# 	 input: input signal (possibly complex-valued)
#        frequency: value of the frequency
#        scale: size of the window
#	 plot: if set to TRUE, plotss the real part of cgt on the graphic
#		device.
#
#       output:
#       -------
#        Routput + i Ioutput: voice gabor transform (complex 1D array)
#
#########################################################################
{
   oldinput <- input
   isize <- length(oldinput)

   tmp <- adjust.length(oldinput)
   input <- tmp$signal
   newsize <- length(input)

   Routput <- numeric(newsize)
   Ioutput <- numeric(newsize)
   dim(input) <- c(newsize,1)

   z <- .C("Svgabor",
            as.single(input),
            Rtmp = as.double(Routput),
            Itmp = as.double(Ioutput),
            as.single(frequency),
            as.integer(newsize),
            as.single(scale))

   Routput <- z$Rtmp
   Ioutput <- z$Itmp
   i <- sqrt(as.complex(-1))
   if(plot==T) {
      plot.ts(Re(z$tmp));
      title("Real part of Gabor transform");
   }

   Routput[1:isize] + Ioutput[1:isize] * i
}



gabor <- function(sigsize, location, frequency, scale)
#########################################################################
#       gabor:   
#       ------
#        Generates a Gabor for given location and frequency
#
#       input:
#       ------
#        sigsize: signal size (dimension of the array)
#        location: location of the wavelet
#        frequency: value of the frequency
#        scale: size of the window
#
#       output:
#       -------
#        z$gabor.r + z$gabor.i * i: gabor (complex 1D array
#          of size sigsize)
#
#########################################################################
{
   gabor.r <- numeric(sigsize)
   gabor.i <- numeric(sigsize)


   z <- .C("gabor_time",
            as.single(frequency),
            as.single(scale),
            as.integer(location),
            gabor.r = as.double(gabor.r),
            gabor.i = as.double(gabor.i),
            as.integer(sigsize))

   i <- sqrt(as.complex(-1))

   z$gabor.r + z$gabor.i * i
}


vecgabor <- function(sigsize, nbnodes, location, frequency, scale)
#########################################################################
#       vecgabor:   
#       --------
#        Generates Gabor functions for given locations and frequencies
#	 on a ridge.
#
#       input:
#       ------
#        sigsize: signal size (dimension of the array)
#        nbnodes: number of ridge samples
#        location: b coordinates of the ridge samples
#        frequency: acoordinates of the ridge samples
#        scale: size of the window
#
#       output:
#       -------
#        z$gabor.r + z$gabor.i * i: 2D array containing the gabor
#           functions located on the ridge
#
#########################################################################
{

   
   gabor.r <- numeric(nbnodes * sigsize)
   gabor.i <- numeric(nbnodes * sigsize)


   z <- .C("vgabor_time",
            as.single(frequency),
            as.single(scale),
            as.integer(location),
            gabor.r = as.double(gabor.r),
            gabor.i = as.double(gabor.i),
            as.integer(sigsize),
            as.integer(nbnodes))

   i <- sqrt(as.complex(-1))

   z$gabor.r + z$gabor.i * i
}











#########################################################################
#      $Log: Pca_Climbers.S,v $
#
#               (c) Copyright  1997
#                          by                                   
#      Author: Rene Carmona, Bruno Torresani, Wen-Liang Hwang   
#                  Princeton University
#                  All right reserved                           
#########################################################################





hescrc <- function(tfrep, tfspec = numeric(dim(tfrep)[2]), grida = 10, gridb = 20,
        bstep = 3, iteration = 10000, rate = .001, seed = -7, 
	nbclimb = 10, flag.int = TRUE, chain= TRUE, flag.temp=FALSE,lineflag=F)
#########################################################################
#  hescrc:   Time-frequency multiple ridge estimation (hessian climbers)
#  ------
# 	use the hessian climber algorithm to evaluate ridges of
#    	   continuous Time-Frequency transform
#
#      input:
#      ------
# 	tfrep: wavelet or Gabor transform
#	tfspec: additional potential (coming from learning the noise)
#       grida : window of a
#       gridb : window of b
#       pct : the percent number of points in window 2*grida and 2*gridb to select histogram
#       count : the maximal number of repetitive selection, if location was chosen before
#       iteration: number of iterations
#       rate: initial value of the temperature
#       seed: initialization for random numbers generator
#       nbclimb: number of crazy climbers
#       bstep: step size for the climber walk in the time direction
#	flag.int: if set to TRUE, computes the integral on the ridge.
#	chain: if set to TRUE, chains the ridges.
#	flag.temp: if set to TRUE, keeps a constant temperature.
#       linfflag: if set to TRUE, the line segments oriented to the
#                 where the climber will move freely at the block is shown in 
#                 image   	
#       
#      output:
#      -------
#       beemap: 2D array containing the (weighted or unweighted)
#               occupation measure (integrated with respect to time)
#       hesmap: 2D array containing number from 1 to 4, denoting the
#               direction a climber will go at some position; where
#               1 moving freely along b, 3 freely along a, 2 along 
#               line a = -b, and 4 along the line a = b. The restricted
#               move is perdendicular to the free move.
#
#########################################################################
{

   tfspectrum <- tfspec 

   d <- dim(tfrep)
   sigsize <- d[1]
   nscale <- d[2]

   sqmodulus <- Re(tfrep*Conj(tfrep))
   for (k in 1:nscale)
     sqmodulus[,k] <- sqmodulus[,k] - tfspectrum[k]
   image(sqmodulus)

  # number of grid size
   nbx <- as.integer(sigsize/gridb)
   nby <- as.integer(nscale/grida)
   if((sigsize/gridb - nbx) > 0) nbx <- nbx + 1
   if((nscale/grida - nby) > 0) nby <- nby + 1

   nbblock <- nbx * nby

  # first two locations at each block stores the lower-left corner of the block
  # followed by the locations of up-right corner and by the negative values
  # of hessian matrix from xx, xy, yx, and yy at the last.

   tst <- matrix(0, 8, nbblock)
   dim(tst) <- c(nbblock * 8,1)

   dim(sqmodulus) <- c(sigsize * nscale, 1)

   z <- .C("Shessianmap",
	as.double(sqmodulus),
        as.integer(sigsize),
	as.integer(nscale),
	ublock = as.integer(nbblock),
	as.integer(gridb),
	as.integer(grida),
	tst =as.single(tst))

  
   tst <- z$tst
   dim(tst) <- c(8, nbblock)

  # actual number of block 	

   ublock <- z$ublock

  # principle component analysis and calculate the direction of each block


  pcamap <- matrix(1,sigsize,nscale)

  # first eigenvector of a block
   eigv1 <- matrix(0,ublock,2)

   # to draw eigenvector in a block 
   if(lineflag) {	
     lng <- min(grida, gridb)
     oldLng <- lng
     linex <- numeric(oldLng)
     liney <- numeric(oldLng)
   }


   for(j in 1:ublock) {
        
	left <- max(1,as.integer(tst[1,j]))
	down <- max(1,as.integer(tst[2,j]))
	right <- min(as.integer(tst[3,j]),sigsize)
	up <- min(as.integer(tst[4,j]),nscale)

	ctst <- tst[5:8,j]
	dim(ctst) <- c(2,2)
        ctst <- t(ctst)
        eig<- eigen(ctst)

	
#	cat("first eig = ",eig$values[1]," ratio in eigenvalues =  ",abs(eig$values[1]/eig$values[2])," \n")

        # eigen vector 1
        eigv1[j,] <- eig$vector[,1]


        # shifted center position of a block
	if(lineflag) {
          centerx <- as.integer((left + right)/2)
          centery <- as.integer((down + up)/2)
	}

        # theta of the eigen vector 1; -pi < theta <= pi

        theta <- atan((eigv1[j,])[2],(eigv1[j,])[1])

        # along x : denote 1 in the hesmap

        if(((theta <= pi/8) && (theta > -pi/8)) ||
	   (theta > 7*pi/8 || theta <= -7*pi/8)) {

	  pcamap[left:(right-1),down:(up-1)] <- 1

	  if(lineflag) {
             Lng <- min(lng,(sigsize-centerx))
             if(Lng != oldLng) {
               linex <- numeric(Lng)
               liney <- numeric(Lng)
             }
             llng <- as.integer(Lng/2)
	     x1 <- max(centerx-llng,1)
             x2 <- min(sigsize,x1 + Lng-1)
             linex[] <- x1:x2
             liney[] <- centery
            lines(linex,liney)
	}
	}

        # along x=y : denoted as 4 in hesmap

        if(((theta > pi/8) && (theta <= 3*pi/8)) ||
	   ((theta > -7*pi/8) && (theta <= -5*pi/8))) {

	  pcamap[left:(right-1),down:(up-1)] <- 4

	  if(lineflag) {
            Lng <- min(lng,sigsize-centerx)
            Lng <- min(Lng,nscale-centery)
            if(Lng != oldLng) {
              linex <- numeric(Lng)
              liney <- numeric(Lng)
            }
            llng <- as.integer(Lng/2)
            x1 <- max(centerx-llng,1)
            x2 <- min(sigsize,x1 + Lng-1)
	    y1 <- max(centery-llng,1)
            y2 <- min(sigsize,y1 + Lng-1)

            linex[] <- x1:x2
            liney[] <- y1:y2
            lines(linex,liney)
           }
	}


        # along y : as 3 in hesmap

        if(((theta > 3*pi/8) && (theta <= 5*pi/8)) ||
	   ((theta > -5*pi/8) && (theta <= -3*pi/8))) {

	  pcamap[left:(right-1),down:(up-1)] <- 3

	  if(lineflag) {
             Lng <- min(lng, nscale-centery)
             if(Lng != oldLng) {
               linex <- numeric(Lng)
               liney <- numeric(Lng)
             }
            llng <- as.integer(Lng/2)
	    y1 <- max(centery-llng,1)
            y2 <- min(sigsize,y1 + Lng-1)

            linex[] <- centerx
            liney[] <- y1:y2
            lines(linex,liney)
          }
	}


        # along x=-y : as 2 in hesmap

        if(((theta > 5*pi/8) && (theta <= 7*pi/8)) ||
	   ((theta > -3*pi/8) && (theta <= -pi/8))) {

	  pcamap[left:(right-1),down:(up-1)] <- 2

	  if(lineflag) {
            Lng <- min(lng,nscale-centery)
            Lng <- min(Lng,sigsize-centerx)
            if(Lng != oldLng) {
               linex <- numeric(Lng)
               liney <- numeric(Lng)
            }
            llng <- as.integer(Lng/2)

	    x1 <- max(centerx-llng,1)
            x2 <- min(sigsize,x1 + Lng-1)
	    y1 <- max(centery-llng,1)
            y2 <- min(sigsize,y1 + Lng-1)

            linex[] <- x1:x2
            liney[] <- y2:y1
            lines(linex,liney)
          }
	}

	  if(lineflag) oldLng <- Lng
   }

#   if(lineflag==F) image(pcamap)	

# From the following on, it is similar to the process of crazy climbers ....

   beemap <- matrix(0,sigsize,nscale) 
   dim(beemap) <- c(nscale * sigsize, 1)
   dim(pcamap) <- c(nscale * sigsize, 1)
     
   z <- .C("Spca_annealing",
        as.double(sqmodulus),
        beemap= as.double(beemap),
	as.integer(pcamap),
        as.single(rate),
	as.integer(sigsize),
        as.integer(nscale),
        as.integer(iteration),
	as.integer(seed),
        as.integer(bstep),
        as.integer(nbclimb),
	as.integer(flag.int),
	as.integer(chain),
	as.integer(flag.temp))

   beemap <- z$beemap
   dim(beemap) <- c(sigsize, nscale)
   dim(pcamap) <- c(sigsize, nscale)
   image(beemap)
   list(beemap = beemap, pcamap = pcamap)
}

   



#########################################################################
#
#               (c) Copyright  1997
#                          by                                   
#      Author: Rene Carmona, Bruno Torresani, Wen-Liang Hwang   
#                  Princeton University
#                  All right reserved                           
#########################################################################




dw <- function(inputdata, maxresoln, scale=FALSE, NW=6, plot=TRUE)
#*********************************************************************#
# dw computes the Daubechies wavelet transform from resolution 1 to
# the given maximum resolution.  If phi.flag is TRUE, it returns a 
# list consisting psi and phi wavelet transforms as two matrices.
#*********************************************************************#
{
  if ( ! ((NW > 1) && (NW <11)) )
    stop("NW has to be between 2 and 10")

  x <- adjust.length(inputdata)
  s <- x$signal
  np <- x$length
  
  check.maxresoln(maxresoln, np)

  phi <- matrix(0, nrow=maxresoln+1, ncol=np)
  psi <- matrix(0, nrow=maxresoln, ncol=np)

  # Convert matrices into vectors
  phi <- t(phi)
  dim(phi) <- c(length(phi), 1)
  psi <- t(psi)
  dim(psi) <- c(length(psi), 1)

  z <- .C("daubechies_wt", 
	phi=as.single(phi),
        psi=as.single(psi),
	as.single(s),
	as.integer(NW),
	as.integer(maxresoln),
	as.integer(np) )

  # Convert vectors into original matrices
  phi <- t(z$phi)
  dim(phi) <- c(np, maxresoln+1)
  psi <- t(z$psi)
  dim(psi) <- c(np, maxresoln)

  if(plot) plotwt(s, psi, phi, maxresoln, scale)

  phi <- matrix(phi[,(maxresoln+1)], ncol=1, nrow=np)
  list( original=s, Wf=psi, Sf=phi, maxresoln=maxresoln, np=np, NW=NW )
}


ddw <- function(inputdata, maxresoln, scale=FALSE, NW=6 )
#*********************************************************************#
# ddw computes the discrete Daubechies wavelet transform from 
# resolution 1 to the given maximum resolution.  
#*********************************************************************#
{
  if ( ! ((NW > 1) && (NW <11)) )
    stop("NW has to be between 2 and 10")

  x <- adjust.length(inputdata)
  s <- x$signal
  np <- x$length	

  check.maxresoln(maxresoln, np)

  phi <- matrix(0, nrow=maxresoln+1, ncol=np)
  psi <- matrix(0, nrow=maxresoln, ncol=np)

  # Convert matrices into vectors
  phi <- t(phi)
  dim(phi) <- c(length(phi), 1)
  psi <- t(psi)
  dim(psi) <- c(length(psi), 1)

  z <- .C("compute_ddwave", 
	phi=as.single(phi),
	psi=as.single(psi),
	as.single(s),
	as.integer(maxresoln),
	as.integer(np),
	as.integer(NW) )

  # Convert vectors into original matrices
  phi <- t(z$phi)
  dim(phi) <- c(np, maxresoln+1)
  psi <- t(z$psi)
  dim(psi) <- c(np, maxresoln)

  plotwt(s, psi, phi, maxresoln, scale)
  
  phi <- matrix(phi[,(maxresoln+1)], ncol=1, nrow=np)
  list( original=s, Wf=psi, Sf=phi, maxresoln=maxresoln, np=np, NW=NW )
}











#########################################################################
#      $Log: MGabor.S,v $
#########################################################################
#
#               (c) Copyright  1997
#                          by                                   
#      Author: Rene Carmona, Bruno Torresani, Wen-Liang Hwang   
#                  Princeton University 
#                  All right reserved                           
#########################################################################





mcgt <- function(input, nvoice, freqstep = (1/nvoice), nscales = 10,
	scalestep = 5, initscale = 0, crit = 0, plot = TRUE,
	tchatche = FALSE)
#########################################################################
#       mcgt:   
#       ----
#        multiwindow continuous Gabor transform function:
# 	  compute the continuous Gabor transform with gaussian windows.
#	  selects the optimal window using L^p norm or entropy
#         optimization
#
#       input:
#       ------
# 	 input: input signal (possibly complex-valued)
#	 nvoice: number of frequency bands
#        freqstep: sampling rate for the frequency axis
#        nscales: number of scales considered
#        scalestep: increment for scale parameter
#        initscale: initial value form the scale parameter
#        crit: criterion for optimization (L^p norm or entropy)
#	 plot: if set to TRUE, displays the modulus of cwt on the graphic
#		device.
#        rchatche: if set to TRUE, prints intermediate results
#
#       output:
#       -------
#        output: optimal continuous (complex) Gabor transform
#
#########################################################################
{
   oldinput <- input
   isize <- length(oldinput)

   tmp <- adjust.length(oldinput)
   input <- tmp$signal
   newsize <- length(input)

   pp <- nvoice
   Routput <- matrix(0,newsize,pp)
   Ioutput <- matrix(0,newsize,pp)
   output <- matrix(0,newsize,pp)
   dim(Routput) <- c(pp * newsize,1)
   dim(Ioutput) <- c(pp * newsize,1)
   dim(input) <- c(newsize,1)


   norm <- 1
   i <- sqrt(as.complex(-1))
   lpoptnorm <- 0
   entoptnorm <- 100000000

   optsca <- 0


#####################
# Loop over scales: #
#####################

   sca <- initscale
   for(k in 1:nscales){

      sca <- sca + scalestep

# Compute Gabor transform
# -----------------------
      z <- .C("Sgabor",
            as.single(input),
            Rtmp = as.double(Routput),
            Itmp = as.double(Ioutput),
            as.integer(nvoice),
	    as.single(freqstep),
            as.integer(newsize),
            as.single(sca))

      Routput <- z$Rtmp
      Ioutput <- z$Itmp
      dim(Routput) <- c(newsize,pp)
      dim(Ioutput) <- c(newsize,pp)

# Compute L^2 norm for normalization
# ----------------------------------
      pexp <- as.double(2)
      z <- .C("Lpnorm",
            ltwonorm = as.double(norm),
	    as.double(pexp),
            as.double(Routput),
            as.double(Ioutput),
            as.integer(newsize),
            as.integer(nvoice))
#    cat("l2 norm=",z$ltwonorm,"\n")

# Normalize
# ---------
      Routput <- Routput/z$ltwonorm
      Ioutput <- Ioutput/z$ltwonorm

      if (crit == 0){
         z <- .C("entropy",
            lpnorm = as.double(norm),
            as.double(Routput),
            as.double(Ioutput),
            as.integer(newsize),
            as.integer(nvoice))
         if(tchatche){
            cat("     scale=",sca,"; entropy=",z$lpnorm,"\n")
         }

         if (z$lpnorm < entoptnorm){
            entoptnorm <- z$lpnorm
            optsca <- sca
            output <- Routput[1:isize,] + Ioutput[1:isize,] * i
         }
      optnorm <- entoptnorm
      }

      else {

# Compute L^p norm
# ----------------
         pexp <- as.double(crit)
         z <- .C("Lpnorm",
            lpnorm = as.double(norm),
	    as.double(pexp),
            as.double(Routput),
            as.double(Ioutput),
            as.integer(newsize),
            as.integer(nvoice))
         if(tchatche){
            cat("     scale=",sca,"; l",pexp," norm=",z$lpnorm,"\n")
         }

         if (z$lpnorm > lpoptnorm){
            lpoptnorm <- z$lpnorm
            optsca <- sca
            output <- Routput[1:isize,] + Ioutput[1:isize,] * i
         }
      optnorm <- lpoptnorm
      }
   }

   cat("   Optimal scale: ",optsca,"\n")

   if(plot) {
      image(Mod(output),xlab="Time",ylab="Frequency")
      title("Gabor Transform Modulus")
   }

   output

}


#########################################################################
#       $Log: Noise.S,v $
# Revision 1.2  1995/04/05  18:56:55  bruno
# *** empty log message ***
#
# Revision 1.1  1995/04/02  01:04:16  bruno
# Initial revision
#
#               (c) Copyright  1997
#                          by                                   
#      Author: Rene Carmona, Bruno Torresani, Wen-Liang Hwang   
#                  Princeton University
#                  All right reserved                           
#########################################################################




tfpct <- function(input,percent=.8,plot=TRUE)
#########################################################################
#       tfpct:   
#       ------
# 	 compute a percentile of time-frequency representation frequency 
#		by frequency
#
#       input:
#       ------
# 	 input: modulus of the continuous wavelet transform (2D array)
#	 percent: value of the percentile to be computed
#        plot: if set to TRUE, displays the values of the energy as a 
#                  function of the scale
#
#       output:
#       -------
#        output: 1D array of size nbscales containing the noise estimate
#
#########################################################################
{

  nscale <- dim(input)[2]

  output <- numeric(nscale)
  for(i in 1:nscale) output[i] <- quantile(input[,i],percent)

  if(plot)plot.ts(output)
  output
}


tfmean <- function(input,plot=TRUE)
#########################################################################
#       tfmean:   
#       -------
# 	 compute the mean of time-frequency representation frequency 
#		by frequency
#
#       input:
#       ------
# 	 input: modulus of the continuous wavelet transform (2D array)
#        plot: if set to TRUE, displays the values of the energy as a 
#                  function of the scale
#
#       output:
#       -------
#        output: 1D array of size nbscales containing the noise estimate
#
#########################################################################
{

  nscale <- dim(input)[2]

  output <- numeric(nscale)
  for(i in 1:nscale) output[i] <- mean(input[,i])

  if(plot) plot.ts(output)
  output
}


tfvar <- function(input,plot=TRUE)
#########################################################################
#       tfvar:   
#       ------
# 	 compute the variance of time-frequency representation frequency 
#		by frequency
#
#       input:
#       ------
# 	 input: modulus of the continuous wavelet transform (2D array)
#        plot: if set to TRUE, displays the values of the energy as a 
#                  function of the scale
#
#       output:
#       -------
#        output: 1D array of size nbscales containing the noise estimate
#
#########################################################################
{

  nscale <- dim(input)[2]

  output <- numeric(nscale)
  for(i in 1:nscale) output[i] <- var(input[,i])

  if(plot) plot.ts(output)
  output
}




















#########################################################################
#      $Log: Pca_Climbers.S,v $
#
#               (c) Copyright  1997
#                          by                                   
#      Author: Rene Carmona, Bruno Torresani, Wen-Liang Hwang   
#                  Princeton University 
#                  All right reserved                           
#########################################################################





pcacrc <- function(tfrep, tfspec = numeric(dim(tfrep)[2]), grida = 10,
	gridb = 20, pct = 0.1, count = 10000, bstep = 3, iteration = 10000,
	rate = .001, seed = -7, nbclimb = 10, flag.int = TRUE, chain= TRUE,
	flag.temp=FALSE,lineflag=F,flag.cwt=F)
#########################################################################
#  pcacrc:   Time-frequency multiple ridge estimation (pca climbers)
#  ------
# 	use the pca climber algorithm to evaluate ridges of
#    	   continuous Time-Frequency transform
#
#      input:
#      ------
# 	tfrep: wavelet or Gabor transform
#	tfspec: additional potential (coming from learning the noise)
#       grida : window of a
#       gridb : window of b
#       pct : the percent number of points in window 2*grida and
#	  2*gridb to select histogram
#       count : the maximal number of repetitive selection, if
#	  location was chosen before
#       iteration: number of iterations
#       rate: initial value of the temperature
#       seed: initialization for random numbers generator
#       nbclimb: number of crazy climbers
#       bstep: step size for the climber walk in the time direction
#	flag.int: if set to TRUE, computes the integral on the ridge.
#	chain: if set to TRUE, chains the ridges.
#	flag.temp: if set to TRUE, keeps a constant temperature.
#	flag.cwt: if set to TRUE, the bstep as well as the block size
#	 of a and b are adjusted according to the scale of wavelet
#	transform(not implemented yet).
#       linfflag: if set to TRUE, the line segments oriented to the
#                 where the climber will move freely at the block is shown in 
#                 image   	
#       
#      output:
#      -------
#       beemap: 2D array containing the (weighted or unweighted)
#               occupation measure (integrated with respect to time)
#       pcamap: 2D array containing number from 1 to 4, denoting the
#               direction a climber will go at some position; where
#               1 moving freely along b, 3 freely along a, 2 along 
#               line a = -b, and 4 along the line a = b. The restricted
#               move is perdendicular to the free move.
#
#########################################################################
{

   tfspectrum <- tfspec 

   d <- dim(tfrep)
   sigsize <- d[1]
   nscale <- d[2]

   sqmodulus <- Re(tfrep*Conj(tfrep))
   for (k in 1:nscale)
     sqmodulus[,k] <- sqmodulus[,k] - tfspectrum[k]
   image(sqmodulus)

  # percentage of points to be selected in a block  
  # considering overlapping 1/2 on each coordinate with neighborhood 
   nbpoint <- as.integer(2* grida * 2 * gridb * pct)
#   nbpoint <- as.integer(grida * gridb * pct)

  # number of grid size
   nbx <- as.integer(sigsize/gridb)
   nby <- as.integer(nscale/grida)
   if((sigsize/gridb - nbx) > 0) nbx <- nbx + 1
   if((nscale/grida - nby) > 0) nby <- nby + 1

   nbblock <- nbx * nby

  # first two locations at each block stores the lower-left corner of the block
  # followed by the locations of up-right corner and by the coordinates of the 
  # selected points in the block. The locations are in the order of (x, y)

   tstsize <- 2 * nbpoint + 4

   tst <- matrix(0, tstsize, nbblock)
   dim(tst) <- c(nbblock * tstsize,1)

  # the map of points selected in all the block
   pointmap <- matrix(0,sigsize,nscale)
   dim(pointmap) <- c(sigsize * nscale, 1)
   dim(sqmodulus) <- c(sigsize * nscale, 1)

   z <- .C("Spointmap",
	as.double(sqmodulus),
        as.integer(sigsize),
	as.integer(nscale),
	as.integer(gridb),
	as.integer(grida),
	as.integer(nbblock),
	as.integer(nbpoint),
        pointmap = as.integer(pointmap),
	tst =as.single(tst),
	as.integer(tstsize),
	as.integer(count),
	as.integer(seed))
  
   pointmap <- z$pointmap
   tst <- z$tst
   dim(pointmap) <- c(sigsize,nscale)
   dim(tst) <- c(tstsize, nbblock)

   

  # principle component analysis and calculate the direction at each block


  pcamap <- matrix(1,sigsize,nscale)

  # first eigenvector of a block
   eigv1 <- matrix(0,nbblock,2)

   # to draw eigenvector in a block 
   if(lineflag) {	
     lng <- min(grida, gridb)
     oldLng <- lng
     linex <- numeric(oldLng)
     liney <- numeric(oldLng)
   }


   for(j in 1:nbblock) {
        
	left <- max(1,as.integer(tst[1,j]))
	down <- max(1,as.integer(tst[2,j]))
	right <- min(as.integer(tst[3,j]),sigsize)
	up <- min(as.integer(tst[4,j]),nscale)

	input <- tst[5:tstsize,j]
	dim(input) <- c(2,nbpoint)
        input <- t(input)
        ctst <- var(input)
        eig<- eigen(ctst)

	# cat(" ratio in eigenvalues =  ",abs(eig$values[1]/eig$values[2])," \n")

        # eigen vector 1
        eigv1[j,] <- eig$vector[,1]


        # shifted center position of a block
	if(lineflag) {
          centerx <- as.integer((left + right)/2)
          centery <- as.integer((down + up)/2)
	}

        # theta of the eigen vector 1; -pi < theta <= pi

        theta <- atan((eigv1[j,])[2],(eigv1[j,])[1])

        # along x : denote 1 in the pcamap

        if(((theta <= pi/8) && (theta > -pi/8)) ||
	   (theta > 7*pi/8 || theta <= -7*pi/8)) {

	  pcamap[left:(right-1),down:(up-1)] <- 1

	  if(lineflag) {
             Lng <- min(lng,(sigsize-centerx))
             if(Lng != oldLng) {
               linex <- numeric(Lng)
               liney <- numeric(Lng)
             }
             llng <- as.integer(Lng/2)
	     x1 <- max(centerx-llng,1)
             x2 <- min(sigsize,x1 + Lng-1)
             linex[] <- x1:x2
             liney[] <- centery
            lines(linex,liney)
	}
	}

        # along x=y : denoted as 2 in pcamap

        if(((theta > pi/8) && (theta <= 3*pi/8)) ||
	   ((theta > -7*pi/8) && (theta <= -5*pi/8))) {

	  pcamap[left:(right-1),down:(up-1)] <- 2

	  if(lineflag) {
            Lng <- min(lng,sigsize-centerx)
            Lng <- min(Lng,nscale-centery)
            if(Lng != oldLng) {
              linex <- numeric(Lng)
              liney <- numeric(Lng)
            }
            llng <- as.integer(Lng/2)
            x1 <- max(centerx-llng,1)
            x2 <- min(sigsize,x1 + Lng-1)
	    y1 <- max(centery-llng,1)
            y2 <- min(sigsize,y1 + Lng-1)

            linex[] <- x1:x2
            liney[] <- y1:y2
            lines(linex,liney)
           }
	}


        # along y : as 3 in pcamap

        if(((theta > 3*pi/8) && (theta <= 5*pi/8)) ||
	   ((theta > -5*pi/8) && (theta <= -3*pi/8))) {

	  pcamap[left:(right-1),down:(up-1)] <- 3

	  if(lineflag) {
             Lng <- min(lng, nscale-centery)
             if(Lng != oldLng) {
               linex <- numeric(Lng)
               liney <- numeric(Lng)
             }
            llng <- as.integer(Lng/2)
	    y1 <- max(centery-llng,1)
            y2 <- min(sigsize,y1 + Lng-1)

            linex[] <- centerx
            liney[] <- y1:y2
            lines(linex,liney)
          }
	}


        # along x=-y : as 4 in pcamap

        if(((theta > 5*pi/8) && (theta <= 7*pi/8)) ||
	   ((theta > -3*pi/8) && (theta <= -pi/8))) {

	  pcamap[left:(right-1),down:(up-1)] <- 4

	  if(lineflag) {
            Lng <- min(lng,nscale-centery)
            Lng <- min(Lng,sigsize-centerx)
            if(Lng != oldLng) {
               linex <- numeric(Lng)
               liney <- numeric(Lng)
            }
            llng <- as.integer(Lng/2)

	    x1 <- max(centerx-llng,1)
            x2 <- min(sigsize,x1 + Lng-1)
	    y1 <- max(centery-llng,1)
            y2 <- min(sigsize,y1 + Lng-1)

            linex[] <- x1:x2
            liney[] <- y2:y1
            lines(linex,liney)
          }
	}

	  if(lineflag) oldLng <- Lng
   }

# From the following on, it is similar to the process of crazy climbers ....

   beemap <- matrix(0,sigsize,nscale) 
   dim(beemap) <- c(nscale * sigsize, 1)
   dim(pcamap) <- c(nscale * sigsize, 1)
     
   z <- .C("Spca_annealing",
        as.double(sqmodulus),
        beemap= as.double(beemap),
	as.integer(pcamap),
        as.single(rate),
	as.integer(sigsize),
        as.integer(nscale),
        as.integer(iteration),
	as.integer(seed),
        as.integer(bstep),
        as.integer(nbclimb),
	as.integer(flag.int),
	as.integer(chain),
	as.integer(flag.temp))

   beemap <- z$beemap
   dim(beemap) <- c(sigsize, nscale)
   dim(pcamap) <- c(sigsize, nscale)
   image(beemap)
   list(beemap = beemap, pcamap = pcamap)
}

   


pcafamily <-function(pcaridge,orientmap,maxchnlng=as.numeric(dim(pcaridge)[1])+10,bstep= 1, nbchain = 100, ptile = 0.05)
#########################################################################
#     pcafamily:
#     ---------
#     chain the ridges obtained by pca climber.
#
#      input:
#      ------
#       pcaridge: ridge found by pca climber
#       orientmap : the first eigen vector direction at each position
#       bstep: maximal length for a gap in a ridge
#       nbchain: maximum number of chains
#       ptile: relative threshold for the ridges 
#       maxchnlng: maximal chain length
#
#      output:
#      ------
#       ordered: ordered map: image containing the ridges
#                (displayed with different colors)
#       chain: 2D array containing the chained ridges, according
#		to the chain data structure:
#		chain[,1]: first point of the ridge
#		chain[,2]: length of the chain
#		chain[,3:(chain[,2]+2)]: values of the ridge
#	nbchain: number of chains.
#
#########################################################################
{
   d <- dim(pcaridge)
   sigsize <- d[1]
   nscale <- d[2]
   threshold <- range(pcaridge)[2] * ptile
   sz <- 2 * (maxchnlng);
   chain <- matrix(-1,nbchain,sz)
   orderedmap <- matrix(0,sigsize,nscale)	
   dim(chain) <- c(nbchain * sz, 1)
   dim(pcaridge) <- c(nscale * sigsize, 1)
   dim(orderedmap) <- c(nscale * sigsize, 1)
   dim(orientmap) <- c(nscale * sigsize, 1)
	
   
   z <- .C("Spca_family",
        as.double(pcaridge),
        as.integer(orientmap),
        orderedmap = as.single(orderedmap),
	chain = as.integer(chain),
	chainnb = as.integer(nbchain),
	as.integer(sigsize),
        as.integer(nscale),
        as.integer(bstep),
	as.single(threshold),
	as.integer(maxchnlng))
	
   orderedmap <- z$orderedmap
   chain <- z$chain

   dim(orderedmap) <- c(sigsize,nscale)   
   dim(chain) <- c(nbchain,sz)   
   nbchain <- z$chainnb
   chain <- chain +1
   chain[,1] <- chain[,1]-1

   list(ordered = orderedmap, chain = chain, nbchain=nbchain)
}   



pcamaxima <- function(beemap,orientmap)
{
   d <- dim(beemap)
   sigsize <- d[1]
   nscale <- d[2]  
   tfmaxima <- matrix(0,sigsize,nscale)
   dim(orientmap) <- c(nscale * sigsize, 1)
   dim(beemap) <- c(nscale * sigsize, 1)   
   dim(tfmaxima) <- c(nscale * sigsize, 1)   

   z <- .C("Stf_pcaridge",
        as.double(beemap),
        tfmaxima = as.double(tfmaxima),
	as.integer(sigsize),
        as.integer(nscale),
        as.integer(orientmap))
	
   tfmaxima <- z$tfmaxima
   dim(tfmaxima) <- c(sigsize,nscale)   
   tfmaxima
}


simplepcarec <- function(siginput, tfinput, beemap, orientmap, bstep = 5,
	ptile = .01, plot = 2)
#########################################################################
#     simplepcarec:
#     -------------
#      Simple reconstruction of a real valued signal from ridges found by 
#      pca climbers.
#
#      input:
#      ------
#      siginput: input signal
#      tfinput: continuous time-frequency transform (output of cwt or cgt)
#      beemap: output of pca climber algorithm
#      orientmap: pca dirction of climber
#      bstep: used for the chaining
#      ptile: 
#      plot: if set to 1, displays the signal, the components, and
#            signal and reconstruction one after another. If set to
#            2, displays the signal, the components, and the
#            reconstruction on the same page. Else, no plot.
#
#      output:
#      -------
#      rec: reconstructed signal
#      ordered: image of the ridges (with different colors)
#      comp: 2D array containing the signals reconstructed from ridges
#
#########################################################################
{
   tmp <- pcafamily(beemap,orientmap,ptile=ptile)
   chain <- tmp$chain
   nbchain <- tmp$nbchain
   rec <- numeric(length(siginput))

   if(plot != F){
     par(mfrow=c(2,1))
     plot.ts(siginput)
     title("Original signal")
     image(tmp$ordered)
     title("Chained Ridges")
   }

   npl(1)

   sigsize <- dim(tfinput)[1]
   nscale <- dim(tfinput)[2]
   rec <- numeric(sigsize)

   tmp3 <- matrix(0+0i,nbchain,sigsize)

   for (j in 1:nbchain){
      for(i in 1:chain[j,1]) {
        cat("The (b,a) at chain", j,"is (",chain[j,2*i+1],chain[j,2*i],")\n")
        tmp3[j,chain[j,2*i+1]] <- tmp3[j,chain[j,2*i+1]]+
                                  2*tfinput[chain[j,2*i+1],chain[j,2*i]]
      }
      rec <- rec + Re(tmp3[j,])
   }

   if(plot == 1){
      par(mfrow=c(2,1))
      par(cex=1.1)
      plot.ts(siginput)
      title("Original signal")
      plot.ts(rec)
      title("Reconstructed signal")
   }
   else if (plot == 2){
      par(mfrow=c(nbchain+2,1))
      par(mar=c(2,4,4,4))
      par(cex=1.1)
      par(err=-1)
      plot.ts(siginput)
      title("Original signal")

      for (j in 1:nbchain)
         plot.ts(Re(tmp3[j,]))
        
      plot.ts(rec)
      title("Reconstructed signal")
   }
   list(rec=rec, ordered=tmp$ordered,chain = tmp$chain, comp=tmp3)
}





#########################################################################
#      $Log: Pca_Climbers.S,v $
#
#               (c) Copyright  1997
#                          by                                   
#      Author: Rene Carmona, Bruno Torresani, Wen-Liang Hwang   
#                  Princeton University 
#                  All right reserved                           
#########################################################################





pcacrc <- function(tfrep, tfspec = numeric(dim(tfrep)[2]), grida = 10,
	gridb = 20, pct = 0.1, count = 10000, bstep = 3, iteration = 10000,
	rate = .001, seed = -7, nbclimb = 10, flag.int = TRUE, chain= TRUE,
	flag.temp=FALSE,lineflag=F,flag.cwt=F,nvoice=0)
#########################################################################
#  pcacrc:   Time-frequency multiple ridge estimation (pca climbers)
#  ------
# 	use the pca climber algorithm to evaluate ridges of
#    	   continuous Time-Frequency transform
#
#      input:
#      ------
# 	tfrep: wavelet or Gabor transform
#	tfspec: additional potential (coming from learning the noise)
#       grida : window of a
#       gridb : window of b
#       pct : the percent number of points in window 2*grida and
#	  2*gridb to select histogram
#       count : the maximal number of repetitive selection, if
#	  location was chosen before
#       iteration: number of iterations
#       rate: initial value of the temperature
#       seed: initialization for random numbers generator
#       nbclimb: number of crazy climbers
#       bstep: step size for the climber walk in the time direction
#	flag.int: if set to TRUE, computes the integral on the ridge.
#	chain: if set to TRUE, chains the ridges.
#	flag.temp: if set to TRUE, keeps a constant temperature.
#	flag.cwt: if set to TRUE, the bstep as well as the block size
#	 of a and b are adjusted according to the scale of wavelet
#	transform(not implemented yet).
#       nvoice: number of voice per octave (cwt)
#       linfflag: if set to TRUE, the line segments oriented to the
#                 where the climber will move freely at the block is shown in 
#                 image   	
#       
#      output:
#      -------
#       beemap: 2D array containing the (weighted or unweighted)
#               occupation measure (integrated with respect to time)
#       pcamap: 2D array containing number from 1 to 4, denoting the
#               direction a climber will go at some position; where
#               1 moving freely along b, 3 freely along a, 2 along 
#               line a = -b, and 4 along the line a = b. The restricted
#               move is perdendicular to the free move.
#
#########################################################################
{

   tfspectrum <- tfspec 

   d <- dim(tfrep)
   sigsize <- d[1]
   nscale <- d[2]

   sqmodulus <- Re(tfrep*Conj(tfrep))
   for (k in 1:nscale)
     sqmodulus[,k] <- sqmodulus[,k] - tfspectrum[k]
   image(sqmodulus)

  # percentage of points to be selected in a block  
  # considering overlapping 1/2 on each coordinate with neighborhood 
   nbpoint <- as.integer(2* grida * 2 * gridb * pct)
#   nbpoint <- as.integer(grida * gridb * pct)

  # number of grid size
   nbx <- as.integer(sigsize/gridb)
   nby <- as.integer(nscale/grida)
   if((sigsize/gridb - nbx) > 0) nbx <- nbx + 1
   if((nscale/grida - nby) > 0) nby <- nby + 1

   nbblock <- nbx * nby

  # first two locations at each block stores the lower-left corner of the block
  # followed by the locations of up-right corner and by the coordinates of the 
  # selected points in the block. The locations are in the order of (x, y)

   tstsize <- 2 * nbpoint + 4

   tst <- matrix(0, tstsize, nbblock)
   dim(tst) <- c(nbblock * tstsize,1)

  # the map of points selected in all the block
   pointmap <- matrix(0,sigsize,nscale)
   dim(pointmap) <- c(sigsize * nscale, 1)
   dim(sqmodulus) <- c(sigsize * nscale, 1)

   z <- .C("Spointmap",
	as.double(sqmodulus),
        as.integer(sigsize),
	as.integer(nscale),
	as.integer(gridb),
	as.integer(grida),
	as.integer(nbblock),
	as.integer(nbpoint),
        pointmap = as.integer(pointmap),
	tst =as.single(tst),
	as.integer(tstsize),
	as.integer(count),
	as.integer(seed),
        as.integer(flag.cwt),
        as.integer(nvoice))
  
   pointmap <- z$pointmap
   tst <- z$tst
   dim(pointmap) <- c(sigsize,nscale)
   dim(tst) <- c(tstsize, nbblock)

   

  # principle component analysis and calculate the direction at each block


  pcamap <- matrix(1,sigsize,nscale)

  # first eigenvector of a block
   eigv1 <- matrix(0,nbblock,2)

   # to draw eigenvector in a block 
   if(lineflag) {	
     lng <- min(grida, gridb)
     oldLng <- lng
     linex <- numeric(oldLng)
     liney <- numeric(oldLng)
   }


   for(j in 1:nbblock) {
        
	left <- max(1,as.integer(tst[1,j]))
	down <- max(1,as.integer(tst[2,j]))
	right <- min(as.integer(tst[3,j]),sigsize)
	up <- min(as.integer(tst[4,j]),nscale)

	input <- tst[5:tstsize,j]
	dim(input) <- c(2,nbpoint)
        input <- t(input)
        ctst <- var(input)
        eig<- eigen(ctst)

	# cat(" ratio in eigenvalues =  ",abs(eig$values[1]/eig$values[2])," \n")

        # eigen vector 1
        eigv1[j,] <- eig$vector[,1]


        # shifted center position of a block
	if(lineflag) {
          centerx <- as.integer((left + right)/2)
          centery <- as.integer((down + up)/2)
	}

        # theta of the eigen vector 1; -pi < theta <= pi

        theta <- atan((eigv1[j,])[2],(eigv1[j,])[1])

        # along x : denote 1 in the pcamap

        if(((theta <= pi/8) && (theta > -pi/8)) ||
	   (theta > 7*pi/8 || theta <= -7*pi/8)) {

	  pcamap[left:(right-1),down:(up-1)] <- 1

	  if(lineflag) {
             Lng <- min(lng,(sigsize-centerx))
             if(Lng != oldLng) {
               linex <- numeric(Lng)
               liney <- numeric(Lng)
             }
             llng <- as.integer(Lng/2)
	     x1 <- max(centerx-llng,1)
             x2 <- min(sigsize,x1 + Lng-1)
             linex[] <- x1:x2
             liney[] <- centery
            lines(linex,liney)
	}
	}

        # along x=y : denoted as 2 in pcamap

        if(((theta > pi/8) && (theta <= 3*pi/8)) ||
	   ((theta > -7*pi/8) && (theta <= -5*pi/8))) {

	  pcamap[left:(right-1),down:(up-1)] <- 2

	  if(lineflag) {
            Lng <- min(lng,sigsize-centerx)
            Lng <- min(Lng,nscale-centery)
            if(Lng != oldLng) {
              linex <- numeric(Lng)
              liney <- numeric(Lng)
            }
            llng <- as.integer(Lng/2)
            x1 <- max(centerx-llng,1)
            x2 <- min(sigsize,x1 + Lng-1)
	    y1 <- max(centery-llng,1)
            y2 <- min(sigsize,y1 + Lng-1)

            linex[] <- x1:x2
            liney[] <- y1:y2
            lines(linex,liney)
           }
	}


        # along y : as 3 in pcamap

        if(((theta > 3*pi/8) && (theta <= 5*pi/8)) ||
	   ((theta > -5*pi/8) && (theta <= -3*pi/8))) {

	  pcamap[left:(right-1),down:(up-1)] <- 3

	  if(lineflag) {
             Lng <- min(lng, nscale-centery)
             if(Lng != oldLng) {
               linex <- numeric(Lng)
               liney <- numeric(Lng)
             }
            llng <- as.integer(Lng/2)
	    y1 <- max(centery-llng,1)
            y2 <- min(sigsize,y1 + Lng-1)

            linex[] <- centerx
            liney[] <- y1:y2
            lines(linex,liney)
          }
	}


        # along x=-y : as 4 in pcamap

        if(((theta > 5*pi/8) && (theta <= 7*pi/8)) ||
	   ((theta > -3*pi/8) && (theta <= -pi/8))) {

	  pcamap[left:(right-1),down:(up-1)] <- 4

	  if(lineflag) {
            Lng <- min(lng,nscale-centery)
            Lng <- min(Lng,sigsize-centerx)
            if(Lng != oldLng) {
               linex <- numeric(Lng)
               liney <- numeric(Lng)
            }
            llng <- as.integer(Lng/2)

	    x1 <- max(centerx-llng,1)
            x2 <- min(sigsize,x1 + Lng-1)
	    y1 <- max(centery-llng,1)
            y2 <- min(sigsize,y1 + Lng-1)

            linex[] <- x1:x2
            liney[] <- y2:y1
            lines(linex,liney)
          }
	}

	  if(lineflag) oldLng <- Lng
   }

# From the following on, it is similar to the process of crazy climbers ....

   beemap <- matrix(0,sigsize,nscale) 
   dim(beemap) <- c(nscale * sigsize, 1)
   dim(pcamap) <- c(nscale * sigsize, 1)
     
   z <- .C("Spca_annealing",
        as.double(sqmodulus),
        beemap= as.double(beemap),
	as.integer(pcamap),
        as.single(rate),
	as.integer(sigsize),
        as.integer(nscale),
        as.integer(iteration),
	as.integer(seed),
        as.integer(bstep),
        as.integer(nbclimb),
	as.integer(flag.int),
	as.integer(chain),
	as.integer(flag.temp))

   beemap <- z$beemap
   dim(beemap) <- c(sigsize, nscale)
   dim(pcamap) <- c(sigsize, nscale)
   image(beemap)
   list(beemap = beemap, pcamap = pcamap)
}

   


pcafamily <-function(pcaridge,orientmap,maxchnlng=as.numeric(dim(pcaridge)[1])+10,bstep= 1, nbchain = 100, ptile = 0.05)
#########################################################################
#     pcafamily:
#     ---------
#     chain the ridges obtained by pca climber.
#
#      input:
#      ------
#       pcaridge: ridge found by pca climber
#       orientmap : the first eigen vector direction at each position
#       bstep: maximal length for a gap in a ridge
#       nbchain: maximum number of chains
#       ptile: relative threshold for the ridges 
#       maxchnlng: maximal chain length
#
#      output:
#      ------
#       ordered: ordered map: image containing the ridges
#                (displayed with different colors)
#       chain: 2D array containing the chained ridges, according
#		to the chain data structure:
#		chain[,1]: first point of the ridge
#		chain[,2]: length of the chain
#		chain[,3:(chain[,2]+2)]: values of the ridge
#	nbchain: number of chains.
#
#########################################################################
{
   d <- dim(pcaridge)
   sigsize <- d[1]
   nscale <- d[2]
   threshold <- range(pcaridge)[2] * ptile
   sz <- 2 * (maxchnlng);
   chain <- matrix(-1,nbchain,sz)
   orderedmap <- matrix(0,sigsize,nscale)	
   dim(chain) <- c(nbchain * sz, 1)
   dim(pcaridge) <- c(nscale * sigsize, 1)
   dim(orderedmap) <- c(nscale * sigsize, 1)
   dim(orientmap) <- c(nscale * sigsize, 1)
	
   
   z <- .C("Spca_family",
        as.double(pcaridge),
        as.integer(orientmap),
        orderedmap = as.single(orderedmap),
	chain = as.integer(chain),
	chainnb = as.integer(nbchain),
	as.integer(sigsize),
        as.integer(nscale),
        as.integer(bstep),
	as.single(threshold),
	as.integer(maxchnlng))
	
   orderedmap <- z$orderedmap
   chain <- z$chain

   dim(orderedmap) <- c(sigsize,nscale)   
   dim(chain) <- c(nbchain,sz)   
   nbchain <- z$chainnb
   chain <- chain +1
   chain[,1] <- chain[,1]-1

   list(ordered = orderedmap, chain = chain, nbchain=nbchain)
}   



pcamaxima <- function(beemap,orientmap)
{
   d <- dim(beemap)
   sigsize <- d[1]
   nscale <- d[2]  
   tfmaxima <- matrix(0,sigsize,nscale)
   dim(orientmap) <- c(nscale * sigsize, 1)
   dim(beemap) <- c(nscale * sigsize, 1)   
   dim(tfmaxima) <- c(nscale * sigsize, 1)   

   z <- .C("Stf_pcaridge",
        as.double(beemap),
        tfmaxima = as.double(tfmaxima),
	as.integer(sigsize),
        as.integer(nscale),
        as.integer(orientmap))
	
   tfmaxima <- z$tfmaxima
   dim(tfmaxima) <- c(sigsize,nscale)   
   tfmaxima
}


simplepcarec <- function(siginput, tfinput, beemap, orientmap, bstep = 5,
	ptile = .01, plot = 2)
#########################################################################
#     simplepcarec:
#     -------------
#      Simple reconstruction of a real valued signal from ridges found by 
#      pca climbers.
#
#      input:
#      ------
#      siginput: input signal
#      tfinput: continuous time-frequency transform (output of cwt or cgt)
#      beemap: output of pca climber algorithm
#      orientmap: pca dirction of climber
#      bstep: used for the chaining
#      ptile: 
#      plot: if set to 1, displays the signal, the components, and
#            signal and reconstruction one after another. If set to
#            2, displays the signal, the components, and the
#            reconstruction on the same page. Else, no plot.
#
#      output:
#      -------
#      rec: reconstructed signal
#      ordered: image of the ridges (with different colors)
#      comp: 2D array containing the signals reconstructed from ridges
#
#########################################################################
{
   tmp <- pcafamily(beemap,orientmap,ptile=ptile)
   chain <- tmp$chain
   nbchain <- tmp$nbchain
   rec <- numeric(length(siginput))

   if(plot != F){
     par(mfrow=c(2,1))
     plot.ts(siginput)
     title("Original signal")
     image(tmp$ordered)
     title("Chained Ridges")
   }

   npl(1)

   sigsize <- dim(tfinput)[1]
   nscale <- dim(tfinput)[2]
   rec <- numeric(sigsize)

   tmp3 <- matrix(0+0i,nbchain,sigsize)

   for (j in 1:nbchain){
      for(i in 1:chain[j,1]) {
        cat("The (b,a) at chain", j,"is (",chain[j,2*i+1],chain[j,2*i],")\n")
        tmp3[j,chain[j,2*i+1]] <- tmp3[j,chain[j,2*i+1]]+
                                  2*tfinput[chain[j,2*i+1],chain[j,2*i]]
      }
      rec <- rec + Re(tmp3[j,])
   }

   if(plot == 1){
      par(mfrow=c(2,1))
      par(cex=1.1)
      plot.ts(siginput)
      title("Original signal")
      plot.ts(rec)
      title("Reconstructed signal")
   }
   else if (plot == 2){
      par(mfrow=c(nbchain+2,1))
      par(mar=c(2,4,4,4))
      par(cex=1.1)
      par(err=-1)
      plot.ts(siginput)
      title("Original signal")

      for (j in 1:nbchain)
         plot.ts(Re(tmp3[j,]))
        
      plot.ts(rec)
      title("Reconstructed signal")
   }
   list(rec=rec, ordered=tmp$ordered,chain = tmp$chain, comp=tmp3)
}





#########################################################################
#               (c) Copyright  1997                             
#                          by                                   
#      Author: Rene Carmona, Bruno Torresani, Wen-Liang Hwang   
#                  Princeton University
#                  All right reserved                           
#########################################################################


#########################################################################
#
#	Functions to reconstruct a signal from the output of
#	 the pca climber algorithm.
#
#########################################################################




pcarec <- function(siginput,inputwt, beemap, orientmap, noct, nvoice, compr,
	maxchnlng=as.numeric(dim(beemap)[1])+10,minnbnodes = 2, 
	w0 = 2*pi, nbchain=100,bstep = 1,ptile =.01,para=5,plot=2,check=F)
#########################################################################
#     pcarec:
#     -------
#      Reconstruction of a real valued signal from ridges found by 
#      pca climbers on a wavelet transform.
#
#      input:
#      ------
#      siginput: input signal
#      inputwt: continuous wavelet transform (output of cwt)
#      beemap: output of pca climber algorithm
#      orientmap: principle direction of pca climber
#      noct: number of octaves (powers of 2)
#      nvoice: number of different scales per octave
#      compr: subsampling rate for the ridge
#      maxchnlng: maxlength of a chain of ridges (a,b)
#      bstep: used for the chaining
#      ptile: 
#      para:
#      plot: plot the original and reconstructed signal in display
#      check: check whether the reconstructed signal keeps the values
#             at ridges
#
#      output:
#      -------
#      rec: reconstructed signal
#      ordered: image of the ridges (with different colors)
#      comp: 2D array containing the signals reconstructed from ridges
#
#########################################################################
{

   tmp <- pcafamily(beemap,orientmap,maxchnlng=maxchnlng,bstep=bstep,nbchain=nbchain,ptile=ptile)
   chain <- tmp$chain
   nbchain <- tmp$nbchain
   ordered <- tmp$ordered
   sigsize <- length(siginput)
   rec <- numeric(sigsize)
   plnb <- 0

   if(plot != F){
     par(mfrow=c(2,1))
     plot.ts(siginput)
     title("Original signal")
     image(tmp$ordered)
     title("Chained Ridges")
   }

   sol <- matrix(0,nbchain,sigsize)

   totnbnodes <- 0
   idx <- numeric(nbchain)	
   p <- 0

   if(check==T) {
      inputskel <- matrix(0+0i,nbchain,sigsize)
      solskel <- matrix(0+0i, nbchain, sigsize)
   }

   for (j in 1:nbchain){
      
      sol[j,] <- 0
      nbnode <- chain[j,1]
      bnode <- numeric(nbnode)
      anode <- numeric(nbnode)
      for(k in 1:nbnode) {
        anode[k] <- chain[j,2*k]
        bnode[k] <- chain[j,2*k+1]
      }
      cat("Chain number",j,"\n")

      tmp2 <- pcaregrec(siginput[min(bnode):max(bnode)],
          inputwt[min(bnode):max(bnode),],
          anode,bnode,compr,noct,nvoice, 
          w0 = w0, para = para,minnbnodes = minnbnodes, check=check);

      if(is.list(tmp2)==T) {
        totnbnodes <- totnbnodes + tmp2$nbnodes
        if((sigsize-min(bnode)) > (length(tmp2$sol)-tmp2$bstart)){
          np <- length(tmp2$sol) - tmp2$bstart
          sol[j,min(bnode):(np+min(bnode))]<-tmp2$sol[tmp2$bstart:length(tmp2$sol)]
        }
        else {
	  np <- sigsize - min(bnode)
          sol[j,min(bnode):sigsize]<-tmp2$sol[tmp2$bstart:(np+tmp2$bstart)]
        }
        if(min(bnode) < tmp2$bstart) {
          np <- min(bnode)-1
          sol[j,1:min(bnode)]<-tmp2$sol[(tmp2$bstart-np):(tmp2$bstart)]
        }          
        else {
          np <- tmp2$bstart-1
          sol[j,(min(bnode)-np):min(bnode)]<-tmp2$sol[1:(tmp2$bstart)]
        }

        if(check==T) {
           bridge <- tmp2$bnode
           aridge <- tmp2$anode
           wtsol <- cwt(sol[j,],noct,nvoice)
           for(k in 1:length(bridge)) solskel[j,k] <- wtsol[bridge[k],aridge[k]]
           for(k in 1:length(bridge)) inputskel[j,k] <- inputwt[bridge[k],aridge[k]]
        }
        rec <- rec+sol[j,]
      }
      plnb <- plnb + 1
      p <- p + 1
      idx[p] <- j
   }

   if(plot == 1){
      par(mfrow=c(2,1))
      par(cex=1.1)
      plot.ts(siginput)
      title("Original signal")
      plot.ts(Re(rec))
	
      title("Reconstructed signal")
   }
   else if (plot == 2){
      par(mfrow=c(plnb+2,1))
      par(mar=c(2,4,4,4))
      par(cex=1.1)
      par(err=-1)
      plot.ts(siginput)
      title("Original signal")

      for (j in 1:p)
        plot.ts(sol[idx[j],]);
      plot.ts(Re(rec))	
      title("Reconstructed signal")
   }

   cat("Total number of ridge samples used: ",totnbnodes,"\n")

   par(mfrow=c(1,1))
   if(check==T) list(rec=rec,ordered=ordered,chain=chain,
                     comp=tmp,inputskel=inputskel,solskel=solskel,lam=tmp2$lam)
   else list(rec=rec, ordered=ordered,chain=chain,comp=tmp)
}

PcaRidgeSampling <- function(anode, bnode, compr)
#########################################################################
#    PcaRidgeSampling:
#    ----------------
#      Given a ridge (a,b)(for the Gabor transform), returns a 
#        (regularly) subsampled version of length nbnodes.
#
#      Input:
#      ------
#      anode: scale coordinates of the (unsampled) ridge
#      bnode: time coordinates of the (unsampled) ridge
#      compr: compression ratio to obtain from the ridge
#
#      Output:
#      -------
#      bnode: time coordinates of the ridge samples (from 1 to nbnodes
#             step by compr)
#      anode: scale coordinates of the ridge samples
#      nbnode: number of node sampled
#
#########################################################################
{
   # number of node at the ridge
   nbnode <- length(anode)

   # compr to be integer
   compr <- as.integer(compr)
   if(compr < 1) compr <- 1

   # sample rate is compr at the ridge
   k <- 1
   count <- 1
   while(k < nbnode) {
     anode[count] <- anode[k]
     bnode[count] <- bnode[k]
     k <- k + compr 
     count <- count + 1
   }

   anode[count] <- anode[nbnode]
   bnode[count] <- bnode[nbnode]
   
   aridge <- numeric(count)
   bridge <- numeric(count)
   aridge <- anode[1:count]
   bridge <- bnode[1:count]
   
   list(anode = aridge, bnode = bridge, nbnodes = count)
}


pcaregrec <- function(siginput,cwtinput,anode,bnode,compr,noct,nvoice,
	w0 = 2*pi, plot = F, para = 5, minnbnodes = 2, check=F)
#########################################################################
#     pcaregrec:
#     ---------
#      Reconstruction of a real valued signal from a (continuous) ridge
#      (uses a regular sampling of the ridge)
#
#      input:
#      -------
#      siginput: input signal
#      cwtinput: Continuous wavelet transform (output of cwt)
#      anode: the scale coordinates of the (unsampled) ridge
#      bnode: the time coordinates of the (unsampled) ridge
#      compr: subsampling rate for the wavelet coefficients (at scale 1)
#      noct: number of octaves (powers of 2)
#      nvoice: number of different scales per octave
#      epsilon: coeff of the Q2 term in reconstruction kernel
#      w0: central frequency of Morlet wavelet
#      fast: if set to TRUE, the kernel is computed using Riemann
#           sums instead of Romberg's quadrature
#      plot: if set to TRUE, displays original and reconstructed signals
#      para: constant for extending the support of reconstructed signal
#            outside the ridge. 
#      check: if set to TRUE, computes the wavelet transform of the
#            reconstructed signal
#      minnbnodes: minimum number of nodes for the reconstruction
#      real: if set to TRUE, only uses constraints on the real part
#            of the transform for the reconstruction.
#
#      output:
#      -------
#      sol: reconstruction from a ridge
#      A: <wavelets,dualwavelets> matrix 
#      lam: coefficients of dual wavelets in reconstructed signal.
#      dualwave: array containing the dual wavelets.
#      Q2: second part of the reconstruction kernel
#      nbnodes: number of nodes used for the reconstruction.
#
##########################################################################
{

#  Generate Sampled ridge
#   tmp <- wRidgeSampling(anode,compr,nvoice)
#   bnode <- tmp$node
#   anode <- tmp$phinode
#   nbnodes <- tmp$nbnodes

   tmp <- PcaRidgeSampling(anode, bnode, compr)
   bnode <- tmp$bnode
   anode <- tmp$anode
   nbnodes <- tmp$nbnodes

   cat("Sampled nodes (b,a): \n")
   for(j in 1:nbnodes) cat("(",bnode[j],anode[j],")")
   cat("\n")

   if(nbnodes < minnbnodes){
      cat(" Chain too small\n")
      NULL
   }
   else {

    a.min <- 2 * 2^(min(anode)/nvoice)
    a.max <- 2 * 2^(max(anode)/nvoice)
    b.min <- min(bnode)
    b.max <- max(bnode)
    b.max <- (b.max-b.min+1) + round(para * a.max)
    b.min <- (b.min-b.min+1) - round(para * a.max)

# location 2-b.min is the place where the min(bnode) is.
# The position at 2-b.min of returned signal is aligned to 
# the position at min(b.node) for the reconstruction

    bnode <- bnode-min(bnode)+2-b.min
    b.inc <- 1
    np <- as.integer((b.max - b.min)/b.inc) +1
    cat("(size:",np,",",nbnodes,"sampled nodes):\n")

    # Generating the Q2 term in reconstruction kernel
    Q2 <- 0
	
    # Generating the Q1 term in reconstruction kernel
    one <- numeric(np)
    one[] <- 1
    Qinv <-  1/one

    tmp2 <- pcaridrec(cwtinput,bnode,anode,noct,nvoice,
              Qinv,np,w0=w0,check=check)

    if(plot == T){
       par(mfrow=c(2,1))
       plot.ts(Re(siginput))
       title("Original signal")
       plot.ts(Re(tmp2$sol))
       title("Reconstructed signal")
    }
    lam <- tmp2$lam
#    if(check==T) {
#      N <- length(lam)/2
#      mlam <- numeric(N)
#      for(j in 1:N)
#         mlam[j] <- Mod(lam[j] + i * lam[N + j])
#      npl(2)
#      plot.ts(lam,xlab="Number",ylab="Lambda Value")
#      title("Lambda Profile")
#      plot.ts(sort(mlam),xlab="Number",ylab="Mod of Lambda")
#    }
    list(sol = tmp2$sol,A = tmp2$A,
          lam = tmp2$lam, dualwave = tmp2$dualwave, 
          morvelets = tmp2$morvelets,pQ2 = Q2, nbnodes=nbnodes,
          bstart=2-b.min, anode=tmp$anode,bnode=tmp$bnode)
   }
	
}


pcaridrec <- function(cwtinput,bnode,anode,noct,voice,
	Qinv, np, w0 = 2*pi, check = F)
#########################################################################
#     pcaridrec:
#     ---------
#      Reconstruction of a real valued signal from a ridge, given
#       the kernel of the bilinear form.
#
#      input:
#      ------
#      cwtinput: Continuous wavelet transform (output of cwt)
#      bnode: time coordinates of the ridge samples
#      anode: scale coordinates of the ridge samples
#      noct: number of octaves (powers of 2)
#      nvoice: number of different scales per octave
#      Qinv: inverse of the reconstruction kernel
#      epsilon: coefficient of the Q2 term in the reconstruction kernel
#      np: number of samples of the reconstructed signal
#      w0: central frequency of Morlet wavelet
#      check: if set to TRUE, computes the wavelet transform of the
#            reconstructed signal
#
#      output:
#      -------
#      sol: reconstruction from a ridge
#      A: <wavelets,dualwavelets> matrix 
#      lam: coefficients of dual wavelets in reconstructed signal.
#      morwave: array of morlet wavelets located on the ridge samples
#      dualwave: array containing the dual wavelets.
#      solskel: wavelet transform of sol, restricted to the ridge
#      inputskel: wavelet transform of signal, restricted to the ridge
#
#########################################################################
{
   # number of sampled nodes in a ridge
   N <- length(bnode)

   aridge <- anode
   bridge <- bnode

   morvelets <- pcamorwave(bridge,aridge,nvoice,np,N)

   cat("morvelets;  ")

   sk <- pcazeroskeleton(cwtinput,Qinv,morvelets,bridge,aridge,N)    

   cat("skeleton.\n")

   solskel <- 0 #not needed if check not done
   inputskel <- 0 #not needed if check not done

   list(sol=sk$sol,A=sk$A,lam=sk$lam,dualwave=sk$dualwave,morvelets=morvelets,
	solskel=solskel,inputskel = inputskel)
}


pcazeroskeleton <- function(cwtinput,Qinv,morvelets,bridge,aridge,N)
#########################################################################
#     pcazeroskeleton:
#     ----------------
#      Computation of the reconstructed signal from the ridge when
#       the epsilon parameter is set to 0.
#
#      input:
#      -------
#      cwtinput: Continuous wavelet transform (output of cwt)
#      Qinv: 1D array (warning: different from skeleton) containing
#            the diagonal of the inverse of reconstruction kernel.
#      morvelets: array of morlet wavelets located on the ridge samples
#      bridge: time coordinates of the ridge samples
#      aridge: scale coordinates of the ridge samples
#      N: number of complex constraints
#
#      output:
#      -------
#      sol: reconstruction from a ridge
#      A: <wavelets,dualwavelets> matrix 
#      lam: coefficients of dual wavelets in reconstructed signal.
#      dualwave: array containing the dual wavelets.
#
##########################################################################
{
   tmp1 <- dim(morvelets)[1]
   constraints2 <- dim(morvelets)[2]
   dualwave <- matrix(0,tmp1,constraints2)
   for (j in 1:tmp1) {
      dualwave[j,] <- morvelets[j,]*Qinv[j]
   }
   A <- t(morvelets) %*% dualwave

   rskel <- numeric(2*N)

   if(is.vector(cwtinput) == T) {
     for(j in 1:N) {
      rskel[j] <- Re(cwtinput[aridge[j]])
      rskel[N+j] <- -Im(cwtinput[aridge[j]])
     }
   }
   else {
     for(j in 1:N) {
      rskel[j] <- Re(cwtinput[bridge[j]-bridge[1]+1,aridge[j]])
      rskel[N+j] <- -Im(cwtinput[bridge[j]-bridge[1]+1,aridge[j]])
     }
   }

   B <- SVD(A)
   d <- B$d
   Invd <- numeric(length(d))
   for(j in 1:(2*N))
     if(d[j] < 1.0e-6) 
       Invd[j] <- 0
    else Invd[j] <- 1/d[j]
   lam <- B$v %*% diag(Invd) %*% t(B$u) %*% rskel

   sol <- dualwave%*%lam
   list(lam=lam,sol=sol,dualwave=dualwave,A=A)
}


pcamorwave <- function(bridge, aridge, nvoice, np, N, w0 = 2*pi)
#########################################################################
#     pcamorwave: ( the same function of morwave )
#     -----------
#      Generation of the wavelets located on the ridge.
#
#      input:
#      ------
#      bridge: time coordinates of the ridge samples
#      aridge: scale coordinates of the ridge samples
#      nvoice: number of different scales per octave
#      np: size of the reconstruction kernel
#      N: number of complex constraints
#      w0: central frequency of Morlet wavelet
#
#      output:
#      -------
#      morvelets: array of morlet wavelets located on the ridge samples
#
#########################################################################
{

   morvelets <- matrix(0,np,2*N)

   aridge <-  2 * 2^((aridge - 1)/nvoice)
   tmp <- vecmorlet(np,N,bridge,aridge,w0 = w0)
   dim(tmp) <- c(np,N)
   morvelets[,1:N] <- Re(tmp[,1:N])
   morvelets[,(N+1):(2*N)] <- Im(tmp[,1:N])

#  Another way of generating morvelets
#   dirac <- 1:np
#   dirac[] <- 0
#   for(k in 1:N)  {
#      dirac[bridge[k]] <- 1.0;
#      scale <- 2 * 2^((aridge[k]-1)/nvoice);
#      cwtdirac <- vwt(dirac,scale);
#      morvelets[,k] <- Re(cwtdirac);
#      morvelets[,k+N] <- Im(cwtdirac);
#      dirac[] <- 0
#    }

#  Still another way of generating morvelets
#   for(k in 1:N)  {
#      scale <- 2 * 2^((aridge[k]-1)/nvoice);
#      tmp <- morlet(np,bridge[k],scale)/sqrt(2*pi);
#      morvelets[,k] <- Re(tmp);
#      morvelets[,k+N] <- Im(tmp);
#    }

    morvelets
}



#########################################################################
#      $Log: Ridge_Annealing.S,v $
#########################################################################
#
#               (c) Copyright  1997
#                          by                                   
#      Author: Rene Carmona, Bruno Torresani, Wen-Liang Hwang   
#                  University of California, Irvine             
#                  All right reserved                           
#########################################################################





corona <- function(tfrep, guess, tfspec = numeric(dim(tfrep)[2]),
	subrate = 1, temprate = 3, mu = 1, lambda = 2*mu,
	iteration = 1000000, seed = -7, stagnant = 20000, costsub = 1,plot=T)
#########################################################################
#       corona:
#       -------
#        ridge extraction using simulated annealing
# 	  compute the continuous wavelet transform ridge from the cwt
#	   modulus, using simulated annealing
#
#       Input:
#       ------
# 	 tfrep: wavelet or Gabor transform.
#	 guess: initial guess for ridge function (1D array).
#        tfspec: estimate for the contribution of the noise to the
#               transform modulus (1D array)
#        subrate: subsampling rate for the transform modulus.
#	 temprate: constant (numerator) in the temperature schedule.
#        lambda: coefficient in front of phi' in the cost function
#        mu: coefficient in front of phi'' in the cost function
#        iteration: maximal number of iterations
#        seed: initialization for random numbers generator
#        stagnant: maximum number of steps without move (for the
#                  stopping criterion)
#        costsub: subsampling of the cost function in cost
#               costsub = 1 means that the whole cost function
#               is returned
#        plot:	 when set(default), some results will be shown on 
#	        the display	
#
#       Output:
#       -------
#        ridge: 1D array (of length sigsize) containing the ridge
#        cost: 1D array containing the cost function
#
#########################################################################
{

  sigsize <- dim(tfrep)[1]
  nscale <- dim(tfrep)[2]
  blocksize <- 1	
  costsize <- as.integer(iteration/blocksize) + 1	
  cost <- 1:costsize

  modulus <- Mod(tfrep)
  tfspectrum <- tfspec

  cost[] <- 0
  phi <- as.integer(guess)
  count <- 0	
  dim(phi) <- c(sigsize,1)	
  dim(modulus) <- c(sigsize * nscale, 1)
   
  smodsize <- as.integer(sigsize/subrate)

  if((sigsize/subrate - as.integer(sigsize/subrate)) > 0)
    smodsize <- as.integer(sigsize/subrate) + 1

  smodulus <- matrix(0,smodsize,nscale)
  dim(smodulus) <- c(smodsize * nscale, 1)


  z <- .C("Smodulus_smoothing",
        as.double(modulus),
        smodulus = as.double(smodulus),
	as.integer(sigsize),
        smodsize = as.integer(smodsize),            
        as.integer(nscale),
        as.integer(subrate))

  smodulus <- z$smodulus
  smodsize <- z$smodsize
  smodulus <- smodulus * smodulus
  dim(smodulus) <- c(smodsize, nscale)
  for (k in 1:nscale)
    smodulus[,k] <- smodulus[,k] - tfspectrum[k]

  dim(smodulus) <- c(smodsize * nscale, 1)
  cat("dimension of smodulus",dim(smodulus),"\n")

   z <- .C("Sridge_annealing",
	cost = as.single(cost),
	as.double(smodulus),
	phi = as.single(phi),
	as.single(lambda),
	as.single(mu),
	as.single(temprate),
	as.integer(sigsize),
	as.integer(nscale),
	as.integer(iteration),
	as.integer(stagnant),
	as.integer(seed),
	nb = as.integer(count),
	as.integer(subrate),
	as.integer(costsub),
        as.integer(smodsize))

  count <- z$nb
  cat("Number of iterations:",(count-1)*costsub,"(stagnant=",stagnant,")\n")
  if(plot == T) {
   image(Mod(tfrep))
   lines(z$phi+1)
  }	
  list(ridge = z$phi+1, cost = z$cost[1:count])
}



  
coronoid <- function(tfrep, guess, tfspec = numeric(dim(tfrep)[2]),
	subrate = 1, temprate = 3, mu = 1, lambda = 2*mu,
	iteration = 1000000, seed = -7, stagnant = 20000,
	costsub = 1,plot=T)
#########################################################################
#       coronoid:   
#       ----------
#        ridge extraction using modified simulated annealing
#	   Modification of Sridge_annealing, the cost function is
#	   replaced with another one, in which the smoothness penalty
#	   is proportional to the transform modulus.
#
#       Input:
#       ------
# 	 tfrep: wavelet or Gabor transform.
#	 guess: initial guess for ridge function (1D array).
#        tfspec: estimate for the contribution of the noise to the
#               transform modulus (1D array)
#        subrate: subsampling rate for the transform modulus.
#	 temprate: constant (numerator) in the temperature schedule.
#        lambda: coefficient in front of phi' in the cost function
#        mu: coefficient in front of phi'' in the cost function
#        iteration: maximal number of iterations
#        seed: initialization for random numbers generator
#        stagnant: maximum number of steps without move (for the
#                  stopping criterion)
#        costsub: subsampling of the cost function in cost
#               costsub=1 means that the whole cost function
#               is returned
#        plot:	 when set(default), some results will be shown on 
#	        the display	
#
#       Output:
#       -------
#        ridge: 1D array (of length sigsize) containing the ridge
#        cost: 1D array containing the cost function
#
#########################################################################
{

  sigsize <- dim(tfrep)[1]
  nscale <- dim(tfrep)[2]
  costsize <- as.integer(iteration/costsub) + 1	
  cost <- 1:costsize

  modulus <- Mod(tfrep)
  tfspectrum <- tfspec

  cost[] <- 0
  phi <- as.integer(guess)
  count <- 0	
  dim(phi) <- c(sigsize,1)	
  dim(modulus) <- c(sigsize * nscale, 1)

  smodsize <- as.integer(sigsize/subrate)

  if((sigsize/subrate - as.integer(sigsize/subrate)) > 0)
    smodsize <- as.integer(sigsize/subrate) + 1

  smodulus <- matrix(0,smodsize,nscale)
  dim(smodulus) <- c(smodsize * nscale, 1)

  z <- .C("Smodulus_smoothing",
        as.double(modulus),
        smodulus = as.double(smodulus),
	as.integer(sigsize),
        smodsize = as.integer(smodsize),            
        as.integer(nscale),
        as.integer(subrate))

  smodulus <- z$smodulus
  smodsize <- z$smodsize
  smodulus <- smodulus * smodulus
  dim(smodulus) <- c(smodsize, nscale)
  for (k in 1:nscale)
    smodulus[,k] <- smodulus[,k] - tfspectrum[k]
  dim(smodulus) <- c(smodsize * nscale, 1)


   z <- .C("Sridge_coronoid",
        cost = as.single(cost),
        as.double(smodulus),
        phi = as.single(phi),
        as.single(lambda),
	as.single(mu),
        as.single(temprate),
	as.integer(sigsize),
        as.integer(nscale),
        as.integer(iteration),
	as.integer(stagnant),
	as.integer(seed),
	nb = as.integer(count),
        as.integer(subrate),
	as.integer(costsub),
        as.integer(smodsize))

  count <- z$nb
  cat("Number of iterations:",(count-1)*costsub,"(stagnant=",stagnant,")\n")
  if(plot == T) {	
    image(Mod(tfrep))
    lines(z$phi+1)
  }
  list(ridge = z$phi+1, cost = z$cost[1:count])
}
  























#########################################################################
#      $Log: Ridge_Icm.S,v $
#########################################################################
#
#               (c) Copyright  1997                             
#                          by                                   
#      Author: Rene Carmona, Bruno Torresani, Wen-Liang Hwang   
#                  Princeton University 
#                  All right reserved                           
#########################################################################




icm <- function(modulus, guess, tfspec = numeric(dim(modulus)[2]),
	subrate = 1, mu = 1, lambda = 2*mu, iteration = 100)
#########################################################################
#       icm:   
#       ----
#        ridge extraction using icm algorithm
#
#       Input:
#       ------
# 	 modulus: modulus of the wavelet or Gabor transform.
#	 guess: initial guess for ridge function (1D array).
#        tfspec: estimate for the contribution of the noise to the
#               transform square modulus (1D array).
#        subrate: subsampling rate for the transform square modulus.
#        lambda: coefficient in front of phi'' in the cost function
#        mu: coefficient in front of phi' in the cost function
#        iteration: maximal number of iterations
#
#       Output:
#       -------
#        ridge: 1D array (of length sigsize) containing the ridge
#        cost: 1D array containing the cost function
#
#########################################################################
{

  sigsize <- dim(modulus)[1]
  nscale <- dim(modulus)[2]
  costsize <- iteration	
  cost <- 1:costsize
  tfspectrum <- tfspec

  cost[] <- 0
  phi <- as.integer(guess)
  count <- 0	
  dim(phi) <- c(sigsize,1)	
  dim(modulus) <- c(sigsize * nscale, 1)

  smodsize <- as.integer(sigsize/subrate)

  if((sigsize/subrate - as.integer(sigsize/subrate)) > 0)
    smodsize <- as.integer(sigsize/subrate) + 1

  smodulus <- matrix(0,smodsize,nscale)
  dim(smodulus) <- c(smodsize * nscale, 1)

  if (subrate != 1){
  z <- .C("Smodulus_smoothing",
        as.double(modulus),
        smodulus = as.double(smodulus),
	as.integer(sigsize),
        smodsize = as.integer(smodsize),            
        as.integer(nscale),
        as.integer(subrate))
  
  smodulus <- z$smodulus
  smodsize <- z$smodsize
  }
  else smodulus <- modulus
  smodulus <- smodulus * smodulus
  dim(smodulus) <- c(smodsize, nscale)
  for (k in 1:nscale)
    smodulus[,k] <- smodulus[,k] - tfspectrum[k]
  dim(smodulus) <- c(smodsize * nscale, 1)

  z <- .C("Sridge_icm",
        cost = as.single(cost),
        as.double(smodulus),
        phi = as.single(phi),
        as.single(lambda),
	as.single(mu),
	as.integer(sigsize),
        as.integer(nscale),
        as.integer(iteration),
	nb = as.integer(count),
        as.integer(subrate),
        as.integer(smodsize))

  count <- z$nb
  cat("Number of iterations:",count,"\n")
  list(ridge = z$phi+1, cost = z$cost[1:count])
}
  




















#########################################################################
#      $Log: Ridge_Recons.S,v $
# Revision 1.2  1995/04/05  18:56:55  bruno
# *** empty log message ***
#
# Revision 1.1  1995/04/02  01:04:16  bruno
# Initial revision
#
#
#               (c) Copyright  1997                             
#                          by     
#      Author: Rene Carmona, Bruno Torresani, Wen-Liang Hwang   
#                  Princeton University 
#                  All right reserved                           
#########################################################################





irregrec <- function(siginput, cwtinput,phi,compr,noct,nvoice,
	epsilon = 0.5, w0 = 2*pi, prob = 0.8, fast = F, plot = F, para = 0,
        hflag = F, check = F, minnbnodes = 2, real = F)
#########################################################################
#     irregrec:
#     ---------
#      Reconstruction of a real valued signal from a (continuous) ridge
#      (uses a irregular sampling of the ridge)
#
#      input:
#      -------
#      siginput: input signal
#      cwtinput: Continuous wavelet transform (output of cwt)
#      phi: (unsampled) ridge
#      compr: subsampling rate for the wavelet coefficients (at scale 1)
#      noct: number of octaves (powers of 2)
#      nvoice: number of different scales per octave
#      epsilon: coeff of the Q2 term in reconstruction kernel
#      w0: central frequency of Morlet wavelet
#      fast: if set to TRUE, the kernel is computed using Riemann
#           sums instead of Romberg's quadrature
#      plot: if set to TRUE, displays original and reconstructed signals
#      para: constant for extending the support of reconstructed signal
#            outside the ridge. 
#      check: if set to TRUE, computes the wavelet transform of the
#            reconstructed signal
#      minnbnodes: minimum number of nodes for the reconstruction
#      real: if set to TRUE, only uses constraints on the real part
#            of the transform for the reconstruction.
#      prob: only keep the lambda's values greater than this percential after regular sampling
#
#      output:
#      -------
#      sol: reconstruction from a ridge
#      A: <wavelets,dualwavelets> matrix 
#      lam: coefficients of dual wavelets in reconstructed signal.
#      dualwave: array containing the dual wavelets.
#      solskel: wavelet transform of sol, restricted to the ridge
#      inputskel: wavelet transform of signal, restricted to the ridge
#      Q2: second part of the reconstruction kernel
#      nbnodes: number of nodes used for the reconstruction.
#
##########################################################################
{

#
#  Generate (regularly) sampled ridge
#


   tmp <- wRidgeSampling(phi,compr,nvoice)

   node <- tmp$node
cat("node at = ", node, "\n")
   phinode <- tmp$phinode
   nbnodes <- tmp$nbnodes
   phinode <- as.integer(phinode)	

   if(nbnodes < minnbnodes){
      cat(" Chain too small\n")
      NULL 
   }

   phi.x.min <- 2 * 2^(phinode[1]/nvoice)
   phi.x.max <- 2 * 2^(phinode[length(node)]/nvoice)

   x.min <- node[1]
   x.max <- node[length(node)]
   x.max <- x.max + round(para * phi.x.max)
   x.min <- x.min - round(para * phi.x.min)

   node <- node + 1 - x.min
   x.inc <- 1
   np <- as.integer((x.max - x.min)/x.inc) +1

    # Generating the Q2 term in reconstruction kernel
    if(epsilon == 0)
       Q2 <- 0
    else {
       if (fast == F)
         Q2 <- rkernel(node, phinode, nvoice, x.min = x.min,x.max = x.max,w0 = w0)
       else
	 Q2 <- fastkernel(node, phinode, nvoice, x.min = x.min,x.max = x.max,w0 = w0)
     }
    cat(" kernel; ")

    # Generating the Q1 term in reconstruction kernel
    if (hflag == T){
       one <- numeric(np)
       one[] <- 1
     }
    else{
       one <- numeric(np)
       one[] <- 1
    }

    if (epsilon !=0 ){
      Q <-  epsilon * Q2
      for(j in 1:np)
         Q[j,j] <- Q[j,j] + one[j]
      Qinv <- solve(Q)
    }
    else{
       Qinv <-  1/one
     }

    tmp2 <- ridrec(cwtinput, node, phinode, noct, nvoice, Qinv,
                            epsilon, np, w0 = w0, check = check, real = real)

# 
# Now perform the irregular sampling 
#
# C and Splus ...

   tmp <- RidgeIrregSampling(phi, node, tmp2$lam, prob)
   node <- tmp$node

  cat("node at  ", node,"\n")

   phinode <- tmp$phinode

  cat("phinode at  ", phinode,"\n")
   nbnodes <- tmp$nbnodes
   phinode <- as.integer(phinode)	

   if(nbnodes < minnbnodes){
      cat(" Chain too small\n")
      NULL 
   }

   phi.x.min <- 2 * 2^(phinode[1]/nvoice)
   phi.x.max <- 2 * 2^(phinode[length(node)]/nvoice)

   x.min <- node[1]
   x.max <- node[length(node)]

   x.max <- x.max + round(para * phi.x.max)
   x.min <- x.min - round(para * phi.x.min)

   node <- node + 1 - x.min
   x.inc <- 1


   np <- as.integer((x.max - x.min)/x.inc) +1
   cat("(size:",np,",",nbnodes,"nodes):")

    # Generating the Q2 term in reconstruction kernel
    if(epsilon == 0)
       Q2 <- 0
    else {
       if (fast == F)
         Q2 <- rkernel(node, phinode, nvoice, x.min = x.min,x.max = x.max,w0 = w0)
       else
	 Q2 <- fastkernel(node, phinode, nvoice, x.min = x.min,x.max = x.max,w0 = w0)
     }
    cat(" kernel; ")

    # Generating the Q1 term in reconstruction kernel
    if (hflag == T){
       one <- numeric(np)
       one[] <- 1
     }
    else{
       one <- numeric(np)
       one[] <- 1
    }

    if (epsilon !=0 ){
      Q <-  epsilon * Q2
      for(j in 1:np)
        Q[j,j] <- Q[j,j] + one[j]
      Qinv <- solve(Q)
    }
    else{
      Qinv <-  1/one
    }


    tmp2 <- ridrec(cwtinput, node, phinode, noct, nvoice, Qinv,
                   epsilon, np, w0 = w0, check = check, real = real)


    if(plot == T){
       par(mfrow=c(2,1))
       plot.ts(Re(siginput))
       title("Original signal")
       plot.ts(Re(tmp2$sol))
       title("Reconstructed signal")
    }


    list(sol = tmp2$sol, A = tmp2$A, lam = tmp2$lam, dualwave = tmp2$dualwave,
	  morvelets = tmp2$morvelets, solskel = tmp2$solskel,
          inputskel = tmp2$inputskel, Q2 = Q2, nbnodes = nbnodes)
}



#########################################################################
#      $Log: Ridge_Recons.S,v $
#########################################################################
#
#               (c) Copyright  1997
#                          by     
#      Author: Rene Carmona, Bruno Torresani, Wen-Liang Hwang   
#                  Princeton University 
#                  All right reserved                           
#########################################################################





regrec <- function(siginput, cwtinput,phi,compr,noct,nvoice,
	epsilon = 0, w0 = 2*pi, fast = F, plot = F, para = 0,
        hflag = F, check = F, minnbnodes = 2, real = F)
#########################################################################
#     regrec:
#     -------
#      Reconstruction of a real valued signal from a (continuous) ridge
#      (uses a regular sampling of the ridge)
#
#      input:
#      -------
#      siginput: input signal
#      cwtinput: Continuous wavelet transform (output of cwt)
#      phi: (unsampled) ridge
#      compr: subsampling rate for the wavelet coefficients (at scale 1)
#      noct: number of octaves (powers of 2)
#      nvoice: number of different scales per octave
#      epsilon: coeff of the Q2 term in reconstruction kernel
#      w0: central frequency of Morlet wavelet
#      fast: if set to TRUE, the kernel is computed using Riemann
#           sums instead of Romberg's quadrature
#      plot: if set to TRUE, displays original and reconstructed signals
#      para: constant for extending the support of reconstructed signal
#            outside the ridge. 
#      check: if set to TRUE, computes the wavelet transform of the
#            reconstructed signal
#      minnbnodes: minimum number of nodes for the reconstruction
#      real: if set to TRUE, only uses constraints on the real part
#            of the transform for the reconstruction.
#
#      output:
#      -------
#      sol: reconstruction from a ridge
#      A: <wavelets,dualwavelets> matrix 
#      lam: coefficients of dual wavelets in reconstructed signal.
#      dualwave: array containing the dual wavelets.
#      morvelets: array containing the wavelets on sampled ridge.
#      solskel: wavelet transform of sol, restricted to the ridge
#      inputskel: wavelet transform of signal, restricted to the ridge
#      Q2: second part of the reconstruction kernel
#      nbnodes: number of nodes used for the reconstruction.
#
##########################################################################
{
#  Generate Sampled ridge

   tmp <- wRidgeSampling(phi,compr,nvoice)

   node <- tmp$node
   phinode <- tmp$phinode
   nbnodes <- tmp$nbnodes
   phinode <- as.integer(phinode)	

   if(nbnodes < minnbnodes){
      cat(" Chain too small\n")
      NULL
   }
   else {
   phi.x.min <- 2 * 2^(phinode[1]/nvoice)
   phi.x.max <- 2 * 2^(phinode[length(node)]/nvoice)

   x.min <- node[1]
   x.max <- node[length(node)]
   x.max <- x.max + round(para * phi.x.max)
   x.min <- x.min - round(para * phi.x.min)

   node <- node + 1 - x.min
   x.inc <- 1
   np <- as.integer((x.max - x.min)/x.inc) +1
   cat("(size:",np,",",nbnodes,"nodes):")

   fast <- F

    # Generating the Q2 term in reconstruction kernel
    if(epsilon == 0)
       Q2 <- 0
    else {
       if (fast == F)
         Q2 <- rkernel(node, phinode, nvoice, x.min = x.min,x.max = x.max,w0 = w0)
       else
	 Q2 <- fastkernel(node, phinode, nvoice, x.min = x.min,x.max = x.max,w0 = w0)
    }
	
    # Generating the Q1 term in reconstruction kernel
    if (hflag == T){
       one <- numeric(np)
       one[] <- 1
     }
    else{
       one <- numeric(np)
       one[] <- 1
    }

     if (epsilon !=0 ){
        Q <-  epsilon * Q2
        for(j in 1:np)
            Q[j,j] <- Q[j,j] + one[j]
           Qinv <- solve(Q)
     }
     else{
        Qinv <-  1/one
      }

   tmp2 <- ridrec(cwtinput,node,phinode,noct,nvoice,
             Qinv,epsilon,np, w0=w0, check=check,real=real)

    if(plot == T){
       par(mfrow=c(2,1))
       plot.ts(Re(siginput))
       title("Original signal")
       plot.ts(Re(tmp2$sol))
       title("Reconstructed signal")
    }
    npl(2)
    lam <- tmp2$lam
    if(plot==T) {
       plot.ts(lam,xlab="Number",ylab="Lambda Value")
       title("Lambda Profile")
     }
    N <- length(lam)/2
    mlam <- numeric(N)
    for(j in 1:N)
       mlam[j] <- Mod(lam[j] + lam[N + j]*(1i))
    if(plot==T) plot.ts(sort(mlam))

    list(sol = tmp2$sol, A = tmp2$A, lam = tmp2$lam, dualwave = tmp2$dualwave,
	  morvelets = tmp2$morvelets, solskel = tmp2$solskel,
          inputskel = tmp2$inputskel, Q2 = Q2, nbnodes = nbnodes)
	}
	
}


ridrec <- function(cwtinput, node, phinode, noct, nvoice,
	Qinv, epsilon, np, w0 = 2*pi, check = F, real = F)
#########################################################################
#     ridrec:
#     ------
#      Reconstruction of a real valued signal from a ridge, given
#       the kernel of the bilinear form.
#
#      input:
#      ------
#      cwtinput: Continuous wavelet transform (output of cwt)
#      node: time coordinates of the ridge samples
#      phinode: scale coordinates of the ridge samples
#      noct: number of octaves (powers of 2)
#      nvoice: number of different scales per octave
#      Qinv: inverse of the reconstruction kernel
#      epsilon: coefficient of the Q2 term in the reconstruction kernel
#      np: number of samples of the reconstructed signal
#      w0: central frequency of Morlet wavelet
#      check: if set to TRUE, computes the wavelet transform of the
#            reconstructed signal
#      real: if set to TRUE, only uses constraints on the real part
#            of the transform for the reconstruction.
#
#      output:
#      -------
#      sol: reconstruction from a ridge
#      A: <wavelets,dualwavelets> matrix 
#      lam: coefficients of dual wavelets in reconstructed signal.
#      dualwave: array containing the dual wavelets.
#      morvelets: array of morlet wavelets located on the ridge samples
#      solskel: wavelet transform of sol, restricted to the ridge
#      inputskel: wavelet transform of signal, restricted to the ridge
#
#########################################################################
{
   N <- length(node)

   aridge <- phinode
   bridge <- node

   if (real == T)
      morvelets <- morwave2(bridge,aridge,nvoice,np,N)
   else
      morvelets <- morwave(bridge,aridge,nvoice,np,N)

   cat("morvelets;  ")


   if( epsilon == 0){
      if (real == T)
         sk <- zeroskeleton2(cwtinput,Qinv,morvelets,bridge,aridge,N)
      else
         sk <- zeroskeleton(cwtinput,Qinv,morvelets,bridge,aridge,N)      
   }
   else { 
      if(real == T) 
         sk <- skeleton2(cwtinput,Qinv,morvelets,bridge,aridge,N)
      else
         sk <- skeleton(cwtinput,Qinv,morvelets,bridge,aridge,N)
   }
   cat("skeleton.\n")

   solskel <- 0 #not needed if check not done
   inputskel <- 0 #not needed if check not done

   if(check == T){
      wtsol <- cwt(Re(sk$sol),noct,nvoice)
      solskel <- complex(N)
      for(j in 1:N) solskel[j] <- wtsol[bridge[j],aridge[j]]

      inputskel <- complex(N)
      for (j in 1:N)
         inputskel[j] <- cwtinput[bridge[j],aridge[j]]
   }

   list(sol=sk$sol,A=sk$A,lam=sk$lam,dualwave=sk$dualwave,morvelets=morvelets,
	solskel=solskel,inputskel = inputskel)
}



morwave <- function(bridge, aridge, nvoice, np, N, w0 = 2*pi)
#########################################################################
#     morwave:
#     -------
#      Generation of the wavelets located on the ridge.
#
#      input:
#      ------
#      bridge: time coordinates of the ridge samples
#      aridge: scale coordinates of the ridge samples
#      nvoice: number of different scales per octave
#      np: size of the reconstruction kernel
#      N: number of complex constraints
#      w0: central frequency of Morlet wavelet
#
#      output:
#      -------
#      morvelets: array of morlet wavelets located on the ridge samples
#
#########################################################################
{

   morvelets <- matrix(0,np,2*N)

   aridge <-  2 * 2^((aridge - 1)/nvoice)
   tmp <- vecmorlet(np,N,bridge,aridge,w0 = w0)
   dim(tmp) <- c(np,N)
   morvelets[,1:N] <- Re(tmp[,1:N])
   morvelets[,(N+1):(2*N)] <- Im(tmp[,1:N])

#  Another way of generating morvelets
#   dirac <- 1:np
#   dirac[] <- 0
#   for(k in 1:N)  {
#      dirac[bridge[k]] <- 1.0;
#      scale <- 2 * 2^((aridge[k]-1)/nvoice);
#      cwtdirac <- vwt(dirac,scale,open=FALSE);
#      morvelets[,k] <- Re(cwtdirac);
#      morvelets[,k+N] <- Im(cwtdirac);
#      dirac[] <- 0
#    }

#  Still another way of generating morvelets
#   for(k in 1:N)  {
#      scale <- 2 * 2^((aridge[k]-1)/nvoice);
#      tmp <- morlet(np,bridge[k],scale)/sqrt(2*pi);
#      morvelets[,k] <- Re(tmp);
#      morvelets[,k+N] <- Im(tmp);
#    }

    morvelets
}



morwave2 <- function(bridge, aridge, nvoice, np, N, w0 = 2*pi)
#########################################################################
#     morwave2:
#     ---------
#      Generation of the real parts of morlet wavelets
#        located on the ridge.
#
#      input:
#      -------
#      bridge: time coordinates of the ridge samples
#      aridge: scale coordinates of the ridge samples
#      nvoice: number of different scales per octave
#      np: size of the reconstruction kernel
#      N: number of real constraints
#      w0: central frequency of Morlet wavelet
#
#      output:
#      -------
#      morvelets: array of morlet wavelets located on the ridge samples
#
#########################################################################
{

   morvelets <- matrix(0,np,N)

   aridge <-  2 * 2^((aridge - 1)/nvoice)
   morvelets <- Re(vecmorlet(np,N,bridge,aridge,w0 = w0))
   dim(morvelets)<- c(np,N)

# changed by Wen 5/22
#
#   aridge <-  2 * 2^((aridge - 1)/nvoice)
#   morvelets <- Re(vecmorlet(np,N,bridge,aridge,w0 = w0))
#   dim(morvelets) <- c(np,N)

   morvelets
}



skeleton <- function(cwtinput, Qinv, morvelets, bridge, aridge, N)
#########################################################################
#     skeleton:
#     --------
#      Computation of the reconstructed signal from the ridge
#
#      input:
#      -------
#      cwtinput: Continuous wavelet transform (output of cwt)
#      Qinv: inverse of the reconstruction kernel (2D array)
#      morvelets: array of morlet wavelets located on the ridge samples
#      bridge: time coordinates of the ridge samples
#      aridge: scale coordinates of the ridge samples
#      N: number of complex constraints
#
#      output:
#      -------
#      sol: reconstruction from a ridge
#      A: <wavelets,dualwavelets> matrix 
#      lam: coefficients of dual wavelets in reconstructed signal.
#      dualwave: array containing the dual wavelets.
#
##########################################################################
{
   dualwave <- Qinv %*% morvelets
   A <- t(morvelets) %*% dualwave

   rskel <- numeric(2*N)
   for(j in 1:N) {
      rskel[j] <- Re(cwtinput[bridge[j]-bridge[1]+1,aridge[j]]);
      rskel[N+j] <- -Im(cwtinput[bridge[j]-bridge[1]+1,aridge[j]])
   }

  B <- SVD(A)
  d <- B$d
  Invd <- numeric(length(d))
  for(j in 1:(2*N))
    if(d[j] < 1.0e-6) 
       Invd[j] <- 0
    else Invd[j] <- 1/d[j]

  lam <- B$v %*% diag(Invd) %*% t(B$u) %*% rskel
   sol <- dualwave%*%lam

   list(lam=lam,sol=sol,dualwave=dualwave,A=A)
}


skeleton2 <- function(cwtinput, Qinv, morvelets, bridge, aridge, N)
#########################################################################
#     skeleton2:
#     ---------
#      Computation of the reconstructed signal from the ridge, in the
#        case of real constraints.
#
#      input:
#      -------
#      cwtinput: Continuous wavelet transform (output of cwt)
#      Qinv: inverse of the reconstruction kernel (2D array)
#      morvelets: array of morlet wavelets located on the ridge samples
#      bridge: time coordinates of the ridge samples
#      aridge: scale coordinates of the ridge samples
#      N: number of real constraints
#
#      output:
#      -------
#      sol: reconstruction from a ridge
#      A: <wavelets,dualwavelets> matrix 
#      lam: coefficients of dual wavelets in reconstructed signal.
#      dualwave: array containing the dual wavelets.
#
##########################################################################
{
   dualwave <- Qinv %*% morvelets
   A <- t(morvelets) %*% dualwave

   rskel <- numeric(N)
   for(j in 1:N) {
      rskel[j] <- Re(cwtinput[bridge[j]-bridge[1]+1,aridge[j]]);
   }
   B <- SVD(A)
   d <- B$d
   Invd <- numeric(length(d))
   for(j in 1:N)
     if(d[j] < 1.0e-6) 
       Invd[j] <- 0
    else Invd[j] <- 1/d[j]

   lam <- B$v %*% diag(Invd) %*% t(B$u) %*% rskel

   sol <- dualwave%*%lam

   list(lam=lam,sol=sol,dualwave=dualwave,A=A)
}



zeroskeleton <- function(cwtinput, Qinv, morvelets, bridge, aridge, N)
#########################################################################
#     zeroskeleton:
#     -------------
#      Computation of the reconstructed signal from the ridge when
#       the epsilon parameter is set to 0.
#
#      input:
#      -------
#      cwtinput: Continuous wavelet transform (output of cwt)
#      Qinv: 1D array (warning: different from skeleton) containing
#            the diagonal of the inverse of reconstruction kernel.
#      morvelets: array of morlet wavelets located on the ridge samples
#      bridge: time coordinates of the ridge samples
#      aridge: scale coordinates of the ridge samples
#      N: number of complex constraints
#
#      output:
#      -------
#      sol: reconstruction from a ridge
#      A: <wavelets,dualwavelets> matrix 
#      lam: coefficients of dual wavelets in reconstructed signal.
#      dualwave: array containing the dual wavelets.
#
##########################################################################
{
   tmp1 <- dim(morvelets)[1]
   tmp2 <- dim(morvelets)[2]
cat(dim(morvelets))
   dualwave <- matrix(0,tmp1,tmp2)
   for (j in 1:tmp1)
      dualwave[j,] <- morvelets[j,]*Qinv[j]
   A <- t(morvelets) %*% dualwave

   rskel <- numeric(2*N)

   for(j in 1:N) {
      rskel[j] <- Re(cwtinput[bridge[j]-bridge[1]+1,aridge[j]]);
      rskel[N+j] <- -Im(cwtinput[bridge[j]-bridge[1]+1,aridge[j]])
   }

   B <- SVD(A)
   d <- B$d
   Invd <- numeric(length(d))
   for(j in 1:(2*N))
     if(d[j] < 1.0e-6) 
       Invd[j] <- 0
    else Invd[j] <- 1/d[j]

   lam <- B$v %*% diag(Invd) %*% t(B$u) %*% rskel

   sol <- dualwave%*%lam

   list(lam=lam,sol=sol,dualwave=dualwave,A=A)
}




zeroskeleton2 <- function(cwtinput, Qinv, morvelets, bridge, aridge, N)
#########################################################################
#     zeroskeleton2:
#     -------------
#      Computation of the reconstructed signal from the ridge when
#       the epsilon parameter is set to 0 and the constraints are real.
#
#      input:
#      -------
#      cwtinput: Continuous wavelet transform (output of cwt)
#      Qinv: 1D array (warning: different than in skeleton) containing
#            the diagonal of the inverse of reconstruction kernel.
#      morvelets: array of morlet wavelets located on the ridge samples
#      bridge: time coordinates of the ridge samples
#      aridge: scale coordinates of the ridge samples
#      N: number of real constraints
#
#      output:
#      -------
#      sol: reconstruction from a ridge
#      A: <wavelets,dualwavelets> matrix 
#      lam: coefficients of dual wavelets in reconstructed signal.
#      dualwave: array containing the dual wavelets.
#
##########################################################################
{
   tmp1 <- dim(morvelets)[1]
   tmp2 <- dim(morvelets)[2]
   dualwave <- matrix(0,tmp1,tmp2)
   for (j in 1:tmp1)
      dualwave[j,] <- morvelets[j,]*Qinv[j]
   A <- t(morvelets) %*% dualwave

   rskel <- numeric(N)

   for(j in 1:N) 
     rskel[j] <- Re(cwtinput[bridge[j]-bridge[1]+1,aridge[j]]);

   B <- SVD(A)
   d <- B$d
   Invd <- numeric(length(d))
   for(j in 1:N)
     if(d[j] < 1.0e-6) 
       Invd[j] <- 0
    else Invd[j] <- 1/d[j]
   lam <- B$v %*% diag(Invd) %*% t(B$u) %*% rskel

   sol <- dualwave%*%lam

   list(lam=lam,sol=sol,dualwave=morvelets,A=A)
}





regrec2 <- function(siginput, cwtinput, phi, nbnodes, noct, nvoice, Q2,
	epsilon = 0.5, w0 = 2*pi , plot = F)
#########################################################################
#     regrec2:
#     -------
#      Reconstruction of a real valued signal from a (continuous) ridge
#      (uses a regular sampling of the ridge), from a precomputed kernel.
#
#      input:
#      -------
#      siginput: input signal
#      cwtinput: Continuous wavelet transform (output of cwt)
#      phi: (unsampled) ridge
#      nbnodes: number of samples on the ridge
#      noct: number of octaves (powers of 2)
#      nvoice: number of different scales per octave
#      Q2: second part of the reconstruction kernel
#      epsilon: coeff of the Q2 term in reconstruction kernel
#      w0: central frequency of Morlet wavelet
#      plot: if set to TRUE, displays original and reconstructed signals
#
#      output: same as ridrec
#      -------
#      sol: reconstruction from a ridge
#      A: <wavelets,dualwavelets> matrix 
#      lam: coefficients of dual wavelets in reconstructed signal.
#      dualwave: array containing the dual wavelets.
#      morvelets: array of morlet wavelets located on the ridge samples
#      solskel: wavelet transform of sol, restricted to the ridge
#      inputskel: wavelet transform of signal, restricted to the ridge
#
#########################################################################
{
   tmp <- RidgeSampling(phi,nbnodes)
   node <- tmp$node
   phinode <- tmp$phinode

   np <- dim(Q2)[1]
   one <- diag(np)
   Q <- one + epsilon * Re(Q2)

   if(epsilon == 0)
      Qinv <- one
   else
      Qinv <- solve(Q)

   tmp2 <- ridrec(cwtinput,node,phinode,noct,nvoice,Qinv,epsilon,np)

   if(plot == T){
      par(mfrow=c(2,1))
      plot.ts(Re(siginput))
      title("Original signal")
      plot.ts(Re(tmp2$sol))
      title("Reconstructed signal")
   }
   tmp2
}




RidgeSampling <- function(phi, nbnodes)
#########################################################################
#    RidgeSampling:
#    --------------
#      Given a ridge phi (for the Gabor transform), returns a 
#        (regularly) subsampled version of length nbnodes.
#
#      Input:
#      ------
#      phi: ridge
#      nbnodes: number of samples.
#
#      Output:
#      -------
#      node: time coordinates of the ridge samples
#      phinode: frequency coordinates of the ridge samples
#
#########################################################################
{
   node <- numeric(nbnodes)
   phinode <- numeric(nbnodes)
   N <- nbnodes
   LL <- length(phi)

   for (i in 1:(N+1)) node[i] <- as.integer(((LL-1)*i+nbnodes-LL+1)/nbnodes)
   for (i in 1:(N+1)) phinode[i] <- phi[node[i]]

   list(node = node,phinode = phinode)
}



wRidgeSampling <- function(phi, compr, nvoice)
#########################################################################
#    wRidgeSampling:
#    ---------------
#      Given a ridge phi, returns a subsampled version of length
#      nbnodes (wavelet case).
#
#      Input:
#      ------
#      phi: ridge
#      compr: subsampling rate for the wavelet coefficients (at scale 1)
#             when compr <= 0, uses all the points in a ridge
#
#      Output:
#      -------
#      node: time coordinates of the ridge samples
#      phinode: frequency coordinates of the ridge samples
#      nbnodes: number of samples.
#
#########################################################################
{


   LL <- length(phi)
   node <- numeric(LL)
   phinode <- numeric(LL)

   LN <- 1
   j <- 1
   node[1] <- 1

# We multiply a number on the phi. This number is obtained experimentally
# If we have a ridge which is linear of scale, then the subsampling is
# adapted with s, however, in the wavelet case, the ridge is linear 
# in log(scale), the subsampling rate is therefore adapted to k*s^l
# l here is the constant 100/15, and k is compr 
   aridge <- 2 * 2^((phi-1)/(nvoice * 100/15))

#   while (node[j] < LL){
#      j <- j + 1
#      LN <- round(compr * aridge[node[j-1]])
#      cat("LN=",LN)
#      node[j] <- node[j-1] + LN
#      cat("; node=",node[j],"\n")
#      phinode[j] <- phi[node[j]]
#      cat("; phinode=",phinode[j],"\n")
#   }

   nbnodes <- 0
   while (node[j] < LL){
      phinode[j] <- phi[node[j]]
      nbnodes <- nbnodes + 1

# To guarantee at least advance by one 
       LN <- max(1,round(compr * aridge[node[j]]))

      j <- j + 1
      node[j] <- node[j-1] + LN
   }


   tmpnode <- numeric(nbnodes)
   tmpphinode <- numeric(nbnodes)
   tmpnode <- node[1:nbnodes]
   tmpphinode <- phinode[1:nbnodes]

   list(node = tmpnode,phinode = tmpphinode, nbnodes = nbnodes)
}




sridrec <- function(tfinput, ridge)
#########################################################################
#     sridrec
#     -------
#      Simple reconstruction of a real valued signal from a ridge,
#      by restriction of the wavelet transform.
#
#      Input:
#      ------
#     tfinput: Continuous time-frequency transform (output of cwt or cgt)
#     ridge: (unsampled) ridge
#
#      Output:
#      -------
#      rec: reconstructed signal
#
#########################################################################
{
   bsize <- dim(tfinput)[1]
   asize <- dim(tfinput)[2]
  
   rec <- numeric(bsize)
   rec[] <- 0 + 0i

   for(i in 1:bsize) {
      if(ridge[i] > 0)
        rec[i] <- rec[i] + 2 * tfinput[i,ridge[i]]
   }
  
   Re(rec)
}




#########################################################################
#
#               (c) Copyright  1997                             
#                          by                                   
#      Author: Rene Carmona, Bruno Torresani, Wen-Liang Hwang   
#                  Princeton University 
#                  All right reserved                           
#########################################################################

#*********************************************************************#
# mnpval
#*********************************************************************#

mnpval <- function(inputdata, maxresoln, wl=128, scale=FALSE)
{
  x <- adjust.length(inputdata)
  s <- x$signal
  np <- x$length
  num.of.windows <- (np/wl - 1) * 4 + 1

  pval <- matrix(0, nrow=maxresoln, ncol=np)
  pval <- t(pval)
  dim(pval) <- c(length(pval), 1)

  z <- .C("normal_pval_compute",
	a=as.single(pval),
	as.single(s),
	as.integer(maxresoln),
	as.integer(np),
	as.integer(num.of.windows),
	as.integer(wl) )

  pval <- t(z$a)
  dim(pval) <- c(np, maxresoln)

  plot.result(pval, s, maxresoln, scale)
  list( original=s, pval=pval, maxresoln=maxresoln, np=np )
}

#*********************************************************************#
# mbpval
#*********************************************************************#

mbpval <- function(inputdata, maxresoln, wl=128, scale=FALSE)
{
  x <- adjust.length(inputdata)
  s <- x$signal
  np <- x$length
  num.of.windows <- (np/wl - 1) * 4 + 1

  pval <- matrix(0, nrow=maxresoln, ncol=np)
  pval <- t(pval)
  dim(pval) <- c(length(pval), 1)

  z <- .C("compute_mallat_bootstrap_pval",
	a=as.single(pval),
	as.single(s),
	as.integer(maxresoln),
	as.integer(np),
	as.integer(num.of.windows),
	as.integer(wl) )

  pval <- t(z$a)
  dim(pval) <- c(np, maxresoln)

  plot.result(pval, s, maxresoln, scale)
  list( original=s, pval=pval, maxresoln=maxresoln, np=np )
}

#*********************************************************************#
# mntrim
#*********************************************************************#

mntrim <- function(extrema, scale=FALSE, prct=.95)
{
  s <- extrema$original
  maxresoln <- extrema$maxresoln
  np <- extrema$np
  sample.size <- 128
  nthresh <- 1:maxresoln

  # Find the threshold for each resoln

  z <- .C("nthresh_compute",
	a=as.single(nthresh),
	as.single(s),
	as.integer(maxresoln),
	as.integer(sample.size),
	as.double(prct) )

  nthresh <- z$a

  trim <- matrix(0, nrow=np, ncol=maxresoln)

  for (j in 1:maxresoln)
  {
    # Keep the extrema if the absolute value of the extrema >= threshold
    temp <- (abs(extrema$extrema[,j]) >= nthresh[j])
    trim[,j] <- temp * extrema$extrema[,j]
  }

  cat("number of extrema left =", sum(trim!=0), "\n")

  plot.result(trim, s, maxresoln, scale)
  list( original=s, extrema=trim, Sf=extrema$Sf, maxresoln=maxresoln, np=np )
}

#*********************************************************************#
# mbtrim
#*********************************************************************#

mbtrim <- function(extrema, scale=FALSE, prct=.95)
{
  s <- extrema$original
  maxresoln <- extrema$maxresoln
  np <- extrema$np
  sample.size <- 128
  bthresh <- 1:maxresoln

  # Find the threshold for each resoln

  z <- .C("bthresh_compute",
	a=as.single(bthresh),
	as.single(s),
	as.integer(maxresoln),
	as.integer(sample.size),
	as.double(prct) )

  bthresh <- z$a

  trim <- matrix(0, nrow=np, ncol=maxresoln)
  for (j in 1:maxresoln)
  {
    # Keep the extrema if the absolute value of the extrema >= threshold
    temp <- (abs(extrema$extrema[,j]) >= bthresh[j])
    trim[,j] <- temp * extrema$extrema[,j]
  }

  cat("number of extrema left =", sum(trim!=0), "\n")

  plot.result(trim, s, maxresoln, scale)
  list( original=s, extrema=trim, Sf=extrema$Sf, maxresoln=maxresoln, np=np )
}










#########################################################################
#      $Log: Skernel.S,v $
# Revision 1.2  1995/04/05  18:56:55  bruno
# *** empty log message ***
#
# Revision 1.1  1995/04/02  01:04:16  bruno
# Initial revision
#
#               (c) Copyright  1997                             
#                          by                                   
#      Author: Rene Carmona, Bruno Torresani, Wen-Liang Hwang   
#                  Princeton University
#                  All right reserved                           
#########################################################################





kernel <- function(node, phinode, nvoice, x.inc = 1, x.min = node[1],
	x.max = node[length(node)], w0 = 2*pi, plot = F)
#########################################################################
#     kernel:   
#     ------
#      Computes the cost from the sample of points on the estimated ridge
#      and the matrix used in the reconstruction of the original signal
#      The output is a lng x lng matrix of complex numbers, lng being the
#      number of points at which the signal is to be reconstructed. 
#      The dependence upon the original signal is only through the 
#      sample (node,phinode) of the ridge. This should be the output
#      of a previously run estimation procedure.
#
#     Input:
#     ------
#      node: values of the variable b for the nodes of the ridge
#      phinode: values of the scale variable a for the nodes of the ridge
#      nvoice: number of scales within 1 octave
#      x.inc: step unit for the computation of the kernel
#      x.min: minimal value of x for the computation of Q2
#      x.max: maximal value of x for the computation of Q2
#      w0: central frequency of the wavelet
#      plot: if set to TRUE, displays the modulus of the matrix of Q2
#
#     Output:
#     -------
#      ker: matrix of the Q2 kernel
#
#########################################################################
{

cat("x.min = ",x.min,"\n")
cat("x.max = ",x.max,"\n")
cat("x.inc = ",x.inc,"\n")
   lng <- as.integer((x.max-x.min)/x.inc)+1
   nbnode <- length(node)

   cat("lng = ",lng,"\n");

   b.start <- x.min - 50
   b.end <- x.max + 50
   
   ker.r <- matrix(0, lng, lng)
   ker.i <- matrix(0, lng, lng)

   dim(ker.r) <- c(lng * lng, 1)
   dim(ker.i) <- c(lng * lng, 1)
   phinode <-  2 * 2^(phinode/nvoice)
  

  z <- .C("kernel",
        ker.r = as.double(ker.r),
        ker.i = as.double(ker.i),
        as.integer(x.min),
        as.integer(x.max),
	as.integer(x.inc),
	as.integer(lng),
        as.double(node),
        as.double(phinode),
        as.integer(nbnode),
        as.double(w0),
        as.double(b.start),
        as.double(b.end))


  ker.r <- z$ker.r
  ker.i <- z$ker.i
  dim(ker.r) <- c(lng, lng)
  dim(ker.i) <- c(lng, lng)
  ker <- matrix(0, lng, lng)
  i <- sqrt(as.complex(-1))
  ker <- ker.r + i * ker.i  

  if(plot == T){
     par(mfrow=c(1,1))
     image(Mod(ker))
     title("Matrix of the reconstructing kernel (modulus)")
  }

  ker
}

rkernel <- function(node, phinode, nvoice, x.inc = 1, x.min = node[1],
	x.max = node[length(node)], w0 = 2*pi, plot = F)
#########################################################################
#     rkernel:   
#     -------
#      Same as kernel, except that a real valued kernel is computed
#      (precisely the real part of the previous kernel); this applies to
#      real signals.
#
#     Input:
#     ------
#      node: values of the variable b for the nodes of the ridge
#      phinode: values of the scale variable a for the nodes of the ridge
#      nvoice: number of scales within 1 octave
#      x.inc: step unit for the computation of the kernel
#      x.min: minimal value of x for the computation of Q2
#      x.max: maximal value of x for the computation of Q2
#      w0: central frequency of the wavelet
#      plot: if set to TRUE, displays the modulus of the matrix of Q2
#
#     Output:
#     -------
#      ker: matrix of the Q2 kernel
#
#########################################################################
{
   lng <- as.integer((x.max-x.min)/x.inc)+1
   nbnode <- length(node)

   b.start <- x.min - 50
   b.end <- x.max + 50

   ker <- matrix(0, lng, lng)
   dim(ker) <- c(lng * lng, 1)
   phinode <-  2 * 2^(phinode/nvoice)

  z <- .C("rkernel",
        ker = as.double(ker),
        as.integer(x.min),
        as.integer(x.max),
	as.integer(x.inc),
	as.integer(lng),
        as.double(node),
        as.double(phinode),
        as.integer(nbnode),
        as.double(w0),
        as.double(b.start),
        as.double(b.end))

  ker <- z$ker
  dim(ker) <- c(lng, lng)

  if(plot == T){
     par(mfrow=c(1,1))
     image(Mod(ker))
     title("Matrix of the Q2 kernel (modulus)")
  }

  ker
}



fastkernel <- function(node, phinode, nvoice, x.inc = 1, x.min = node[1],
	x.max = node[length(node)], w0 = 2*pi, plot = F)
#########################################################################
#     fastkernel:   
#     -----------
#	   Same as kernel, except that the kernel is computed
#	     using Riemann sums instead of Romberg integration.
#
#     Input:
#     ------
#      node: values of the variable b for the nodes of the ridge
#      phinode: values of the scale variable a for the nodes of the ridge
#      nvoice: number of scales within 1 octave
#      x.inc: step unit for the computation of the kernel
#      x.min: minimal value of x for the computation of Q2
#      x.max: maximal value of x for the computation of Q2
#      w0: central frequency of the wavelet
#      plot: if set to TRUE, displays the modulus of the matrix of Q2
#
#     Output:
#     -------
#      ker: matrix of the Q2 kernel
#
#########################################################################
{
   lng <- as.integer((x.max-x.min)/x.inc)+1
   nbnode <- length(node)

   b.start <- x.min - 50
   b.end <- x.max + 50

   ker.r <- matrix(0, lng, lng)
   ker.i <- matrix(0, lng, lng)
   dim(ker.r) <- c(lng * lng, 1)
   dim(ker.i) <- c(lng * lng, 1)
   phinode <-  2 * 2^(phinode/nvoice)

  z <- .C("fastkernel",
        ker.r = as.double(ker.r),
        ker.i = as.double(ker.i),
        as.integer(x.min),
        as.integer(x.max),
	as.integer(x.inc),
	as.integer(lng),
        as.double(node),
        as.double(phinode),
        as.integer(nbnode),
        as.double(w0),
        as.double(b.start),
        as.double(b.end))


  ker.r <- z$ker.r
  ker.i <- z$ker.i
  dim(ker.r) <- c(lng, lng)
  dim(ker.i) <- c(lng, lng)
  ker <- matrix(0, lng, lng)
  i <- sqrt(as.complex(-1))
  ker <- ker.r + i * ker.i  

  if(plot == T){
     par(mfrow=c(1,1))
     image(Mod(ker))
     title("Matrix of the reconstructing kernel (modulus)")
  }

  ker
}



zerokernel <- function(x.inc = 1, x.min,x.max)
#########################################################################
#     zerokernel:   
#     -----------
#	Returns a zero kernel.
#
#     Input:
#     ------
#      x.inc: step unit for the computation of the kernel
#      x.min: minimal value of x for the computation of Q2
#      x.max: maximal value of x for the computation of Q2
#
#     Output:
#     -------
#      ker: matrix of the Q2 kernel
#
#########################################################################
{
   lng <- as.integer((x.max-x.min)/x.inc)+1

   ker <- matrix(0, lng, lng)

   ker
}




gkernel <- function(node, phinode, freqstep, scale, x.inc = 1,
	x.min = node[1], x.max = node[length(node)], plot = F)
#########################################################################
#     gkernel:   
#     -------
#	Same as kernel, for the case of Gabor transform.
#
#     input:
#     ------
#      node: values of the variable b for the nodes of the ridge
#      phinode: values of the frequency variable for the ridge nodes
#      freqstep: sampling rate for the frequency axis
#      scale: size of the window
#      x.inc: step unit for the computation of the kernel
#      x.min: minimal value of x for the computation of Q2
#      x.max: maximal value of x for the computation of Q2
#      plot: if set to TRUE, displays the modulus of the matrix of Q2
#
#     output:
#     -------
#      ker: matrix of the Q2 kernel
#
#########################################################################
{
   lng <- as.integer((x.max-x.min)/x.inc)+1
   nbnode <- length(node)
   b.start <- node[1] - 50
   b.end <- node[length(node)] + 50
   ker <- matrix(0, lng, lng)
   dim(ker) <- c(lng * lng, 1)

   phinode <- phinode * freqstep


  z <- .C("gkernel",
        ker = as.double(ker),
        as.integer(x.min),
        as.integer(x.max),
	as.integer(x.inc),
	as.integer(lng),
        as.double(node),
        as.double(phinode),
        as.integer(nbnode),
        as.double(scale),
        as.double(b.start),
        as.double(b.end))


  ker <- z$ker
  dim(ker) <- c(lng, lng)

  if(plot == T){
     par(mfrow=c(1,1))
     image(Mod(ker))
     title("Matrix of the reconstructing kernel (modulus)")
  }

  ker
}



fastgkernel <- function(node, phinode, freqstep, scale, x.inc = 1,
	x.min = node[1], x.max = node[length(node)], plot = F)
#########################################################################
#     fastgkernel:   
#     ------------
#	Same as gkernel, except that the kernel is computed
#	  using Riemann sums instead of Romberg integration.
#
#
#     input:
#     ------
#      node: values of the variable b for the nodes of the ridge
#      phinode: values of the frequency variable for the ridge nodes
#      freqstep: sampling rate for the frequency axis
#      scale: size of the window
#      x.inc: step unit for the computation of the kernel
#      x.min: minimal value of x for the computation of Q2
#      x.max: maximal value of x for the computation of Q2
#      plot: if set to TRUE, displays the modulus of the matrix of Q2
#
#     output:
#     -------
#      ker: matrix of the Q2 kernel
#
#########################################################################
{
   lng <- as.integer((x.max-x.min)/x.inc)+1
   nbnode <- length(node)

   b.start <- x.min - 50
   b.end <- x.max + 50

   ker <- matrix(0, lng, lng)
   dim(ker) <- c(lng * lng, 1)
   phinode <-  phinode*freqstep


  z <- .C("fastgkernel",
        ker = as.double(ker),
        as.integer(x.min),
        as.integer(x.max),
	as.integer(x.inc),
	as.integer(lng),
        as.double(node),
        as.double(phinode),
        as.integer(nbnode),
        as.double(scale),
        as.double(b.start),
        as.double(b.end))


  ker <- z$ker
  dim(ker) <- c(lng, lng)

  if(plot == T){
     par(mfrow=c(1,1))
     image(Mod(ker))
     title("Matrix of the reconstructing kernel (modulus)")
  }

  ker
}









#########################################################################
#      $Log: Snake_Annealing.S,v $
# Revision 1.2  1995/04/05  18:56:55  bruno
# *** empty log message ***
#
# Revision 1.1  1995/04/02  01:04:16  bruno
# Initial revision
#
#               (c) Copyright  1997
#                          by                                   
#      Author: Rene Carmona, Bruno Torresani, Wen-Liang Hwang   
#                  Princeton University
#                  All right reserved                           
#########################################################################



snake <- function(tfrep, guessA, guessB, snakesize = length(guessB),
	tfspec = numeric(dim(modulus)[2]), subrate = 1, temprate = 3,
	muA = 1, muB = muA, lambdaB = 2*muB, lambdaA = 2*muA,
	iteration = 1000000, seed = -7, costsub = 1, stagnant = 20000,plot=T)
#########################################################################
#     snake:   
#     ------
#      ridge extraction using snakes and simulated annealing
#
#       input:
#       ------
# 	 tfrep: the wavelet or Gabor transform.
#	 guessA: initial guess for ridge function (frequency coordinate).
#	 guessB: initial guess for ridge function (time coordinate).
#        tfspec: estimate for the contribution of the noise to the
#               transform modulus (1D array)
#        subrate: subsampling rate for the transform modulus.
#	 temprate: constant (numerator) in the temperature schedule.
#	 muA, muB: coeff of phi' and rho' in cost function
#	 lambdaA,lambdaB: same thing for 2nd derivatives
#        iteration: maximal number of iterations
#        seed: initialization for random numbers generator
#        costsub: subsampling of the cost function in cost
#               costsub=1 means that the whole cost function
#               is returned
#        stagnant: maximum number of steps without move (for the
#                  stopping criterion)
#        plot :	 when set (by default), certain results will be
#	         displayed	
#
#       output:
#       ------
#        A: 1D array containing the frequency coordinate of the ridge
#        B: 1D array containing the time coordinate of the ridge
#        cost: 1D array containing the cost function
#
#########################################################################
{

  sigsize <- dim(tfrep)[1]
  nscale <- dim(tfrep)[2]
  costsize <- as.integer(iteration/costsub) + 1	
  cost <- 1:costsize
  modulus <- Mod(tfrep)	
  tfspectrum <- tfspec

  cost[] <- 0
  phi <- as.integer(guessA)
  rho <- as.integer(guessB)
  count <- 0	
  dim(phi) <- c(length(phi),1)	
  dim(rho) <- c(length(rho),1)	
  if(plot== T)  image(modulus)
  dim(modulus) <- c(sigsize * nscale, 1)


  smodsize <- as.integer(sigsize/subrate)
  if((sigsize/subrate - as.integer(sigsize/subrate)) > 0)
    smodsize <- as.integer(sigsize/subrate) + 1

  smodulus <- matrix(0,smodsize,nscale)       
  dim(smodulus) <- c(smodsize * nscale, 1)

  z <- .C("Smodulus_smoothing",
        as.double(modulus),
        smodulus = as.double(smodulus),
	as.integer(sigsize),
        smodsize = as.integer(smodsize),            
        as.integer(nscale),
        as.integer(subrate))

  smodulus <- z$smodulus
  smodsize <- z$smodsize
  smodulus <- smodulus * smodulus
  dim(smodulus) <- c(smodsize, nscale)
  for (k in 1:nscale)
    smodulus[,k] <- smodulus[,k] - tfspectrum[k]
  dim(smodulus) <- c(smodsize * nscale, 1)


   z <- .C("Ssnake_annealing",
        cost = as.single(cost),
        as.double(smodulus),
        phi = as.single(phi),
        rho = as.single(rho),
        as.single(lambdaA),
	as.single(muA),
        as.single(lambdaB),
	as.single(muB),
        as.single(temprate),
	as.integer(sigsize),
	as.integer(snakesize),
        as.integer(nscale),
        as.integer(iteration),
	as.integer(stagnant),
	as.integer(seed),
	nb = as.integer(count),
        as.integer(subrate),
	as.integer(costsub),
        as.integer(smodsize))

  count <- z$nb
  cat("Number of moves:",count,"(",stagnant," still steps)\n")
  if(plot==T) lines(z$rho+1, z$phi+1)
  list(A = z$phi+1, B = z$rho+1, cost = z$cost[1:count])

}
  


   
snakoid <- function(modulus, guessA, guessB, snakesize = length(guessB),
	tfspec = numeric(dim(modulus)[2]), subrate = 1, temprate = 3,
	muA = 1, muB = muA, lambdaB = 2*muB, lambdaA = 2*muA,
	iteration = 1000000, seed = -7, costsub = 1, stagnant = 20000,plot=T)
#########################################################################
#     snakoid:   
#     ----------
#      ridge extraction using snakes and simulated annealing
#
#       input:
#       ------
# 	 modulus: modulus of the wavelet or Gabor transform.
#	 guessA: initial guess for ridge function (frequency coordinate).
#	 guessB: initial guess for ridge function (time coordinate).
#        tfspec: estimate for the contribution of the noise to the
#               transform modulus (1D array)
#        subrate: subsampling rate for the transform modulus.
#	 temprate: constant (numerator) in the temperature schedule.
#	 muA, muB: coeff of phi' and rho' in cost function
#	 lambdaA,lambdaB: same thing for 2nd derivatives
#        iteration: maximal number of iterations
#        seed: initialization for random numbers generator
#        costsub: subsampling of the cost function in cost
#               costsub=1 means that the whole cost function
#               is returned
#        stagnant: maximum number of steps without move (for the
#                  stopping criterion)
#        plot:	 when set(default), some results will be displayed	
#
#       output:
#       ------
#        A: 1D array containing the frequency coordinate of the ridge
#        B: 1D array containing the time coordinate of the ridge
#        cost: 1D array containing the cost function
#
#########################################################################
{

  sigsize <- dim(modulus)[1]
  nscale <- dim(modulus)[2]
  costsize <- as.integer(iteration/costsub) + 1	
  cost <- 1:costsize
  tfspectrum <- tfspec 

  cost[] <- 0
  phi <- guessA
  rho <- guessB
  count <- 0	
  dim(phi) <- c(length(phi),1)	
  dim(rho) <- c(length(rho),1)
  if(plot== T)  image(modulus)	
  dim(modulus) <- c(sigsize * nscale, 1)


  smodsize <- as.integer(sigsize/subrate)
  if((sigsize/subrate - as.integer(sigsize/subrate)) > 0)
    smodsize <- as.integer(sigsize/subrate) + 1

  smodulus <- matrix(0,smodsize,nscale)       
  dim(smodulus) <- c(smodsize * nscale, 1)

  z <- .C("Smodulus_smoothing",
        as.double(modulus),
        smodulus = as.double(smodulus),
	as.integer(sigsize),
        smodsize = as.integer(smodsize),            
        as.integer(nscale),
        as.single(subrate))

  smodulus <- z$smodulus
  smodsize <- z$smodsize
  smodulus <- smodulus * smodulus
  dim(smodulus) <- c(smodsize, nscale)
  for (k in 1:nscale)
    smodulus[,k] <- smodulus[,k] - tfspectrum[k]
  dim(smodulus) <- c(smodsize * nscale, 1)


   z <- .C("Ssnakenoid_annealing",
        cost = as.single(cost),
        as.double(smodulus),
        phi = as.single(phi),
        rho = as.single(rho),
        as.single(lambda),
	as.single(mu),
        as.single(lambda2),
	as.single(mu2),
        as.single(temprate),
	as.integer(sigsize),
	as.integer(snakesize),
        as.integer(nscale),
        as.integer(iteration),
	as.integer(stagnant),
	as.integer(seed),
	nb = as.integer(count),
        as.integer(subrate),
	as.integer(costsub),
        as.integer(smodsize))

  count <- z$nb
  cat("Number of moves:",count,"(",stagnant," still steps)\n")
  if(plot== T)  lines(z$rho+1, z$phi+1) 
  list(A = z$phi+1, B = z$rho+1, cost = z$cost[1:count])

}
  


snakeview <- function(modulus,snake)
#########################################################################
#      snakeview:
#      ----------
#       Restrict time-frequency transform to a snake
#
#       input:
#       ------
#	 modulus: modulus of the wavelet (or Gabor) transform
#        snake: B and A components of a snake
#
#       output:
#       -------
#        img: 2D array containing the restriction of the transform
#             modulus to the snake
#
#########################################################################
{
   A <- snake$A
   B <- snake$B
   snakesize <- length(A)
   sigsize <- dim(modulus)[1]
   nscale <- dim(modulus)[2]

   img <- matrix(0,sigsize,nscale)
   for(i in 1:snakesize)
     img[B[i],A[i]] <- modulus[B[i],A[i]]

   img
}

     









#########################################################################
#       $Log: TF_Maxima.S,v $
#
#               (c) Copyright  1997
#                          by                                   
#      Author: Rene Carmona, Bruno Torresani, Wen-Liang Hwang   
#                  Princeton University 
#                  All right reserved                           
#########################################################################





tfgmax <- function(input, plot = TRUE)
#########################################################################
#  tfgmax:  
#  -------
# 	 Continuous time-frequency transform global maxima:
#          compute the continuous wavelet transform global maxima (for
#    	   fixed position)
#
#      input:
#      ------
# 	input: continuous time-frequency transform (2D array)
#	plot: if set to TRUE, displays the maxima of cwt on the graphic
#          device.
#
#      output:
#      -------
#       output: values of the maxima (1D array)
#       pos: positions of the maxima (1D array)
#       
#########################################################################
{

  sigsize <- dim(input)[1]
  pp <- dim(input)[2]
  input1 <- input
	
  output <- matrix(0,sigsize,pp)
  dim(input1) <- c(sigsize * pp,1)
  dim(output) <- c(sigsize * pp,1)
  posvector <- 1:sigsize
  posvector[] <- 0

  z <- .C("Scwt_gmax",
           as.double(input1),
           output = as.double(output),
           as.integer(sigsize),
           as.integer(pp),
           pos = as.integer(posvector))

  output <- z$output
  pos <- z$pos

  dim(output) <- c(sigsize, pp)
  if(plot)image(output)
  list(output=output, pos= pos)
}



tflmax <- function(input, plot = TRUE)
#########################################################################
#  tflmax:
#  -------
# 	continuous time-frequency transform local maxima:
#        compute the time-frequency transform local maxima (for
#    	 fixed position)
#
#      input:
#      ------
# 	input: continuous time-frequency transform (2D array)
#	plot: if set to TRUE, displays the maxima of cwt on the graphic
#          device.
#
#      output:
#      -------
#       output: values of the maxima (2D array)
#
#########################################################################
{

  sigsize <- dim(input)[1]
  pp <- dim(input)[2]
  input1 <- input

  output <- matrix(0,sigsize,pp)
  dim(input1) <- c(sigsize * pp,1)
  dim(output) <- c(sigsize * pp,1)

  z <- .C("Scwt_mridge",
           as.double(input1),
           output = as.double(output),
           as.integer(sigsize),
           as.integer(pp))

  output <- z$output
  dim(output) <- c(sigsize, pp)
  if(plot)image(output)
  output
}



cleanph <- function(tfrep, thresh = .01, plot = TRUE)
#########################################################################
#  cleanph:
#  --------
# 	sets to zero the phase of time-frequency transform when
#        modulus is below a certain value.
#
#      input:
#      ------
# 	tfrep: continuous time-frequency transform (2D array)
#       thresh: (relative) threshold.
#	plot: if set to TRUE, displays the maxima of cwt on the graphic
#          device.
#
#      output:
#      -------
#       output: thresholded phase (2D array)
#
#########################################################################
{

  thrmod1 <- Mod(tfrep)
  thrmod2 <- Mod(tfrep)
  limit <- range(thrmod1)[2]*thresh
  thrmod1 <- (Mod(tfrep) > limit)
  thrmod2 <- (Mod(tfrep) <= limit)

  output <- thrmod1*Arg(tfrep) - pi * thrmod2

  if(plot)image(output)
  output
}

























#########################################################################
#      $Log: WV.S,v $
#########################################################################
#
#               (c) Copyright  1998
#                          by                                   
#      Author: Rene Carmona, Bruno Torresani, Wen-Liang Hwang   
#                  Princeton University 
#                  All right reserved                           
#########################################################################





WV <- function(input, nvoice, freqstep = (1/nvoice), plot = TRUE)
#########################################################################
#       WV:   
#       ---
#        Wigner-Ville function
# 	  compute the Wigner-Ville transform, without any smoothing
#
#       input:
#       ------
# 	 input: input signal (possibly complex-valued)
#	 nvoice: number of frequency bands
#        freqstep: sampling rate for the frequency axis
#	 plot: if set to TRUE, displays the modulus of cwt on the graphic
#		device.
#
#       output:
#       -------
#        output: (complex) Wigner-Ville transform
#
#########################################################################
{
   oldinput <- input
   isize <- length(oldinput)

   tmp <- adjust.length(oldinput)
   input <- tmp$signal
   newsize <- length(input)

#   pp <- nvoice
   pp <- newsize
   Routput <- matrix(0,newsize,pp)
   Ioutput <- matrix(0,newsize,pp)
   output <- matrix(0,newsize,pp)
   dim(Routput) <- c(pp * newsize,1)
   dim(Ioutput) <- c(pp * newsize,1)
   dim(input) <- c(newsize,1)


   z <- .C("WV",
            as.single(input),
            Rtmp = as.double(Routput),
            Itmp = as.double(Ioutput),
            as.integer(nvoice),
	    as.single(freqstep),
            as.integer(newsize))

   Routput <- z$Rtmp
   Ioutput <- z$Itmp
   dim(Routput) <- c(newsize,pp)
   dim(Ioutput) <- c(newsize,pp)

       
   i <- sqrt(as.complex(-1))
   output <- Routput[1:isize,] + Ioutput[1:isize,] * i
   if(plot) {
      image(Mod(output[,1:isize/2]),xlab="Time",ylab="Frequency")
      title("Wigner-Ville Transform Modulus")
   }
   output
}












#########################################################################
#    $log: extrema.S,v $
#########################################################################
#    (c) Copyright 1997
#               by
#   Author: Rene Carmona, Bruno Torresani, Wen L. Hwang, Andrea Wang
#              Princeton University 
#              All right reserved
########################################################################




ext <- function(wt, scale=FALSE, plot=TRUE)
#****************************************************************
# ext
# ---
#   compute the extrema of the dyadic wavelet transform.
#
# input
# -----
#   wt: wavelet transform
#   scale: a flag indicating if the extrema at each 
#	   resolution will be plotted at the same sacle.
#
# output
# ------
#   original: original signal
#   extrema: extrema representation
#   Sf: coarse resolution of signal
#   maxresoln: number of decomposition
#   np: size of signal
#****************************************************************
{
  s <- wt$original
  maxresoln <- wt$maxresoln
  np <- wt$np

  extrema <- matrix(0, nrow=maxresoln, ncol=np)
  extrema <- t(extrema)
  dim(extrema) <- c(length(extrema), 1)

  t(wt$Wf)   # Note: transposed wt is not assigned to wt
  dim(wt$Wf) <- c(length(wt$Wf), 1)

  z <- .C("modulus_maxima", 
	a=as.single(extrema),
	as.single(wt$Wf),
	as.integer(maxresoln), 
	as.integer(np) )

  extrema <- t(z$a)
  dim(extrema) <- c(np, maxresoln)
  
  cat("number of extrema =", sum(extrema!=0), "\n")

  if(plot)  plot.result(extrema, s, maxresoln, scale)

    list(original=s,extrema=extrema,Sf=wt$Sf,maxresoln=maxresoln,np=np)
}












#########################################################################
#      $Log: gRidge_Recons.S,v $
# Revision 1.2  1995/04/05  18:56:55  bruno
# *** empty log message ***
#
# Revision 1.1  1995/04/02  01:04:16  bruno
# Initial revision
#
#
#               (c) Copyright  1997                             
#                          by                                   
#      Author: Rene Carmona, Bruno Torresani, Wen-Liang Hwang   
#                  University of California, Irvine             
#                  All right reserved                           
#########################################################################




girregrec <- function(siginput, gtinput, phi, nbnodes, nvoice,
	freqstep, scale, epsilon = 0.5,fast = F, prob= 0.8, plot = F,
	para = 0, hflag = F, real = F, check = F)
#########################################################################
#     girregrec:
#     ---------
#      Reconstruction of a real valued signal from a (continuous) 
#      Gabor ridge (uses a irregular sampling of the ridge)
#
#      input:
#      -------
#      siginput: input signal
#      gtinput: Continuous gabor transform (output of cgt)
#      phi: (unsampled) ridge
#      nbnodes: number of nodes used for the reconstruction.
#      nvoice: number of different scales per octave
#      freqstep: sampling rate for the frequency axis
#      epsilon: coeff of the Q2 term in reconstruction kernel
#      fast: if set to TRUE, the kernel is computed using Riemann
#           sums instead of Romberg's quadrature
#      plot: if set to TRUE, displays original and reconstructed signals
#      para: no comment !
#      prob: sampling (irregularly) at the places of the ridge where their lambdas 
#            are more than prob in lambda distribution.
#      hflag: 
#      real: if set to TRUE, only uses constraints on the real part
#            of the transform for the reconstruction.
#      check: if set to TRUE, computes the wavelet transform of the
#            reconstructed signal
#      minnbnodes: minimum number of nodes for the reconstruction
#
#      output:
#      -------
#      sol: reconstruction from a ridge
#      A: <wavelets,dualwavelets> matrix 
#      lam: coefficients of dual wavelets in reconstructed signal.
#      dualwave: array containing the dual wavelets.
#      solskel: wavelet transform of sol, restricted to the ridge
#      inputskel: wavelet transform of signal, restricted to the ridge
#      Q2: second part of the reconstruction kernel
#########################################################################
{

#
# first performing regular sampling
#

   tmp <- RidgeSampling(phi,nbnodes)

   node <- tmp$node
   phinode <- tmp$phinode

   phi.x.min <- scale
   phi.x.max <- scale
   
   x.min <- node[1]
   x.max <- node[length(node)]
   x.max <- x.max + round(para * phi.x.max)
   x.min <- x.min - round(para * phi.x.min)

   node <- node - x.min + 1
   
   x.inc <- 1
   np <- as.integer((x.max-x.min)/x.inc)+1

   cat(" (np:",np,")")

   if(epsilon == 0)
       Q2 <- 0
   else {
      if (fast == F)
	Q2 <- gkernel(node,phinode,freqstep,scale,x.min = x.min, x.max = x.max)
      else
	Q2 <- fastgkernel(node,phinode,freqstep,scale,x.min = x.min,
		x.max = x.max)
      }
   cat(" kernel;")

# Generating the Q1 term in reconstruction kernel
   if (hflag == T)
      one <- gsampleOne(node,scale,np)
   else{
      one <- numeric(np)
      one[] <- 1
   }

   if (epsilon !=0 ){
     Q <-  epsilon * Q2
     for(j in 1:np)
        Q[j,j] <- Q[j,j] + one[j]
     Qinv <- solve(Q)
    }
    else{
       Qinv <-  1/one
    }

   tmp2 <- gridrec(gtinput,node,phinode,nvoice,
	freqstep,scale,Qinv,epsilon,np, real = real, check = check)

# 
# Now perform the irregular sampling 
#  

   tmp <- RidgeIrregSampling(phinode, node, tmp2$lam, prob)

   node <- tmp$node
   phinode <- tmp$phinode

   phi.x.min <- scale
   phi.x.max <- scale
   
   x.min <- node[1]
   x.max <- node[length(node)]
   x.max <- x.max + round(para * phi.x.max)
   x.min <- x.min - round(para * phi.x.min)

   node <- node - x.min + 1
   
   x.inc <- 1
   np <- as.integer((x.max-x.min)/x.inc)+1

   cat(" (np:",np,")")

   if(epsilon == 0)
       Q2 <- 0
   else {
      if (fast == F)
	Q2 <- gkernel(node,phinode,freqstep,scale,x.min = x.min, x.max = x.max)
      else
	Q2 <- fastgkernel(node,phinode,freqstep,scale,x.min = x.min,
		x.max = x.max)
      }
   cat(" kernel;")

# Generating the Q1 term in reconstruction kernel
   if (hflag == T)
      one <- gsampleOne(node,scale,np)
   else{
      one <- numeric(np)
      one[] <- 1
   }


   if (epsilon !=0 ){
      Q <-  epsilon * Q2
      for(j in 1:np)
         Q[j,j] <- Q[j,j] + one[j]
      Qinv <- solve(Q)
   }
   else{
      Qinv <-  1/one
    }

   tmp2 <- gridrec(gtinput,node,phinode,nvoice,
	freqstep,scale,Qinv,epsilon,np, real = real, check = check)


   if(plot == T){
      par(mfrow=c(2,1))
      plot.ts(Re(siginput))
      title("Original signal")
      plot.ts(Re(tmp2$sol))
      title("Reconstructed signal")
   }

   list(sol=tmp2$sol,A=tmp2$A,lam=tmp2$lam,dualwave=tmp2$dualwave,
	gaborets=tmp2$gaborets, solskel=tmp2$solskel,
	inputskel = tmp2$inputskel, Q2 = Q2)

}




RidgeIrregSampling <- function(phinode, node, lam, prob)
#########################################################################
#    RidgeSamplingData:
#    -----------------
#      Given a ridge phi (for the Gabor transform), and lam returns a 
#        (data-driven) subsampled version of length nbnodes.
#
#      Input:
#      ------
#      phinode: ridge at node
#      lam : vector of 2 * number of regular sampled node
#      prob : the percentage of lam to be preserved
#
#      Output:
#      -------
#      node: time coordinates of the ridge sampled by lambda
#      phinode: frequency coordinates of the ridge sampled by lambda
#
#########################################################################
{
   nbnodes <- length(node)   
   mlam <- numeric(nbnodes)
   for(j in 1:nbnodes) 
     mlam[j] <- Mod(lam[j]  + i *  lam[nbnodes + j])

   pct <- quantile(mlam,prob)

   count <- sum(mlam > pct)

   newnode <- numeric(count)
   newphinode <- numeric(count)

   count <- 0
   for(j in 1:nbnodes) {
     if(mlam[j] > pct) {
        newnode[count+1] <- node[j]
        newphinode[count+1] <- phinode[j]             
        count <- count + 1
     }
   }

   list(node = newnode,phinode = newphinode, nbnodes = count)
}





#########################################################################
#      $Log: gRidge_Recons.S,v $
# Revision 1.2  1995/04/05  18:56:55  bruno
# *** empty log message ***
#
# Revision 1.1  1995/04/02  01:04:16  bruno
# Initial revision
#
#
#               (c) Copyright  1997                             
#                          by                                   
#      Author: Rene Carmona, Bruno Torresani, Wen-Liang Hwang   
#                  Princeton University 
#                  All right reserved                           
#########################################################################





gregrec <- function(siginput, gtinput, phi, nbnodes, nvoice,
	freqstep, scale, epsilon = 0, fast = F, plot = F,
	para = 0, hflag = F, real = F, check = F)
#########################################################################
#     gregrec:
#     -------
#      Reconstruction of a real valued signal from a (continuous) 
#      Gabor ridge (uses a regular sampling of the ridge)
#
#      input:
#      -------
#      siginput: input signal
#      gtinput: Continuous gabor transform (output of cgt)
#      phi: (unsampled) ridge
#      nbnodes: number of nodes used for the reconstruction.
#      nvoice: number of different scales per octave
#      freqstep: sampling rate for the frequency axis
#      epsilon: coeff of the Q2 term in reconstruction kernel
#      fast: if set to TRUE, the kernel is computed using Riemann
#           sums instead of Romberg's quadrature
#      plot: if set to TRUE, displays original and reconstructed signals
#      para: scale parameter for extrapolating the ridges.
#      hflag:if set to T, uses $Q_1$ as first term in the kernel.
#      real: if set to TRUE, only uses constraints on the real part
#            of the transform for the reconstruction.
#      check: if set to TRUE, computes the wavelet transform of the
#            reconstructed signal
#
#      output:
#      -------
#      sol: reconstruction from a ridge
#      A: <gaborlets,dualgaborlets> matrix 
#      lam: coefficients of dual wavelets in reconstructed signal.
#      dualwave: array containing the dual wavelets.
#      gaborets: array containing the wavelets on sampled ridge.
#      solskel: Gabor transform of sol, restricted to the ridge
#      inputskel: Gabor transform of signal, restricted to the ridge
#      Q2: second part of the reconstruction kernel
#########################################################################
{
   tmp <- RidgeSampling(phi,nbnodes)
   node <- tmp$node
   phinode <- tmp$phinode

   phi.x.min <- scale
   phi.x.max <- scale
   
   x.min <- node[1]
   x.max <- node[length(node)]
   x.max <- x.max + round(para * phi.x.max)
   x.min <- x.min - round(para * phi.x.min)

   node <- node - x.min + 1
   
   x.inc <- 1
   np <- as.integer((x.max-x.min)/x.inc)+1

   cat(" (np:",np,")")

   if(epsilon == 0)
       Q2 <- 0
   else {
      if (fast == F)
	Q2 <- gkernel(node,phinode,freqstep,scale,x.min = x.min, x.max = x.max)
      else
	Q2 <- fastgkernel(node,phinode,freqstep,scale,x.min = x.min,
		x.max = x.max)
      }
   cat(" kernel;")

# Generating the Q1 term in reconstruction kernel
   if (hflag == T)
      one <- gsampleOne(node,scale,np)
   else{
      one <- numeric(np)
      one[] <- 1
   }


   if (epsilon !=0 ){
     Q <-  epsilon * Q2
       for(j in 1:np)
            Q[j,j] <- Q[j,j] + one[j]
         Qinv <- solve(Q)
      }
      else{
         Qinv <-  1/one
	 }

   tmp2 <- gridrec(gtinput,node,phinode,nvoice,
	freqstep,scale,Qinv,epsilon,np, real = real, check = check)

   if(plot == T){
      par(mfrow=c(2,1))
      plot.ts(Re(siginput))
      title("Original signal")
      plot.ts(Re(tmp2$sol))
      title("Reconstructed signal")
   }

npl(2)
lam <- tmp2$lam
plot.ts(lam,xlab="Number",ylab="Lambda Value")
title("Lambda Profile")
N <- length(lam)/2
mlam <- numeric(N)
for(j in 1:N)
   mlam[j] <- Mod(lam[j] + i * lam[N + j])
plot.ts(sort(mlam))

   list(sol=tmp2$sol,A=tmp2$A,lam=tmp2$lam,dualwave=tmp2$dualwave,
	gaborets=tmp2$gaborets, solskel=tmp2$solskel,
	inputskel = tmp2$inputskel, Q2 = Q2)

}




gridrec <- function(gtinput, node, phinode, nvoice,
	freqstep, scale, Qinv, epsilon, np, real = F, check = F)
#########################################################################
#     gridrec:
#     --------
#      Reconstruction of a real valued signal from a gabor ridge.
#
#      input:
#      ------
#      gtinput: Continuous gabor transform (output of cgt)
#      node: time coordinates of the ridge samples.
#      phinode: frequency coordinates of the ridge samples.
#      nvoice: number of different frequencies.
#      freqstep: sampling rate for the frequency axis.
#      scale: scale of the window.
#      Qinv: inverse of the reconstruction kernel
#      epsilon: coefficient of the Q2 term in the reconstruction kernel
#      np: number of samples of the reconstructed signal
#      real: if set to TRUE, only uses constraints on the real part
#            of the transform for the reconstruction.
#      check: if set to TRUE, computes the Gabor transform of the
#            reconstructed signal.
#
#      output:
#      -------
#      sol: reconstruction from a ridge
#      A: <gaborlets,dualgaborlets> matrix.
#      lam: coefficients of dual wavelets in reconstructed signal.
#      dualwave: array containing the dual gaborlets.
#      gaborets: array of gaborlets located on the ridge samples
#      solskel: Gabor transform of sol, restricted to the ridge
#      inputskel: Gabor transform of signal, restricted to the ridge
#
#########################################################################
{
   N <- length(node)


   omegaridge <- phinode
   bridge <- node
   if(real == T)
      gaborets <- gwave2(bridge,omegaridge,nvoice,freqstep,scale,np,N)
   else
      gaborets <- gwave(bridge,omegaridge,nvoice,freqstep,scale,np,N)

   cat("gaborets;")

   if(epsilon == 0){
      if(real == T)
         sk <- zeroskeleton2(gtinput,Qinv,gaborets,bridge,omegaridge,N)
      else
         sk <- zeroskeleton(gtinput,Qinv,gaborets,bridge,omegaridge,N)
      }

   else
      sk <- skeleton(gtinput,Qinv,gaborets,bridge,omegaridge,N)

   cat("skeleton.\n")

   solskel <- 0 #not needed if check not done
   inputskel <- 0 #not needed if check not done

   if(check == T){
      solskel <- complex(N)
      inputskel <- complex(N)
      gtsol <- cgt(sk$sol,nvoice,freqstep,plot=F)
      for(j in 1:N) solskel[j] <- gtsol[bridge[j],omegaridge[j]]

      for (j in 1:N)
         inputskel[j] <- gtinput[bridge[j],omegaridge[j]]
   }

   list(sol=sk$sol,A=sk$A,lam=sk$lam,dualwave=sk$dualwave,gaborets=gaborets,
	solskel=solskel,inputskel = inputskel)
}


gwave <- function(bridge,omegaridge,nvoice,freqstep,scale,np,N)
#########################################################################
#     gwave:
#     ------
#      Generation of the gaborets located on the ridge
#
#      input:
#      ------
#      bridge: time coordinates of the ridge samples
#      omegaridge: frequency coordinates of the ridge samples
#      nvoice: number of different scales per octave
#      freqstep: sampling rate for the frequency axis
#      scale: scale of the window
#      np: size of the reconstruction kernel
#      N: number of complex constraints
#
#      output:
#      -------
#      gaborets: array of morlet wavelets located on the ridge samples
#
#########################################################################
{

   gaborets <- matrix(0,np,2*N)

    omegaridge <- omegaridge * freqstep
    tmp <- vecgabor(np,N,bridge,omegaridge,scale)
    dim(tmp) <- c(np,N)
    gaborets[,1:N] <- Re(tmp[,1:N])
    gaborets[,(N+1):(2*N)] <- Im(tmp[,1:N])

# Alternative way of generating gaborets
#   for(k in 1:N)  {
#      omega <- omegaridge[k]*freqstep;
#      tmp <- gabor(np,bridge[k],omega,scale);
#      gaborets[,k] <- Re(tmp);
#      gaborets[,k+N] <- Im(tmp);
#    }


# Still another way of generating gaborets
#
#   dirac <- 1:np
#   dirac[] <- 0
#   for(k in 1:N)  {
#      dirac[bridge[k]] <- 1.0;
#      omega <- omegaridge[k]*freqstep;
#      gtdirac <- vgt(dirac,omega,scale,open=FALSE,plot=F);
#      gaborets[,k] <- Re(gtdirac);
#      gaborets[,k+N] <- Im(gtdirac);
#      dirac[] <- 0
#    }

    gaborets
}



gwave2 <- function(bridge,omegaridge,nvoice,freqstep,scale,np,N)
#########################################################################
#     gwave2:
#     -------
#      Generation of the real parts of gaborets located on the ridge
#
#      input:
#      ------
#      bridge: time coordinates of the ridge samples
#      omegaridge: frequency coordinates of the ridge samples
#      nvoice: number of different scales per octave
#      freqstep: sampling rate for the frequency axis
#      scale: scale of the window
#      np: size of the reconstruction kernel
#      N: number of complex constraints
#
#      output:
#      -------
#      gaborets: array of morlet wavelets located on the ridge samples
#
#########################################################################
{
    omegaridge <- omegaridge * freqstep
    gaborets <- Re(vecgabor(np,N,bridge,omegaridge,scale))
    dim(gaborets) <- c(np,N)


#   for(k in 1:N)  {
#      omega <- omegaridge[k]*freqstep;
#      tmp <- gabor(np,bridge[k],omega,scale);
#      gaborets[,k] <- Re(tmp);
#    }
    gaborets
}




gsampleOne <- function(node,scale,np)
#########################################################################
#     gsampleOne:
#     -----------
#
#       Generate a ``sampled identity"
#
#     input:
#     ------
#      node: location of the reconstruction gabor functions
#      scale: scale of the gabor functions
#      np: size of the reconstructed signal
#
#
#     output:
#     -------
#      dia: diagonal of the ``sampled'' Q1 term (1D vector)
#
#########################################################################
{
   dia <- numeric(np)
   N <- length(node)

   normalization <- sqrt(pi) * scale

   for(j in 1:np) {
     tmp1 <- (j-node)/scale
     tmp1 <- exp(-(tmp1 * tmp1))

     dia[j] <- sum(tmp1) #/normalization
    }
   dia

#   for(j in 1:np) {
#     tmp <- 0
#     for(k in 1:N) {
#       b <- node[k]
#       tmp1 <- (j-b)/scale
#       tmp1 <- exp(-(tmp1 * tmp1))
#       tmp <- tmp + tmp1
#     }
#     dia[j] <- tmp #/normalization
#   }
   dia
}






#########################################################################
#    (c) Copyright 1997
#               by
#   Author: Rene Carmona, Bruno Torresani, Wen L. Hwang, Andrea Wang
#              Princeton University 
#              All right reserved
########################################################################



mrecons <- function(extrema,filtername="Gaussian1",readflag= TRUE)
#****************************************************************
#  mrecons
#  --------
#    Reconstruct from dyadic extrema. The reconstructed signal
#    preserves locations and values at extrema. This is Carmona's
#    method of extrema reconstruction.	
# 
#  input
#  -----
#    extrema : the extrema representation 
#    filtername : filters (Gaussain1 stands filters corresponding to
#	Mallat and Zhong's wavelet. And, the Haar stands for the
#       filters for Haar basis.	
#    readflag: if set to TRUE, read kernel from a precomputed file.
#
#  output
#  ------
#    f: signal reconstructed from the wavelet transform extrema
#    g: f plus mean of f
#    h: f plus coarse scale of f
#****************************************************************
{
  pure <- extrema$original
  x <- adjust.length(pure)
  pure.data <- x$signal

  data <- extrema$original
  maxresoln <- extrema$maxresoln
  np <- extrema$np

  if (2^maxresoln > np)
    stop("The support of the coarsest scale filter is bigger than signal size")
  if ((readflag && (maxresoln > 9)) || (np > 4096))
    stop("Can not use readflag")

  dim(extrema$extrema) <- c(length(extrema$extrema),1)

  f <- 1:np

  z <- .C("extrema_reconst",
	as.character(filtername),
	a=as.single(f),	
	as.single(extrema$extrema),
	as.integer(maxresoln),
	as.integer(np),
        as.integer(readflag))

  f <- z$a 

  par(oma=c(2,0,0,0))

  g <- f + mean(pure.data)
  h <- f + extrema$Sf

  lim <- c(min(f, g, h, pure.data), max(f, g, h, pure.data))
  par(mfrow=c(3,1))

  plot.ts( f, pure.data, bty="n", ylim=lim )
  mtext("Reconstruction", cex=1.2, side=1, line=4 )

  plot.ts( g, pure.data, bty="n", ylim=lim )
  mtext("Reconstruction + Mean", cex=1.2, side=1, line=4 )

  plot.ts( h, pure.data, bty="n", ylim=lim )
  mtext("Reconstruction + Sf", cex=1.2, side=1, line=4 )

  par(mfrow=c(1,1))
  par(mar=c(5.1, 4.1, 4.1, 2.1))

  list( f= f, g= g, h = h )
}


dwinverse <- function(wt,filtername="Gaussian1")
#****************************************************************
#  dwinverse
#  ---------
#    Invert the dyadic wavelet transform.
# 
#  input
#  -----
#    wt: the wavelet transform 
#    filtername : filters used. ("Gaussian1" stands for the filters
#       corresponds to those of Mallat and Zhong's wavlet. And
#       "Haar" stands for the filters of Haar basis.	
#
#  output
#  ------
#    f: the reconstructed signal
#****************************************************************
{

  Sf <- wt$Sf
  Sf <- c(Sf)	
  dim(Sf) <- c(length(Sf),1)
  Wf <- wt$Wf
  Wf <- c(Wf)	
  dim(Wf) <- c(length(Wf),1)
  np <- wt$np
  maxresoln <- wt$maxresoln
  fback <- 1:np
  dim(fback) <- c(length(fback),1)

  z <- .C("Sinverse_wavelet_transform",
	a=as.single(fback),
        as.single(Sf),
        as.single(Wf),
	as.integer(maxresoln),
	as.integer(np), 
        as.character(filtername) )

  f <- z$a 
  f
}




#########################################################################
#    $log: mw.S,v $
#########################################################################
#
#    (c) Copyright 1997
#               by
#   Author: Rene Carmona, Bruno Torresani, Wen L. Hwang, Andrea Wang
#              Princeton University 
#              All right reserved
########################################################################




check.maxresoln <- function( maxresoln, np )
#*********************************************************************#
# check.maxresoln
# ---------------
# stop when the size of 2^maxresoln is no less than signal size
#
# input 
# -----
# maxresoln: number of decomposition
# np: signal size
#
# output
# ------
#*********************************************************************#
{
  if ( 2^(maxresoln+1) > np )
    stop("maxresoln is too large for the given signal")
}


mw <- function(inputdata,maxresoln,filtername="Gaussian1",scale=FALSE,
	       plot=TRUE)
#*********************************************************************#
# mw
# --
#   mw computes the wavelet decomposition by Mallat and Zhong's wavelet.
#
# input
# -----
#   inputdata: either a text file or an S object of input data.
#   maxresoln: number of decomposition
#   filename: name of filter (Gaussian1 stands for Mallat and Zhong's filter ;
#             Haar stands for the Haar basis).	
#   scale: when set, the wavelet transform at each scale will be plotted 
#          with the same scale.
#   plot : indicate if the wavelet transform at each scale will be plotted.
#          		
# output
# ------
#   original: original signal
#   Wf: wavelet transform of signal
#   Sf: signal at a lower resolution
#   maxresoln: number of decomposition
#   np: size of signal
#*********************************************************************#
{
  x <- adjust.length(inputdata)
  s <- x$signal
  np <- x$length

  Sf <- matrix(0, nrow=(maxresoln+1), ncol=np)
  Wf <- matrix(0, nrow=maxresoln, ncol=np)  

  # Convert a matrix maxresoln by np matrix into a vector
  Sf <- t(Sf)
  dim(Sf) <- c(length(Sf),1)
  Wf <- t(Wf)
  dim(Wf) <- c(length(Wf),1)

  y <- .C("Sf_compute", 
	Sf=as.single(Sf), 
	as.single(s), 
	as.integer(maxresoln), 
	as.integer(np),
	as.character(filtername))

  z <- .C("Wf_compute",
	Wf=as.single(Wf),
	as.single(y$Sf),
	as.integer(maxresoln),
	as.integer(np),
	as.character(filtername))

  # Convert the vectors into the original matrix
  Sf <- t(y$Sf)

  Sf <- Sf[(maxresoln*np+1): ((maxresoln+1)*np)]

  Wf <- t(z$Wf)
  dim(Wf) <- c(np, maxresoln)

  if(plot)
    plotwt(s, Wf, Sf, maxresoln, scale)

  list( original=s, Wf=Wf, Sf=Sf, maxresoln=maxresoln, np=np )
}








#########################################################################
#    (c) Copyright 1997
#               by
#   Author: Rene Carmona, Bruno Torresani, Wen L. Hwang, Andrea Wang
#              Princeton University
#              All right reserved
########################################################################


plotwt <- function(original, psi, phi, maxresoln, scale=F, yaxtype="s")
#*********************************************************************#
# plotwt
# ------
# Plot wavelet transform
# Output the original signal, wavelet transform Wf starting resoln 1,
# and Sf at the last resoln  
# Used by mw and dw
#*********************************************************************#
{
  par.orig <- par()  # Save the original plot parameters

  par(mfrow = c((maxresoln+2), 1))

  if ( !scale )
  {
    plot.ts(original,bty="n",yaxt=yaxtype)
    for (j in 1:maxresoln)
      plot.ts(psi[,j],bty="n",yaxt=yaxtype)
    plot.ts(phi, bty="n",yaxt=yaxtype)
  }
  else
  {
    limit <- c(min(original, psi, phi), max(original, psi, phi))
    plot.ts(original, ylim=limit, bty="n",yaxt=yaxtype)
    for (j in 1:maxresoln)
      plot.ts(psi[,j], ylim=limit, bty="n",yaxt=yaxtype)
    plot.ts(phi, ylim=limit, bty="n",yaxt=yaxtype)
  }

  par(par.orig)  
}


plot.result <- function(result, original, maxresoln, scale=F, yaxtype="s")
#*********************************************************************#
# Plot.result
# -----------
#    function for the output of the following S functions:
# 	dnpval, dbpval, region, dntrim
#	mnpval,	ext, mntrim, localvar
#*********************************************************************#

{
  par.orig <- par()  # Save the original plot parameters
  par(mar=c(2.0, 2.0, 0.5, 0.5))
  par(oma=c(1,0,0,0)) # Change 1 to a larger number to allow a larger
		      # bottom margin 

  par(mfrow = c(maxresoln+1, 1))

  if ( !scale )
  {
    plot.ts(original, bty="n",yaxt=yaxtype)
    for (j in 1:maxresoln)
      plot.ts(result[,j], bty="n",yaxt=yaxtype)
  }
  else
  {
    ymin <- min(result)
    ymax <- max(result)
    plot.ts(original, bty="n")
    for (j in 1:maxresoln)
      plot.ts(result[,j], ylim=c(ymin, ymax), bty="n",yaxt=yaxtype)
  }
 
  par(par.orig)
  par(mfrow=c(1,1))
}


wpl <- function(dwtrans)
#*********************************************************************#
# wpl
# ---
# Plot dyadic wavelet transform (output of mw).
#*********************************************************************#
{
  maxresoln <- dwtrans$maxresoln
  plotwt(dwtrans$original, dwtrans$Wf, dwtrans$Sf,maxresoln)
  cat("")
}


epl <- function(dwext)
#*********************************************************************#
# epl
# ---
# Plot wavelet transform extrema (output of ext).
#*********************************************************************#
{
  maxresoln <- dwext$maxresoln
  plot.result(dwext$extrema, dwext$original, maxresoln)
  cat("")
}




############################################################
# Functions processing radar images and their example. 
############################################################

######################################
# processing Radar images from V. Chen
######################################

#b727 <- matrix(scan("b727s.dat"),ncol=512,byrow=T)
#b727 <- matrix(scan("b727r.dat"),ncol=256,byrow=T)

#ncol <- dim(b727)[2]/2
#b727r<- b727[,-1+2*(1:ncol)]   
#b727i<- b727[,2*(1:ncol)]       

#B727 <- t(Conj(b727r + i*b727i))
#nrow <- dim(B727)[1]
#ncol <- dim(B727)[2]


#shift fft 1D data
#-----------------

fftshift <- function(x) {
    lng <- length(x)
    y <- numeric(lng)
    y[1:(lng/2)] <- x[(lng/2+1):(lng)]
    y[(lng/2+1):lng] <- x[1:(lng/2)]
    y
}


#* The following tow S-routine concludes our Radar experiments
#-------------------------------------------------------------
#b727 <- matrix(scan("b727s.dat"),ncol=512,byrow=T)
#b727 <- matrix(scan("b727r.dat"),ncol=256,byrow=T)

#Generating the original image
#-----------------------------

#e.g. B727 <- showRadar(b727)

showRadar<- function(x,plot=TRUE) {
   ncol <- dim(x)[2]/2
   xr<- x[,-1+2*(1:ncol)]   
   xi<- x[,2*(1:ncol)]       

   y <- t(Conj(xr + i*xi))

   oy <- t(apply(apply(y,2,fftshift),1,fftshift))
   oy <- apply(oy,2,fft)
   oy <- t(apply(apply(oy,2,fftshift),1,fftshift))
   if(plot) image(Mod(t(oy)))
   oy
}

#Enhancing Radar Image by Gabor Transform (V.Chen)
#------------------------------------------------


#gtime <- 128
#scale <- 50

#e.g. B727 <- cgtRadar(b727,gtime,scale)
#e.g. 
#     b727 <- matrix(scan("b727s.dat"),ncol=512,byrow=T)
#     b727 <- matrix(scan("b727r.dat"),ncol=256,byrow=T)
#     ncol <- dim(b727)[2]/2
#     b727r<- b727[,-1+2*(1:ncol)]   
#     b727i<- b727[,2*(1:ncol)]       

#     b727 <- t(Conj(b727r + i*b727i))
#     B727 <- cgtRadar(b727,gtime,scale,flag=F)
#     image(t(apply(Mod(B727$cgtout[,2,]),1,fftshift))) 

cgtRadar <- function(x,gtime,scale,flag=TRUE,plot=TRUE) {

  y <- x
  
  if(flag) {
    ncol <- dim(x)[2]/2
    xr<- x[,-1+2*(1:ncol)]   
    xi<- x[,2*(1:ncol)]       
    
    y <- t(Conj(xr + i*xi))
  }
  
  nrow <- dim(y)[1]
  ncol <- dim(y)[2]
  cgty <- array(0+0i,c(ncol,nrow,gtime))
  cgtyy <- array(0+0i,c(ncol,nrow,gtime))
  for(k in 1:ncol) cgty[k,,] <- cgt(Re(y[,k]),gtime,2/gtime,scale,plot=F) +
    i*(cgt(Im(y[,k]),gtime,2/gtime,scale,plot=F))
  for(k in 1:nrow) cgtyy[,k,] <- t(apply(cgty[,k,],1,fftshift))
  
  oy <- t(apply(apply(Mod(cgty),c(1,3),mean),1,fftshift))
  if(plot)  image(oy)
  
  list(output=oy,cgtout=cgtyy)
}

# Generate X, Y, Z for the command : cloud

cloudXYZ <- function(dim1, dim2,dim3) {
  II <- rep(1:dim1,dim2*dim3)
  dim(II) <- c(dim1,dim2,dim3)
  
  KK <- rep(1:dim3,dim2)
  dim(KK) <- c(dim3,dim2)
  KK <- t(KK)
  KK <- rep(KK,dim1)
  dim(KK) <- c(dim2 * dim3, dim1)
  KK <- t(KK)
  dim(KK) <- c(dim1, dim2, dim3)
  
  JJ <- rep(1:dim2,dim1)
  dim(JJ) <- c(dim2,dim1)
  JJ <- t(JJ)
  JJ <- rep(JJ,dim3)
  dim(JJ) <- c(dim1,dim2,dim3)
 
  list(II=II, JJ=JJ, KK=KK)
}

# Commands show points :

#cgtout <- B727$cgtout
#tt <-  cloudXYZ(dim(cgtout)[1], dim(cgtout)[2], dim(cgtout)[3])
#range(abs(cgtout))
#SS <- (abs(cgtout) > 4.0)
#XX <- tt$II[SS]
#YY <- tt$JJ[SS]
#ZZ <- tt$KK[SS]
#XX[length(XX)+1] <- 0
#XX[length(XX)+1] <- dim1
#YY[length(YY)+1] <- 0
#YY[length(YY)+1] <- dim2
#ZZ[length(ZZ)+1] <- 0
#ZZ[length(ZZ)+1] <- dim3
#cloud(ZZ~XX*YY)


#showpoints <- function(cgtout,low,high) {
# X <- 1:crossrange
# XX <- X %*% t(rep(1,crossrange))
# YY <- t(XX)
# freq <- dim(cgtout)[2]
# ZZ <- matrix(0,crossrange,crossrange)
# modmax <- max(Mod(cgtout))
# for(j in 1:crossrange)
#  for(k in 1:crossrange)
#    ZZ[j,k] <- maxpos(Mod(cgtout[j,,k]),low*modmax,high*modmax)
# cloud(ZZ~XX*YY)
#}


#maxpos <- function(x,low,high) {
#  pos <- 1       
#  lng <- length(x)
#  maxval <- max(x)
#  if((maxval >= low) & (maxval <= high)) { 
#   for(j in 1:lng) {
#     if(x[j] >= maxval) {
#       pos <- j
#     }
#   }
#  }
#  pos
#}



#showpoints <- function(cgtout,low,high) {
# crossrange <- dim(cgtout)[1]
# gabortime <- dim(cgtout)[3]
# freq <- dim(cgtout)[2]
# modmax <- max(Mod(cgtout))
# tt <-persp(matrix(1,crossrange,gabortime),zlim=1:freq)
# TT <- (Mod(cgtout) > low * modmax) * (Mod(cgtout) <= high * modmax)
# for(k in 1:gabortime) {
#    for(j in 1:freq) {
#      for(i in 1:crossrange) { 
#        if(TT[i,j,k]==T) 
#          points(perspp(i,k,j,tt))
#      }
#    }
#  }
# tt
#}



#shift fft 1D data
#-----------------

fftshift <- function(x) {
    lng <- length(x)
    y <- numeric(lng)
    y[1:(lng/2)] <- x[(lng/2+1):(lng)]
    y[(lng/2+1):lng] <- x[1:(lng/2)]
    y
}


#b727 <- matrix(scan("b727s.dat"),ncol=512,byrow=T)
#b727 <- matrix(scan("b727r.dat"),ncol=256,byrow=T)

#ncol <- dim(b727)[2]/2
#b727r<- b727[,-1+2*(1:ncol)]   
#b727i<- b727[,2*(1:ncol)]       

#B727 <- t(Conj(b727r + i*b727i))
#nrow <- dim(B727)[1]
#ncol <- dim(B727)[2]

#Generating the original image
#-----------------------------
#* apply fftshift twice to swap the center of the matrix


#o727 <- t(apply(apply(B727,2,fftshift),1,fftshift))
#o727 <- apply(o727,2,fft)
#o727 <- t(apply(apply(o727,2,fftshift),1,fftshift))



#*** only for display to compare with matlab output
#image(Mod(t(o727)))

#Applying Gabor Transform to ...
#-------------------------------
#gtime <- 100
#scale <- 50

#BB727 <- B727
#cgt727 <- array(0+0i,c(ncol,nrow,gtime))
#for(k in 1:ncol) cgt727[k,,] <-
#cgt(Re(BB727[,k]),gtime,2/gtime,scale,plot=F)+i*(cgt(Im(BB727[,k]),gtime,2/gtime,scale,plot=F))

#image(apply(apply(Mod(cgt727),c(1,3),mean),1,fftshift))


#* The following tow S-routine concludes our Radar experiments
#-------------------------------------------------------------
#b727 <- matrix(scan("b727s.dat"),ncol=512,byrow=T)
#b727 <- matrix(scan("b727r.dat"),ncol=256,byrow=T)

#Generating the original image
#-----------------------------

#e.g. B727 <- showRadar(b727)

showRadar<- function(x) {
   ncol <- dim(x)[2]/2
   xr<- x[,-1+2*(1:ncol)]   
   xi<- x[,2*(1:ncol)]       

   y <- t(Conj(xr + i*xi))

   oy <- t(apply(apply(y,2,fftshift),1,fftshift))
   oy <- apply(oy,2,fft)
   oy <- t(apply(apply(oy,2,fftshift),1,fftshift))
   image(Mod(t(oy)))
   oy
}

#Enhancing Radar Image by Gabor Transform
#----------------------------------------


#gtime <- 128
#scale <- 50

#e.g. B727 <- cgtRadar(b727,gtime,scale)
#e.g. 
#     ncol <- dim(b727)[2]/2
#     b727r<- b727[,-1+2*(1:ncol)]   
#     b727i<- b727[,2*(1:ncol)]       

#     b727 <- t(Conj(b727r + i*b727i))
#     B727 <- cgtRadar(b727,gtime,scale,flag=F)

cgtRadar <- function(x,gtime,scale,flag=TRUE) {

   y <- x

   if(flag) {
     ncol <- dim(x)[2]/2
     xr<- x[,-1+2*(1:ncol)]   
     xi<- x[,2*(1:ncol)]       

     y <- t(Conj(xr + i*xi))
   }

   nrow <- dim(y)[1]
   ncol <- dim(y)[2]
   cgty <- array(0+0i,c(ncol,nrow,gtime))
   for(k in 1:ncol) cgty[k,,] <- cgt(Re(y[,k]),gtime,2/gtime,scale,plot=F)+i*(cgt(Im(y[,k]),gtime,2/gtime,scale,plot=F))

   oy <- apply(apply(Mod(cgty),c(1,3),mean),1,fftshift)
   for(k in 1:nrow) cgty[,k,] <- apply(cgty[,k,],1,fftshift)
   image(oy)
   list(output=oy,cgtout=cgty)
}








# The commands perform the test data of the robust of the
# reconstructions from ridges in wavelet transform
# The command used to generate run of testing
# 128 run for each db in c(25,20,10,5,0)
# chbat signal with 5 noct,10 voice, 5 compr 
#tt <- RunRec(chbat,5,10,5,128,c(25,20,10,5,0))

#plot.ts(HOWAREYOU)
#cgtHOWAREYOU <- cgt(HOWAREYOU,70,0.01,100)
#clHOWAREYOU <- crc(Mod(cgtHOWAREYOU),nbclimb=1000)
#cfHOWAREYOU <- cfamily(clHOWAREYOU,ptile=0.001)
#image(cfHOWAREYOU$ordered > 0) 
#HOWord <- cfHOWAREYOU$ordered > 0 

robustrec <- function(x,nvoice,freqstep,scale,db,ptile) {

   lng <- length(x)
   nx <- x
   wn <- rnorm(lng,0,sqrt(var(nx)/InverseDB(db))) 
   nx <- nx + wn

   cgtnx <- cgt(nx,nvoice,freqstep,scale,plot=F)
   crcnx <- crc(Mod(cgtnx), nbclimb = 1000)
   cfnx <- cfamily(crcnx,,ptile=ptile)
   nxordered <- (cfnx$ordered > 0)
     
   nxordered
}

# db is a vector of DB
# number of experiments (nrun) for each db (ndb) and each radius (nr)
# vdb <- c(20,15,10,5,0,-5)
# pvec <- c(0.001,0.001,0.001,0.01,0.01,0.02)
# pr <- c(0,1)
# tt <- RunRec(HOWord,HOWAREYOU,70,0.01,100,1,vdb,pvec,pr)
# MM <- apply(tt,c(1,2),mean)

RunRec <- function(HOWord,x,nvoice,freqstep,scale,nrun,vdb,pvec,pr) {
  ndb <- length(vdb)
  nr <- length(pr)
  Run <- array(0,c(ndb,nr,nrun))

  for(ii in 1:nrun) {
cat("The number of run:",ii,"\n")

     for(jj in 1:ndb) {
        db <- vdb[jj]
        ptile <- pvec[jj]
        tt <- robustrec(x,nvoice,freqstep,scale,db,ptile)
        for(kk in 1:nr) {
          rr <- pr[kk]
          Run[jj,kk,ii] <- RidgeDist(tt,HOWord,rr)
        }
      }
  }
  Run
}

#wn <- rnorm(2048,0,sqrt(var(CLICKS)/InverseDB(1)))

#10*log10(var(CLICKS)/var(wn))

#plot.ts(HOWAREYOU)
#cgtHOWAREYOU <- cgt(HOWAREYOU,70,0.01,100)
#clHOWAREYOU <- crc(Mod(cgtHOWAREYOU),nbclimb=1000)
#cfHOWAREYOU <- cfamily(clHOWAREYOU,ptile=0.001)
#image(cfHOWAREYOU$ordered > 0) 
#HOWord <- cfHOWAREYOU$ordered > 0 




band <- function(x,r) {

xx <- x

lng <- length(xx)
yy <- logical(lng)
x.rts <- as.rts(xx)
for(k in 1:r) {
  y.rts <- lag(x.rts,-k) 
  zz <- as.logical(y.rts)
  yy[(1+(k)):lng] <- zz[1:(lng-(k))]
  xx <- yy | xx
  y.rts <- lag(x.rts,k)
  zz <- as.logical(y.rts)
  yy[1:(lng-(k))] <- zz[(1+(k)):lng]
  yy[(lng-(k)+1):lng] <- F
  xx <- yy | xx
}
xx
}

Sausage <- function(A,r)
{
 if( abs(r) > 0 ) {
   r <- abs(r)
   B <- band(c(A),r)
   dim(B) <- c(dim(A)[1],dim(A)[2])

   C <- band(c(t(B)),r)
   #C <- band(c(t(A)),r)

   dim(C) <- c(dim(t(A))[1],dim(t(A))[2])
   B <- B | t(C)
 }
 else {
  B <- A
 }
 B
}

# Measure the distance of the two set M1 and M2: [0,1]
# (#(M1 - (M1&Sausage(M2))) +  #(M2 - (M2&Sausage(M1))))/(#M1+#M2)

RidgeDist <- function(m1,m2,r) {
  ms1 <- Sausage(m1,r)
  ms2 <- Sausage(m2,r)
  m3 <- m1 & ms2
  m4 <- m1 & (!m3)
  m5 <- m2 & ms1
  m6 <- m2 & (!m5)
  d <- (sum(m4) + sum(m6))/(sum(m1)+sum(m2))
  d
}

# confidence level:

low <- 0.05
high <- 0.95

confident <- function(x,low,high) {
 ndb <- dim(x)[1]
 nr <- dim(x)[2]
 nrun <- dim(x)[3]
 for(jj in 1:ndb) {
   for(kk in 1:nr) {
      tmp <- x[jj,kk,]
      ttmp <- quantile(c(tmp),c(low,high))
cat("<dB r> = ",jj,kk,"<low high>=",ttmp,"\n")
   }
  }
}
 


#**************************************************************#
#              (c) Copyright  1997                             #
#                         by                                   #
#     Author: Rene Carmona, Bruno Torresani, Wen L. Hwang      #
#                 Princeton University                         #
#                 All right reserved                           #
#**************************************************************#


# compute the svd (single value decomposition)
# The Input/Output of the command is the same as the
# commands svd in Splus

SVD <- function(a)  
{
  m <- dim(a)[1]
  n <- dim(a)[2]

# diagonal elements
  w <- numeric(n)
  dim(w) <- c(n,1)	 

  v <- matrix(0,n,n)
  v <- c(v)
  dim(v) <- c(length(v),1)
  a <- c(a)
  dim(a) <- c(length(a),1)

  z <- .C("Ssvdecomp",
	u= as.double(a),
	as.integer(m),
	as.integer(n),
        w = as.double(w),
        v = as.double(v))

  u <- z$u
  dim(u) <- c(n,n)	
  w <- z$w
  dim(w) <- c(n,1)
  v <- z$v
  dim(v) <- c(n,n)
	
  list(d=w,u=u,v=v)
}





.First.lib <- function(lib, pkg) library.dynam("Rwave",pkg,lib)
