From 534969d42ca5168d645678345cd21242fe41f389 Mon Sep 17 00:00:00 2001 From: Phil Burk Date: Tue, 30 Dec 2014 16:53:03 -0800 Subject: Initial commit of code. --- src/com/softsynth/math/FourierMath.java | 254 ++++++++++++++++++++++++++++++++ 1 file changed, 254 insertions(+) create mode 100644 src/com/softsynth/math/FourierMath.java (limited to 'src/com/softsynth/math/FourierMath.java') diff --git a/src/com/softsynth/math/FourierMath.java b/src/com/softsynth/math/FourierMath.java new file mode 100644 index 0000000..82d4ec9 --- /dev/null +++ b/src/com/softsynth/math/FourierMath.java @@ -0,0 +1,254 @@ +/* + * Copyright 2009 Phil Burk, Mobileer Inc + * + * 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. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +package com.softsynth.math; + +//Simple Fast Fourier Transform. +public class FourierMath { + static private final int MAX_SIZE_LOG_2 = 16; + static BitReverseTable[] reverseTables = new BitReverseTable[MAX_SIZE_LOG_2]; + static DoubleSineTable[] sineTables = new DoubleSineTable[MAX_SIZE_LOG_2]; + static FloatSineTable[] floatSineTables = new FloatSineTable[MAX_SIZE_LOG_2]; + + private static class DoubleSineTable { + double[] sineValues; + + DoubleSineTable(int numBits) { + int len = 1 << numBits; + sineValues = new double[1 << numBits]; + for (int i = 0; i < len; i++) { + sineValues[i] = Math.sin((i * Math.PI * 2.0) / len); + } + } + } + + private static double[] getDoubleSineTable(int n) { + DoubleSineTable sineTable = sineTables[n]; + if (sineTable == null) { + sineTable = new DoubleSineTable(n); + sineTables[n] = sineTable; + } + return sineTable.sineValues; + } + + private static class FloatSineTable { + float[] sineValues; + + FloatSineTable(int numBits) { + int len = 1 << numBits; + sineValues = new float[1 << numBits]; + for (int i = 0; i < len; i++) { + sineValues[i] = (float) Math.sin((i * Math.PI * 2.0) / len); + } + } + } + + private static float[] getFloatSineTable(int n) { + FloatSineTable sineTable = floatSineTables[n]; + if (sineTable == null) { + sineTable = new FloatSineTable(n); + floatSineTables[n] = sineTable; + } + return sineTable.sineValues; + } + + private static class BitReverseTable { + int[] reversedBits; + + BitReverseTable(int numBits) { + reversedBits = new int[1 << numBits]; + for (int i = 0; i < reversedBits.length; i++) { + reversedBits[i] = reverseBits(i, numBits); + } + } + + static int reverseBits(int index, int numBits) { + int i, rev; + + for (i = rev = 0; i < numBits; i++) { + rev = (rev << 1) | (index & 1); + index >>= 1; + } + + return rev; + } + } + + private static int[] getReverseTable(int n) { + BitReverseTable reverseTable = reverseTables[n]; + if (reverseTable == null) { + reverseTable = new BitReverseTable(n); + reverseTables[n] = reverseTable; + } + return reverseTable.reversedBits; + } + + /** + * Calculate the amplitude of the sine wave associated with each bin of a complex FFT result. + * + * @param ar + * @param ai + * @param magnitudes + */ + public static void calculateMagnitudes(double ar[], double ai[], double[] magnitudes) { + for (int i = 0; i < magnitudes.length; ++i) { + magnitudes[i] = Math.sqrt((ar[i] * ar[i]) + (ai[i] * ai[i])); + } + } + + /** + * Calculate the amplitude of the sine wave associated with each bin of a complex FFT result. + * + * @param ar + * @param ai + * @param magnitudes + */ + public static void calculateMagnitudes(float ar[], float ai[], float[] magnitudes) { + for (int i = 0; i < magnitudes.length; ++i) { + magnitudes[i] = (float) Math.sqrt((ar[i] * ar[i]) + (ai[i] * ai[i])); + } + } + + public static void transform(int sign, int n, double ar[], double ai[]) { + double scale = Math.sqrt(1.0 / n); + + int numBits = FourierMath.numBits(n); + int[] reverseTable = getReverseTable(numBits); + double[] sineTable = getDoubleSineTable(numBits); + int mask = n - 1; + int cosineOffset = n / 4; // phase offset between cos and sin + + int i, j; + for (i = 0; i < n; i++) { + j = reverseTable[i]; + if (j >= i) { + double tempr = ar[j] * scale; + double tempi = ai[j] * scale; + ar[j] = ar[i] * scale; + ai[j] = ai[i] * scale; + ar[i] = tempr; + ai[i] = tempi; + } + } + + int mmax, stride; + int numerator = sign * n; + for (mmax = 1, stride = 2 * mmax; mmax < n; mmax = stride, stride = 2 * mmax) { + int phase = 0; + int phaseIncrement = numerator / (2 * mmax); + for (int m = 0; m < mmax; ++m) { + double wr = sineTable[(phase + cosineOffset) & mask]; + double wi = sineTable[phase]; + + for (i = m; i < n; i += stride) { + j = i + mmax; + double tr = (wr * ar[j]) - (wi * ai[j]); + double ti = (wr * ai[j]) + (wi * ar[j]); + ar[j] = ar[i] - tr; + ai[j] = ai[i] - ti; + ar[i] += tr; + ai[i] += ti; + } + + phase = (phase + phaseIncrement) & mask; + } + mmax = stride; + } + } + + public static void transform(int sign, int n, float ar[], float ai[]) { + float scale = (float) Math.sqrt(1.0 / n); + + int numBits = FourierMath.numBits(n); + int[] reverseTable = getReverseTable(numBits); + float[] sineTable = getFloatSineTable(numBits); + int mask = n - 1; + int cosineOffset = n / 4; // phase offset between cos and sin + + int i, j; + for (i = 0; i < n; i++) { + j = reverseTable[i]; + if (j >= i) { + float tempr = ar[j] * scale; + float tempi = ai[j] * scale; + ar[j] = ar[i] * scale; + ai[j] = ai[i] * scale; + ar[i] = tempr; + ai[i] = tempi; + } + } + + int mmax, stride; + int numerator = sign * n; + for (mmax = 1, stride = 2 * mmax; mmax < n; mmax = stride, stride = 2 * mmax) { + int phase = 0; + int phaseIncrement = numerator / (2 * mmax); + for (int m = 0; m < mmax; ++m) { + float wr = sineTable[(phase + cosineOffset) & mask]; + float wi = sineTable[phase]; + + for (i = m; i < n; i += stride) { + j = i + mmax; + float tr = (wr * ar[j]) - (wi * ai[j]); + float ti = (wr * ai[j]) + (wi * ar[j]); + ar[j] = ar[i] - tr; + ai[j] = ai[i] - ti; + ar[i] += tr; + ai[i] += ti; + } + + phase = (phase + phaseIncrement) & mask; + } + mmax = stride; + } + } + + /** + * Calculate log2(n) + * + * @param powerOf2 must be a power of two, for example 512 or 1024 + * @return for example, 9 for an input value of 512 + */ + public static int numBits(int powerOf2) { + int i; + assert ((powerOf2 & (powerOf2 - 1)) == 0); // is it a power of 2? + for (i = -1; powerOf2 > 0; powerOf2 = powerOf2 >> 1, i++) + ; + return i; + } + + /** + * Calculate an FFT in place, modifying the input arrays. + * + * @param n + * @param ar + * @param ai + */ + public static void fft(int n, double ar[], double ai[]) { + transform(1, n, ar, ai); // TODO -1 or 1 + } + + /** + * Calculate an inverse FFT in place, modifying the input arrays. + * + * @param n + * @param ar + * @param ai + */ + public static void ifft(int n, double ar[], double ai[]) { + transform(-1, n, ar, ai); // TODO -1 or 1 + } +} -- cgit v1.2.3