doku newon / DokuFFTPACK

widely known free FFT library.

Dependents:   Peach_AudioChannelDividerAndCompensator

bluestein.c@0:d7215946c769, 2015-10-11 (annotated)

Committer:
dokunewon
Date:
Sun Oct 11 07:20:35 2015 +0000
Revision:
0:d7215946c769
First version.

Who changed what in which revision?

UserRevisionLine numberNew contents of line
dokunewon 0:d7215946c769 1 /*
dokunewon 0:d7215946c769 2 * This file is part of libfftpack.
dokunewon 0:d7215946c769 3 *
dokunewon 0:d7215946c769 4 * libfftpack is free software; you can redistribute it and/or modify
dokunewon 0:d7215946c769 5 * it under the terms of the GNU General Public License as published by
dokunewon 0:d7215946c769 6 * the Free Software Foundation; either version 2 of the License, or
dokunewon 0:d7215946c769 7 * (at your option) any later version.
dokunewon 0:d7215946c769 8 *
dokunewon 0:d7215946c769 9 * libfftpack is distributed in the hope that it will be useful,
dokunewon 0:d7215946c769 10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
dokunewon 0:d7215946c769 11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
dokunewon 0:d7215946c769 12 * GNU General Public License for more details.
dokunewon 0:d7215946c769 13 *
dokunewon 0:d7215946c769 14 * You should have received a copy of the GNU General Public License
dokunewon 0:d7215946c769 15 * along with libfftpack; if not, write to the Free Software
dokunewon 0:d7215946c769 16 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
dokunewon 0:d7215946c769 17 */
dokunewon 0:d7215946c769 18
dokunewon 0:d7215946c769 19 /*
dokunewon 0:d7215946c769 20 * libfftpack is being developed at the Max-Planck-Institut fuer Astrophysik
dokunewon 0:d7215946c769 21 * and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt
dokunewon 0:d7215946c769 22 * (DLR).
dokunewon 0:d7215946c769 23 */
dokunewon 0:d7215946c769 24
dokunewon 0:d7215946c769 25 /*
dokunewon 0:d7215946c769 26 * Copyright (C) 2005, 2006, 2007, 2008 Max-Planck-Society
dokunewon 0:d7215946c769 27 * \author Martin Reinecke
dokunewon 0:d7215946c769 28 */
dokunewon 0:d7215946c769 29
dokunewon 0:d7215946c769 30 #include <math.h>
dokunewon 0:d7215946c769 31 #include <stdlib.h>
dokunewon 0:d7215946c769 32 #include "mytype.h"
dokunewon 0:d7215946c769 33 #include "fftpack.h"
dokunewon 0:d7215946c769 34 #include "bluestein.h"
dokunewon 0:d7215946c769 35
dokunewon 0:d7215946c769 36 /* returns the sum of all prime factors of n */
dokunewon 0:d7215946c769 37 size_t prime_factor_sum (size_t n)
dokunewon 0:d7215946c769 38 {
dokunewon 0:d7215946c769 39 size_t result=0,x,limit,tmp;
dokunewon 0:d7215946c769 40 while (((tmp=(n>>1))<<1)==n)
dokunewon 0:d7215946c769 41 { result+=2; n=tmp; }
dokunewon 0:d7215946c769 42
dokunewon 0:d7215946c769 43 limit=(size_t)sqrt(n+0.01);
dokunewon 0:d7215946c769 44 for (x=3; x<=limit; x+=2)
dokunewon 0:d7215946c769 45 while ((tmp=(n/x))*x==n)
dokunewon 0:d7215946c769 46 {
dokunewon 0:d7215946c769 47 result+=x;
dokunewon 0:d7215946c769 48 n=tmp;
dokunewon 0:d7215946c769 49 limit=(size_t)sqrt(n+0.01);
dokunewon 0:d7215946c769 50 }
dokunewon 0:d7215946c769 51 if (n>1) result+=n;
dokunewon 0:d7215946c769 52
dokunewon 0:d7215946c769 53 return result;
dokunewon 0:d7215946c769 54 }
dokunewon 0:d7215946c769 55
dokunewon 0:d7215946c769 56 /* returns the smallest composite of 2, 3 and 5 which is >= n */
dokunewon 0:d7215946c769 57 static size_t good_size(size_t n)
dokunewon 0:d7215946c769 58 {
dokunewon 0:d7215946c769 59 size_t maxfactors=1, i, j, k, f2=1, f3, f5, bestfac, guessfac;
dokunewon 0:d7215946c769 60 while ((n>>maxfactors)>0)
dokunewon 0:d7215946c769 61 ++maxfactors;
dokunewon 0:d7215946c769 62 bestfac=1<<maxfactors;
dokunewon 0:d7215946c769 63
dokunewon 0:d7215946c769 64 for (i=0; i<maxfactors; ++i)
dokunewon 0:d7215946c769 65 {
dokunewon 0:d7215946c769 66 f3=1;
dokunewon 0:d7215946c769 67 for (j=0; j<maxfactors-i; ++j)
dokunewon 0:d7215946c769 68 {
dokunewon 0:d7215946c769 69 if (f2*f3>bestfac) break;
dokunewon 0:d7215946c769 70 f5=1;
dokunewon 0:d7215946c769 71 for (k=0; k<maxfactors-i-j; ++k)
dokunewon 0:d7215946c769 72 {
dokunewon 0:d7215946c769 73 guessfac = f2*f3*f5;
dokunewon 0:d7215946c769 74 if (guessfac>=bestfac) break;
dokunewon 0:d7215946c769 75 if ((guessfac>=n) && (guessfac<bestfac))
dokunewon 0:d7215946c769 76 bestfac=guessfac;
dokunewon 0:d7215946c769 77 f5*=5;
dokunewon 0:d7215946c769 78 }
dokunewon 0:d7215946c769 79 f3*=3;
dokunewon 0:d7215946c769 80 }
dokunewon 0:d7215946c769 81 f2*=2;
dokunewon 0:d7215946c769 82 }
dokunewon 0:d7215946c769 83 return bestfac;
dokunewon 0:d7215946c769 84 }
dokunewon 0:d7215946c769 85
dokunewon 0:d7215946c769 86 void bluestein_i (size_t n, FLOAT **tstorage)
dokunewon 0:d7215946c769 87 {
dokunewon 0:d7215946c769 88 static const FLOAT pi=3.14159265358979323846;
dokunewon 0:d7215946c769 89 size_t n2=good_size(n*2-1);
dokunewon 0:d7215946c769 90 size_t m, coeff;
dokunewon 0:d7215946c769 91 FLOAT angle, xn2;
dokunewon 0:d7215946c769 92 FLOAT *bk, *bkf, *work;
dokunewon 0:d7215946c769 93 FLOAT pibyn=pi/n;
dokunewon 0:d7215946c769 94 *tstorage = RALLOC(FLOAT,2+2*n+8*n2+16);
dokunewon 0:d7215946c769 95 ((size_t *)(*tstorage))[0]=n2;
dokunewon 0:d7215946c769 96 bk = *tstorage+2;
dokunewon 0:d7215946c769 97 bkf = *tstorage+2+2*n;
dokunewon 0:d7215946c769 98 work= *tstorage+2+2*(n+n2);
dokunewon 0:d7215946c769 99
dokunewon 0:d7215946c769 100 /* initialize b_k */
dokunewon 0:d7215946c769 101 bk[0] = 1;
dokunewon 0:d7215946c769 102 bk[1] = 0;
dokunewon 0:d7215946c769 103
dokunewon 0:d7215946c769 104 coeff=0;
dokunewon 0:d7215946c769 105 for (m=1; m<n; ++m)
dokunewon 0:d7215946c769 106 {
dokunewon 0:d7215946c769 107 coeff+=2*m-1;
dokunewon 0:d7215946c769 108 if (coeff>=2*n) coeff-=2*n;
dokunewon 0:d7215946c769 109 angle = pibyn*coeff;
dokunewon 0:d7215946c769 110 bk[2*m] = cos(angle);
dokunewon 0:d7215946c769 111 bk[2*m+1] = sin(angle);
dokunewon 0:d7215946c769 112 }
dokunewon 0:d7215946c769 113
dokunewon 0:d7215946c769 114 /* initialize the zero-padded, Fourier transformed b_k. Add normalisation. */
dokunewon 0:d7215946c769 115 xn2 = 1./n2;
dokunewon 0:d7215946c769 116 bkf[0] = bk[0]*xn2;
dokunewon 0:d7215946c769 117 bkf[1] = bk[1]*xn2;
dokunewon 0:d7215946c769 118 for (m=2; m<2*n; m+=2)
dokunewon 0:d7215946c769 119 {
dokunewon 0:d7215946c769 120 bkf[m] = bkf[2*n2-m] = bk[m] *xn2;
dokunewon 0:d7215946c769 121 bkf[m+1] = bkf[2*n2-m+1] = bk[m+1] *xn2;
dokunewon 0:d7215946c769 122 }
dokunewon 0:d7215946c769 123 for (m=2*n;m<=(2*n2-2*n+1);++m)
dokunewon 0:d7215946c769 124 bkf[m]=0.;
dokunewon 0:d7215946c769 125 cffti (n2,work);
dokunewon 0:d7215946c769 126 cfftf (n2,bkf,work);
dokunewon 0:d7215946c769 127 }
dokunewon 0:d7215946c769 128
dokunewon 0:d7215946c769 129 void bluestein (size_t n, FLOAT *data, FLOAT *tstorage, int isign)
dokunewon 0:d7215946c769 130 {
dokunewon 0:d7215946c769 131 size_t n2=*((size_t *)tstorage);
dokunewon 0:d7215946c769 132 size_t m;
dokunewon 0:d7215946c769 133 FLOAT *bk, *bkf, *akf, *work;
dokunewon 0:d7215946c769 134 bk = tstorage+2;
dokunewon 0:d7215946c769 135 bkf = tstorage+2+2*n;
dokunewon 0:d7215946c769 136 work= tstorage+2+2*(n+n2);
dokunewon 0:d7215946c769 137 akf = tstorage+2+2*n+6*n2+16;
dokunewon 0:d7215946c769 138
dokunewon 0:d7215946c769 139 /* initialize a_k and FFT it */
dokunewon 0:d7215946c769 140 if (isign>0)
dokunewon 0:d7215946c769 141 for (m=0; m<2*n; m+=2)
dokunewon 0:d7215946c769 142 {
dokunewon 0:d7215946c769 143 akf[m] = data[m]*bk[m] - data[m+1]*bk[m+1];
dokunewon 0:d7215946c769 144 akf[m+1] = data[m]*bk[m+1] + data[m+1]*bk[m];
dokunewon 0:d7215946c769 145 }
dokunewon 0:d7215946c769 146 else
dokunewon 0:d7215946c769 147 for (m=0; m<2*n; m+=2)
dokunewon 0:d7215946c769 148 {
dokunewon 0:d7215946c769 149 akf[m] = data[m]*bk[m] + data[m+1]*bk[m+1];
dokunewon 0:d7215946c769 150 akf[m+1] =-data[m]*bk[m+1] + data[m+1]*bk[m];
dokunewon 0:d7215946c769 151 }
dokunewon 0:d7215946c769 152 for (m=2*n; m<2*n2; ++m)
dokunewon 0:d7215946c769 153 akf[m]=0;
dokunewon 0:d7215946c769 154
dokunewon 0:d7215946c769 155 cfftf (n2,akf,work);
dokunewon 0:d7215946c769 156
dokunewon 0:d7215946c769 157 /* do the convolution */
dokunewon 0:d7215946c769 158 if (isign>0)
dokunewon 0:d7215946c769 159 for (m=0; m<2*n2; m+=2)
dokunewon 0:d7215946c769 160 {
dokunewon 0:d7215946c769 161 FLOAT im = -akf[m]*bkf[m+1] + akf[m+1]*bkf[m];
dokunewon 0:d7215946c769 162 akf[m ] = akf[m]*bkf[m] + akf[m+1]*bkf[m+1];
dokunewon 0:d7215946c769 163 akf[m+1] = im;
dokunewon 0:d7215946c769 164 }
dokunewon 0:d7215946c769 165 else
dokunewon 0:d7215946c769 166 for (m=0; m<2*n2; m+=2)
dokunewon 0:d7215946c769 167 {
dokunewon 0:d7215946c769 168 FLOAT im = akf[m]*bkf[m+1] + akf[m+1]*bkf[m];
dokunewon 0:d7215946c769 169 akf[m ] = akf[m]*bkf[m] - akf[m+1]*bkf[m+1];
dokunewon 0:d7215946c769 170 akf[m+1] = im;
dokunewon 0:d7215946c769 171 }
dokunewon 0:d7215946c769 172
dokunewon 0:d7215946c769 173
dokunewon 0:d7215946c769 174 /* inverse FFT */
dokunewon 0:d7215946c769 175 cfftb (n2,akf,work);
dokunewon 0:d7215946c769 176
dokunewon 0:d7215946c769 177 /* multiply by b_k* */
dokunewon 0:d7215946c769 178 if (isign>0)
dokunewon 0:d7215946c769 179 for (m=0; m<2*n; m+=2)
dokunewon 0:d7215946c769 180 {
dokunewon 0:d7215946c769 181 data[m] = bk[m] *akf[m] - bk[m+1]*akf[m+1];
dokunewon 0:d7215946c769 182 data[m+1] = bk[m+1]*akf[m] + bk[m] *akf[m+1];
dokunewon 0:d7215946c769 183 }
dokunewon 0:d7215946c769 184 else
dokunewon 0:d7215946c769 185 for (m=0; m<2*n; m+=2)
dokunewon 0:d7215946c769 186 {
dokunewon 0:d7215946c769 187 data[m] = bk[m] *akf[m] + bk[m+1]*akf[m+1];
dokunewon 0:d7215946c769 188 data[m+1] =-bk[m+1]*akf[m] + bk[m] *akf[m+1];
dokunewon 0:d7215946c769 189 }
dokunewon 0:d7215946c769 190 }