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