bluestein.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, 2006, 2007 Max-Planck-Society
00027  *  \author Martin Reinecke
00028  */
00029 
00030 #include <math.h>
00031 #include <stdlib.h>
00032 #include "fftpack.h"
00033 #include "bluestein.h"
00034 
00035 /* returns the sum of all prime factors of n */
00036 int prime_factor_sum (int n)
00037   {
00038   int result=0,x,limit,tmp;
00039   while (((tmp=(n>>1))<<1)==n)
00040     { result+=2; n=tmp; }
00041 
00042   limit=sqrt(n+0.01);
00043   for (x=3; x<=limit; x+=2)
00044   while ((tmp=(n/x))*x==n)
00045     {
00046     result+=x;
00047     n=tmp;
00048     limit=sqrt(n+0.01);
00049     }
00050   if (n>1) result+=n;
00051 
00052   return result;
00053   }
00054 
00055 /* returns the smallest composite of 2, 3 and 5 which is >= n */
00056 static int good_size(int n)
00057   {
00058   int maxfactors=1, i, j, k, f2=1, f3, f5, bestfac, guessfac;
00059   while ((n>>maxfactors)>0)
00060     ++maxfactors;
00061   bestfac=1<<maxfactors;
00062 
00063   for (i=0; i<maxfactors; ++i)
00064     {
00065     f3=1;
00066     for (j=0; j<maxfactors-i; ++j)
00067       {
00068       if (f2*f3>bestfac) break;
00069       f5=1;
00070       for (k=0; k<maxfactors-i-j; ++k)
00071         {
00072         guessfac = f2*f3*f5;
00073         if (guessfac>=bestfac) break;
00074         if ((guessfac>=n) && (guessfac<bestfac))
00075           bestfac=guessfac;
00076         f5*=5;
00077         }
00078       f3*=3;
00079       }
00080     f2*=2;
00081     }
00082   return bestfac;
00083   }
00084 
00085 void bluestein_i (int n, double **tstorage)
00086   {
00087   static const double pi=3.14159265358979323846;
00088   int n2=good_size(n*2-1);
00089   int m, coeff;
00090   double angle, xn2;
00091   double *bk, *bkf, *work;
00092   double pibyn=pi/n;
00093   *tstorage = (double *)malloc (sizeof(double)*(1+2*n+8*n2+15));
00094   ((int *)(*tstorage))[0]=n2;
00095   bk  = *tstorage+1;
00096   bkf = *tstorage+1+2*n;
00097   work= *tstorage+1+2*(n+n2);
00098 
00099 /* initialize b_k */
00100   bk[0] = 1;
00101   bk[1] = 0;
00102 
00103   coeff=0;
00104   for (m=1; m<n; ++m)
00105     {
00106     coeff+=2*m-1;
00107     if (coeff>=2*n) coeff-=2*n;
00108     angle = pibyn*coeff;
00109     bk[2*m] = cos(angle);
00110     bk[2*m+1] = sin(angle);
00111     }
00112 
00113 /* initialize the zero-padded, Fourier transformed b_k. Add normalisation.  */
00114   xn2 = 1./n2;
00115   bkf[0] = bk[0]*xn2;
00116   bkf[1] = bk[1]*xn2;
00117   for (m=2; m<2*n; m+=2)
00118     {
00119     bkf[m]   = bkf[2*n2-m]   = bk[m]   *xn2;
00120     bkf[m+1] = bkf[2*n2-m+1] = bk[m+1] *xn2;
00121     }
00122   for (m=2*n;m<=(2*n2-2*n+1);++m)
00123     bkf[m]=0.;
00124   cffti (n2,work);
00125   cfftf (n2,bkf,work);
00126   }
00127 
00128 void bluestein (int n, double *data, double *tstorage, int isign)
00129   {
00130   int n2=*((int *)tstorage);
00131   int m;
00132   double *bk, *bkf, *akf, *work;
00133   bk  = tstorage+1;
00134   bkf = tstorage+1+2*n;
00135   work= tstorage+1+2*(n+n2);
00136   akf = tstorage+1+2*n+6*n2+15;  
00137 
00138 /* initialize a_k and FFT it */
00139   if (isign>0)
00140     for (m=0; m<2*n; m+=2)
00141       {
00142       akf[m]   = data[m]*bk[m]   - data[m+1]*bk[m+1];
00143       akf[m+1] = data[m]*bk[m+1] + data[m+1]*bk[m];
00144       }
00145   else
00146     for (m=0; m<2*n; m+=2)
00147       {
00148       akf[m]   = data[m]*bk[m]   + data[m+1]*bk[m+1];
00149       akf[m+1] =-data[m]*bk[m+1] + data[m+1]*bk[m];
00150       }
00151   for (m=2*n; m<2*n2; ++m)
00152     akf[m]=0;
00153 
00154   cfftf (n2,akf,work);
00155 
00156 /* do the convolution */
00157   if (isign>0)
00158     for (m=0; m<2*n2; m+=2)
00159       {
00160       double im = -akf[m]*bkf[m+1] + akf[m+1]*bkf[m];
00161       akf[m  ]  =  akf[m]*bkf[m]   + akf[m+1]*bkf[m+1];
00162       akf[m+1]  = im;
00163       }
00164   else
00165     for (m=0; m<2*n2; m+=2)
00166       {
00167       double im = akf[m]*bkf[m+1] + akf[m+1]*bkf[m];
00168       akf[m  ]  = akf[m]*bkf[m]   - akf[m+1]*bkf[m+1];
00169       akf[m+1]  = im;
00170       }
00171 
00172 
00173 /* inverse FFT */
00174   cfftb (n2,akf,work);
00175 
00176 /* multiply by b_k* */
00177   if (isign>0)
00178     for (m=0; m<2*n; m+=2)
00179       {
00180       data[m]   = bk[m]  *akf[m] - bk[m+1]*akf[m+1];
00181       data[m+1] = bk[m+1]*akf[m] + bk[m]  *akf[m+1];
00182       }
00183   else
00184     for (m=0; m<2*n; m+=2)
00185       {
00186       data[m]   = bk[m]  *akf[m] + bk[m+1]*akf[m+1];
00187       data[m+1] =-bk[m+1]*akf[m] + bk[m]  *akf[m+1];
00188       }
00189   }