plectrum

ConstantQ.java (12653B)

---

 /*
 *      _______                       _____   _____ _____  
 *     |__   __|                     |  __ \ / ____|  __ \ 
 *        | | __ _ _ __ ___  ___  ___| |  | | (___ | |__) |
 *        | |/ _` | '__/ __|/ _ \/ __| |  | |\___ \|  ___/ 
 *        | | (_| | |  \__ \ (_) \__ \ |__| |____) | |     
 *        |_|\__,_|_|  |___/\___/|___/_____/|_____/|_|     
 *                                                         
 * -------------------------------------------------------------
 *
 * TarsosDSP is developed by Joren Six at IPEM, University Ghent
 *  
 * -------------------------------------------------------------
 *
 *  Info: http://0110.be/tag/TarsosDSP
 *  Github: https://github.com/JorenSix/TarsosDSP
 *  Releases: http://0110.be/releases/TarsosDSP/
 *  
 *  TarsosDSP includes modified source code by various authors,
 *  for credits and info, see README.
 * 
 */
 
 /*
 *      _______                       _____   _____ _____  
 *     |__   __|                     |  __ \ / ____|  __ \ 
 *        | | __ _ _ __ ___  ___  ___| |  | | (___ | |__) |
 *        | |/ _` | '__/ __|/ _ \/ __| |  | |\___ \|  ___/ 
 *        | | (_| | |  \__ \ (_) \__ \ |__| |____) | |     
 *        |_|\__,_|_|  |___/\___/|___/_____/|_____/|_|     
 *                                                         
 * -----------------------------------------------------------
 *
 *  TarsosDSP is developed by Joren Six at 
 *  The Royal Academy of Fine Arts & Royal Conservatory,
 *  University College Ghent,
 *  Hoogpoort 64, 9000 Ghent - Belgium
 *  
 *  http://tarsos.0110.be/tag/TarsosDSP
 *  https://github.com/JorenSix/TarsosDSP
 *  http://tarsos.0110.be/releases/TarsosDSP/
 * 
 */
 /* 
  * Copyright (c) 2006, Karl Helgason
  * 
  * 2007/1/8 modified by p.j.leonard
  * 
  * All rights reserved.
  * 
  * Redistribution and use in source and binary forms, with or without
  * modification, are permitted provided that the following conditions
  * are met:
  * 
  *    1. Redistributions of source code must retain the above copyright
  *       notice, this list of conditions and the following disclaimer.
  *    2. Redistributions in binary form must reproduce the above
  *       copyright notice, this list of conditions and the following
  *       disclaimer in the documentation and/or other materials
  *       provided with the distribution.
  *    3. The name of the author may not be used to endorse or promote
  *       products derived from this software without specific prior
  *       written permission.
  * 
  * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
  * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
  * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY
  * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
  * GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
  * IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
  * OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN
  * IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  */
 
 package be.tarsos.dsp;
 
 import be.tarsos.dsp.util.fft.FFT;
 
 /**
  * Implementation of the Constant Q Transform.<br> References:
  * <p>
  * Judith C. Brown, <a
  * href="http://www.wellesley.edu/Physics/brown/pubs/cq1stPaper.pdf">
  * Calculation of a constant Q spectral transform</a>, J. Acoust. Soc. Am.,
  * 89(1): 425-434, 1991.
  * </p>
  * <p>
  * Judith C. Brown and Miller S. Puckette, <a
  * href="http://www.wellesley.edu/Physics/brown/pubs/effalgV92P2698-P2701.pdf"
  * >An efficient algorithm for the calculation of a constant Q transform</a>, J.
  * Acoust. Soc. Am., Vol. 92, No. 5, November 1992
  * </p>
  * <p>
  * Benjamin Blankertz, <a href=
  * "http://wwwmath1.uni-muenster.de/logik/org/staff/blankertz/constQ/constQ.pdf"
  * >The Constant Q Transform</a>
  * </p>
  * 
  * 
  * @author Joren Six
  * @author Karl Helgason
  * @author P.J Leonard
  */
 public class ConstantQ implements AudioProcessor {
 
 
 	/**
 	 * The minimum frequency, in Hertz. The Constant-Q factors are calculated
 	 * starting from this frequency.
 	 */
 	private final float minimumFrequency ;
 
 	/**
 	 * The maximum frequency in Hertz.
 	 */
 	private final float maximumFreqency;
 
 	/**
 	 * The length of the underlying FFT.
 	 */
 	private int fftLength;
 
 	/**
 	 * Lists the start of each frequency bin, in Hertz. 
 	 */
 	private final float[] frequencies;
 
 	private final float[][] qKernel;
 
 	private final int[][] qKernel_indexes;
 	
 	/**
 	 * The array with constant q coefficients. If you for
 	 * example are interested in coefficients between 256 and 1024 Hz
 	 * (2^8 and 2^10 Hz) and you requested 12 bins per octave, you
 	 * will need 12 bins/octave * 2 octaves * 2 entries/bin = 48
 	 * places in the output buffer. The coefficient needs two entries
 	 * in the output buffer since they are complex numbers.
 	 */
 	private final float[] coefficients;
 	
 	/**
 	 * The output buffer with constant q magnitudes. If you for example are
 	 * interested in coefficients between 256 and 1024 Hz (2^8 and 2^10 Hz) and
 	 * you requested 12 bins per octave, you will need 12 bins/octave * 2
 	 * octaves = 24 places in the output buffer.
 	 */
 	private final float[] magnitudes;
 	
 	/**
 	 * The number of bins per octave.
 	 */
 	private int binsPerOctave;
 
 	/**
 	 * The underlying FFT object.
 	 */
 	private FFT fft;
 
 
 	public ConstantQ(float sampleRate, float minFreq, float maxFreq,float binsPerOctave) {
 		this(sampleRate,minFreq,maxFreq,binsPerOctave,0.001f,1.0f);
 	}
 
 	
 	public ConstantQ(float sampleRate, float minFreq, float maxFreq,float binsPerOctave, float threshold,float spread) {
 		this.minimumFrequency = minFreq;
 		this.maximumFreqency = maxFreq;
 		this.binsPerOctave = (int) binsPerOctave;
 		
 		// Calculate Constant Q		
 		double q = 1.0 / (Math.pow(2, 1.0 / binsPerOctave) - 1.0) / spread;
 
 		// Calculate number of output bins
 		int numberOfBins = (int) Math.ceil(binsPerOctave * Math.log(maximumFreqency / minimumFrequency) / Math.log(2));
 		
 		// Initialize the coefficients array (complex number so 2 x number of bins)
 		coefficients = new float[numberOfBins*2];
 		
 		// Initialize the magnitudes array
 		magnitudes = new float[numberOfBins];
 
 		
 		// Calculate the minimum length of the FFT to support the minimum
 		// frequency
 		float calc_fftlen = (float) Math.ceil(q * sampleRate / minimumFrequency);
 		
 		// No need to use power of 2 FFT length.
 		fftLength = (int) calc_fftlen; 
 		
 		//System.out.println(fftLength);
 		//The FFT length needs to be a power of two for performance reasons:
 		fftLength = (int) Math.pow(2, Math.ceil(Math.log(calc_fftlen) / Math.log(2)));
 		
 		// Create FFT object
 		fft = new FFT(fftLength);
 		qKernel = new float[numberOfBins][];
 		qKernel_indexes = new int[numberOfBins][];
 		frequencies = new float[numberOfBins];
 
 		// Calculate Constant Q kernels
 		 float[] temp = new float[fftLength*2];
 		 float[] ctemp = new float[fftLength*2];
 		 int[] cindexes = new int[fftLength];
 		    for (int i = 0; i < numberOfBins; i++) {
 		    	float[] sKernel = temp;
 		    	// Calculate the frequency of current bin
 		    	frequencies[i] = (float) (minimumFrequency * Math.pow(2, i/binsPerOctave ));
 		    	
 		    	// Calculate length of window
 		    	int len = (int)Math.min(Math.ceil( q * sampleRate / frequencies[i]), fftLength);
 		    	
 		    	for (int j = 0; j < len; j++) {
 		    		
 		    		double window = -.5*Math.cos(2.*Math.PI*(double)j/(double)len)+.5;; // Hanning Window
 		    		// double window = -.46*Math.cos(2.*Math.PI*(double)j/(double)len)+.54; // Hamming Window
 
 		    		window /= len;
 		    		
 		    		// Calculate kernel
 		    	    double x = 2*Math.PI * q * (double)j/(double)len;
 		    		sKernel[j*2] = (float) (window * Math.cos(x));
 		    		sKernel[j*2+1] = (float) (window * Math.sin(x));	
 				}
 		    	for (int j = len*2; j < fftLength*2; j++) {
 		    		sKernel[j] = 0;
 		    	}
 
 			// Perform FFT on kernel
 			fft.complexForwardTransform(sKernel);
 
 			// Remove all zeros from kernel to improve performance
 			float[] cKernel = ctemp;
 
 			int k = 0;
 			for (int j = 0, j2 = sKernel.length - 2; j < sKernel.length/2; j+=2,j2-=2)
 	    	{
 	    		double absval = Math.sqrt(sKernel[j]*sKernel[j] + sKernel[j+1]*sKernel[j+1]);
 	    	    absval += Math.sqrt(sKernel[j2]*sKernel[j2] + sKernel[j2+1]*sKernel[j2+1]);	    	    
 	    		if(absval > threshold)
 	    		{
 	    			cindexes[k] = j;
 	    			cKernel[2*k] = sKernel[j] + sKernel[j2];
 	    			cKernel[2*k + 1] = sKernel[j + 1] + sKernel[j2 + 1];
 	    			k++;
 	    		}	    		
 	    	}
 
 			sKernel = new float[k * 2];
 			int[] indexes = new int[k];
 
 			for (int j = 0; j < k * 2; j++)
 				sKernel[j] = cKernel[j];
 			for (int j = 0; j < k; j++)
 				indexes[j] = cindexes[j];
 
 			// Normalize fft output
 			for (int j = 0; j < sKernel.length; j++)
 				sKernel[j] /= fftLength;
 
 			// Perform complex conjugate on sKernel
 			for (int j = 1; j < sKernel.length; j += 2)
 				sKernel[j] = -sKernel[j];
 			
 			for (int j = 0; j < sKernel.length; j ++)
 				sKernel[j] = -sKernel[j];
 
 			qKernel_indexes[i] = indexes;
 			qKernel[i] = sKernel;
 		}
 	}
 
 	/**
 	 * Take an input buffer with audio and calculate the constant Q
 	 * coefficients.
 	 * 
 	 * @param inputBuffer
 	 *            The input buffer with audio.
 	 * 
 	 *            
 	 */
 	public void calculate(float[] inputBuffer) {
 		fft.forwardTransform(inputBuffer);
 		for (int i = 0; i < qKernel.length; i++) {
 			float[] kernel = qKernel[i];
 			int[] indexes = qKernel_indexes[i];
 			float t_r = 0;
 			float t_i = 0;
 			for (int j = 0, l = 0; j < kernel.length; j += 2, l++) {
 				int jj = indexes[l];
 				float b_r = inputBuffer[jj];
 				float b_i = inputBuffer[jj + 1];
 				float k_r = kernel[j];
 				float k_i = kernel[j + 1];
 				// COMPLEX: T += B * K
 				t_r += b_r * k_r - b_i * k_i;
 				t_i += b_r * k_i + b_i * k_r;
 			}
 			coefficients[i * 2] = t_r;
 			coefficients[i * 2 + 1] = t_i;
 		}
 	}
 	
 	/**
 	 * Take an input buffer with audio and calculate the constant Q magnitudes.
 	 * @param inputBuffer The input buffer with audio.
 	 */
 	public void calculateMagintudes(float[] inputBuffer) {
 		calculate(inputBuffer);
 		for(int i = 0 ; i < magnitudes.length ; i++){
 			magnitudes[i] = (float) Math.sqrt(coefficients[i*2] * coefficients[i*2] + coefficients[i*2+1] * coefficients[i*2+1]); 
 		}
 	}
 	
 
 	@Override
 	public boolean process(AudioEvent audioEvent) {
 		float[] audioBuffer = audioEvent.getFloatBuffer().clone();
 		if(audioBuffer.length != getFFTlength()){
 			throw new IllegalArgumentException(String.format("The length of the fft (%d) should be the same as the length of the audio buffer (%d)",getFFTlength(),audioBuffer.length));
 		}
 		calculateMagintudes(audioBuffer);
 		return true;
 	}
 
 	@Override
 	public void processingFinished() {
 		// Do nothing.
 	}
 	
 	//----GETTERS
 
 	/**
 	 * @return The list of starting frequencies for each band. In Hertz.
 	 */
 	public float[] getFreqencies() {
 		return frequencies;
 	}
 	
 	/**
 	 * Returns the Constant Q magnitudes calculated for the previous audio
 	 * buffer. Beware: the array is reused for performance reasons. If your need
 	 * to cache your results, please copy the array.
 	 * @return The output buffer with constant q magnitudes. If you for example are
 	 * interested in coefficients between 256 and 1024 Hz (2^8 and 2^10 Hz) and
 	 * you requested 12 bins per octave, you will need 12 bins/octave * 2
 	 * octaves = 24 places in the output buffer.
 	 */
 	public float[] getMagnitudes() {
 		return magnitudes;
 	}
 
 	
 	/**
 	 * Return the Constant Q coefficients calculated for the previous audio
 	 * buffer. Beware: the array is reused for performance reasons. If your need
 	 * to cache your results, please copy the array.
 	 * 
 	 * @return The array with constant q coefficients. If you for example are
 	 *         interested in coefficients between 256 and 1024 Hz (2^8 and 2^10
 	 *         Hz) and you requested 12 bins per octave, you will need 12
 	 *         bins/octave * 2 octaves * 2 entries/bin = 48 places in the output
 	 *         buffer. The coefficient needs two entries in the output buffer
 	 *         since they are complex numbers.
 	 */
 	public float[] getCoefficients() {
 		return coefficients;
 	}
 
 	/**
 	 * @return The number of coefficients, output bands.
 	 */
 	public int getNumberOfOutputBands() {
 		return frequencies.length;
 	}
 
 	/**
 	 * @return The required length the FFT.
 	 */
 	public int getFFTlength() {
 		return fftLength;
 	}
 	
 	/**
 	 * @return the number of bins every octave.
 	 */
 	public int getBinsPerOctave(){
 		return binsPerOctave;
 	}
 }