c@225: /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
c@225: 
c@225: /*
c@225:     QM DSP Library
c@225: 
c@225:     Centre for Digital Music, Queen Mary, University of London.
c@309:     This file 2005-2006 Christian Landone.
c@309: 
mathieu@321:     Modifications:
mathieu@321: 
mathieu@321:     - delta threshold
mathieu@321:     Description: add delta threshold used as offset in the smoothed
mathieu@321:     detection function
mathieu@321:     Author: Mathieu Barthet
mathieu@321:     Date: June 2010
mathieu@321: 
c@309:     This program is free software; you can redistribute it and/or
c@309:     modify it under the terms of the GNU General Public License as
c@309:     published by the Free Software Foundation; either version 2 of the
c@309:     License, or (at your option) any later version.  See the file
c@309:     COPYING included with this distribution for more information.
c@225: */
c@225: 
c@225: #include "DFProcess.h"
c@241: #include "maths/MathUtilities.h"
c@225: 
c@272: #include <cstring>
c@434: #include <algorithm>
c@434: 
c@272: 
c@225: //////////////////////////////////////////////////////////////////////
c@225: // Construction/Destruction
c@225: //////////////////////////////////////////////////////////////////////
c@225: 
cannam@503: DFProcess::DFProcess( DFProcConfig config )
c@225: {
c@225:     filtSrc = NULL;
cannam@483:     filtDst = NULL;     
c@225:     m_filtScratchIn = NULL;
c@225:     m_filtScratchOut = NULL;
c@225: 
c@225:     m_FFOrd = 0;
c@225: 
cannam@503:     initialise( config );
c@225: }
c@225: 
c@225: DFProcess::~DFProcess()
c@225: {
c@225:     deInitialise();
c@225: }
c@225: 
cannam@503: void DFProcess::initialise( DFProcConfig config )
c@225: {
cannam@503:     m_length = config.length;
cannam@503:     m_winPre = config.winPre;
cannam@503:     m_winPost = config.winPost;
cannam@503:     m_alphaNormParam = config.AlphaNormParam;
c@225: 
cannam@503:     m_isMedianPositive = config.isMedianPositive;
c@225: 
c@225:     filtSrc = new double[ m_length ];
c@225:     filtDst = new double[ m_length ];
c@225: 
c@417:     Filter::Parameters params;
cannam@483:     params.a = std::vector<double>
cannam@504:         (config.LPACoeffs, config.LPACoeffs + config.LPOrd + 1);
cannam@483:     params.b = std::vector<double>
cannam@504:         (config.LPBCoeffs, config.LPBCoeffs + config.LPOrd + 1);
c@417:     
c@417:     m_FiltFilt = new FiltFilt(params);
cannam@483:         
mathieu@321:     //add delta threshold
cannam@503:     m_delta = config.delta;
c@225: }
c@225: 
c@225: void DFProcess::deInitialise()
c@225: {
c@225:     delete [] filtSrc;
c@225:     delete [] filtDst;
c@225:     delete [] m_filtScratchIn;
c@225:     delete [] m_filtScratchOut;
c@225:     delete m_FiltFilt;
c@225: }
c@225: 
c@225: void DFProcess::process(double *src, double* dst)
c@225: {
c@283:     if (m_length == 0) return;
c@283: 
c@225:     removeDCNormalize( src, filtSrc );
c@225: 
c@225:     m_FiltFilt->process( filtSrc, filtDst, m_length );
c@225: 
c@225:     medianFilter( filtDst, dst );
c@225: }
c@225: 
c@225: 
c@225: void DFProcess::medianFilter(double *src, double *dst)
c@225: {
c@299:     int i,k,j,l;
c@299:     int index = 0;
c@225: 
c@225:     double val = 0;
c@225: 
c@225:     double* y = new double[ m_winPost + m_winPre + 1];
c@225:     memset( y, 0, sizeof( double ) * ( m_winPost + m_winPre + 1) );
c@225: 
c@225:     double* scratch = new double[ m_length ];
c@225: 
cannam@483:     for( i = 0; i < m_winPre; i++) {
cannam@483:         
cannam@483:         if (index >= m_length) {
cannam@483:             break;
cannam@483:         }
c@299: 
cannam@483:         k = i + m_winPost + 1;
c@225: 
cannam@483:         for( j = 0; j < k; j++) {
cannam@483:             y[ j ] = src[ j ];
cannam@483:         }
cannam@483:         scratch[ index ] = MathUtilities::median( y, k );
cannam@483:         index++;
c@225:     }
c@225: 
cannam@483:     for(  i = 0; i + m_winPost + m_winPre < m_length; i ++) {
cannam@483:         
cannam@483:         if (index >= m_length) {
cannam@483:             break;
cannam@483:         }
cannam@483:                          
cannam@483:         l = 0;
cannam@483:         for(  j  = i; j < ( i + m_winPost + m_winPre + 1); j++) {
cannam@483:             y[ l ] = src[ j ];
cannam@483:             l++;
cannam@483:         }
c@299: 
cannam@483:         scratch[index] = MathUtilities::median( y, (m_winPost + m_winPre + 1 ));
cannam@483:         index++;
c@225:     }
c@225: 
cannam@483:     for( i = std::max( m_length - m_winPost, 1); i < m_length; i++) {
cannam@483:         
cannam@483:         if (index >= m_length) {
cannam@483:             break;
cannam@483:         }
c@299: 
cannam@483:         k = std::max( i - m_winPre, 1);
c@225: 
cannam@483:         l = 0;
cannam@483:         for( j = k; j < m_length; j++) {
cannam@483:             y[ l ] = src[ j ];
cannam@483:             l++;
cannam@483:         }
cannam@483:                 
cannam@483:         scratch[index] = MathUtilities::median( y, l);
cannam@483:         index++;
c@225:     }
c@225: 
cannam@483:     for( i = 0; i < m_length; i++ ) {
cannam@483:         //add a delta threshold used as an offset when computing the smoothed detection function
cannam@483:         //(helps to discard noise when detecting peaks) 
cannam@483:         val = src[ i ] - scratch[ i ] - m_delta;
cannam@483:                 
cannam@483:         if( m_isMedianPositive ) {
cannam@483:             if( val > 0 ) {
cannam@483:                 dst[ i ]  = val;
cannam@483:             } else {
cannam@483:                 dst[ i ]  = 0;
cannam@483:             }
cannam@483:         } else {
cannam@483:             dst[ i ]  = val;
cannam@483:         }
c@225:     }
cannam@483:         
c@225:     delete [] y;
c@225:     delete [] scratch;
c@225: }
c@225: 
c@225: 
c@225: void DFProcess::removeDCNormalize( double *src, double*dst )
c@225: {
c@225:     double DFmax = 0;
c@225:     double DFMin = 0;
c@225:     double DFAlphaNorm = 0;
c@225: 
c@225:     MathUtilities::getFrameMinMax( src, m_length, &DFMin, &DFmax );
c@225: 
c@225:     MathUtilities::getAlphaNorm( src, m_length, m_alphaNormParam, &DFAlphaNorm );
c@225: 
cannam@483:     for (int i = 0; i < m_length; i++) {
cannam@483:         dst[ i ] = ( src[ i ] - DFMin ) / DFAlphaNorm; 
c@225:     }
c@225: }