Chris@2: /*
Chris@2: 	Copyright (C) 2001, 2006 by Simon Dixon
Chris@2: 
Chris@2: 	This program is free software; you can redistribute it and/or modify
Chris@2: 	it under the terms of the GNU General Public License as published by
Chris@2: 	the Free Software Foundation; either version 2 of the License, or
Chris@2: 	(at your option) any later version.
Chris@2: 
Chris@2: 	This program is distributed in the hope that it will be useful,
Chris@2: 	but WITHOUT ANY WARRANTY; without even the implied warranty of
Chris@2: 	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
Chris@2: 	GNU General Public License for more details.
Chris@2: 
Chris@2: 	You should have received a copy of the GNU General Public License along
Chris@2: 	with this program (the file gpl.txt); if not, download it from
Chris@2: 	http://www.gnu.org/licenses/gpl.txt or write to the
Chris@2: 	Free Software Foundation, Inc.,
Chris@2: 	51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
Chris@2: */
Chris@2: 
Chris@2: package at.ofai.music.audio;
Chris@2: 
Chris@2: /** Class for computing a windowed fast Fourier transform.
Chris@2:  *  Implements some of the window functions for the STFT from
Chris@2:  *  Harris (1978), Proc. IEEE, 66, 1, 51-83.
Chris@2:  */
Chris@2: public class FFT {
Chris@2: 
Chris@2: 	/** used in {@link FFT#fft(double[], double[], int)} to specify
Chris@2: 	 *  a forward Fourier transform */
Chris@2: 	public static final int FORWARD = -1;
Chris@2: 	/** used in {@link FFT#fft(double[], double[], int)} to specify
Chris@2: 	 *  an inverse Fourier transform */
Chris@2: 	public static final int REVERSE = 1;
Chris@2: 	/** used in {@link FFT#makeWindow(int,int,int)} to specify a
Chris@2: 	 *  rectangular window function */
Chris@2: 	public static final int RECT = 0;
Chris@2: 	/** used in {@link FFT#makeWindow(int,int,int)} to specify a
Chris@2: 	 *  Hamming window function */
Chris@2: 	public static final int HAMMING = 1;
Chris@2: 	/** used in {@link FFT#makeWindow(int,int,int)} to specify a
Chris@2: 	 *  61-dB 3-sample Blackman-Harris window function */
Chris@2: 	public static final int BH3 = 2;
Chris@2: 	/** used in {@link FFT#makeWindow(int,int,int)} to specify a
Chris@2: 	 *  74-dB 4-sample Blackman-Harris window function */
Chris@2: 	public static final int BH4 = 3;
Chris@2: 	/** used in {@link FFT#makeWindow(int,int,int)} to specify a
Chris@2: 	 *  minimum 3-sample Blackman-Harris window function */
Chris@2: 	public static final int BH3MIN = 4;
Chris@2: 	/** used in {@link FFT#makeWindow(int,int,int)} to specify a
Chris@2: 	 *  minimum 4-sample Blackman-Harris window function */
Chris@2: 	public static final int BH4MIN = 5;
Chris@2: 	/** used in {@link FFT#makeWindow(int,int,int)} to specify a
Chris@2: 	 *  Gaussian window function */
Chris@2: 	public static final int GAUSS = 6; 
Chris@2: 	static final double twoPI = 2 * Math.PI;
Chris@2: 
Chris@2: 	/** The FFT method. Calculation is inline, for complex data stored
Chris@2: 	 *  in 2 separate arrays. Length of input data must be a power of two.
Chris@2: 	 *  @param re        the real part of the complex input and output data
Chris@2: 	 *  @param im        the imaginary part of the complex input and output data
Chris@2: 	 *  @param direction the direction of the Fourier transform (FORWARD or
Chris@2: 	 *  REVERSE)
Chris@2: 	 *  @throws IllegalArgumentException if the length of the input data is
Chris@2: 	 *  not a power of 2
Chris@2: 	 */
Chris@2: 	public static void fft(double re[], double im[], int direction) {
Chris@2: 		int n = re.length;
Chris@2: 		int bits = (int)Math.rint(Math.log(n) / Math.log(2));
Chris@2: 		if (n != (1 << bits))
Chris@2: 			throw new IllegalArgumentException("FFT data must be power of 2");
Chris@2: 		int localN;
Chris@2: 		int j = 0;
Chris@2: 		for (int i = 0; i < n-1; i++) {
Chris@2: 			if (i < j) {
Chris@2: 				double temp = re[j];
Chris@2: 				re[j] = re[i];
Chris@2: 				re[i] = temp;
Chris@2: 				temp = im[j];
Chris@2: 				im[j] = im[i];
Chris@2: 				im[i] = temp;
Chris@2: 			}
Chris@2: 			int k = n / 2;
Chris@2: 			while ((k >= 1) &&  (k - 1 < j)) {
Chris@2: 				j = j - k;
Chris@2: 				k = k / 2;
Chris@2: 			}
Chris@2: 			j = j + k;
Chris@2: 		}
Chris@2: 		for(int m = 1; m <= bits; m++) {
Chris@2: 			localN = 1 << m;
Chris@2: 			double Wjk_r = 1;
Chris@2: 			double Wjk_i = 0;
Chris@2: 			double theta = twoPI / localN;
Chris@2: 			double Wj_r = Math.cos(theta);
Chris@2: 			double Wj_i = direction * Math.sin(theta);
Chris@2: 			int nby2 = localN / 2;
Chris@2: 			for (j = 0; j < nby2; j++) {
Chris@2: 				for (int k = j; k < n; k += localN) {
Chris@2: 					int id = k + nby2;
Chris@2: 					double tempr = Wjk_r * re[id] - Wjk_i * im[id];
Chris@2: 					double tempi = Wjk_r * im[id] + Wjk_i * re[id];
Chris@2: 					re[id] = re[k] - tempr;
Chris@2: 					im[id] = im[k] - tempi;
Chris@2: 					re[k] += tempr;
Chris@2: 					im[k] += tempi;
Chris@2: 				}
Chris@2: 				double wtemp = Wjk_r;
Chris@2: 				Wjk_r = Wj_r * Wjk_r  - Wj_i * Wjk_i;
Chris@2: 				Wjk_i = Wj_r * Wjk_i  + Wj_i * wtemp;
Chris@2: 			}
Chris@2: 		}
Chris@2: 	} // fft()
Chris@2: 
Chris@2: 	/** Computes the power spectrum of a real sequence (in place).
Chris@2: 	 *  @param re the real input and output data; length must be a power of 2
Chris@2: 	 */
Chris@2: 	public static void powerFFT(double[] re) {
Chris@2: 		double[] im = new double[re.length];
Chris@2: 		fft(re, im, FORWARD);
Chris@2: 		for (int i = 0; i < re.length; i++)
Chris@2: 			re[i] = re[i] * re[i] + im[i] * im[i];
Chris@2: 	} // powerFFT()	
Chris@2: 
Chris@2: 	/** Converts a real power sequence from to magnitude representation,
Chris@2: 	 *  by computing the square root of each value.
Chris@2: 	 *  @param re the real input (power) and output (magnitude) data; length
Chris@2: 	 *  must be a power of 2
Chris@2: 	 */
Chris@2: 	public static void toMagnitude(double[] re) {
Chris@2: 		for (int i = 0; i < re.length; i++)
Chris@2: 			re[i] = Math.sqrt(re[i]);
Chris@2: 	} // toMagnitude()
Chris@2: 	
Chris@2: 	/** Computes the magnitude spectrum of a real sequence (in place).
Chris@2: 	 *  @param re the real input and output data; length must be a power of 2
Chris@2: 	 */
Chris@2: 	public static void magnitudeFFT(double[] re) {
Chris@2: 		powerFFT(re);
Chris@2: 		toMagnitude(re);
Chris@2: 	} // magnitudeFFT()
Chris@2: 
Chris@2: 	/** Computes a complex (or real if im[] == {0,...}) FFT and converts
Chris@2: 	 *  the results to polar coordinates (power and phase). Both arrays
Chris@2: 	 *  must be the same length, which is a power of 2.
Chris@2: 	 *  @param re the real part of the input data and the power of the output
Chris@2: 	 *  data
Chris@2: 	 *  @param im the imaginary part of the input data and the phase of the
Chris@2: 	 *  output data
Chris@2: 	 */
Chris@2: 	public static void powerPhaseFFT(double[] re, double[] im) {
Chris@2: 		fft(re, im, FORWARD);
Chris@2: 		for (int i = 0; i < re.length; i++) {
Chris@2: 			double pow = re[i] * re[i] + im[i] * im[i];
Chris@2: 			im[i] = Math.atan2(im[i], re[i]);
Chris@2: 			re[i] = pow;
Chris@2: 		}
Chris@2: 	} // powerPhaseFFT()
Chris@2: 	
Chris@2: 	/** Inline computation of the inverse FFT given spectral input data
Chris@2: 	 *  in polar coordinates (power and phase).
Chris@2: 	 *  Both arrays must be the same length, which is a power of 2.
Chris@2: 	 *  @param pow the power of the spectral input data (and real part of the
Chris@2: 	 *  output data)
Chris@2: 	 *  @param ph the phase of the spectral input data (and the imaginary part
Chris@2: 	 *  of the output data)
Chris@2: 	 */
Chris@2: 	public static void powerPhaseIFFT(double[] pow, double[] ph) {
Chris@2: 		toMagnitude(pow);
Chris@2: 		for (int i = 0; i < pow.length; i++) {
Chris@2: 			double re = pow[i] * Math.cos(ph[i]);
Chris@2: 			ph[i] = pow[i] * Math.sin(ph[i]);
Chris@2: 			pow[i] = re;
Chris@2: 		}
Chris@2: 		fft(pow, ph, REVERSE);
Chris@2: 	} // powerPhaseIFFT()
Chris@2: 	
Chris@2: 	/** Computes a complex (or real if im[] == {0,...}) FFT and converts
Chris@2: 	 *  the results to polar coordinates (magnitude and phase). Both arrays
Chris@2: 	 *  must be the same length, which is a power of 2.
Chris@2: 	 *  @param re the real part of the input data and the magnitude of the
Chris@2: 	 *  output data
Chris@2: 	 *  @param im the imaginary part of the input data and the phase of the
Chris@2: 	 *  output data
Chris@2: 	 */
Chris@2: 	public static void magnitudePhaseFFT(double[] re, double[] im) {
Chris@2: 		powerPhaseFFT(re, im);
Chris@2: 		toMagnitude(re);
Chris@2: 	} // magnitudePhaseFFT()
Chris@2: 
Chris@2: 
Chris@2: 	/** Fill an array with the values of a standard Hamming window function
Chris@2: 	 *  @param data the array to be filled
Chris@2: 	 *  @param size the number of non zero values; if the array is larger than
Chris@2: 	 *  this, it is zero-padded symmetrically at both ends 
Chris@2: 	 */
Chris@2: 	static void hamming(double[] data, int size) {
Chris@2: 		int start = (data.length - size) / 2;
Chris@2: 		int stop = (data.length + size) / 2;
Chris@2: 		double scale = 1.0 / (double)size / 0.54;
Chris@2: 		double factor = twoPI / (double)size;
Chris@2: 		for (int i = 0; start < stop; start++, i++)
Chris@2: 			data[i] = scale * (25.0/46.0 - 21.0/46.0 * Math.cos(factor * i));
Chris@2: 	} // hamming()
Chris@2: 
Chris@2: 	/** Fill an array with the values of a minimum 4-sample Blackman-Harris
Chris@2: 	 *  window function
Chris@2: 	 *  @param data the array to be filled
Chris@2: 	 *  @param size the number of non zero values; if the array is larger than
Chris@2: 	 *  this, it is zero-padded symmetrically at both ends 
Chris@2: 	 */
Chris@2: 	static void blackmanHarris4sMin(double[] data, int size) {
Chris@2: 		int start = (data.length - size) / 2;
Chris@2: 		int stop = (data.length + size) / 2;
Chris@2: 		double scale = 1.0 / (double)size / 0.36;
Chris@2: 		for (int i = 0; start < stop; start++, i++)
Chris@2: 			data[i] = scale * ( 0.35875 -
Chris@2: 								0.48829 * Math.cos(twoPI * i / size) +
Chris@2: 								0.14128 * Math.cos(2 * twoPI * i / size) -
Chris@2: 								0.01168 * Math.cos(3 * twoPI * i / size));
Chris@2: 	} // blackmanHarris4sMin()
Chris@2: 
Chris@2: 	/** Fill an array with the values of a 74-dB 4-sample Blackman-Harris
Chris@2: 	 *  window function
Chris@2: 	 *  @param data the array to be filled
Chris@2: 	 *  @param size the number of non zero values; if the array is larger than
Chris@2: 	 *  this, it is zero-padded symmetrically at both ends 
Chris@2: 	 */
Chris@2: 	static void blackmanHarris4s(double[] data, int size) {
Chris@2: 		int start = (data.length - size) / 2;
Chris@2: 		int stop = (data.length + size) / 2;
Chris@2: 		double scale = 1.0 / (double)size / 0.4;
Chris@2: 		for (int i = 0; start < stop; start++, i++)
Chris@2: 			data[i] = scale * ( 0.40217 -
Chris@2: 								0.49703 * Math.cos(twoPI * i / size) +
Chris@2: 								0.09392 * Math.cos(2 * twoPI * i / size) -
Chris@2: 								0.00183 * Math.cos(3 * twoPI * i / size));
Chris@2: 	} // blackmanHarris4s()
Chris@2: 
Chris@2: 	/** Fill an array with the values of a minimum 3-sample Blackman-Harris
Chris@2: 	 *  window function
Chris@2: 	 *  @param data the array to be filled
Chris@2: 	 *  @param size the number of non zero values; if the array is larger than
Chris@2: 	 *  this, it is zero-padded symmetrically at both ends 
Chris@2: 	 */
Chris@2: 	static void blackmanHarris3sMin(double[] data, int size) {
Chris@2: 		int start = (data.length - size) / 2;
Chris@2: 		int stop = (data.length + size) / 2;
Chris@2: 		double scale = 1.0 / (double) size / 0.42;
Chris@2: 		for (int i = 0; start < stop; start++, i++)
Chris@2: 			data[i] = scale * ( 0.42323 -
Chris@2: 								0.49755 * Math.cos(twoPI * i / size) +
Chris@2: 								0.07922 * Math.cos(2 * twoPI * i / size));
Chris@2: 	} // blackmanHarris3sMin()
Chris@2: 
Chris@2: 	/** Fill an array with the values of a 61-dB 3-sample Blackman-Harris
Chris@2: 	 *  window function
Chris@2: 	 *  @param data the array to be filled
Chris@2: 	 *  @param size the number of non zero values; if the array is larger than
Chris@2: 	 *  this, it is zero-padded symmetrically at both ends 
Chris@2: 	 */
Chris@2: 	static void blackmanHarris3s(double[] data, int size) {
Chris@2: 		int start = (data.length - size) / 2;
Chris@2: 		int stop = (data.length + size) / 2;
Chris@2: 		double scale = 1.0 / (double) size / 0.45;
Chris@2: 		for (int i = 0; start < stop; start++, i++)
Chris@2: 			data[i] = scale * ( 0.44959 -
Chris@2: 								0.49364 * Math.cos(twoPI * i / size) +
Chris@2: 								0.05677 * Math.cos(2 * twoPI * i / size));
Chris@2: 	} // blackmanHarris3s()
Chris@2: 
Chris@2: 	/** Fill an array with the values of a Gaussian window function
Chris@2: 	 *  @param data the array to be filled
Chris@2: 	 *  @param size the number of non zero values; if the array is larger than
Chris@2: 	 *  this, it is zero-padded symmetrically at both ends 
Chris@2: 	 */
Chris@2: 	static void gauss(double[] data, int size) { // ?? between 61/3 and 74/4 BHW
Chris@2: 		int start = (data.length - size) / 2;
Chris@2: 		int stop = (data.length + size) / 2;
Chris@2: 		double delta = 5.0 / size;
Chris@2: 		double x = (1 - size) / 2.0 * delta;
Chris@2: 		double c = -Math.PI * Math.exp(1.0) / 10.0;
Chris@2: 		double sum = 0;
Chris@2: 		for (int i = start; i < stop; i++) {
Chris@2: 			data[i] = Math.exp(c * x * x);
Chris@2: 			x += delta;
Chris@2: 			sum += data[i];
Chris@2: 		}
Chris@2: 		for (int i = start; i < stop; i++)
Chris@2: 			data[i] /= sum;
Chris@2: 	} // gauss()
Chris@2: 
Chris@2: 	/** Fill an array with the values of a rectangular window function
Chris@2: 	 *  @param data the array to be filled
Chris@2: 	 *  @param size the number of non zero values; if the array is larger than
Chris@2: 	 *  this, it is zero-padded symmetrically at both ends 
Chris@2: 	 */
Chris@2: 	static void rectangle(double[] data, int size) {
Chris@2: 		int start = (data.length - size) / 2;
Chris@2: 		int stop = (data.length + size) / 2;
Chris@2: 		for (int i = start; i < stop; i++)
Chris@2: 			data[i] = 1.0 / (double) size;
Chris@2: 	} // rectangle()
Chris@2: 
Chris@2: 	/** Returns an array of values of a normalised smooth window function,
Chris@2: 	 *  as used for performing a short time Fourier transform (STFT).
Chris@2: 	 *  All functions are normalised by length and coherent gain.
Chris@2: 	 *  More information on characteristics of these functions can be found
Chris@2: 	 *  in F.J. Harris (1978), On the Use of Windows for Harmonic Analysis
Chris@2: 	 *  with the Discrete Fourier Transform, <em>Proceedings of the IEEE</em>,
Chris@2: 	 *  66, 1, 51-83.
Chris@2: 	 *  @param choice  the choice of window function, one of the constants
Chris@2: 	 *  defined above
Chris@2: 	 *  @param size    the size of the returned array
Chris@2: 	 *  @param support the number of non-zero values in the array
Chris@2: 	 *  @return the array containing the values of the window function
Chris@2: 	 */
Chris@2: 	public static double[] makeWindow(int choice, int size, int support) {
Chris@2: 		double[] data = new double[size];
Chris@2: 		if (support > size)
Chris@2: 			support = size;
Chris@2: 		switch (choice) {
Chris@2: 			case RECT:		rectangle(data, support);			break;
Chris@2: 			case HAMMING:	hamming(data, support);				break;
Chris@2: 			case BH3:		blackmanHarris3s(data, support);	break;
Chris@2: 			case BH4:		blackmanHarris4s(data, support);	break;
Chris@2: 			case BH3MIN:	blackmanHarris3sMin(data, support);	break;
Chris@2: 			case BH4MIN:	blackmanHarris4sMin(data, support);	break;
Chris@2: 			case GAUSS:		gauss(data, support);				break;
Chris@2: 			default:		rectangle(data, support);			break;
Chris@2: 		}
Chris@2: 		return data;
Chris@2: 	} // makeWindow()
Chris@2: 
Chris@2: 	/** Applies a window function to an array of data, storing the result in
Chris@2: 	 *  the data array.
Chris@2: 	 *  Performs a dot product of the data and window arrays. 
Chris@2: 	 *  @param data   the array of input data, also used for output
Chris@2: 	 *  @param window the values of the window function to be applied to data
Chris@2: 	 */
Chris@2: 	public static void applyWindow(double[] data, double[] window) {
Chris@2: 		for (int i = 0; i < data.length; i++)
Chris@2: 			data[i] *= window[i];
Chris@2: 	} // applyWindow()
Chris@2: 
Chris@2: 	/** Unit test of the FFT class.
Chris@2: 	 *  Performs a forward and inverse FFT on a 1MB array of random values
Chris@2: 	 *  and checks how closely the values are preserved.
Chris@2: 	 *  @param args ignored
Chris@2: 	 */
Chris@2: 	public static void main(String[] args) {
Chris@2: 		final int SZ = 1024 * 1024;
Chris@2: 		double[] r1 = new double[SZ];
Chris@2: 		double[] i1 = new double[SZ];
Chris@2: 		double[] r2 = new double[SZ];
Chris@2: 		double[] i2 = new double[SZ];
Chris@2: 		for (int j = 0; j < SZ; j++) {
Chris@2: 			r1[j] = r2[j] = Math.random();
Chris@2: 			i1[j] = i2[j] = Math.random();
Chris@2: 		}
Chris@2: 		System.out.println("start");
Chris@2: 		fft(r2, i2, FORWARD);
Chris@2: 		System.out.println("reverse");
Chris@2: 		fft(r2, i2, REVERSE);
Chris@2: 		System.out.println("result");
Chris@2: 		double err = 0;
Chris@2: 		for (int j = 0; j < SZ; j++)
Chris@2: 			err += Math.abs(r1[j] - r2[j] / SZ) + Math.abs(i1[j] - i2[j] / SZ);
Chris@2: 		System.out.printf( "Err: %12.10f   Av: %12.10f\n", err, err / SZ);
Chris@2: 	} // main()
Chris@2: 
Chris@2: } // class FFT