#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