Mercurial > hg > js-dsp-test

annotate fft/fft.js/test.html @ 40:223f770b5341 kissfft-double tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Try a double-precision kissfft
author Chris Cannam
date Wed, 07 Sep 2016 10:40:32 +0100
parents 66f9fd5ac611
children
rev   line source
Chris@25 1 <!DOCTYPE html>
Chris@25 2
Chris@25 3 <head>
Chris@25 4 <title>fft.js</title>
Chris@25 5 <script src='lib/real.js'></script>
Chris@25 6 <script src='lib/complex.js'></script>
Chris@25 7 <script>
Chris@25 8
Chris@25 9 function kahanDiff(x, xOffset, xStride, y, yOffset, yStride, n, run) {
Chris@25 10 var sum = 0.0, compensation = 0.0
Chris@25 11
Chris@25 12 for (var i = 0; i < x.length / run / xStride; i++) {
Chris@25 13 for (var j = 0; j < run; j++) {
Chris@25 14 var v = Math.abs(x[xOffset + run * xStride * i + j] - y[yOffset + run * yStride * i + j] / n) - compensation
Chris@25 15
Chris@25 16 var t = sum + v
Chris@25 17
Chris@25 18 compensation = (t - sum) - v
Chris@25 19
Chris@25 20 sum = t
Chris@25 21
Chris@25 22 if (isNaN(sum)) {
Chris@25 23 debugger
Chris@25 24 }
Chris@25 25 }
Chris@25 26 }
Chris@25 27
Chris@25 28 return sum
Chris@25 29 }
Chris@25 30
Chris@25 31 function kahanStrideDiff(x, xOffset, xStride, run) {
Chris@25 32 var sum = 0.0, compensation = 0.0
Chris@25 33
Chris@25 34 for (var i = 0; i < x.length / run / xStride; i++) {
Chris@25 35 for (var j = run; j < run * xStride; j++) {
Chris@25 36 var v = Math.abs(x[xOffset + run * xStride * i + j]) - compensation
Chris@25 37
Chris@25 38 var t = sum + v
Chris@25 39
Chris@25 40 compensation = (t - sum) - v
Chris@25 41
Chris@25 42 sum = t
Chris@25 43
Chris@25 44 if (isNaN(sum)) {
Chris@25 45 debugger
Chris@25 46 }
Chris@25 47 }
Chris@25 48 }
Chris@25 49
Chris@25 50 return sum
Chris@25 51 }
Chris@25 52
Chris@25 53 function kahanOffsetDiff(x, xOffset) {
Chris@25 54 var sum = 0.0, compensation = 0.0
Chris@25 55
Chris@25 56 for (var i = 0; i < xOffset; i++) {
Chris@25 57 var v = Math.abs(x[i]) - compensation
Chris@25 58
Chris@25 59 var t = sum + v
Chris@25 60
Chris@25 61 compensation = (t - sum) - v
Chris@25 62
Chris@25 63 sum = t
Chris@25 64
Chris@25 65 if (isNaN(sum)) {
Chris@25 66 debugger
Chris@25 67 }
Chris@25 68 }
Chris@25 69
Chris@25 70 return sum
Chris@25 71 }
Chris@25 72
Chris@25 73 function testComplex(n) {
Chris@25 74 var i = new Float64Array(2 * n), o1 = new Float64Array(2 * n), o2 = new Float64Array(2 * n)
Chris@25 75 var fft = new FFT.complex(n, false), ifft = new FFT.complex(n, true)
Chris@25 76
Chris@25 77 for (var j = 0; j < (2 * n); j++) {
Chris@25 78 i[j] = Math.random()
Chris@25 79 }
Chris@25 80
Chris@25 81 fft.process(o1, 0, 1, i, 0, 1)
Chris@25 82 ifft.process(o2, 0, 1, o1, 0, 1)
Chris@25 83
Chris@25 84 return kahanDiff(i, 0, 1, o2, 0, 1, n, 2)
Chris@25 85 }
Chris@25 86
Chris@25 87 function testComplexStride(n) {
Chris@25 88 var iStride = ~~(12 * Math.random()) + 1, o1Stride = ~~(16 * Math.random()) + 1, o2Stride = ~~(5 * Math.random()) + 1
Chris@25 89
Chris@25 90 var i = new Float64Array(2 * iStride * n), o1 = new Float64Array(2 * o1Stride * n), o2 = new Float64Array(2 * o2Stride * n)
Chris@25 91 var fft = new FFT.complex(n, false), ifft = new FFT.complex(n, true)
Chris@25 92
Chris@25 93 for (var j = 0; j < n; j++) {
Chris@25 94 i[2 * j * iStride + 0] = Math.random()
Chris@25 95 i[2 * j * iStride + 1] = Math.random()
Chris@25 96 }
Chris@25 97
Chris@25 98 fft.process(o1, 0, o1Stride, i, 0, iStride)
Chris@25 99 ifft.process(o2, 0, o2Stride, o1, 0, o1Stride)
Chris@25 100
Chris@25 101 return [kahanDiff(i, 0, iStride, o2, 0, o2Stride, n, 2), kahanStrideDiff(o2, 0, o2Stride, 2)]
Chris@25 102 }
Chris@25 103
Chris@25 104 function testRealToComplex(n) {
Chris@25 105 var i = new Float64Array(n), o1 = new Float64Array(2 * n), o2 = new Float64Array(2 * n)
Chris@25 106 var fft = new FFT.complex(n, false), ifft = new FFT.complex(n, true)
Chris@25 107
Chris@25 108 for (var j = 0; j < n; j++) {
Chris@25 109 i[j] = Math.random()
Chris@25 110 }
Chris@25 111
Chris@25 112 fft.simple(o1, i, 'real')
Chris@25 113 ifft.simple(o2, o1)
Chris@25 114
Chris@25 115 return [kahanDiff(i, 0, 1, o2, 0, 2, n, 1), kahanStrideDiff(o2, 2, 1)]
Chris@25 116 }
Chris@25 117
Chris@25 118 function testRealToComplexWithOffset(n) {
Chris@25 119 var iOffset = ~~(12 * Math.random()), o1Offset = ~~(12 * Math.random()), o2Offset = ~~(12 * Math.random())
Chris@25 120
Chris@25 121 var i = new Float64Array(iOffset + n), o1 = new Float64Array(o1Offset + 2 * n), o2 = new Float64Array(o2Offset + 2 * n)
Chris@25 122 var fft = new FFT.complex(n, false), ifft = new FFT.complex(n, true)
Chris@25 123
Chris@25 124 for (var j = 0; j < n; j++) {
Chris@25 125 i[iOffset + j] = Math.random()
Chris@25 126 }
Chris@25 127
Chris@25 128 fft.process(o1, o1Offset, 1, i, iOffset, 1, 'real')
Chris@25 129 ifft.process(o2, o2Offset, 1, o1, o1Offset, 1)
Chris@25 130
Chris@25 131 return [kahanDiff(i, 1, o2, 2, n, 1), kahanStrideDiff(o2, 2, 1), kahanOffsetDiff(o1, o1Offset), kahanOffsetDiff(o2, o2Offset)]
Chris@25 132 }
Chris@25 133
Chris@25 134
Chris@25 135 var count = 0
Chris@25 136 for (var i = 2; i < 100; i++) {
Chris@25 137 var v = testComplex(i)
Chris@25 138 if (isNaN(v) || v > 1e-12) {
Chris@25 139 console.log('Complex', i, v), count++
Chris@25 140 }
Chris@25 141 }
Chris@25 142
Chris@25 143 for (var i = 2; i < 100; i++) {
Chris@25 144 var v = testComplexStride(i)
Chris@25 145 if (isNaN(v[0]) || v[0] > 1e-12 || v[1] > 0) {
Chris@25 146 console.log('Complex Strided', i, v), count++
Chris@25 147 }
Chris@25 148 }
Chris@25 149
Chris@25 150 for (var i = 2; i < 100; i++) {
Chris@25 151 var v = testRealToComplex(i)
Chris@25 152 if (isNaN(v[0]) || isNaN(v[1]) || v[0] > 1e-12 || v[1] > 1e-12) {
Chris@25 153 console.log('Real to Complex', i, v), count++
Chris@25 154 }
Chris@25 155 }
Chris@25 156
Chris@25 157 for (var i = 2; i < 100; i++) {
Chris@25 158 var v = testRealToComplexWithOffset(i)
Chris@25 159 if (isNaN(v[0]) || isNaN(v[1]) || isNaN(v[2]) || isNaN(v[3]) || v[0] > 1e-12 || v[1] > 1e-12 || v[2] > 0 || v[3] > 0) {
Chris@25 160 console.log('Real to Complex Offset', i, v), count++
Chris@25 161 }
Chris@25 162 }
Chris@25 163
Chris@25 164 var f0 = new Float64Array(32), r0 = new Float64Array(64), fft = new FFT.complex(32, false); f0[0] = 1, f0[31] = 1
Chris@25 165 fft.simple(r0, f0, 'real')
Chris@25 166
Chris@25 167 console.log(r0)
Chris@25 168
Chris@25 169 console.log("Failcount", count)
Chris@25 170 </script>
Chris@25 171 </head>