signal_functions.c

LastWave

/*..........................................................................*/
/*                                                                          */
/*      L a s t W a v e    P a c k a g e 'signal' 2.1                       */
/*                                                                          */
/*      Copyright (C) 1998-2003 Emmanuel Bacry.                             */
/*      email : [e-mail]                              */
/*                                                                          */
/*..........................................................................*/
/*                                                                          */
/*      This program is a free software, you can redistribute it and/or     */
/*      modify it under the terms of the GNU General Public License as      */
/*      published by the Free Software Foundation; either version 2 of the  */
/*      License, or (at your option) any later version                      */
/*                                                                          */
/*      This program is distributed in the hope that it will be useful,     */
/*      but WITHOUT ANY WARRANTY; without even the implied warranty of      */
/*      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       */
/*      GNU General Public License for more details.                        */
/*                                                                          */
/*      You should have received a copy of the GNU General Public License   */
/*      along with this program (in a file named COPYRIGHT);                */
/*      if not, write to the Free Software Foundation, Inc.,                */
/*      59 Temple Place, Suite 330, Boston, MA  02111-1307  USA             */
/*                                                                          */
/*..........................................................................*/




/****************************************************************************/
/*                                                                          */
/*  signal_functions.c                                                      */
/*                   Miscellaneous useful functions on signals				*/
/*                                                                          */
/****************************************************************************/



#include "lastwave.h"
#include "signals.h"


/*
 * Update a signal by adding a value at its end 
 */
 
#define INCRALLOC 5000

void AppendFloat2Signal(SIGNAL sig,LWFLOAT value)
{
  SIGNAL newSig;
  LWFLOAT *f;
  int i;
  
  if (sig->size+1 > sig->sizeMallocY) {
    
  
    newSig = NewSignal();
    SizeSignal(newSig,sig->size+INCRALLOC,YSIG);
    CopySig(sig,newSig);
    f = sig->Y;
    i = sig->sizeMallocY;
    sig->Y = newSig->Y;
    sig->sizeMallocY = newSig->sizeMallocY;
    newSig->Y = f;    
    newSig->sizeMallocY = i;
    f = sig->X;
    i = sig->sizeMallocX;
    sig->X = newSig->X;
    sig->sizeMallocX = newSig->sizeMallocX;
    newSig->X = f;    
    newSig->sizeMallocX = i;
    DeleteSignal(newSig);
  }
  
  sig->size++;
  sig->Y[sig->size-1]=value;
}    
    
/*********************************/
/* Print the values of a signal  */
/*********************************/

/*?????????? void C_Print(char **argv)
{
  SIGNAL signal;
  int j;

  argv = ParseArgv(argv,tSIGNALI,&signal,0);

  for (j = 0; j < signal->size; j++) {
    Printf("% .2e  %.8g\n", XSig(signal,j),signal->Y[j]);
  }
}*/

/*********************************/
/*     Put a signal to 0         */
/*********************************/

void ZeroSig(SIGNAL sig)
{
  int i;
  
  for(i=0;i<sig->size;i++) sig->Y[i] = 0;
  
}

/* 
 * Compute the Min and the max of a signal 
 * If *pxMin < *pxMax the yMin and yMax are computed only between these values
 */
void MinMaxSig(SIGNAL signal,LWFLOAT *pxMin,LWFLOAT *pxMax,LWFLOAT *pyMin,LWFLOAT *pyMax,
               int *piyMin,int *piyMax,int flagCausal)
{
  int i;
  int iMin,iMax;

  if(signal->size == 0) 
    Errorf("MinMaxSig() : signal is empty signal");

  if (flagCausal == YES) {
    iMin = signal->firstp;
    iMax = signal->lastp;
  }
  else {
    iMin = 0;
    iMax = signal->size-1;
  }
  
  if (*pxMin < *pxMax) {
    i = ISig(signal,*pxMin);
    iMin = MAX(iMin,i);
    i = ISig(signal,*pxMax);
    iMax = MIN(iMax,i);
  }
  
  for (i = iMin, (*pxMin) = (*pxMax) = XSig(signal,iMin),
              (*piyMin) = (*piyMax) = iMin, 
  			  (*pyMin) = (*pyMax) = signal->Y[iMin]; i <= iMax; i++) 
  	{
  		(*pxMin) = MIN(XSig(signal,i),(*pxMin));
  		if (signal->Y[i] < *pyMin) {
  		  (*pyMin) = signal->Y[i];
  		  *piyMin = i;
  		}
  		if (signal->Y[i] > *pyMax) {
  		  (*pyMax) = signal->Y[i];
  		  *piyMax = i;
  		}
  		(*pxMax) = MAX(XSig(signal,i),(*pxMax));
	}
	
}


/*****************************************************************/
/* This procedure thresholds a signal                            */
/*   flagX == YES ==> the threshold is on the X's                */
/*   flagY == YES ==> the threshold is on the Y's                */
/*   flagMin == YES ==> a min value is specified in 'min'        */
/*   flagMax == YES ==> a max value is specified in 'max'        */
/*****************************************************************/

void ThreshSig(SIGNAL in, SIGNAL out, int flagX, int flagY,
                  int flagMin, LWFLOAT min, int flagMax,LWFLOAT max)
{
  int i,j;
  LWFLOAT x,y;
  
  if (in == out) Errorf("ThreshSig() : the input and output signals must be different");
  
  if (flagY == YES) SizeSignal(out,in->size,XYSIG);
  else SizeSignal(out,in->size,in->type);
  
  out->firstp = out->lastp = -1;
  out->x0 = in->x0;
  out->dx = in->dx;
  out->size = 0;
  
  for (i=0,j=0;i<in->size;i++) {
    x = XSig(in,i);
    y = in->Y[i];
   
    if (flagX == YES && ((flagMin == YES && x < min) || (flagMax == YES && x > max))) continue;
    if (flagY == YES && ((flagMin == YES && y < min) || (flagMax == YES && y > max))) continue;

    if (out->firstp == -1 && i >= in->firstp) out->firstp = j;
    if (i <= in->lastp) out->lastp = j;
    
    if (out->type == XYSIG) out->X[j] = x;
    else if (j == 0) out->x0 = x;
    else if (j == 1) out->dx = x - out->x0;
      
    out->Y[j++] = y;
    out->size++; 
  }
  
  if (out->firstp == -1 || out->lastp == -1) {
    out->firstp = 1;
    out->lastp = 0;
  }
}


/* The corresponding command */
void C_Thresh(char **argv)
{
  SIGNAL in,out;
  int flagX,flagY,flagMin,flagMax;
  LWFLOAT min,max;
  char *arg1,*arg2;
  char opt;
    
  argv = ParseArgv(argv,tSIGNALI,&in,tSIGNAL,&out,-1);
  if (in == out) Errorf("Signals must be different");
  
  flagX = flagY = NO;
  while(opt = ParseOption(&argv)) { 
   switch(opt) {
   case 'x': 
     argv = ParseArgv(argv,tWORD,&arg1,tWORD,&arg2,-1);
     if (flagY == YES) Errorf("Can't set both option '-x' and '-y'");
     flagX = YES;
     if (!strcmp(arg1,"*")) flagMin = NO;
     else {
       flagMin = YES;
       ParseFloat(arg1,&min);
     }
     if (!strcmp(arg2,"*")) flagMax = NO;
     else {
       flagMax = YES;
       ParseFloat(arg2,&max);
     }
     break;
   case 'y': 
     argv = ParseArgv(argv,tWORD,&arg1,tWORD,&arg2,-1);
     if (flagX == YES) Errorf("Can't set both option '-x' and '-y'");
     flagY = YES;
     if (!strcmp(arg1,"*")) flagMin = NO;
     else {
       flagMin = YES;
       ParseFloat(arg1,&min);
     }
     if (!strcmp(arg2,"*")) flagMax = NO;
     else {
       flagMax = YES;
       ParseFloat(arg2,&max);
     }
     break;
    default:
     ErrorOption(opt);
   }
 }    
 NoMoreArgs(argv);

 if (flagMin == NO && flagMax == NO) Errorf("You have to specify either the Min or the Max");
   
 ThreshSig(in,out,flagX,flagY,flagMin,min,flagMax,max);
}
 

/****************************************/
/*    Sort a signal according           */
/****************************************/

/* Sorting two arrays of LWFLOAT        */
/* The sorting is made according to x */
/* The array y can be == NULL         */
/* 'n' is the size of the arrays      */

void SortArrays(LWFLOAT *x, LWFLOAT *y, int n)
{
  int l,j,ir,i;
  LWFLOAT xx,yy;
		
  l = ( n>>1)+1;
  ir = n;
  for (;;) {
    if(l>1) {
      xx = x[--l-1];
      if (y != NULL) yy = y[l-1];
    }
    else {
      xx = x[ir-1];
      x[ir-1] = x[0];
      if (y != NULL) {
        yy = y[ir-1];
        y[ir-1] = y[0];
      }
      if(--ir ==1) { 
      	x[0] = xx; 
      	if (y != NULL) y[0] = yy;
      	return; 
      }
    }
    i=l;
    j=l<<1;
    while( j<= ir) {
      if( j<ir && x[j-1]< x[j]) ++j;
      if(xx < x[j-1]) {
        x[i-1] = x[j-1];
        if (y != NULL) y[i-1] = y[j-1];
        j+=(i=j);
      }
      else j=ir+1;
    }
    x[i-1] = xx;
    if (y != NULL) y[i-1] = yy;
  }
}

/* Sort a signal either according to X or to Y */
void SortSig(SIGNAL signal)
{
  LWFLOAT *x,*y;
  int n;
	
  if (signal->size == 0) Errorf("SortSig() : Signal must be of size != 0");
  if (signal->size == 1) return;

  /* Do we need to sort ? */
  if (signal->type == XYSIG) {
    for (n=0;n<signal->size-1;n++) {
      if (signal->X[n] > signal->X[n+1]) break;
    }
    if (n == signal->size-1) return;
  }
  	    	
  /* Let's sort */	    	
  n = signal->size;
  switch (signal->type) {
    case XYSIG: 
      x = signal->X; 
      y = signal->Y; 
      break;
    case YSIG: 
      x = signal->Y; 
      if (signal->type == XYSIG) y = signal->X; 
      else y = NULL;
      break;
    default:
      Errorf("SortSig() : (Weired bug) Unknown signal type");
  }
  
  SortArrays(x,y,n);
}

/* The corresponding command */
void C_Sort(char **argv)
{
  SIGNAL signal;
  
  argv = ParseArgv(argv,tSIGNALI,&signal,0);
  
  SortSig(signal);
}


/******************************************/
/*    Make one XY signal from two signals */
/******************************************/

/* flag == YES  ==> out->X = in1->Y and out->Y = in2->Y            */
/* flag == NO   ==> 'in1->X' is put in 'in2' and 'in1->Y' in 'out' */

void MergeSig(SIGNAL in1, SIGNAL in2, SIGNAL out, int flag)
{
  int i;
  
  if (in1 == out || in2 == out)
    Errorf("MergeSig() : the output signal must be different from the input signals");
  
  if (flag == NO) {
    SizeSignal(out,in1->size,YSIG);
    SizeSignal(in2,in1->size,YSIG);
    
    for(i=0;i<in1->size;i++) {
      in2->Y[i] = XSig(in1,i);
      out->Y[i] = in1->Y[i];
    }
  } 
  else {
    SizeSignal(out,in1->size,XYSIG);  
    for(i=0;i<in1->size;i++) {
      out->X[i] = in1->Y[i];
      out->Y[i] = in2->Y[i];
    }
    SortSig(out);
  }
}



/*************************************
 *
 * Extract a signal from another one 
 *
 *************************************/

void ExtractSig(SIGNAL sig,SIGNAL sigOut, int borderType, int firstPoint,int newSize)
{
  SIGNAL sigOut1;
  int i,j,f;
  
  /* Some Checkings */
  if (sig->type == XYSIG && (firstPoint < 0 || firstPoint+newSize > sig->size))
    Errorf("ExtractSig() : Sorry I do not know how to manage border effects of a XY signal !");
    
  /* Set the sigOut1 */
  if (sigOut == sig) sigOut1 = NewSignal();
  else sigOut1 = sigOut;
    
  /* Init the sigOut1 */
  SizeSignal(sigOut1,newSize,sig->type);
  sigOut1->dx = sig->dx;  
  sigOut1->x0 = sig->x0 + firstPoint*sig->dx;
  sigOut1->firstp = sig->firstp - firstPoint;
  if (sigOut1->firstp <=0) sigOut1->firstp = 0;
  sigOut1->lastp = sig->lastp - firstPoint;
  if (sigOut1->lastp >= sigOut1->size) sigOut1->lastp = sigOut1->size-1;

  /* The main loop in case of XY sig */
  if (sig->type == XYSIG) {
    for (i= firstPoint,j=0;j<newSize;i++,j++) {
      sigOut1->Y[j] = sig->Y[i];
      sigOut1->X[j] = sig->X[i];
    }
  }

  /* The main loop in case of Y sig */
  else {  
  switch( borderType) {

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

    case BorderPad0 :
      for (i= firstPoint,j=0;j<newSize;i++,j++) {
        if (i<0) sigOut1->Y[j] = 0;
        else if (i>=sig->size) sigOut1->Y[j] = 0;
        else sigOut1->Y[j] = sig->Y[i];
      }