|
| 1 | +//see https://github.com/g200kg/Fft-asm.js |
| 2 | + |
| 3 | + |
| 4 | +module.exports = FftModule; |
| 5 | + |
| 6 | +function FftModule(sz,asm) { |
| 7 | + var i,j,k; |
| 8 | + this.sz=sz; |
| 9 | + this.bufsz=sz*32; |
| 10 | + if(this.bufsz<4096) |
| 11 | + this.bufsz=4096; |
| 12 | + this.heap=new ArrayBuffer(this.bufsz); |
| 13 | + this.foreign=new ArrayBuffer(4096); |
| 14 | + this.flt64=new Float64Array(this.heap); |
| 15 | + if(asm) |
| 16 | + this.fftasm=FftModuleAsm(window,this.foreign,this.heap); |
| 17 | + else |
| 18 | + this.fftasm=FftModuleNoAsm(window,this.foreign,this.heap); |
| 19 | + this.fftasm.setup(sz); |
| 20 | + var t,th=Math.PI; |
| 21 | + for(i=1,j=0;i<sz;i<<=1) { |
| 22 | + t=0.0; |
| 23 | + for(k=0;k<i;++k,j+=2) { |
| 24 | + t+=th; |
| 25 | + this.flt64[sz*2+j]=Math.cos(t); |
| 26 | + this.flt64[sz*2+j+1]=Math.sin(t); |
| 27 | + } |
| 28 | + th*=0.5; |
| 29 | + } |
| 30 | + this.bitrev=new Array(sz); |
| 31 | + this.bitrev[0]=0; |
| 32 | + this.bitrev[sz-1]=sz-1; |
| 33 | + for(j=1,i=0;j<sz-1;++j) { |
| 34 | + for(var k=sz>>1;k>(i^=k);k>>=1) |
| 35 | + ; |
| 36 | + this.bitrev[j]=i; |
| 37 | + } |
| 38 | + this.fft=function(real,imag,normalize,asm) { |
| 39 | + var i; |
| 40 | + for(i=0;i<this.sz;++i) { |
| 41 | + this.flt64[this.bitrev[i]]=real[i]; |
| 42 | + this.flt64[this.sz+this.bitrev[i]]=imag[i]; |
| 43 | + } |
| 44 | + this.fftasm.fft(normalize); |
| 45 | + for(i=0;i<this.sz;++i) { |
| 46 | + real[i]=this.flt64[i]; |
| 47 | + imag[i]=this.flt64[this.sz+i]; |
| 48 | + } |
| 49 | + } |
| 50 | + this.fftmag=function(real,imag) { |
| 51 | + var i; |
| 52 | + for(i=0;i<this.sz;++i) { |
| 53 | + this.flt64[this.bitrev[i]]=real[i]; |
| 54 | + this.flt64[this.sz+this.bitrev[i]]=imag[i]; |
| 55 | + } |
| 56 | + this.fftasm.fft(1); |
| 57 | + this.fftasm.mag(); |
| 58 | + for(i=0;i<this.sz;++i) { |
| 59 | + real[i]=this.flt64[i]; |
| 60 | + } |
| 61 | + } |
| 62 | +} |
| 63 | + |
| 64 | + |
| 65 | + |
| 66 | +// |
| 67 | +//FFT Module without "use asm" for comparison |
| 68 | +// |
| 69 | +function FftModuleNoAsm(stdlib, foreign, heap) { |
| 70 | + //"use asm"; |
| 71 | + var FLOAT64 = new stdlib.Float64Array(heap); |
| 72 | + var sqrt=stdlib.Math.sqrt; |
| 73 | + var sz=0; |
| 74 | + function setup(n) { |
| 75 | + n=n|0; |
| 76 | + sz=n; |
| 77 | + } |
| 78 | + function fft(normalize) { |
| 79 | + normalize=normalize|0; |
| 80 | + var m=0, mh=0, i=0, j=0, k=0; |
| 81 | + var wr=0.0, wi=0.0, xr=0.0, xi=0.0, w=0.0, r=0.0; |
| 82 | + var tcnt=0; |
| 83 | + for(mh=1;(m=mh<<1)<=(sz|0);mh=m) { |
| 84 | + for(i=0;(i|0)<(mh|0);i=(i+1)|0) { |
| 85 | + wr=+FLOAT64[(sz*16+tcnt)>>3]; |
| 86 | + wi=+FLOAT64[(sz*16+tcnt+8)>>3]; |
| 87 | + tcnt=(tcnt+16)|0; |
| 88 | + for(j=i;(j|0)<(sz|0);j=(j+m)|0) { |
| 89 | + k=(j+mh)|0; |
| 90 | + xr=wr*FLOAT64[k<<3>>3]-wi*FLOAT64[(sz+k)<<3>>3]; |
| 91 | + xi=wr*FLOAT64[(sz+k)<<3>>3]+wi*FLOAT64[k<<3>>3]; |
| 92 | + FLOAT64[k<<3>>3] = +FLOAT64[j<<3>>3]-xr; |
| 93 | + FLOAT64[(sz+k)<<3>>3] = +FLOAT64[(sz+j)<<3>>3]-xi; |
| 94 | + FLOAT64[j<<3>>3] = +FLOAT64[j<<3>>3] + xr; |
| 95 | + FLOAT64[(sz+j)<<3>>3] = +FLOAT64[(sz+j)<<3>>3] + xi; |
| 96 | + } |
| 97 | + } |
| 98 | + } |
| 99 | + if(normalize) { |
| 100 | + r = +(1.0/+(sz|0)); |
| 101 | + for(i=0;(i|0)<(sz|0);i=(i+1)|0) { |
| 102 | + FLOAT64[i<<3>>3]=+FLOAT64[i<<3>>3]*r; |
| 103 | + FLOAT64[(sz+i)<<3>>3]=+FLOAT64[(sz+i)<<3>>3]*r; |
| 104 | + } |
| 105 | + } |
| 106 | + } |
| 107 | + function mag() { |
| 108 | + var i=0,j=0; |
| 109 | + for(i=0;(i|0)<(sz|0);i=(i+1)|0) { |
| 110 | + j=(sz*2-i-1)|0; |
| 111 | + FLOAT64[i<<3>>3]=sqrt(FLOAT64[i<<3>>3]*FLOAT64[i<<3>>3]+FLOAT64[j<<3>>3]*FLOAT64[j<<3>>3]); |
| 112 | + } |
| 113 | + } |
| 114 | + return { |
| 115 | + setup:setup, |
| 116 | + fft:fft, |
| 117 | + mag:mag |
| 118 | + }; |
| 119 | +} |
| 120 | + |
| 121 | + |
| 122 | + |
| 123 | +// |
| 124 | +//FFT Module with "use asm" |
| 125 | +// |
| 126 | +function FftModuleAsm(stdlib, foreign, heap) { |
| 127 | + "use asm"; |
| 128 | + var FLOAT64 = new stdlib.Float64Array(heap); |
| 129 | + var sqrt=stdlib.Math.sqrt; |
| 130 | + var sz=0; |
| 131 | + function setup(n) { |
| 132 | + n=n|0; |
| 133 | + sz=n; |
| 134 | + } |
| 135 | + function fft(normalize) { |
| 136 | + normalize=normalize|0; |
| 137 | + var m=0, mh=0, i=0, j=0, k=0; |
| 138 | + var wr=0.0, wi=0.0, xr=0.0, xi=0.0, w=0.0, r=0.0; |
| 139 | + var tcnt=0; |
| 140 | + for(mh=1;(m=mh<<1)<=(sz|0);mh=m) { |
| 141 | + for(i=0;(i|0)<(mh|0);i=(i+1)|0) { |
| 142 | + wr=+FLOAT64[(sz*16+tcnt)>>3]; |
| 143 | + wi=+FLOAT64[(sz*16+tcnt+8)>>3]; |
| 144 | + tcnt=(tcnt+16)|0; |
| 145 | + for(j=i;(j|0)<(sz|0);j=(j+m)|0) { |
| 146 | + k=(j+mh)|0; |
| 147 | + xr=wr*FLOAT64[k<<3>>3]-wi*FLOAT64[(sz+k)<<3>>3]; |
| 148 | + xi=wr*FLOAT64[(sz+k)<<3>>3]+wi*FLOAT64[k<<3>>3]; |
| 149 | + FLOAT64[k<<3>>3] = +FLOAT64[j<<3>>3]-xr; |
| 150 | + FLOAT64[(sz+k)<<3>>3] = +FLOAT64[(sz+j)<<3>>3]-xi; |
| 151 | + FLOAT64[j<<3>>3] = +FLOAT64[j<<3>>3] + xr; |
| 152 | + FLOAT64[(sz+j)<<3>>3] = +FLOAT64[(sz+j)<<3>>3] + xi; |
| 153 | + } |
| 154 | + } |
| 155 | + } |
| 156 | + if(normalize) { |
| 157 | + r = +(1.0/+(sz|0)); |
| 158 | + for(i=0;(i|0)<(sz|0);i=(i+1)|0) { |
| 159 | + FLOAT64[i<<3>>3]=+FLOAT64[i<<3>>3]*r; |
| 160 | + FLOAT64[(sz+i)<<3>>3]=+FLOAT64[(sz+i)<<3>>3]*r; |
| 161 | + } |
| 162 | + } |
| 163 | + } |
| 164 | + function mag() { |
| 165 | + var i=0,j=0; |
| 166 | + for(i=0;(i|0)<(sz|0);i=(i+1)|0) { |
| 167 | + j=(sz*2-i-1)|0; |
| 168 | + FLOAT64[i<<3>>3]=sqrt(FLOAT64[i<<3>>3]*FLOAT64[i<<3>>3]+FLOAT64[j<<3>>3]*FLOAT64[j<<3>>3]); |
| 169 | + } |
| 170 | + } |
| 171 | + return { |
| 172 | + setup:setup, |
| 173 | + fft:fft, |
| 174 | + mag:mag |
| 175 | + }; |
| 176 | +} |
0 commit comments