signal_functions.c

LastWave

      break;
      
    case BorderPeriodic :
      if (firstPoint < 0) f = ((-firstPoint)/sig->size+1)*sig->size+firstPoint;
      else f = firstPoint;
      for (i= f,j=0;j<newSize;i++,j++) {
        sigOut1->Y[j] = sig->Y[i%sig->size];
      }
      break;

    case BorderMirror :
      if (firstPoint < 0) f = 2*((-firstPoint)/sig->size+1)*sig->size+firstPoint;
      else f = firstPoint;
      for (i= f,j=0;j<newSize;i++,j++) {
        if ((i/sig->size)%2 == 0) sigOut1->Y[j] = sig->Y[i%sig->size];
        else sigOut1->Y[j] = sig->Y[sig->size-1-i%sig->size];
      }
      break;
          
    default: Errorf("ExtractSig() : Bad Border type");
  }
  }
      
  
  /* Then set the sigOut */
  if (sigOut1 != sigOut) {
    CopySig(sigOut1,sigOut);
    DeleteSignal(sigOut1);
  }
  
}


/*********************/
/* Padd a signal     */
/*********************/

void PaddSig(SIGNAL sig,SIGNAL sigOut, int borderType,int newSize)
{
  SIGNAL sig1;
  int i,iMin,iMax,n,j;
  int sizen,size;
  
  /* Compute the new size */
  size = sig->size;
  if (newSize == -1) {
    if (size != 1) newSize = 1 << (1 + (int) (log((double) size-1)/log(2.))); 
    else newSize = 1;
  }
     
  /* If the same size don't do anything */
  if (newSize == size) return;
  
  /* Let's fill the new signal with the old one */
  sig1 = NewSignal();
  SizeSignal(sig1,newSize,YSIG);
  sig1->dx = sig->dx;
  iMin = (newSize-size)/2;
  iMax = iMin+size-1;
  for(i=iMin;i<=iMax;i++) sig1->Y[i] = sig->Y[i-iMin];
  sig1->x0 = sig->x0-iMin*sig->dx;
  
  /* Let's fill everything else */
  switch (borderType) {
    case BorderPad : 
      for(i=0;i<iMin;i++) sig1->Y[i] = sig->Y[0];
      for(i=iMax+1;i<newSize;i++) sig1->Y[i] = sig->Y[size-1];
      break;

    case BorderPad0 : 
      for(i=0;i<iMin;i++) sig1->Y[i] = 0;
      for(i=iMax+1;i<newSize;i++) sig1->Y[i] = 0;
      break;
      
    case BorderPeriodic : 
      n = newSize/size+1;
      sizen = n*size;
      for(i=0;i<iMin;i++) sig1->Y[i] = sig->Y[(i-iMin+sizen)%size];
      for(i=iMax+1;i<newSize;i++) sig1->Y[i] = sig->Y[(i-iMin+sizen)%size];
      break;

    case BorderMirror : 
      n = newSize/(2*size) +1;
      sizen = n*2*size;
      for(i=0;i<iMin;i++) {
        j = (i-iMin+sizen)%(2*size);
        if (j<size) sig1->Y[i] = sig->Y[j];
        else sig1->Y[i] = sig->Y[2*size-j-1];
      }
      for(i=iMax+1;i<newSize;i++) {
        j = (i-iMin+sizen)%(2*size);
        if (j<size) sig1->Y[i] = sig->Y[j];
        else sig1->Y[i] = sig->Y[2*size-j-1];
      }
      break;
    
    default: Errorf("PaddSig() : Bad Border type");
  }
      
  
  CopySig(sig1,sigOut);	
  
  DeleteSignal(sig1);
}

/* The corresponding command */
void C_Padd(char **argv)
{
  SIGNAL signal,signalOut;
  char *borderName;
  int borderType;
  int size;
  char opt;
  
  argv = ParseArgv(argv,tSIGNALI,&signal,tSIGNAL_,NULL,&signalOut,-1);
  if (signalOut == NULL) signalOut = signal;
  
  if (*argv != NULL || **argv != '-') argv = ParseArgv(argv,tWORD_,"*bconst",&borderName,-1);
  
  if(!strcmp(borderName,"pad") || !strcmp(borderName,"*bconst")) borderType = BorderPad;
  else if(!strcmp(borderName,"pad0") || !strcmp(borderName,"*b0")) borderType = BorderPad0;
  else if(!strcmp(borderName,"mir") || !strcmp(borderName,"*bmirror")) borderType = BorderMirror;
  else if(!strcmp(borderName,"per") || !strcmp(borderName,"*bperiodic")) borderType = BorderPeriodic;
  else Errorf("Undefined border effect: %s",borderName);

  size = -1;
  while(opt = ParseOption(&argv))
    {
      switch(opt) 
        {
       case 's': /*final size*/
	      argv = ParseArgv(argv,tINT,&size,-1);
	      if (size<=signal->size) Errorf("Size is not big enough!");
	      break;      
        default:
          ErrorOption(opt);
        }
    }
  NoMoreArgs(argv);
 
  PaddSig(signal,signalOut,borderType,size);
}

/**************************************/
/*    Compute statistical values      */
/**************************************/

/* Get moment of order 'n' (and returns the mean) */
LWFLOAT GetNthMoment(SIGNAL signal, int n, LWFLOAT *pNthMoment,int flagCausal, int flagCentered)
{
  int i,j;
  int iMin,iMax;
  LWFLOAT m,mn,s,f;
  
  if (n< 0) Errorf("GetNthMoment() : 'n' should be positive");
  
  if (n==0) {
    *pNthMoment = 1;
    return(1);
  }
  
  /* Border effects */
  if (flagCausal == YES) {
    iMin = signal->firstp;
    iMax = signal->lastp;
  }
  else {
    iMin = 0;
    iMax = signal->size-1;
  }
    
  /* Compute the mean */
  m = 0;
  if (flagCentered || n == 1) {
    for (i = iMin;i<=iMax;i++) m += signal->Y[i];
    m /= (iMax-iMin+1);
  }
  
  /* If n == 1 we are done ! */
  if (n == 1) {
    *pNthMoment = m;
    return(m);
  }
  
  /* Compute the Nth moment */
  mn = 0;
  for (i = iMin;i<=iMax;i++) {
    f = signal->Y[i]-m;
    s = 1;
    for (j=0;j<n;j++) s *= f;
    mn += s;
  }
  mn /= (iMax-iMin+1);

  *pNthMoment = mn;
  return(m);
}

/* Get absolute moment of order 'f1' (and returns the mean) */
LWFLOAT GetAbsMoment(SIGNAL signal, LWFLOAT f1, LWFLOAT *pMoment,int flagCausal, int flagCentered)
{
  int i;
  int iMin,iMax;
  LWFLOAT m,mn;
  
  if (f1==0) {
    *pMoment = 1;
    return(1);
  }
  
  /* Border effects */
  if (flagCausal == YES) {
    iMin = signal->firstp;
    iMax = signal->lastp;
  }
  else {
    iMin = 0;
    iMax = signal->size-1;
  }
    
  /* Compute the mean */
  m = 0;
  if (flagCentered) {
    for (i = iMin;i<=iMax;i++) m += signal->Y[i];
    m /= (iMax-iMin+1);
  }
  
  /* Compute the Nth moment */
  mn = 0;
  for (i = iMin;i<=iMax;i++) {
    mn += pow(fabs(signal->Y[i]-m),f1);
  }
  mn /= (iMax-iMin+1);

  *pMoment = mn;
  return(m);
}

/* Correlation (if signal2 == NULL the correlation is computed between X and Y) */
LWFLOAT GetCorrelation(SIGNAL signal1,SIGNAL signal2,int flagCausal)
{  
  int i;
  LWFLOAT y,y2,xy,x,x2;
  int iMin,iMax;
  SIGNAL s1,s2;
  char type1,type2;
  
  /* Border effects */
  if (flagCausal == YES) {
    iMin = signal1->firstp;
    iMax = signal1->lastp;
  }
  else {
    iMin = 0;
    iMax = signal1->size-1;
  }
  
  /* The values */
  if (signal2 == NULL) {
    s1 = s2 = signal1;
    type1 = 'X';
    type2 = 'Y';
  }
  else {
    s1 = signal1;
    s2 = signal2;
    type1 = type2 = 'Y';
  }
  
  for (i = iMin,x=y=xy=x2=y2=0.;i<=iMax;i++) {
    xy += XYSig(s1,i,type1)*XYSig(s2,i,type2);
    x += XYSig(s1,i,type1);
    x2 += (XYSig(s1,i,type1)*XYSig(s1,i,type1));
    y += XYSig(s2,i,type2);
    y2 += (XYSig(s2,i,type2)*XYSig(s2,i,type2));
  }
  xy /= (iMax-iMin+1);
  x /= (iMax-iMin+1);
  x2 /= (iMax-iMin+1);
  y /= (iMax-iMin+1);
  y2 /= (iMax-iMin+1);
  
  return((xy-x*y)*(xy-x*y)/((y2-y*y)*(x2-x*x)));
}

/****************************************/
/* Fit a signal with a straight line    */
/****************************************/

/*
 * signal       :   the signal to fit 
 * pA,pB        :   the equation line is y = (*pA)x+(*pB)
 * pSigA,pSigB  :   variance on (*pA) and (*pB)
 * iMin,iMax    :   the fit should be performed between iMin and iMax
 */
 
void LineFitSig(SIGNAL signal,LWFLOAT *pA,LWFLOAT *pSigA,LWFLOAT *pB,LWFLOAT *pSigB,int iMin,int iMax) 
{
  int i;
  LWFLOAT chi2,tmp;
  LWFLOAT t,sxoss,sx=0.0,sy=0.0,st2=0.0,ss,sigdat;

  iMin = MAX(0,iMin);
  iMax = MIN(signal->size-1,iMax);
  if (iMax <= 0 || iMin >= iMax) {
    iMin = 0;
    iMax = signal->size-1;
  }
  
  *pA = 0.0;
  for (i=iMin;i<iMax+1;i++) {
       sx += XSig(signal,i);
       sy += signal->Y[i];
  }
  ss = iMax-iMin+1;
  sxoss = sx/ss;
  for (i=iMin;i<iMax+1;i++){
       t = XSig(signal,i) - sxoss;
       st2 += t*t;
       *pA += t*signal->Y[i];
  }
  *pA /= st2;
  *pB = (sy-sx*(*pA))/ss;
  *pSigB = sqrt((1.0 + sx*sx/(ss*st2))/ss);
  *pSigA =  sqrt(1.0/st2);
  chi2 = 0.0;
  for (i=iMin;i<iMax+1;i++){
   tmp = signal->Y[i]-(*pB)-(*pA)*XSig(signal,i);
   chi2 += tmp*tmp;
  }
  sigdat = sqrt(chi2/(iMax-iMin-1));
  *pSigB *= sigdat;
  *pSigA *= sigdat;
}

/* Get Lp Norm (using the dx field) */
LWFLOAT GetLpNormSig(SIGNAL signal, LWFLOAT p,int flagCausal)
{
  int i;
  int iMin,iMax;
  LWFLOAT lp;
  
  /* Border effects */
  if (flagCausal == YES) {
    iMin = signal->firstp;
    iMax = signal->lastp;
  }
  else {
    iMin = 0;
    iMax = signal->size-1;
  }
     
  /* Compute the Lp Norm */
  lp = 0;
  for (i = iMin;i<=iMax;i++) {
    lp += pow(fabs(signal->Y[i]),p);
  }
  lp *= signal->dx;
  lp = pow(lp,1.0/p);

  return(lp);
}


/*
 * Command for computing some statistical values of a signal 
 */
 
void C_Stats(char **argv)
{
  char *action;
  SIGNAL signal,signal1;
  LWFLOAT f,m,f1;
  char opt;
  LWFLOAT xMin,xMax,yMin,yMax,a,b,sigA,sigB,p;
  int iMin,iMax;
  int flagCausal,flagAbs,flagCentered;
  LISTV lv;
    
  argv = ParseArgv(argv,tWORD,&action,-1);

  /* 'mean' action */
  if (!strcmp(action,"mean")) {
    argv = ParseArgv(argv,tSIGNALI,&signal,-1);
    flagCausal = NO;
    while(opt = ParseOption(&argv)) { 
      switch(opt) {
        case 'c': flagCausal = YES; break;
        default: ErrorOption(opt);
      }
    }    
    NoMoreArgs(argv);
    GetNthMoment(signal,1,&f,flagCausal,NO);
    SetResultFloat(f);
  }

  /* 'var' action */
  else if (!strcmp(action,"var")) {
    argv = ParseArgv(argv,tSIGNALI,&signal,-1);
    flagCausal = NO;
    while(opt = ParseOption(&argv)) { 
      switch(opt) {
        case 'c': flagCausal = YES; break;
        default: ErrorOption(opt);
      }
    }    
    NoMoreArgs(argv);
    GetNthMoment(signal,2,&f,flagCausal,YES);
    SetResultFloat(f);
  }

  /* 'skew' action */
  else if (!strcmp(action,"skew")) {
    argv = ParseArgv(argv,tSIGNALI,&signal,-1);
    flagCausal = NO;
    while(opt = ParseOption(&argv)) { 
      switch(opt) {
        case 'c': flagCausal = YES; break;
        default: ErrorOption(opt);
      }
    }    
    NoMoreArgs(argv);
    GetNthMoment(signal,3,&f,flagCausal,YES);
    SetResultFloat(f);
  }

  /* 'kurt' action */
  else if (!strcmp(action,"kurt")) {
    argv = ParseArgv(argv,tSIGNALI,&signal,-1);
    flagCausal = NO;
    while(opt = ParseOption(&argv)) { 
      switch(opt) {
        case 'c': flagCausal = YES; break;
        default: ErrorOption(opt);
      }
    }    
    NoMoreArgs(argv);
    GetNthMoment(signal,4,&f,flagCausal,YES);
    GetNthMoment(signal,2,&f1,flagCausal,YES);
    SetResultFloat(f/(f1*f1));
  }
    
  /* 'nth' action */
  else if (!strcmp(action,"nth")) {
    argv = ParseArgv(argv,tSIGNALI,&signal,tFLOAT,&f1,-1);
    flagCausal = NO;
    flagAbs = NO;