annotate src/ext/kissfft/kissfft.hh @ 196:da283326bcd3 tip master

Update plugin versions in RDF
author Chris Cannam <cannam@all-day-breakfast.com>
date Fri, 28 Feb 2020 09:43:02 +0000
parents 5ed6e970541b
children
rev   line source
c@174 1 #ifndef KISSFFT_CLASS_HH
c@174 2 #define KISSFFT_CLASS_HH
c@174 3 #include <complex>
c@174 4 #include <vector>
c@174 5
c@174 6 namespace kissfft_utils {
c@174 7
c@174 8 template <typename T_scalar>
c@174 9 struct traits
c@174 10 {
c@174 11 typedef T_scalar scalar_type;
c@174 12 typedef std::complex<scalar_type> cpx_type;
c@174 13 void fill_twiddles( std::complex<T_scalar> * dst ,int nfft,bool inverse)
c@174 14 {
c@174 15 T_scalar phinc = (inverse?2:-2)* acos( (T_scalar) -1) / nfft;
c@174 16 for (int i=0;i<nfft;++i)
c@174 17 dst[i] = exp( std::complex<T_scalar>(0,i*phinc) );
c@174 18 }
c@174 19
c@174 20 void prepare(
c@174 21 std::vector< std::complex<T_scalar> > & dst,
c@174 22 int nfft,bool inverse,
c@174 23 std::vector<int> & stageRadix,
c@174 24 std::vector<int> & stageRemainder )
c@174 25 {
c@174 26 _twiddles.resize(nfft);
c@174 27 fill_twiddles( &_twiddles[0],nfft,inverse);
c@174 28 dst = _twiddles;
c@174 29
c@174 30 //factorize
c@174 31 //start factoring out 4's, then 2's, then 3,5,7,9,...
c@174 32 int n= nfft;
c@174 33 int p=4;
c@174 34 do {
c@174 35 while (n % p) {
c@174 36 switch (p) {
c@174 37 case 4: p = 2; break;
c@174 38 case 2: p = 3; break;
c@174 39 default: p += 2; break;
c@174 40 }
c@174 41 if (p*p>n)
c@174 42 p=n;// no more factors
c@174 43 }
c@174 44 n /= p;
c@174 45 stageRadix.push_back(p);
c@174 46 stageRemainder.push_back(n);
c@174 47 }while(n>1);
c@174 48 }
c@174 49 std::vector<cpx_type> _twiddles;
c@174 50
c@174 51
c@174 52 const cpx_type twiddle(int i) { return _twiddles[i]; }
c@174 53 };
c@174 54
c@174 55 }
c@174 56
c@174 57 template <typename T_Scalar,
c@174 58 typename T_traits=kissfft_utils::traits<T_Scalar>
c@174 59 >
c@174 60 class kissfft
c@174 61 {
c@174 62 public:
c@174 63 typedef T_traits traits_type;
c@174 64 typedef typename traits_type::scalar_type scalar_type;
c@174 65 typedef typename traits_type::cpx_type cpx_type;
c@174 66
c@174 67 kissfft(int nfft,bool inverse,const traits_type & traits=traits_type() )
c@174 68 :_nfft(nfft),_inverse(inverse),_traits(traits)
c@174 69 {
c@174 70 _traits.prepare(_twiddles, _nfft,_inverse ,_stageRadix, _stageRemainder);
c@174 71 }
c@174 72
c@174 73 void transform(const cpx_type * src , cpx_type * dst)
c@174 74 {
c@174 75 kf_work(0, dst, src, 1,1);
c@174 76 }
c@174 77
c@174 78 private:
c@174 79 void kf_work( int stage,cpx_type * Fout, const cpx_type * f, size_t fstride,size_t in_stride)
c@174 80 {
c@174 81 int p = _stageRadix[stage];
c@174 82 int m = _stageRemainder[stage];
c@174 83 cpx_type * Fout_beg = Fout;
c@174 84 cpx_type * Fout_end = Fout + p*m;
c@174 85
c@174 86 if (m==1) {
c@174 87 do{
c@174 88 *Fout = *f;
c@174 89 f += fstride*in_stride;
c@174 90 }while(++Fout != Fout_end );
c@174 91 }else{
c@174 92 do{
c@174 93 // recursive call:
c@174 94 // DFT of size m*p performed by doing
c@174 95 // p instances of smaller DFTs of size m,
c@174 96 // each one takes a decimated version of the input
c@174 97 kf_work(stage+1, Fout , f, fstride*p,in_stride);
c@174 98 f += fstride*in_stride;
c@174 99 }while( (Fout += m) != Fout_end );
c@174 100 }
c@174 101
c@174 102 Fout=Fout_beg;
c@174 103
c@174 104 // recombine the p smaller DFTs
c@174 105 switch (p) {
c@174 106 case 2: kf_bfly2(Fout,fstride,m); break;
c@174 107 case 3: kf_bfly3(Fout,fstride,m); break;
c@174 108 case 4: kf_bfly4(Fout,fstride,m); break;
c@174 109 case 5: kf_bfly5(Fout,fstride,m); break;
c@174 110 default: kf_bfly_generic(Fout,fstride,m,p); break;
c@174 111 }
c@174 112 }
c@174 113
c@174 114 // these were #define macros in the original kiss_fft
c@174 115 void C_ADD( cpx_type & c,const cpx_type & a,const cpx_type & b) { c=a+b;}
c@174 116 void C_MUL( cpx_type & c,const cpx_type & a,const cpx_type & b) { c=a*b;}
c@174 117 void C_SUB( cpx_type & c,const cpx_type & a,const cpx_type & b) { c=a-b;}
c@174 118 void C_ADDTO( cpx_type & c,const cpx_type & a) { c+=a;}
c@174 119 void C_FIXDIV( cpx_type & ,int ) {} // NO-OP for float types
c@174 120 scalar_type S_MUL( const scalar_type & a,const scalar_type & b) { return a*b;}
c@174 121 scalar_type HALF_OF( const scalar_type & a) { return a*.5;}
c@174 122 void C_MULBYSCALAR(cpx_type & c,const scalar_type & a) {c*=a;}
c@174 123
c@174 124 void kf_bfly2( cpx_type * Fout, const size_t fstride, int m)
c@174 125 {
c@174 126 for (int k=0;k<m;++k) {
c@174 127 cpx_type t = Fout[m+k] * _traits.twiddle(k*fstride);
c@174 128 Fout[m+k] = Fout[k] - t;
c@174 129 Fout[k] += t;
c@174 130 }
c@174 131 }
c@174 132
c@174 133 void kf_bfly4( cpx_type * Fout, const size_t fstride, const size_t m)
c@174 134 {
c@174 135 cpx_type scratch[7];
c@174 136 int negative_if_inverse = _inverse * -2 +1;
c@174 137 for (size_t k=0;k<m;++k) {
c@174 138 scratch[0] = Fout[k+m] * _traits.twiddle(k*fstride);
c@174 139 scratch[1] = Fout[k+2*m] * _traits.twiddle(k*fstride*2);
c@174 140 scratch[2] = Fout[k+3*m] * _traits.twiddle(k*fstride*3);
c@174 141 scratch[5] = Fout[k] - scratch[1];
c@174 142
c@174 143 Fout[k] += scratch[1];
c@174 144 scratch[3] = scratch[0] + scratch[2];
c@174 145 scratch[4] = scratch[0] - scratch[2];
c@174 146 scratch[4] = cpx_type( scratch[4].imag()*negative_if_inverse , -scratch[4].real()* negative_if_inverse );
c@174 147
c@174 148 Fout[k+2*m] = Fout[k] - scratch[3];
c@174 149 Fout[k] += scratch[3];
c@174 150 Fout[k+m] = scratch[5] + scratch[4];
c@174 151 Fout[k+3*m] = scratch[5] - scratch[4];
c@174 152 }
c@174 153 }
c@174 154
c@174 155 void kf_bfly3( cpx_type * Fout, const size_t fstride, const size_t m)
c@174 156 {
c@174 157 size_t k=m;
c@174 158 const size_t m2 = 2*m;
c@174 159 cpx_type *tw1,*tw2;
c@174 160 cpx_type scratch[5];
c@174 161 cpx_type epi3;
c@174 162 epi3 = _twiddles[fstride*m];
c@174 163
c@174 164 tw1=tw2=&_twiddles[0];
c@174 165
c@174 166 do{
c@174 167 C_FIXDIV(*Fout,3); C_FIXDIV(Fout[m],3); C_FIXDIV(Fout[m2],3);
c@174 168
c@174 169 C_MUL(scratch[1],Fout[m] , *tw1);
c@174 170 C_MUL(scratch[2],Fout[m2] , *tw2);
c@174 171
c@174 172 C_ADD(scratch[3],scratch[1],scratch[2]);
c@174 173 C_SUB(scratch[0],scratch[1],scratch[2]);
c@174 174 tw1 += fstride;
c@174 175 tw2 += fstride*2;
c@174 176
c@174 177 Fout[m] = cpx_type( Fout->real() - HALF_OF(scratch[3].real() ) , Fout->imag() - HALF_OF(scratch[3].imag() ) );
c@174 178
c@174 179 C_MULBYSCALAR( scratch[0] , epi3.imag() );
c@174 180
c@174 181 C_ADDTO(*Fout,scratch[3]);
c@174 182
c@174 183 Fout[m2] = cpx_type( Fout[m].real() + scratch[0].imag() , Fout[m].imag() - scratch[0].real() );
c@174 184
c@174 185 C_ADDTO( Fout[m] , cpx_type( -scratch[0].imag(),scratch[0].real() ) );
c@174 186 ++Fout;
c@174 187 }while(--k);
c@174 188 }
c@174 189
c@174 190 void kf_bfly5( cpx_type * Fout, const size_t fstride, const size_t m)
c@174 191 {
c@174 192 cpx_type *Fout0,*Fout1,*Fout2,*Fout3,*Fout4;
c@174 193 size_t u;
c@174 194 cpx_type scratch[13];
c@174 195 cpx_type * twiddles = &_twiddles[0];
c@174 196 cpx_type *tw;
c@174 197 cpx_type ya,yb;
c@174 198 ya = twiddles[fstride*m];
c@174 199 yb = twiddles[fstride*2*m];
c@174 200
c@174 201 Fout0=Fout;
c@174 202 Fout1=Fout0+m;
c@174 203 Fout2=Fout0+2*m;
c@174 204 Fout3=Fout0+3*m;
c@174 205 Fout4=Fout0+4*m;
c@174 206
c@174 207 tw=twiddles;
c@174 208 for ( u=0; u<m; ++u ) {
c@174 209 C_FIXDIV( *Fout0,5); C_FIXDIV( *Fout1,5); C_FIXDIV( *Fout2,5); C_FIXDIV( *Fout3,5); C_FIXDIV( *Fout4,5);
c@174 210 scratch[0] = *Fout0;
c@174 211
c@174 212 C_MUL(scratch[1] ,*Fout1, tw[u*fstride]);
c@174 213 C_MUL(scratch[2] ,*Fout2, tw[2*u*fstride]);
c@174 214 C_MUL(scratch[3] ,*Fout3, tw[3*u*fstride]);
c@174 215 C_MUL(scratch[4] ,*Fout4, tw[4*u*fstride]);
c@174 216
c@174 217 C_ADD( scratch[7],scratch[1],scratch[4]);
c@174 218 C_SUB( scratch[10],scratch[1],scratch[4]);
c@174 219 C_ADD( scratch[8],scratch[2],scratch[3]);
c@174 220 C_SUB( scratch[9],scratch[2],scratch[3]);
c@174 221
c@174 222 C_ADDTO( *Fout0, scratch[7]);
c@174 223 C_ADDTO( *Fout0, scratch[8]);
c@174 224
c@174 225 scratch[5] = scratch[0] + cpx_type(
c@174 226 S_MUL(scratch[7].real(),ya.real() ) + S_MUL(scratch[8].real() ,yb.real() ),
c@174 227 S_MUL(scratch[7].imag(),ya.real()) + S_MUL(scratch[8].imag(),yb.real())
c@174 228 );
c@174 229
c@174 230 scratch[6] = cpx_type(
c@174 231 S_MUL(scratch[10].imag(),ya.imag()) + S_MUL(scratch[9].imag(),yb.imag()),
c@174 232 -S_MUL(scratch[10].real(),ya.imag()) - S_MUL(scratch[9].real(),yb.imag())
c@174 233 );
c@174 234
c@174 235 C_SUB(*Fout1,scratch[5],scratch[6]);
c@174 236 C_ADD(*Fout4,scratch[5],scratch[6]);
c@174 237
c@174 238 scratch[11] = scratch[0] +
c@174 239 cpx_type(
c@174 240 S_MUL(scratch[7].real(),yb.real()) + S_MUL(scratch[8].real(),ya.real()),
c@174 241 S_MUL(scratch[7].imag(),yb.real()) + S_MUL(scratch[8].imag(),ya.real())
c@174 242 );
c@174 243
c@174 244 scratch[12] = cpx_type(
c@174 245 -S_MUL(scratch[10].imag(),yb.imag()) + S_MUL(scratch[9].imag(),ya.imag()),
c@174 246 S_MUL(scratch[10].real(),yb.imag()) - S_MUL(scratch[9].real(),ya.imag())
c@174 247 );
c@174 248
c@174 249 C_ADD(*Fout2,scratch[11],scratch[12]);
c@174 250 C_SUB(*Fout3,scratch[11],scratch[12]);
c@174 251
c@174 252 ++Fout0;++Fout1;++Fout2;++Fout3;++Fout4;
c@174 253 }
c@174 254 }
c@174 255
c@174 256 /* perform the butterfly for one stage of a mixed radix FFT */
c@174 257 void kf_bfly_generic(
c@174 258 cpx_type * Fout,
c@174 259 const size_t fstride,
c@174 260 int m,
c@174 261 int p
c@174 262 )
c@174 263 {
c@174 264 int u,k,q1,q;
c@174 265 cpx_type * twiddles = &_twiddles[0];
c@174 266 cpx_type t;
c@174 267 int Norig = _nfft;
c@174 268 cpx_type scratchbuf[p];
c@174 269
c@174 270 for ( u=0; u<m; ++u ) {
c@174 271 k=u;
c@174 272 for ( q1=0 ; q1<p ; ++q1 ) {
c@174 273 scratchbuf[q1] = Fout[ k ];
c@174 274 C_FIXDIV(scratchbuf[q1],p);
c@174 275 k += m;
c@174 276 }
c@174 277
c@174 278 k=u;
c@174 279 for ( q1=0 ; q1<p ; ++q1 ) {
c@174 280 int twidx=0;
c@174 281 Fout[ k ] = scratchbuf[0];
c@174 282 for (q=1;q<p;++q ) {
c@174 283 twidx += fstride * k;
c@174 284 if (twidx>=Norig) twidx-=Norig;
c@174 285 C_MUL(t,scratchbuf[q] , twiddles[twidx] );
c@174 286 C_ADDTO( Fout[ k ] ,t);
c@174 287 }
c@174 288 k += m;
c@174 289 }
c@174 290 }
c@174 291 }
c@174 292
c@174 293 int _nfft;
c@174 294 bool _inverse;
c@174 295 std::vector<cpx_type> _twiddles;
c@174 296 std::vector<int> _stageRadix;
c@174 297 std::vector<int> _stageRemainder;
c@174 298 traits_type _traits;
c@174 299 };
c@174 300 #endif