convol.c

LastWave

           for (n=0; n < size_conv; n++)
	       {
                   i             = (n+n1)  < 0 ? 2*((-(n+n1))/signal1->size+1)*(signal1->size) + (n+n1): (n+n1); 
                   div           = i % (2*signal1->size - 1 );
		   tempout->Y[n] = (0 <= div && div < signal1->size) ? signal1->Y[div] : signal1->Y[2*signal1->size - 1 - div];   
	       }
	   rapide (tempout, signal2, out1);
           //WriteSigFile(out1,"totos",NO, "y",NO);
           DeleteSignal(tempout);           
           memcpy(out->Y, out1->Y + signal2->size - 1, size*sizeof(LWFLOAT));
           DeleteSignal(out1);           
           out->x0     = x_begin_conv;
           out->firstp = MAX(0,-n1+signal1->firstp);
           out->lastp  = out->size-1-MAX(0,n2-signal1->lastp);
           
           break;           
	   }
}
  

void  DirectConvolution (SIGNAL signal1,SIGNAL signal2, SIGNAL signalOut, int borderType, LWFLOAT x_begin_conv, LWFLOAT x_end_conv)
{ 
    int     size_conv, begin, end, a, n1, n2;
    int     k, n, nmk, div;
    LWFLOAT   sum, step, minsupport, maxsupport;
   
 
   step          = signal1->dx;
   signalOut->dx =  step;  
   minsupport    = (x_begin_conv - (signal2->x0 + (signal2->size - 1)*step) - signal1->x0)/step;
   n1            = FLOAT2INT(minsupport);  
   maxsupport    =  (x_end_conv   - (signal2->x0 + (signal2->size - 1)*step)  -signal1->x0)/step;
   n2            = FLOAT2INT(maxsupport);  
   //Printf("n1lente=%d, n2lente=%d \n", n1, n2);

   switch (borderType) {

       case  Border0:
           size_conv         = n2 - n1 + 1 ;
           SizeSignal(signalOut, size_conv, YSIG)    ;
           signalOut->firstp = MAX(0,-n1+signal1->firstp);
           signalOut->lastp  = signalOut->size-1-MAX(0,n2+signal2->size-1-signal1->lastp);
           signalOut->x0     = x_begin_conv;            
           for (n=0; n < size_conv; n++)
	     {
	       sum   = 0;
               a     = n + signal2->size - 1 + n1;  
               begin =   a - signal1->size + 1;
               begin = MAX(begin,0);
               end   = MIN (signal2->size - 1, a);
               for (nmk = a-begin, k = begin; nmk >= a-end && k<=end ; nmk--, k++)
	       	       sum += signal2->Y[k]*signal1->Y[nmk];

               signalOut->Y[n]=sum;
	     }
	   break;

       case  BorderCon:
	   //size_conv         = (x_end_conv - x_begin_conv)/step + 1;
           size_conv         = n2 - n1 + 1 ;
           SizeSignal(signalOut, size_conv, YSIG) ;
           signalOut->firstp = MAX(0,-n1+signal1->firstp);
           signalOut->lastp  = signalOut->size-1-MAX(0,n2+signal2->size-1-signal1->lastp);
           signalOut->x0     =  x_begin_conv;
           for (n=0; n < size_conv; n++)
	     {
	       sum=0;
               a = n + signal2->size - 1 + n1;
               for (nmk = a, k=0 ; nmk >= a - signal2->size + 1 && k < signal2->size; nmk--, k++)
		   sum +=  (nmk<0) ? (signal2->Y[k] * signal1->Y[0]) : ((nmk>=signal1->size) ? (signal2->Y[k] * signal1->Y[signal1->size-1]) : (signal2->Y[k]*signal1->Y[nmk])); 
               signalOut->Y[n] = sum;
	     }
	   break;
   
        case  BorderPer:      
	   //size_conv  = (x_end_conv - x_begin_conv)/step + 1;
           size_conv         = n2 - n1 + 1 ;
           SizeSignal(signalOut, size_conv, YSIG);
           signalOut->firstp = MAX(0,-n1+signal1->firstp);
           signalOut->lastp  = signalOut->size-1-MAX(0,n2+signal2->size-1-signal1->lastp);
           signalOut->x0     = x_begin_conv;
           for (n=0; n < size_conv; n++)
	     {
	       sum=0;
               a = n + signal2->size - 1 + n1;
	       for (k=0; k < signal2->size; k++)
		 {
                 
                   nmk  = (a-k) < 0 ? ((-(a-k))/signal1->size+1)*signal1->size + a-k : a-k;
                   sum += signal2->Y[k] * signal1->Y[nmk%signal1->size]; 
		 }
               signalOut->Y[n] = sum;
	     }
 	   break;

        case  BorderMir:
	  // sans rajout de point
	   //size_conv         = (x_end_conv - x_begin_conv)/step + 1;
           size_conv         = n2 - n1 + 1 ;
           SizeSignal(signalOut, size_conv, YSIG);
           signalOut->x0     =  x_begin_conv;
           signalOut->firstp = MAX(0,-n1+signal1->firstp);
           signalOut->lastp  = signalOut->size-1-MAX(0,n2+signal2->size-1-signal1->lastp);
           for (n=0; n < size_conv; n++)
	     {
	       sum=0;
               a  =  n + signal2->size - 1 + n1;
               for (k=0 ; k< signal2->size; k++)
		 {
                   nmk = (a-k) < 0 ? 2*((-(a-k))/signal1->size+1)*(signal1->size) + (a-k): a-k;
                   div  = (signal1->size == 1) ? 0 : nmk%(2*signal1->size - 2);
		   sum += (0 <= div && div < signal1->size) ?  signal2->Y[k] * signal1->Y[div] : signal2->Y[k] * signal1->Y[2*signal1->size-2-div];   
		 }
             signalOut->Y[n] = sum;
	     }
	break;

         case  BorderMir1:
           // point en plus aux extremites
           //size_conv         = (x_end_conv - x_begin_conv)/step + 1;
           size_conv         = n2 - n1 + 1 ;
           SizeSignal(signalOut, size_conv, YSIG);
           signalOut->x0     = x_begin_conv;
           signalOut->firstp = MAX(0,-n1+signal1->firstp);
           signalOut->lastp  = signalOut->size-1-MAX(0,n2+signal2->size-1-signal1->lastp);
           for (n=0; n < size_conv; n++)
	     {
	       sum=0;
               a = n + signal2->size - 1 + n1;
               for (k=0 ;  k<signal2->size;  k++)
		 {
                   nmk  = (a-k) < 0 ? 2*((-(a-k))/signal1->size+1)*signal1->size + (a-k): a-k;
                   div  = nmk % (2*signal1->size - 1);
		   sum += (0 <= div && div < signal1->size) ? signal2->Y[k] * signal1->Y[div] : signal2->Y[k] * signal1->Y[2*signal1->size - 1 -div];   
		 }
             signalOut->Y[n] = sum;
	     }
	break;
	   }
 
/* end switch */
}

  
void C_Conv (char ** argv)
{
  int     size_conv;
  int     borderType, fastconvolution;
  LWFLOAT   begin_conv, end_conv, step;
  SIGNAL  signal, filter, out; 
  char    * borderName, opt;
  char **argv1;
  LWFLOAT t;
   
   
  argv = ParseArgv(argv, tSIGNALI, &signal, tSIGNALI, &filter, tSIGNAL_, NULL, &out, tSTR, &borderName,-1);
  
  if (!strcmp(borderName,"bconst")) borderType = BorderCon;
  else if(!strcmp(borderName,"b0")) borderType = Border0;
  else if(!strcmp(borderName,"bmirror")) borderType = BorderMir;
  else if(!strcmp(borderName,"bperiodic"))  borderType = BorderPer;
  else if(!strcmp(borderName,"bmirror1")) borderType = BorderMir1;
  else Errorf("Undefined border effect: %s", borderName);

  step  = signal->dx;
  if (step != filter->dx) Errorf("Both signal and filter must have the same dx");

  if (signal==out || filter==out) Errorf("Output signal must be different from input and filter");

  
  if (borderType==BorderCon || borderType==Border0) {  
    begin_conv = signal->x0 + filter->x0;
    size_conv  = signal->size + filter->size - 1;
    end_conv   = begin_conv + (size_conv - 1) * step;
  }
  else {
    begin_conv = signal->x0;
    if (borderType==BorderMir) size_conv  = 2*signal->size - 1;
    else if (borderType==BorderMir1) size_conv  = (signal->size == 1) ? 1 : 2*signal->size;
    else size_conv  = signal->size;
    end_conv   = begin_conv + (size_conv - 1) * step;
  }
 
  fastconvolution = 0;
  while (opt = ParseOption(&argv)) {        
    switch(opt) {
    
    case 'x':
	  argv1 = ParseArgv(argv, tFLOAT_, 0., &begin_conv, tFLOAT_, 0., &end_conv,-1);
      if (argv1 != argv+2) {
	    if (borderType==BorderCon || borderType==Border0 || borderType==BorderPer) {  
          begin_conv = signal->x0;
		  size_conv  = signal->size;
		  end_conv   = begin_conv + (size_conv - 1) * step;
	    }
        if (borderType==BorderMir || borderType==BorderMir1) {  
	      begin_conv = signal->x0;
		  size_conv  = signal->size;
		  end_conv   = begin_conv + (size_conv - 1) * step;
	    }      
	  }
	  else {
	    argv = argv1;
	    if (begin_conv>=end_conv) Errorf("xMin should be greater than xMax");
	  }
      break;

	case 'f':
	  fastconvolution = 1;
	  break;
  
    default:
      ErrorOption(opt);
    } 
  }
  NoMoreArgs(argv); 
  
  t = MyTime();
  if (fastconvolution == 1) FFTConvolution (signal, filter, out, borderType, begin_conv, end_conv);
  else DirectConvolution (signal, filter, out, borderType, begin_conv, end_conv);
  t = MyTime()-t;
  
  SetResultFloat(t); 
}