/*
*
      
_______
                       
_____
   
_____ _____
  

*
     
|__
   
__|
                     
|
  
__ \ / ____|__ \
 

*
        
| | __ _ _ __ ___
  
______| || | (___ | |__) |
*
        
| |/ _` | '__/ __|/ _ \/ __| |
  
| |\___ \|___/
 

*
        
| | (_| | |
  
\__ \ (_) \__ \ |__| |____) | |
     

*
        
|_|\__,_|_|
  
|___/\___/|___/_____/|_____/|_|
     

*
                                                         

* -------------------------------------------------------------
*
* 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;
	
}
}