fix scaling in FFT

3 files changed, 133 insertions, 24 deletions

diff --git a/src/com/jsyn/JSyn.java b/src/com/jsyn/JSyn.java
index c6e35b2..bc6f813 100644
--- a/src/com/jsyn/JSyn.java
+++ b/src/com/jsyn/JSyn.java
@@ -54,8 +54,8 @@ public class JSyn {
     // Update these for every release.
     private final static int VERSION_MAJOR = 16;
     private final static int VERSION_MINOR = 7;
-    private final static int VERSION_REVISION = 4;
-    public final static int BUILD_NUMBER = 457;
+    private final static int VERSION_REVISION = 5;
+    public final static int BUILD_NUMBER = 458;
     private final static long BUILD_TIME = new GregorianCalendar(2014, GregorianCalendar.DECEMBER,
             25).getTime().getTime();
 
diff --git a/src/com/softsynth/math/FourierMath.java b/src/com/softsynth/math/FourierMath.java
index 82d4ec9..d133d7f 100644
--- a/src/com/softsynth/math/FourierMath.java
+++ b/src/com/softsynth/math/FourierMath.java
@@ -123,7 +123,7 @@ public class FourierMath {
     }
 
     public static void transform(int sign, int n, double ar[], double ai[]) {
-        double scale = Math.sqrt(1.0 / n);
+        double scale = (sign > 0) ? (2.0 / n) : (0.5);
 
         int numBits = FourierMath.numBits(n);
         int[] reverseTable = getReverseTable(numBits);
@@ -150,7 +150,7 @@ public class FourierMath {
             int phase = 0;
             int phaseIncrement = numerator / (2 * mmax);
             for (int m = 0; m < mmax; ++m) {
-                double wr = sineTable[(phase + cosineOffset) & mask];
+                double wr = sineTable[(phase + cosineOffset) & mask]; // cosine
                 double wi = sineTable[phase];
 
                 for (i = m; i < n; i += stride) {
@@ -170,7 +170,7 @@ public class FourierMath {
     }
 
     public static void transform(int sign, int n, float ar[], float ai[]) {
-        float scale = (float) Math.sqrt(1.0 / n);
+        float scale = (sign > 0) ? (2.0f / n) : (0.5f);
 
         int numBits = FourierMath.numBits(n);
         int[] reverseTable = getReverseTable(numBits);
@@ -197,7 +197,7 @@ public class FourierMath {
             int phase = 0;
             int phaseIncrement = numerator / (2 * mmax);
             for (int m = 0; m < mmax; ++m) {
-                float wr = sineTable[(phase + cosineOffset) & mask];
+                float wr = sineTable[(phase + cosineOffset) & mask]; // cosine
                 float wi = sineTable[phase];
 
                 for (i = m; i < n; i += stride) {
diff --git a/tests/com/jsyn/util/TestFFT.java b/tests/com/jsyn/util/TestFFT.java
index f7fcce6..c79157d 100644
--- a/tests/com/jsyn/util/TestFFT.java
+++ b/tests/com/jsyn/util/TestFFT.java
@@ -4,9 +4,9 @@
  * Licensed under the Apache License, Version 2.0 (the "License");
  * you may not use this file except in compliance with the License.
  * You may obtain a copy of the License at
- * 
+ *
  *     http://www.apache.org/licenses/LICENSE-2.0
- * 
+ *
  * Unless required by applicable law or agreed to in writing, software
  * distributed under the License is distributed on an "AS IS" BASIS,
  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
@@ -32,7 +32,7 @@ public class TestFFT extends TestCase {
         super.tearDown();
     }
 
-    public void checkSingleSine(int size, int bin) {
+    public void checkSingleSineDouble(int size, int bin) {
         double[] ar = new double[size];
         double[] ai = new double[size];
         double[] magnitudes = new double[size];
@@ -43,9 +43,62 @@ public class TestFFT extends TestCase {
         FourierMath.transform(1, size, ar, ai);
         FourierMath.calculateMagnitudes(ar, ai, magnitudes);
 
-        assertTrue(magnitudes[bin - 1] < 0.001);
-        assertTrue(magnitudes[bin] > 0.5);
-        assertTrue(magnitudes[bin + 1] < 0.001);
+        assertEquals("magnitude", 0.0, magnitudes[bin-1], 0.000001);
+        assertEquals("magnitude", amplitude, magnitudes[bin], 0.000001);
+        assertEquals("magnitude", 0.0, magnitudes[bin+1], 0.000001);
+ /*
+        for (int i = 0; i < magnitudes.length; i++) {
+            System.out.printf("%d = %9.7f\n", i, magnitudes[i]);
+        }
+*/
+
+    }
+
+    public void checkSingleSineFloat(int size, int bin) {
+        float[] ar = new float[size];
+        float[] ai = new float[size];
+        float[] magnitudes = new float[size];
+
+        double amplitude = 1.0;
+        addSineWave(size, bin, ar, amplitude);
+
+        FourierMath.transform(1, size, ar, ai);
+        FourierMath.calculateMagnitudes(ar, ai, magnitudes);
+
+        assertEquals("magnitude", 0.0f, magnitudes[bin-1], 0.000001);
+        assertEquals("magnitude", amplitude, magnitudes[bin], 0.000001);
+        assertEquals("magnitude", 0.0f, magnitudes[bin+1], 0.000001);
+/*
+        for (int i = 0; i < magnitudes.length; i++) {
+            System.out.printf("%d = %9.7f\n", i, magnitudes[i]);
+        }
+*/
+    }
+
+    public void checkMultipleSine(int size, int[] bins, double[] amplitudes) {
+        double[] ar = new double[size];
+        double[] ai = new double[size];
+        double[] magnitudes = new double[size];
+
+        for(int i = 0; i<bins.length; i++) {
+            addSineWave(size, bins[i], ar, amplitudes[i]);
+        }
+
+        FourierMath.transform(1, size, ar, ai);
+        FourierMath.calculateMagnitudes(ar, ai, magnitudes);
+
+        for(int bin = 0; bin<size; bin++) {
+            System.out.printf("%d = %9.7f\n", bin, magnitudes[bin]);
+
+            double amplitude = 0.0;
+            for(int i = 0; i<bins.length; i++) {
+                if ((bin == bins[i]) || (bin == (size - bins[i]))) {
+                    amplitude = amplitudes[i];
+                    break;
+                }
+            }
+            assertEquals("magnitude", amplitude, magnitudes[bin], 0.000001);
+        }
 
     }
 
@@ -54,19 +107,45 @@ public class TestFFT extends TestCase {
         double phaseIncrement = 2.0 * Math.PI * bin / size;
         for (int i = 0; i < size; i++) {
             ar[i] += Math.sin(phase) * amplitude;
-            // System.out.println( i + " = " + ar[i] );
             phase += phaseIncrement;
         }
     }
+    private void addSineWave(int size, int bin, float[] ar, double amplitude) {
+        double phase = 0.0;
+        double phaseIncrement = 2.0 * Math.PI * bin / size;
+        for (int i = 0; i < size; i++) {
+            ar[i] += (float) (Math.sin(phase) * amplitude);
+            phase += phaseIncrement;
+        }
+    }
+
+    public void testSinglesDouble() {
+        checkSingleSineDouble(32, 1);
+        checkSingleSineDouble(32, 4);
+        checkSingleSineDouble(64, 5);
+        checkSingleSineDouble(256, 3);
+    }
+
+    public void testSinglesFloat() {
+        checkSingleSineFloat(32, 1);
+        checkSingleSineFloat(32, 4);
+        checkSingleSineFloat(64, 5);
+        checkSingleSineFloat(256, 3);
+    }
 
-    public void testSingles() {
-        checkSingleSine(32, 1);
-        checkSingleSine(32, 2);
-        checkSingleSine(64, 5);
-        checkSingleSine(256, 3);
+    public void testMultipleSines32() {
+        int[] bins = { 1, 5 };
+        double[] amplitudes = { 1.0, 2.0 };
+        checkMultipleSine(32, bins, amplitudes);
     }
 
-    public void checkInverseFFT(int size, int bin) {
+    public void testMultipleSines64() {
+        int[] bins = { 2, 4, 7 };
+        double[] amplitudes = { 1.0, 0.3, 0.5 };
+        checkMultipleSine(64, bins, amplitudes);
+    }
+
+    public void checkInverseFftDouble(int size, int bin) {
         double[] ar1 = new double[size];
         double[] ai1 = new double[size];
         double[] ar2 = new double[size];
@@ -89,10 +168,40 @@ public class TestFFT extends TestCase {
         }
     }
 
-    public void testInverse() {
-        checkInverseFFT(32, 1);
-        checkInverseFFT(32, 2);
-        checkInverseFFT(128, 17);
-        checkInverseFFT(512, 23);
+    public void checkInverseFftFloat(int size, int bin) {
+        float[] ar1 = new float[size];
+        float[] ai1 = new float[size];
+        float[] ar2 = new float[size];
+        float[] ai2 = new float[size];
+
+        double amplitude = 1.0;
+        addSineWave(size, bin, ar1, amplitude);
+
+        // Save a copy of the source.
+        System.arraycopy(ar1, 0, ar2, 0, size);
+        System.arraycopy(ai1, 0, ai2, 0, size);
+
+        FourierMath.transform(1, size, ar1, ai1); // FFT
+
+        FourierMath.transform(-1, size, ar1, ai1); // IFFT
+
+        for (int i = 0; i < size; i++) {
+            assertEquals(ar2[i], ar1[i], 0.00001);
+            assertEquals(ai2[i], ai1[i], 0.00001);
+        }
+    }
+
+    public void testInverseDouble() {
+        checkInverseFftDouble(32, 1);
+        checkInverseFftDouble(32, 2);
+        checkInverseFftDouble(128, 17);
+        checkInverseFftDouble(512, 23);
+    }
+
+    public void testInverseFloat() {
+        checkInverseFftFloat(32, 1);
+        checkInverseFftFloat(32, 2);
+        checkInverseFftFloat(128, 17);
+        checkInverseFftFloat(512, 23);
     }
 }