signal_functions.c

LastWave

    flagCentered = NO;
    while(opt = ParseOption(&argv)) { 
      switch(opt) {
        case 'a': flagAbs = YES; break;
        case 'c': flagCausal = YES; break;
        case 'C': flagCentered = YES; break;
        default: ErrorOption(opt);
      }
    }    
    NoMoreArgs(argv);
   if (!flagAbs) {
     if (f1 != (int) f1 || f1 < 0) Errorf("<n> should be a positive integer unless you set option '-a'");
     m = GetNthMoment(signal,(int) f1,&f,flagCausal,flagCentered);
   }
   else m = GetAbsMoment(signal,f1,&f,flagCausal,flagCentered);
   SetResultFloat(f);
  }

  /* 'minmax' action */
  else if (!strcmp(action,"minmax")) {
    argv = ParseArgv(argv,tSIGNALI,&signal,-1);
    flagCausal = NO;
    while(opt = ParseOption(&argv)) { 
      switch(opt) {
        case 'c': flagCausal = YES; break;
        default: ErrorOption(opt);
      }
    }    
    NoMoreArgs(argv);
    xMin = 1;
    xMax = -1;
    MinMaxSig(signal,&xMin,&xMax,&yMin,&yMax,&iMin,&iMax,flagCausal);
    lv = TNewListv();
    AppendInt2Listv(lv,iMin);
    AppendInt2Listv(lv,iMax);
    SetResultValue(lv);
  }

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

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

  /* 'fit' action */
  else if (!strcmp(action,"fit")) {
    argv = ParseArgv(argv,tSIGNALI,&signal,-1);
    xMax = -1;
    xMin = 1;
    while(opt = ParseOption(&argv)) { 
      switch(opt) {
        case 'x': argv = ParseArgv(argv,tFLOAT,&xMin,tFLOAT,&xMax,-1); break;
        default: ErrorOption(opt);
      }
    }    
    NoMoreArgs(argv);  
    if (xMin >= xMax) {iMin = -1; iMax = -1;}
    else {
      iMin = ISig(signal,xMin);
      iMax = ISig(signal,xMax);
    }
    
    LineFitSig(signal,&a,&sigA,&b,&sigB,iMin,iMax);

   lv = TNewListv();
   AppendFloat2Listv(lv,a);
   AppendFloat2Listv(lv,sigA);
   AppendFloat2Listv(lv,b);
   AppendFloat2Listv(lv,sigB);
   AppendFloat2Listv(lv,iMin);
   AppendFloat2Listv(lv,iMax);
   SetResultValue(lv);
  }
  
  /* 'print' action */
  else if (!strcmp(action,"print")) {
    argv = ParseArgv(argv,tSIGNALI,&signal,-1);
    flagCausal = NO;
    while(opt = ParseOption(&argv)) { 
      switch(opt) {
        case 'c': flagCausal = YES; break;
        default: ErrorOption(opt);
      }
    }    
    NoMoreArgs(argv);
    m = GetLpNormSig(signal,2,flagCausal);
    Printf("L2 norm   : %.8g\n",m);
    m = GetNthMoment(signal,2,&f,flagCausal,YES);
    Printf("Mean      : %.8g\n",m);
    Printf("Variance  : %.8g\n",f);
    m = GetNthMoment(signal,3,&f,flagCausal,YES);
    Printf("Skewness  : %.8g\n",f);
    m = GetNthMoment(signal,4,&f,flagCausal,YES);
    m = GetNthMoment(signal,2,&f1,flagCausal,YES);
    Printf("Kurtosis  : %.8g\n",f-3*f1*f1);
    xMin = 1;
    xMax = -1;
    MinMaxSig(signal,&xMin,&xMax,&yMin,&yMax,&iMin,&iMax,flagCausal);
    Printf("Minimum at x=%.8g (index = %d) is y = %.8g\n",XSig(signal,iMin),iMin,yMin); 
    Printf("Maximum at x=%.8g (index = %d) is y = %.8g\n",XSig(signal,iMax),iMax,yMax); 
  }

  else Errorf("Unknow action '%s'",action);    	
}


 
/**************************************/
/* Compute the 'n' branches histogram */
/* of the signal Y's.                 */
/**************************************/

/*
 * input    : the signal to make the histogram on 
 * output   : the signal to put the histogram into
 * n        : the number of branches of the histogram
 * ymin,ymax  : the yMin and yMax of the histogram (not used if min > max)
 * xmin,xmax  : the xMin and xMax between which the histogram is made
 * weight   : signal of weights (or NULL). It must have the same size as input 
 */
 
void HistoSig(SIGNAL input, SIGNAL output, int n,LWFLOAT xmin, LWFLOAT xmax,LWFLOAT ymin, LWFLOAT ymax, SIGNAL weight,int flagCausal)
{
  int i,index;
  LWFLOAT min,max;
  int iMin,iMax;
  LWFLOAT theWeight;  

  if (output == input) 
    Errorf("HistoSig() : the output signal must be different from the input signal");

  SizeSignal(output,n,YSIG); 
  ZeroSig(output);
  
  if (ymin >= ymax)
    if( xmin >= xmax) 
      MinMaxSig(input,&xmin,&xmax,&ymin,&ymax,&iMin,&iMax,flagCausal);
    else 
      MinMaxSig(input,&min,&max,&ymin,&ymax,&iMin,&iMax,flagCausal);
  else
    if( xmin >= xmax) 
      MinMaxSig(input,&xmin,&xmax,&min,&max,&iMin,&iMax,flagCausal);
 
  if (ymin == ymax) {ymin -= .5;ymax += .5;}

  if(flagCausal == YES)
    {
      iMin = input->firstp;
      iMax = input->lastp;
    }
  else
    {
      iMin = 0;
      iMax = input->size - 1;
    }

  for (i=iMin;i<=iMax;i++)
    {
      if(XSig(input,i)>xmax || XSig(input,i) < xmin) 	continue;
      index = (int)((input->Y[i]-ymin)*n/(ymax-ymin));
      if (index < 0 || index > n) continue;
      if (index == n) index = n-1;
      
      if (weight == NULL) theWeight = 1;
      else theWeight = weight->Y[i];
      
      output->Y[index] += theWeight;
    }
  
  output->x0 = ymin + (ymax-ymin)/(2*n);
  output->dx = (ymax-ymin)/n;
}


/* The corresponding command */
void C_Histo (char **argv)
{
  SIGNAL input,output;
  int n,flagCausal=NO;
  LWFLOAT xmin,xmax;
  LWFLOAT ymin,ymax;
  SIGNAL weight;
  char opt;
  
  argv = ParseArgv(argv,tSIGNALI,&input,tSIGNAL,&output,tINT,&n,-1);

  xmin = 1;
  xmax = -1;
  ymin = 1;
  ymax = -1;
  weight = NULL;
  while(opt = ParseOption(&argv)) { 

    switch(opt) 
      {
      case 'x':
	argv = ParseArgv(argv,tFLOAT,&xmin,tFLOAT,&xmax,-1);
	if (xmin >= xmax) Errorf("xmin should be smaller than xmax");
	break;
      case 'y':
	argv = ParseArgv(argv,tFLOAT,&ymin,tFLOAT,&ymax,-1);
	if (ymin >= ymax) Errorf("ymin should be smaller than ymax");
	break;
      case 'w':
	argv = ParseArgv(argv,tSIGNAL,&weight,-1);
	if (weight->size != input->size) 
	  Errorf("The weight signal should be of the same size as the input signal");
	break;
      case 'c':
	flagCausal=YES;
	break;
      default:
	ErrorOption(opt);
      }
  }    
  NoMoreArgs(argv);
  
  if (input == output)
    Errorf("Input and Output must be different\n");
  
  HistoSig(input, output, n,xmin,xmax,ymin,ymax,weight,flagCausal);  
}



/**********************************************
 *
 * Convolution of a signal with a filter 
 *
 **********************************************/

void ConvSig(SIGNAL in, SIGNAL filter, SIGNAL out, int borderType,int method)
{
  int firstp, lastp;
  int filterOrigin;
  
  if (out == in) Errorf("ConvSig() : The input signal and the output signal must be different");
  if (out == filter) Errorf("ConvSig() : The filter signal and the output signal must be different");

  SizeSignal(out,in->size,YSIG);
  filterOrigin = (int) ((0-filter->x0)/filter->dx+.5);
  
  if (filterOrigin < 0 || filterOrigin > filter->size-1)
    Errorf("ConvSig() : The origin of the filter should be represented in the filter signal"); 

  cv_sig_init(CV_RC_FORM,in->Y,NULL,in->size);
  cv_flt_init_n(CV_RC_FORM,filter->size,filterOrigin,0,0,filter->Y,NULL);
  cv_set_method(method);
  cv_compute(borderType,out->Y,&firstp,&lastp);

  out->firstp = in->firstp+firstp;
  out->lastp = in->lastp - (in->size-1-lastp);
  
  out->dx = in->dx;
  out->x0 = in->x0;
}
  

void C_OldConv(char **argv)
{
  SIGNAL in,out,filter;
  char *borderName;
  int borderType;
  char opt;
  LWFLOAT time;
  int method = CV_UNDEFINED;

  argv = ParseArgv(argv,tSIGNALI,&in,tSIGNALI,&filter,tSIGNAL,&out,-1);
  
  borderName = "pad";
  if (*argv!=NULL && **argv != '-') argv = ParseArgv(argv,tSTR_,"pad",&borderName,-1);

  if(!strcmp(borderName,"pad")) borderType = BorderPad;
  else if(!strcmp(borderName,"mir")) borderType = BorderMirror;
  else if(!strcmp(borderName,"per")) borderType = BorderPeriodic;
  else if(!strcmp(borderName,"pad0")) borderType = BorderPad0;
  else Errorf("Undefined border effect: %s",borderName);


  while(opt = ParseOption(&argv))
    {
      switch(opt) 
        {
        case 'd': 
          method = CV_DI;
          break;
        case 'm': 
          method = CV_MP;
          break;
        case 'f': 
          method = CV_FT;
          break;          
        default:
          ErrorOption(opt);
        }
    }
  NoMoreArgs(argv);


  time = MyTime();
  
  ConvSig(in,filter,out,borderType,method);

  SetResultFloat(MyTime()-time);                 
}

/* A correlation function */
void CorrSig(SIGNAL in1, SIGNAL in2, SIGNAL out1, LWFLOAT dxmin, LWFLOAT dxmax, char flagNormalized, char flagCausal)
{
  int imin,imax,i,di;
  int jmin,jmax,j;
  int first1,last1;
  int first2,last2;
  int N;
  LWFLOAT mean1,mean2,var1,var2;
  SIGNAL out;
  
  if (in1->dx != in2->dx) Errorf("CorrSig() : Sorry the signals must have the same dx");
  if (out1 == in1 || out1 == in2) out = NewSignal();
  else out = out1;

  imin = (int) (dxmin/in1->dx);
  imax = (int) (dxmax/in1->dx);
  di = (int) ((in2->x0-in1->x0)/in1->dx);

  SizeSignal(out,imax-imin+1,YSIG);
  out->dx = in1->dx;
  out->x0 = imin;
  ZeroSig(out);
  
  if (imin > imax) {
    if (out != out1) DeleteSignal(out);
    Errorf("CorrSig() : Sorry, dxmin should be smaller than dxmax");
  }
    
  if (flagCausal) {
    first1 = in1->firstp;
    last1 = in1->lastp;
    first2 = in2->firstp;
    last2 = in2->lastp;
  }
  else {
    first1 = 0;
    last1 = in1->size-1;
    first2 = 0;
    last2 = in2->size-1;
  }

  N = MIN(last1-first1+1,last2-first2+1);
  
  mean1 = 0;
  var1 = 0;
  for (i = first1; i<=last1;i++) {
    mean1 += in1->Y[i];    
    var1 += in1->Y[i]*in1->Y[i];
  }
  mean1 /= last1-first1+1;
  var1 /= last1-first1;
  var1 -= mean1*mean1;
    
  mean2 = 0;
  var2 = 0;
  for (i = first2; i<=last2;i++) {
    mean2 += in2->Y[i];    
    var2 += in2->Y[i]*in2->Y[i];
  }
  mean2 /= last2-first2+1;
  var2 /= last2-first2;
  var2 -= mean2*mean2;

  if (!flagNormalized) var1 = var2 = 1;
  
  for (i = imin; i<= imax;i++) {
    
    jmin = MAX(first1,first2+i+di);
    jmax = MIN(last1,last2+i+di);
    
    for (j=jmin;j<=jmax;j++) out->Y[i-imin] += (in1->Y[j]-mean1)*(in2->Y[j-i-di]-mean2);

/* NOT CONSISTENT!    out->Y[i-imin] /= (jmax-jmin+1); */

    out->Y[i-imin] /= N; 
    out->Y[i-imin] /= sqrt(var1*var2);

  }        
  
  if (out != out1) {
    CopySig(out,out1);
    DeleteSignal(out);
  }
}

void C_CorrSig(char **argv)
{
  SIGNAL in1,in2,out;
  LWFLOAT dxmin, dxmax;
  char flagCausal,flagNormalized;
  char opt;
  
  argv = ParseArgv(argv,tSIGNALI,&in1,tSIGNALI,&in2,tSIGNAL,&out,tFLOAT,&dxmin,tFLOAT,&dxmax,-1);
  
  if (dxmin >= dxmax) Errorf("Sorry, dxmin should smaller than dxmax");
  
  flagCausal = YES;
  flagNormalized = YES;
  while(opt = ParseOption(&argv)) {
    switch(opt) {
    case 'c':
      flagCausal = NO;
      break;
    case 'n':
      flagNormalized = NO;
      break;
    default: 
      ErrorOption(opt);
    }
  }    
  NoMoreArgs(argv);

  CorrSig(in1,in2,out,dxmin,dxmax,flagNormalized,flagCausal);
}