c@225: /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
c@225: /*
c@225:     QM DSP Library
c@225: 
c@225:     Centre for Digital Music, Queen Mary, University of London.
c@225:     This file copyright 2005-2006 Christian Landone.
c@225:     All rights reserved.
c@225: */
c@225: 
c@225: #include "ConstantQ.h"
c@225: #include "dsp/transforms/FFT.h"
c@225: 
c@225: //---------------------------------------------------------------------------
c@225: // nextpow2 returns the smallest integer n such that 2^n >= x.
c@225: static double nextpow2(double x) {
c@225:     double y = ceil(log(x)/log(2.0));
c@225:     return(y);
c@225: }
c@225: 
c@225: static double squaredModule(const double & xx, const double & yy) {
c@225:     return xx*xx + yy*yy;
c@225: }
c@225: 
c@225: //----------------------------------------------------------------------------
c@225: 
c@225: ConstantQ::ConstantQ( CQConfig Config ) 
c@225: {
c@225:     initialise( Config );
c@225: }
c@225: 
c@225: ConstantQ::~ConstantQ()
c@225: {
c@225:     deInitialise();
c@225: }
c@225: 
c@225: //----------------------------------------------------------------------------
c@225: void ConstantQ::sparsekernel()
c@225: {
c@225:     //generates spectral kernel matrix (upside down?)
c@225:     // initialise temporal kernel with zeros, twice length to deal w. complex numbers
c@225: 
c@225:     double* hammingWindowRe = new double [ m_FFTLength ];
c@225:     double* hammingWindowIm = new double [ m_FFTLength ];
c@225:     double* transfHammingWindowRe = new double [ m_FFTLength ];
c@225:     double* transfHammingWindowIm = new double [ m_FFTLength ];
c@225: 
c@225:     for (unsigned u=0; u < m_FFTLength; u++) 
c@225:     {
c@225: 	hammingWindowRe[u] = 0;
c@225: 	hammingWindowIm[u] = 0;
c@225:     }
c@225: 
c@225: 
c@225:     // Here, fftleng*2 is a guess of the number of sparse cells in the matrix
c@225:     // The matrix K x fftlength but the non-zero cells are an antialiased
c@225:     // square root function. So mostly is a line, with some grey point.
c@225:     m_sparseKernelIs.reserve( m_FFTLength*2 );
c@225:     m_sparseKernelJs.reserve( m_FFTLength*2 );
c@225:     m_sparseKernelRealValues.reserve( m_FFTLength*2 );
c@225:     m_sparseKernelImagValues.reserve( m_FFTLength*2 );
c@225: 	
c@225:     // for each bin value K, calculate temporal kernel, take its fft to
c@225:     //calculate the spectral kernel then threshold it to make it sparse and 
c@225:     //add it to the sparse kernels matrix
c@225:     double squareThreshold = m_CQThresh * m_CQThresh;
c@225: 
c@225:     FFT m_FFT;
c@225: 	
c@225:     for (unsigned k = m_uK; k--; ) 
c@225:     {
c@225: 	// Computing a hamming window
c@225: 	const unsigned hammingLength = (int) ceil( m_dQ * m_FS / ( m_FMin * pow(2,((double)(k))/(double)m_BPO)));
c@225: 	for (unsigned i=0; i<hammingLength; i++) 
c@225: 	{
c@225: 	    const double angle = 2*PI*m_dQ*i/hammingLength;
c@225: 	    const double real = cos(angle);
c@225: 	    const double imag = sin(angle);
c@225: 	    const double absol = hamming(hammingLength, i)/hammingLength;
c@225: 	    hammingWindowRe[ i ] = absol*real;
c@225: 	    hammingWindowIm[ i ] = absol*imag;
c@225: 	}
c@225: 
c@225: 	//do fft of hammingWindow
c@225: 	m_FFT.process( m_FFTLength, 0, hammingWindowRe, hammingWindowIm, transfHammingWindowRe, transfHammingWindowIm );
c@225: 
c@225: 		
c@225: 	for (unsigned j=0; j<( m_FFTLength ); j++) 
c@225: 	{
c@225: 	    // perform thresholding
c@225: 	    const double squaredBin = squaredModule( transfHammingWindowRe[ j ], transfHammingWindowIm[ j ]);
c@225: 	    if (squaredBin <= squareThreshold) continue;
c@225: 		
c@225: 	    // Insert non-zero position indexes, doubled because they are floats
c@225: 	    m_sparseKernelIs.push_back(j);
c@225: 	    m_sparseKernelJs.push_back(k);
c@225: 
c@225: 	    // take conjugate, normalise and add to array sparkernel
c@225: 	    m_sparseKernelRealValues.push_back( transfHammingWindowRe[ j ]/m_FFTLength);
c@225: 	    m_sparseKernelImagValues.push_back(-transfHammingWindowIm[ j ]/m_FFTLength);
c@225: 	}
c@225: 
c@225:     }
c@225: 
c@225:     delete [] hammingWindowRe;
c@225:     delete [] hammingWindowIm;
c@225:     delete [] transfHammingWindowRe;
c@225:     delete [] transfHammingWindowIm;
c@225: 
c@225: }
c@225: 
c@225: //-----------------------------------------------------------------------------
c@225: double* ConstantQ::process( double* fftdata )
c@225: {
c@225:     for (unsigned row=0; row<2*m_uK; row++) 
c@225:     {
c@225: 	m_CQdata[ row ] = 0;
c@225: 	m_CQdata[ row+1 ] = 0;
c@225:     }
c@225:     const unsigned *fftbin = &(m_sparseKernelIs[0]);
c@225:     const unsigned *cqbin  = &(m_sparseKernelJs[0]);
c@225:     const double   *real   = &(m_sparseKernelRealValues[0]);
c@225:     const double   *imag   = &(m_sparseKernelImagValues[0]);
c@225:     const unsigned int sparseCells = m_sparseKernelRealValues.size();
c@225: 	
c@225:     for (unsigned i = 0; i<sparseCells; i++)
c@225:     {
c@225: 	const unsigned row = cqbin[i];
c@225: 	const unsigned col = fftbin[i];
c@225: 	const double & r1  = real[i];
c@225: 	const double & i1  = imag[i];
c@225: 	const double & r2  = fftdata[ (2*m_FFTLength) - 2*col];
c@225: 	const double & i2  = fftdata[ (2*m_FFTLength) - 2*col+1];
c@225: 	// add the multiplication
c@225: 	m_CQdata[ 2*row  ] += (r1*r2 - i1*i2);
c@225: 	m_CQdata[ 2*row+1] += (r1*i2 + i1*r2);
c@225:     }
c@225: 
c@225:     return m_CQdata;
c@225: }
c@225: 
c@225: 
c@225: void ConstantQ::initialise( CQConfig Config )
c@225: {
c@225:     m_FS = Config.FS;
c@225:     m_FMin = Config.min;		// min freq
c@225:     m_FMax = Config.max;		// max freq
c@225:     m_BPO = Config.BPO;		// bins per octave
c@225:     m_CQThresh = Config.CQThresh;// ConstantQ threshold for kernel generation
c@225: 
c@225:     m_dQ = 1/(pow(2,(1/(double)m_BPO))-1);	// Work out Q value for Filter bank
c@225:     m_uK = (unsigned int) ceil(m_BPO * log(m_FMax/m_FMin)/log(2.0));	// No. of constant Q bins
c@225: 
c@225:     // work out length of fft required for this constant Q Filter bank
c@225:     m_FFTLength = (int) pow(2, nextpow2(ceil( m_dQ*m_FS/m_FMin )));
c@225: 
c@225:     m_hop = m_FFTLength/8; // <------ hop size is window length divided by 32
c@225: 
c@225:     // allocate memory for cqdata
c@225:     m_CQdata = new double [2*m_uK];
c@225: }
c@225: 
c@225: void ConstantQ::deInitialise()
c@225: {
c@225:     delete [] m_CQdata;
c@225: }
c@225: 
c@225: void ConstantQ::process(double *FFTRe, double* FFTIm, double *CQRe, double *CQIm)
c@225: {
c@225:     for (unsigned row=0; row<m_uK; row++) 
c@225:     {
c@225: 	CQRe[ row ] = 0;
c@225: 	CQIm[ row ] = 0;
c@225:     }
c@225: 
c@225:     const unsigned *fftbin = &(m_sparseKernelIs[0]);
c@225:     const unsigned *cqbin  = &(m_sparseKernelJs[0]);
c@225:     const double   *real   = &(m_sparseKernelRealValues[0]);
c@225:     const double   *imag   = &(m_sparseKernelImagValues[0]);
c@225:     const unsigned int sparseCells = m_sparseKernelRealValues.size();
c@225: 	
c@225:     for (unsigned i = 0; i<sparseCells; i++)
c@225:     {
c@225: 	const unsigned row = cqbin[i];
c@225: 	const unsigned col = fftbin[i];
c@225: 	const double & r1  = real[i];
c@225: 	const double & i1  = imag[i];
c@225: 	const double & r2  = FFTRe[ m_FFTLength- col];
c@225: 	const double & i2  = FFTIm[ m_FFTLength - col];
c@225: 	// add the multiplication
c@225: 	CQRe[ row ] += (r1*r2 - i1*i2);
c@225: 	CQIm[ row ] += (r1*i2 + i1*r2);
c@225:     }
c@225: }