source: Sophya/trunk/Poubelle/archediab.old/archediab.sources/c/fourier.c@ 2007

#include "math.h"

#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr

void four1(float *data,int nn,int isign);

double souslin(double *y,int n);                /* soustrait la droite de y[0] à y[n-1]  */
void filtre(double *y,int n,double f);
void spectre(double *y,int n);
void spectreRMS(double *y,int n);


/*  --------  filtre passe - bas du tableau double y[] -----------*/

/*double souslin(double *y,int n)               /* soustrait la droite de y[0] à y[n-1]  
{double pas;int i;
pas=(y[n-1]-y[0])/(double)(n-1);
for (i=0;i<n;i++) y[i]=y[i]-y[0]-pas*(double)i;
return(pas);
}
*/

double souslin(double *y,int n)         /* soustrait la pente et la moyenne     */
{double pas,mm=0;int i;
pas=(y[n-1]-y[0])/(double)(n-1);
for (i=0;i<n;i++) {y[i]=y[i]-pas*(double)i;mm+=y[i];}
mm/=(double)n;
for (i=0;i<n;i++) y[i]=y[i]-mm;
return(pas);
}

void filtre(double *y,int n,double f)
{
float *v,u;int i;double pas;
v=(float *)y;
pas=souslin(y,n);
/*
printf("y0 ... = %g %g %g %g %g\n",y[0],y[1],y[2],y[3],y[4]);
printf("filtre avec bande=%g  : pente=%g \n",f,pas);
*/
for (i=0;i<n;i++)               v[i]=(float)y[i];
for (i=n-1;i>=0;i--)    {v[2*i+1]=v[i]/(float)n;v[2*i+2]=0;}

four1(v,n,1);
for (i=(int)f;i<n/2;i++)                        
                {u=((double)i-f)/f;
                u=1/(1+2*u*u);
                v[2*i+1]=v[2*i+1]*u;
                v[2*i+2]=v[2*i+2]*u;
                v[2*n-2*i-1]=v[2*i+1]*u;
                v[2*n-2*i]=v[2*i+2]*u;
                }
four1(v,n,-1);
for (i=0;i<n;i++)               v[i]=v[2*i+1];
for (i=n-1;i>=0;i--)    y[i]= (double)v[i];
for (i=0;i<n;i++)       y[i]=y[i] + pas * (double)(i-n/2);
}



/* -------spectre de fréquence du tableau double y[]-----------*/

void spectreRMS(double *y,int n)
{float *v;int i;float ym=0;
v=(float *)y;
// souslin(y,n);
for (i=0;i<n;i++)               {v[i]=(float)y[i];ym=ym+v[i];}
ym=ym/(float)n;
for (i=n-1;i>=0;i--)            {v[2*i+1]=(v[i]-ym)/(float)n;v[2*i+2]=0;}
four1(v,n,1);
for (i=0;i<n;i++)               v[i]=sqrt(2*(v[2*i+1]*v[2*i+1]+v[2*i+2]*v[2*i+2]));
for (i=n-1;i>=0;i--)    y[i]=(double)v[i];
}




void spectre(double *y,int n)
{float *v;int i;float ym=0;
v=(float *)y;
 souslin(y,n);
for (i=0;i<n;i++)               {v[i]=(float)y[i];ym=ym+v[i];}
ym=ym/(float)n;
for (i=n-1;i>=0;i--)            {v[2*i+1]=(v[i]-ym)/(float)n;v[2*i+2]=0;}
four1(v,n,1);
for (i=0;i<n;i++)               v[i]=v[2*i+1]*v[2*i+1]+v[2*i+2]*v[2*i+2];
for (i=n-1;i>=0;i--)    y[i]=(double)v[i];
}



/* ----------------- transformée de Fourier complexe ---------*/

void four1(float *data,int nn,int isign)   /* data tableau de 1 à 2*nn  */
/* isign=1 --> Tf directe   data = real-imag-real-imag ........ */

                /* isign=-1 --> Tf inverse    */
{
        int n,mmax,m,j,istep,i;
        double wtemp,wr,wpr,wpi,wi,theta;
        float tempr,tempi;

        n=nn << 1;
        j=1;
        for (i=1;i<n;i+=2) {
        
                if (j > i) {
                        SWAP(data[j],data[i]);
                        SWAP(data[j+1],data[i+1]);
                }
                m=n >> 1;
                while (m >= 2 && j > m) {
                        j -= m;
                        m >>= 1;
                }
                j += m;
        }
        mmax=2;
        while (n > mmax) {
                istep=2*mmax;
                theta=6.28318530717959/(isign*mmax);
                wtemp=sin(0.5*theta);
                wpr = -2.0*wtemp*wtemp;
                wpi=sin(theta);
                wr=1.0;
                wi=0.0;
                for (m=1;m<mmax;m+=2) {
                        for (i=m;i<=n;i+=istep) {
                                j=i+mmax;
                                tempr=wr*data[j]-wi*data[j+1];
                                tempi=wr*data[j+1]+wi*data[j];
                                data[j]=data[i]-tempr;
                                data[j+1]=data[i+1]-tempi;
                                data[i] += tempr;
                                data[i+1] += tempi;
                        }
                        wr=(wtemp=wr)*wpr-wi*wpi+wr;
                        wi=wi*wpr+wtemp*wpi+wi;
                }
                mmax=istep;
        }
}

#undef SWAP