dconvolution.c

LastWave

  else {
    l1 = r1 = scale / 2;
    if(filtshift == 0) {
      symsize = sigsize;
      type = 0;
    } 
    else {
      type = 1;
      symsize = sigsize;
      if (parity == 0) sigsize -= sigsize % 2;
      else sigsize -= (1-sigsize%2);
    }
  }

 
  /**************************************************************/
  /** Put the fields of the output signal unless there is      **/
  /** something to remember                                    **/
  /**************************************************************/

  if (flag_remember == YES) {
    if (signal_result->sizeMallocY < sigsize) {
      temp_signal = NewSignal();
      CopySig(signal_result,temp_signal);
      SizeSignal(signal_result,sigsize+200,YSIG);
      CopySig(temp_signal,signal_result);
      DeleteSignal(temp_signal);
    }
    signal_result->size = sigsize;
  }
  else {
    SizeSignal(signal_result,sigsize,YSIG);
    signal_result->dx = signal_input->dx;
    signal_result->x0 = signal_input->x0;
    if (scale==1 && filt->shift==1)
      signal_result->x0 -= 0.5*signal_input->dx;
    signal_result->firstp = signal_input->firstp + halfsizefilt - filtshift1;
    signal_result->param = nbinterfilt + signal_input->param;
  };

  signal_result->lastp = signal_input->lastp - halfsizefilt;
  result = signal_result->Y;



  /**************************************************************/
  /**           The convolution                                **/
  /**************************************************************/
  
  switch(type) {
    
    case 0 :
      for(j = start_index; j < sigsize; j += factor) {
         sum = filter[0] * SYMREF0(j,symsize,input);
         for(k = 1, left = j-l1, right = j+r1;k < filtsize ;
	         k++, left -= scale, right += scale)
	         sum += filter[k] *
	             (-SYMREF0(per(sigsize2,left),symsize,input) + 
	              SYMREF0(per(sigsize2,right),symsize,input));
	     result[j] = sum/fnorm;
	  }
      break;
      
    default :    
      for(j = start_index; j < sigsize; j += factor) {
         sum = filter[0] * SYMREF1(j,symsize,input);
         for(k = 1, left = j-l1, right = j+r1;k < filtsize ;
	         k++, left -= scale, right += scale)
	        sum += filter[k] *
	             (-SYMREF1(per(sigsize2,left),symsize,input) + 
	             SYMREF1(per(sigsize2,right),symsize,input));
         result[j] = sum/fnorm;
      }
  }
}


/**************************************************************/
/* Convolution for the symetrical case	(reconstruction)      */
/**************************************************************/

void conv_symr(SIGNAL signal_input, SIGNAL signal_result, FILTER filt, 
               int factor, int scale, int parity)
{
  LWFLOAT *input = signal_input->Y;
  int sigsize = signal_input->size;
  LWFLOAT *result;
  int filtsize = filt->size;
  LWFLOAT *filter = filt->Y;
  int filtshift = filt->shift;  int j, k, left, right, l1, r1;
  LWFLOAT sum;
  int halfsizefilt, nbinterfilt, filtshift1;
  int symsize,type;
  int sigsize2;
     
    
  if (parity == 0) sigsize2 = 2*(sigsize-sigsize%2);
  else sigsize2 = 2*(sigsize-(1-sigsize%2));   
    
  
  /**************************************************************/
  /** Compute :                                                **/
  /**           halfsizefilt : half size of the scaled filter  **/
  /**           nbinterfilt  : number of interval of length 1  **/
  /**                          of the support of the scaled    **/
  /**                          filter                          **/
  /**           filtshift1   : shift of the scaled filter      **/
  /**                          (according to integer values    **/
  /**************************************************************/
  
  if (scale == 1)
    {
      filtshift1 = filtshift;
      halfsizefilt = filtsize-1;
      nbinterfilt=2*halfsizefilt-filtshift;
    }
  else {
    filtshift1 = 0;
    halfsizefilt = scale*(filtsize-1) - filtshift*scale/2;
    nbinterfilt = 2*halfsizefilt;
  }
  

  if(filtshift == 0) {
    l1 = r1 = scale;
    type = 0;
    symsize = sigsize;
  } 
  else {
    if(scale == 1) {
      l1 = 0;
      r1 = 1;
      type = 1;
      symsize = sigsize;
      sigsize--;
    } 
    else {
      l1 = r1 = scale / 2;
      type = 1;
      symsize = sigsize;
    }
  }

  SizeSignal(signal_result,sigsize,YSIG);
  result = signal_result->Y;
  signal_result->dx = signal_input->dx;
  signal_result->x0 = signal_input->x0;

  switch(type) {
    
    case 0 :
      for(j = 0; j < sigsize; j += factor) {
         sum = filter[0] * SYMREF0(j,symsize,input);
         for(k = 1, left = j-l1, right = j+r1;k < filtsize ;
	         k++, left -= scale, right += scale)
	        sum += filter[k] *
	               (SYMREF0(per(sigsize2,left),symsize,input) + 
	                SYMREF0(per(sigsize2,right),symsize,input));
	     result[j] = sum;
	  }
      break;
      
    default :    
      for(j = 0; j < sigsize; j += factor) {
         sum = filter[0] * SYMREF1(j,symsize,input);
         for(k = 1, left = j-l1, right = j+r1;k < filtsize ;
	         k++, left -= scale, right += scale)
	        sum += filter[k] *
	               (SYMREF1(per(sigsize2,left),symsize,input) + 
	                SYMREF1(per(sigsize2,right),symsize,input));
         result[j] = sum;
      }
  }
}


/**************************************************************/
/* Convolution for the asymetrical case	(reconstruction)      */
/* fnorm is the factor renormalization signal_input should be */
/* multiplied by                                              */
/**************************************************************/

void conv_asyr(SIGNAL signal_input, SIGNAL signal_result, FILTER filt, 
               int factor, int scale, LWFLOAT fnorm, int parity)
{
  LWFLOAT *input = signal_input->Y;
  int sigsize = signal_input->size;
  LWFLOAT *result;
  int filtsize = filt->size;
  LWFLOAT *filter = filt->Y;
  int filtshift = filt->shift;
  int j, k, left, right, l1, r1;
  LWFLOAT sum;
  int halfsizefilt, nbinterfilt, filtshift1;
  int symsize,type;
  
  int sigsize2;
     
    
  if (parity == 0) sigsize2 = 2*(sigsize-sigsize%2);
  else sigsize2 = 2*(sigsize-(1-sigsize%2));
   
    
  
  /**************************************************************/
  /** Compute :                                                **/
  /**           halfsizefilt : half size of the scaled filter  **/
  /**           nbinterfilt  : number of interval of length 1  **/
  /**                          of the support of the scaled    **/
  /**                          filter                          **/
  /**           filtshift1   : shift of the scaled filter      **/
  /**                          (according to integer values    **/
  /**************************************************************/
  
  if (scale == 1)
    {
      filtshift1 = filtshift;
      halfsizefilt = filtsize-1;
      nbinterfilt=2*halfsizefilt-filtshift;
    }
  else {
    filtshift1 = 0;
    halfsizefilt = scale*(filtsize-1) - filtshift*scale/2;
    nbinterfilt = 2*halfsizefilt;
  }
  
  
  if(scale == 1) {
    l1 = 0;
    r1 = 1;
    type = 1;
    symsize = sigsize;
  } 
  else {
    l1 = r1 = scale / 2;
    if(filtshift == 0) {
    type = 0;
    symsize = sigsize;
    } 
    else {
      type = 1;
      symsize = sigsize;
      sigsize++;
    }
  }

  SizeSignal(signal_result,sigsize,YSIG);
  result = signal_result->Y;
  signal_result->dx = signal_input->dx;
  signal_result->x0 = signal_input->x0;

  switch(type) {
    
    case 0 :
      for(j = 0; j < sigsize; j += factor) {
         sum = filter[0] * ASYMREF0(j,symsize,input);
         for(k = 1, left = j-l1, right = j+r1;k < filtsize ;
	         k++, left -= scale, right += scale)
	        sum += filter[k] *
	               (ASYMREF0(per(sigsize2,left),symsize,input) - 
	                ASYMREF0(per(sigsize2,right),symsize,input));
	     result[j] = sum*fnorm;
	  }
      break;
      
    default :    
      for(j = 0; j < sigsize; j += factor) {
         sum = filter[0] * ASYMREF1(j,symsize,input);
         for(k = 1, left = j-l1, right = j+r1;k < filtsize ;
	         k++, left -= scale, right += scale)
	        sum += filter[k] *
	               (ASYMREF1(per(sigsize2,left),symsize,input) - 
	               ASYMREF1(per(sigsize2,right),symsize,input));
         result[j] = sum*fnorm;
      }
  }
}