/**************************************************************************
* fft.c
*
* makecint -A -dl fft.sl -c fft.c
**************************************************************************/
#include <math.h>
#include <stdio.h>
#define FFT_MAX 4096
#define G__FFTSL
/* #define G__STDFFT */
/*******************************************************************/
/* dft_cal.c 1996 1/8 M.Goto */
/*******************************************************************/
int dft_cal(double *x,double *re,double *im,int ndat)
{
int i,j,m;
double *s,*c;
double *rei,*imi;
double a,b;
double sums,sumc;
fprintf(stderr,"original data points=%d DFT data points=%d\n",ndat,ndat);
s = (double*)malloc(sizeof(double)*(ndat));
c = (double*)malloc(sizeof(double)*(ndat));
rei = (double*)malloc(sizeof(double)*(ndat));
imi = (double*)malloc(sizeof(double)*(ndat));
/* prepare sin,cos table and copy input data */
a=0 ;
b=3.141592*2/ndat;
for(i=0;i<ndat;i++) {
s[i] = sin(a);
c[i] = cos(a);
/* printf("%d a=%g cos=%g sin=%g\n",i,a,c[i],s[i]); */
a += b;
rei[i] = re[i];
imi[i] = im[i];
}
/* DFT calculation */
for(i=0;i<ndat;i++) {
sums=sumc=0.0;
for(j=0;j<ndat;j++) {
m = (j*i)%ndat;
sumc += c[m]*rei[j] + s[m]*imi[j];
sums += c[m]*imi[j] - s[m]*rei[j];
}
re[i] = sumc;
im[i] = sums;
}
free((void*)imi);
free((void*)rei);
free((void*)c);
free((void*)s);
return(0);
}
/*******************************************************************/
/* fft_cal.c 1987 1/8 M.Goto */
/*******************************************************************/
int fft_cal(double *x,double *re,double *im,int ndat)
{
int g,h,i,j,k,l,m,n,p,q;
/* double a,b,tmp,*pd,s[FFT_MAX],c[FFT_MAX]; */
double a,b,tmp,*s,*c;
/* ndat ----- must be n=2**m */
n=2;
m=1;
while (ndat > n) {
n *= 2 ;
m++ ;
}
if(ndat!=n) {
dft_cal(x,re,im,ndat);
return(0);
}
fprintf(stderr
,"original data points=%d FFT data points=%d m=%d n=2**m \n"
, ndat , n , m );
/* ndat=n/2; */
s = (double*)malloc(sizeof(double)*(n/2+1));
c = (double*)malloc(sizeof(double)*(n/2+1));
/* initialize sin cos table */
a=0 ;
b=3.141592*2/n;
for(i=0;i<=n/2;i++) {
s[i] = sin(a);
c[i] = cos(a);
/* c[i] = sqrt(1-s[i]*s[i]); */
a += b;
}
/* Butterfly */
l=n;
h=1;
for (g=1;g<=m;g++) {
l /= 2 ;
k=0 ;
for (q=1;q<=h;q++) {
p=0 ;
for (i=k;i<=l+k-1;i++) {
j=i+l ;
a = re[i] - re[j];
b = im[i] - im[j];
re[i] = re[i] + re[j];
im[i] = im[i] + im[j];
if (p == 0) {
re[j] = a;
im[j] = b;
}
else {
re[j] = a*c[p] + b*s[p];
im[j] = b*c[p] - a*s[p];
}
p=p+h;
}
k=k+l+l;
}
h=h+h;
}
free((void*)c);
free((void*)s);
/* Bit Reversal */
for (i=1;i<=n-1;i++) {
q=i;
k=n/2;
j=0;
for (l=1;l<=m;l++) {
j=j+k*(q%2); /* j=j+k*mod(q,2) */
q /= 2; /* q=int(q/2) */
k /= 2;
}
if (i < j) {
tmp= re[i];
re[i]=re[j];
re[j]=tmp;
tmp=im[i];
im[i]=im[j];
im[j]=tmp;
}
}
/* compen */
tmp= re[n/2-1];
re[n/2-1] = re[n/2+1] ;
re[n/2+1] = tmp;
tmp= im[n/2-1];
im[n/2-1]=im[n/2+1];
im[n/2+1]=tmp;
return(0);
}
/*******************************************************************/
/* fft.c 1987 1/8 M.Goto */
/*******************************************************************/
int fft(double *x,double *re,double *im,int ndat)
{
int i;
double dT,T,dfreq;
fft_cal(x,re,im,ndat);
dT= x[1]-x[0];
T = ndat * dT ;
dfreq = 1 / T ;
for (i=0;i<ndat;i++) {
x[i] = dfreq * i ;
#ifndef G__STDFFT
re[i] /= ndat;
im[i] /= ndat;
#endif
}
return(0);
}
/*******************************************************************/
/* ifft.c 1987 1/8 M.Goto */
/*******************************************************************/
int ifft(double *x,double *re,double *im,int ndat)
{
int i,j;
double dT,freq,dfreq,tmp;
fft_cal(x,re,im,ndat);
dfreq= x[1]-x[0];
freq = ndat * dfreq ;
dT = 1 / freq ;
for(i=0;i<ndat;i++) {
x[i] = dT * i;
}
/**************** Time reversal *********************/
for (i=1;i<ndat/2-1;i++) {
j= ndat - i ;
tmp = re[i];
re[i]=re[j];
#ifndef G__STDFFT
re[j]=tmp;
#else
re[j]=tmp/ndat;
#endif
tmp = im[i];
im[i]=im[j];
#ifndef G__STDFFT
im[j]=tmp;
#else
im[j]=tmp/ndat;
#endif
}
return(0);
}
/*******************************************************************/
/* spectrum.c 1987 1/8 M.Goto */
/*******************************************************************/
int spectrum(double *x,double *y,int ndat)
{
int i;
double *im;
im = (double*)calloc(ndat,sizeof(double));
fft(x,y,im,ndat);
#ifndef G__STDFFT
for (i=0;i<ndat;i++) y[i] = sqrt(y[i]*y[i] + im[i]*im[i]) * 1.41421356;
#else
double k = 1.41421356/ndat;
for (i=0;i<ndat;i++) y[i] = sqrt(y[i]*y[i] + im[i]*im[i]) * k;
#endif
free((void*)im);
return(0);
}
/*******************************************************************/
/* colleration.c 1987 1/8 M.Goto */
/*******************************************************************/
int colleration(double *x,double *y,int ndat)
{
/* auto colleration */
int i;
double *im;
double tmp;
spectrum(x,y,ndat);
im = (double*)calloc(ndat,sizeof(double));
ifft(x,y,im,ndat);
tmp=y[0];
for (i=0;i<ndat;i++) y[i] = y[i]/tmp;
return(0);
}
/*******************************************************************/
/* cepstrum.c 1987 1/8 M.Goto */
/*******************************************************************/
int cepstrum(double *x,double *y,int ndat)
{
int i;
spectrum(x,y,ndat);
for (i=0;i<ndat;i++) y[i] = log(y[i]);
spectrum(x,y,ndat);
return(0);
}