ls_fft.c

00001 /*
00002  *  This file is part of libfftpack.
00003  *
00004  *  libfftpack is free software; you can redistribute it and/or modify
00005  *  it under the terms of the GNU General Public License as published by
00006  *  the Free Software Foundation; either version 2 of the License, or
00007  *  (at your option) any later version.
00008  *
00009  *  libfftpack is distributed in the hope that it will be useful,
00010  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
00011  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00012  *  GNU General Public License for more details.
00013  *
00014  *  You should have received a copy of the GNU General Public License
00015  *  along with libfftpack; if not, write to the Free Software
00016  *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
00017  */
00018 
00019 /*
00020  *  libfftpack is being developed at the Max-Planck-Institut fuer Astrophysik
00021  *  and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt
00022  *  (DLR).
00023  */
00024 
00025 /*
00026  *  Copyright (C) 2005 Max-Planck-Society
00027  *  \author Martin Reinecke
00028  */
00029 
00030 #include <stdlib.h>
00031 #include <math.h>
00032 #include <string.h>
00033 #include "bluestein.h"
00034 #include "fftpack.h"
00035 #include "ls_fft.h"
00036 
00037 complex_plan make_complex_plan (int length)
00038   {
00039   complex_plan plan = (complex_plan) malloc(sizeof(complex_plan_i));
00040   int pfsum = prime_factor_sum(length);
00041   double comp1 = length*pfsum;
00042   double comp2 = 2*3*length*log(3.*length);
00043   plan->length=length;
00044   plan->bluestein = (comp2<comp1);
00045   if (plan->bluestein)
00046     bluestein_i (length,&(plan->work));
00047   else
00048     {
00049     plan->work=(double *)malloc((4*length+15)*sizeof(double));
00050     cffti(length, plan->work);
00051     }
00052   return plan;
00053   }
00054 
00055 void kill_complex_plan (complex_plan plan)
00056   {
00057   free(plan->work);
00058   free(plan);
00059   }
00060 
00061 void complex_plan_forward (complex_plan plan, double *data)
00062   {
00063   if (plan->bluestein)
00064     bluestein (plan->length, data, plan->work, -1);
00065   else
00066     cfftf (plan->length, data, plan->work);
00067   }
00068 
00069 void complex_plan_backward (complex_plan plan, double *data)
00070   {
00071   if (plan->bluestein)
00072     bluestein (plan->length, data, plan->work, 1);
00073   else
00074     cfftb (plan->length, data, plan->work);
00075   }
00076 
00077 
00078 real_plan make_real_plan (int length)
00079   {
00080   real_plan plan = (real_plan) malloc(sizeof(real_plan_i));
00081   int pfsum = prime_factor_sum(length);
00082   double comp1 = .5*length*pfsum;
00083   double comp2 = 2*3*length*log(3.*length);
00084   plan->length=length;
00085   plan->bluestein = (comp2<comp1);
00086   if (plan->bluestein)
00087     bluestein_i (length,&(plan->work));
00088   else
00089     {
00090     plan->work=(double *)malloc((2*length+15)*sizeof(double));
00091     rffti(length, plan->work);
00092     }
00093   return plan;
00094   }
00095 
00096 void kill_real_plan (real_plan plan)
00097   {
00098   free(plan->work);
00099   free(plan);
00100   }
00101 
00102 void real_plan_forward_fftpack (real_plan plan, double *data)
00103   {
00104   if (plan->bluestein)
00105     {
00106     int m;
00107     int n=plan->length;
00108     double *tmp = (double *)malloc(2*n*sizeof(double));
00109     for (m=0; m<n; ++m)
00110       {
00111       tmp[2*m] = data[m];
00112       tmp[2*m+1] = 0.;
00113       }
00114     bluestein(n,tmp,plan->work,-1);
00115     data[0] = tmp[0];
00116     memcpy (data+1, tmp+2, (n-1)*sizeof(double));
00117     free (tmp);
00118     }
00119   else
00120     rfftf (plan->length, data, plan->work);
00121   }
00122 
00123 static void fftpack2halfcomplex (double *data, int n)
00124   {
00125   int m;
00126   double *tmp = (double *)malloc(n*sizeof(double));
00127   tmp[0]=data[0];
00128   for (m=1; m<(n+1)/2; ++m)
00129     {
00130     tmp[m]=data[2*m-1];
00131     tmp[n-m]=data[2*m];
00132     }
00133   if (!(n&1))
00134     tmp[n/2]=data[n-1];
00135   memcpy (data,tmp,n*sizeof(double));
00136   free(tmp);
00137   }
00138 
00139 static void halfcomplex2fftpack (double *data, int n)
00140   {
00141   int m;
00142   double *tmp = (double *)malloc(n*sizeof(double));
00143   tmp[0]=data[0];
00144   for (m=1; m<(n+1)/2; ++m)
00145     {
00146     tmp[2*m-1]=data[m];
00147     tmp[2*m]=data[n-m];
00148     }
00149   if (!(n&1))
00150     tmp[n-1]=data[n/2];
00151   memcpy (data,tmp,n*sizeof(double));
00152   free(tmp);
00153   }
00154 
00155 void real_plan_forward_fftw (real_plan plan, double *data)
00156   {
00157   real_plan_forward_fftpack (plan, data);
00158   fftpack2halfcomplex (data,plan->length);
00159   }
00160 
00161 void real_plan_backward_fftpack (real_plan plan, double *data)
00162   {
00163   if (plan->bluestein)
00164     {
00165     int m;
00166     int n=plan->length;
00167     double *tmp = (double *)malloc(2*n*sizeof(double));
00168     tmp[0]=data[0];
00169     tmp[1]=0.;
00170     memcpy (tmp+2,data+1, (n-1)*sizeof(double));
00171     if ((n&1)==0) tmp[n+1]=0.;
00172     for (m=2; m<n; m+=2)
00173       {
00174       tmp[2*n-m]=tmp[m];
00175       tmp[2*n-m+1]=-tmp[m+1];
00176       }
00177     bluestein (n, tmp, plan->work, 1);
00178     for (m=0; m<n; ++m)
00179       data[m] = tmp[2*m];
00180     free (tmp);
00181     }
00182   else
00183     rfftb (plan->length, data, plan->work);
00184   }
00185 
00186 void real_plan_backward_fftw (real_plan plan, double *data)
00187   {
00188   halfcomplex2fftpack (data,plan->length);
00189   real_plan_backward_fftpack (plan, data);
00190   }
00191 
00192 void real_plan_forward_c (real_plan plan, double *data)
00193   {
00194   int m;
00195   int n=plan->length;
00196 
00197   if (plan->bluestein)
00198     {
00199     for (m=1; m<2*n; m+=2)
00200       data[m]=0;
00201     bluestein (plan->length, data, plan->work, -1);
00202     data[1]=0;
00203     for (m=2; m<n; m+=2)
00204       {
00205       double avg;
00206       avg = 0.5*(data[2*n-m]+data[m]);
00207       data[2*n-m] = data[m] = avg;
00208       avg = 0.5*(data[2*n-m+1]-data[m+1]);
00209       data[2*n-m+1] = avg;
00210       data[m+1] = -avg;
00211       }
00212     if ((n&1)==0) data[n+1] = 0.;
00213     }
00214   else
00215     {
00216     for (m=0; m<n; ++m) data[m+1] = data[2*m];
00217     rfftf (n, data+1, plan->work);
00218     data[0] = data[1];
00219     data[1] = 0;
00220     for (m=2; m<n; m+=2)
00221       {
00222       data[2*n-m]   =  data[m];
00223       data[2*n-m+1] = -data[m+1];
00224       }
00225     if ((n&1)==0) data[n+1] = 0.;
00226     }
00227   }
00228 
00229 void real_plan_backward_c (real_plan plan, double *data)
00230   {
00231   int m;
00232   int n=plan->length;
00233 
00234   if (plan->bluestein)
00235     {
00236     data[1]=0;
00237     for (m=2; m<n; m+=2)
00238       {
00239       double avg;
00240       avg = 0.5*(data[2*n-m]+data[m]);
00241       data[2*n-m] = data[m] = avg;
00242       avg = 0.5*(data[2*n-m+1]-data[m+1]);
00243       data[2*n-m+1] = avg;
00244       data[m+1] = -avg;
00245       }
00246     if ((n&1)==0) data[n+1] = 0.;
00247     bluestein (plan->length, data, plan->work, 1);
00248     for (m=1; m<2*n; m+=2)
00249       data[m]=0;
00250     }
00251   else
00252     {
00253     data[1] = data[0];
00254     rfftb (n, data+1, plan->work);
00255     for (m=n-1; m>=0; --m)
00256       {
00257       data[2*m]   = data[m+1];
00258       data[2*m+1] = 0.;
00259       }
00260     }
00261   }