/*..........................................................................*//*                                                                          *//*      L a s t W a v e    P a c k a g e 'signal' 2.1                       *//*                                                                          *//*      Copyright (C) 1998-2002 Xavier Suraud, Emmanuel Bacry.              *//*      email : lastwave@cmap.polytechnique.fr                              *//*                                                                          *//*..........................................................................*//*                                                                          *//*      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             *//*                                                                          *//*..........................................................................*/#include "lastwave.h"#include "signals.h"#define FLOAT2INT(X) ( ((X) >=0) ? (int)((X) + .5) : (int)((X) - .5))void extraction (SIGNAL signal, SIGNAL out, int indice, int size2)  {   int k, cut;          cut = ((indice + 1)*size2 - 1) - (signal->size - 1);   cut =  cut > 0 ? cut : 0;    memcpy(out->Y, (signal->Y)+(indice*size2), (size2 - cut)*sizeof(LWFLOAT));   for (k = size2 - cut; k< 2*size2; k++)	out->Y[k] = 0;  }void multiplication (SIGNAL signal1Real, SIGNAL signal1Im, SIGNAL signal2Real, SIGNAL signal2Im)  {      int     n;     LWFLOAT   real1;     if (signal1Real->size != signal2Real->size || signal1Im->size !=  signal2Im->size)       Errorf("C_multiplication: les tailles doivent etres egales");          for (n=0; n < signal1Real->size; n++)        {	real1             = signal1Real->Y[n];        signal1Real->Y[n] = real1 * signal2Real->Y[n] - signal1Im->Y[n] * signal2Im->Y[n];        signal1Im->Y[n]   = real1 * signal2Im->Y[n]   + signal1Im->Y[n] * signal2Real->Y[n];        }  }      void rapide (SIGNAL signal1, SIGNAL signal2, SIGNAL out)  {    int     size2;    int     begin, end, u, v, indice, n;    double  R;    SIGNAL  b, FfilterIm, FfilterReal, Fa_rReal, Fa_rIm;      n           = (int) ((signal2->size == 1) ? 1 : ceil(log(signal2->size)/log(2)));      size2       =  ((signal2->size < 1<<n) && (signal2->size > (1<<(n-1)))) ? (1<<n): (1<<n-1);               SizeSignal(out, signal1->size + signal2->size, YSIG);     R           = ceil((double)signal1->size / (double)size2);     b           = NewSignal();     SizeSignal(b, 2*size2, YSIG);     b->dx       = signal1->dx;     extraction (signal2, b, 0, size2);     FfilterReal = NewSignal();     FfilterIm   = NewSignal();     Fft(b, NULL, FfilterReal, FfilterIm, 1, NO);     // WriteSigFile(FfilterReal,"totof1",NO, "y",YES);     // WriteSigFile(FfilterIm,  "totof2",NO, "y",YES);               Fa_rReal  = NewSignal();     Fa_rIm    = NewSignal();             for (n = 0; n < signal1->size + signal2->size ; n ++)            out->Y[n] = 0;     for (indice = 0; indice < R; indice ++)	{           extraction (signal1, b, indice, size2);           Fft (b, NULL, Fa_rReal, Fa_rIm, 1, NO);           multiplication (Fa_rReal, Fa_rIm, FfilterReal, FfilterIm);           Fft (Fa_rReal, Fa_rIm, b, NULL, -1, NO);           //WriteSigFile(b,"totos",NO, "y",NO);           begin = indice * size2;           u     = (2 + indice) * size2 - 1;           v     = out->size - 1;           end = MIN(u,v);	   for (n=begin; n<=end; n++)                   out->Y[n] += b->Y[n-indice*size2];	}               DeleteSignal(Fa_rReal);     DeleteSignal(Fa_rIm);     DeleteSignal(b);     DeleteSignal(FfilterReal);     DeleteSignal(FfilterIm);       }    void  FFTConvolution (SIGNAL signal1, SIGNAL signal2, SIGNAL out, int borderType, LWFLOAT x_begin_conv, LWFLOAT x_end_conv){    int     size_conv, size;    int     n, n1, n2 ,i ,div;    LWFLOAT   step, endfilter;    LWFLOAT   minsupport, maxsupport;        SIGNAL  tempout, out1;   tempout    = NewSignal();   out1       = NewSignal();   step       = signal1->dx;   endfilter  = (LWFLOAT)(signal2->size - 1)*step;      /* le support  total du signal qu'il faudra balayer ???*/   minsupport = (x_begin_conv - signal1->x0 - (signal2->x0 + endfilter))/step;   maxsupport = (x_end_conv   - signal1->x0 - signal2->x0)/step;   n1         = FLOAT2INT(minsupport);   n2         = FLOAT2INT(maxsupport);   size       = FLOAT2INT((x_end_conv - x_begin_conv)/step + 1);   //Printf("n1rapide=%d pour %g , n2rapide=%d\n", n1, step, n2);   switch (borderType){       case  Border0:    	   //size        =  (x_end_conv - x_begin_conv)/step + 1;              SizeSignal(out, size, YSIG);           out->dx     = step;           size_conv   = n2 - n1 + 1;                   SizeSignal(tempout, size_conv , YSIG);           tempout->dx = step;             for (n=0, i=n1 ; n < size_conv && i < n1 + size_conv; n++, i++)	       tempout->Y[n] = i<0 ? 0: i>=signal1->size ? 0 : signal1->Y[i];           rapide (tempout, signal2, out1);            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;       case  BorderCon:	   //size      = (x_end_conv - x_begin_conv)/step + 1;           SizeSignal(out, size, YSIG);           out->dx     = step;           size_conv   = n2 - n1 + 1;            SizeSignal(tempout, size_conv, YSIG);           //Printf ("size = %d, size_conv= %d\n\n", size, size_conv);           tempout->dx = step;           for (n=0, i=n1; n < size_conv && i < n1 + size_conv; n++, i++)	       tempout->Y[n] =  i<0 ? signal1->Y[0]: i>=signal1->size ? signal1->Y[signal1->size-1] : signal1->Y[i];           rapide (tempout, signal2, out1);           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);           //WriteSigFile(out,"totos",NO, "y",NO);       	   break;            case  BorderPer:      	   //size       = (x_end_conv - x_begin_conv)/step + 1;           SizeSignal(out, size, YSIG);           out->dx     = step;           size_conv   = n2 - n1 + 1;           SizeSignal(tempout, size_conv, YSIG);           tempout->dx = step;            for (n=0; n < size_conv; n++) 	     {                  i              = (n+n1) < 0 ? ((-(n+n1))/signal1->size+1)*signal1->size + (n+n1): n+n1;                  tempout->Y[n]  = signal1->Y[i%signal1->size];	     }           rapide (tempout, signal2, out1);           DeleteSignal(tempout);            //WriteSigFile(tempout,"totos",NO, "y",NO);                      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;       case  BorderMir:            // sans rajout de point	 //size        = (x_end_conv - x_begin_conv)/step + 1;           SizeSignal(out, size, YSIG);           out->dx     = step;           size_conv   = n2 - n1 + 1;           SizeSignal(tempout, size_conv, YSIG);           tempout->dx =  step;           for (n=0; n < size_conv; n++)	     {                   i             = (n + n1) < 0 ? 2*((-(n + n1))/signal1->size+1)*(signal1->size) + (n + n1) : n+ n1;                   div           = signal1->size == 1 ? 0 : i % (2*signal1->size - 2);		   tempout->Y[n] = (0 <= div && div < signal1->size) ? signal1->Y[div] : signal1->Y[2*signal1->size - 2 - 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;        case  BorderMir1:           //  point en plus aux extremites	   //size      = (x_end_conv - x_begin_conv)/step + 1;           SizeSignal(out, size, YSIG);           out->dx   = step;           size_conv = n2 - n1 + 1;	   //Printf ("size = %d, size_conv= %d\n\n", size, size_conv);           SizeSignal(tempout, size_conv, YSIG);