libs/qm-dsp/dsp/onsets/PeakPicking.cpp

   1 /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
   2
   3 /*
   4     QM DSP Library
   5
   6     Centre for Digital Music, Queen Mary, University of London.
   7     This file 2005-2006 Christian Landone.
   8
   9     Modifications:
  10
  11     - delta threshold
  12     Description: add delta threshold used as offset in the smoothed
  13     detection function
  14     Author: Mathieu Barthet
  15     Date: June 2010
  16
  17     This program is free software; you can redistribute it and/or
  18     modify it under the terms of the GNU General Public License as
  19     published by the Free Software Foundation; either version 2 of the
  20     License, or (at your option) any later version.  See the file
  21     COPYING included with this distribution for more information.
  22 */
  23
  24 #include "PeakPicking.h"
  25 #include "maths/Polyfit.h"
  26
  27 #include <iostream>
  28 #include <cstring>
  29
  30
  31 //////////////////////////////////////////////////////////////////////
  32 // Construction/Destruction
  33 //////////////////////////////////////////////////////////////////////
  34
  35 PeakPicking::PeakPicking( PPickParams Config )
  36 {
  37     m_workBuffer = NULL;
  38     initialise( Config );
  39 }
  40
  41 PeakPicking::~PeakPicking()
  42 {
  43     deInitialise();
  44 }
  45
  46 void PeakPicking::initialise( PPickParams Config )
  47 {
  48     m_DFLength = Config.length ;
  49     Qfilta = Config.QuadThresh.a ;
  50     Qfiltb = Config.QuadThresh.b ;
  51     Qfiltc = Config.QuadThresh.c ;
  52
  53     m_DFProcessingParams.length = m_DFLength;
  54     m_DFProcessingParams.LPOrd = Config.LPOrd;
  55     m_DFProcessingParams.LPACoeffs = Config.LPACoeffs;
  56     m_DFProcessingParams.LPBCoeffs = Config.LPBCoeffs;
  57     m_DFProcessingParams.winPre  = Config.WinT.pre;
  58     m_DFProcessingParams.winPost = Config.WinT.post;
  59     m_DFProcessingParams.AlphaNormParam = Config.alpha;
  60     m_DFProcessingParams.isMedianPositive = false;
  61     m_DFProcessingParams.delta = Config.delta; //add the delta threshold as an adjustable parameter
  62
  63     m_DFSmoothing = new DFProcess( m_DFProcessingParams );
  64
  65     m_workBuffer = new double[ m_DFLength ];
  66     memset( m_workBuffer, 0, sizeof(double)*m_DFLength);
  67 }
  68
  69 void PeakPicking::deInitialise()
  70 {
  71     delete [] m_workBuffer;
  72     delete m_DFSmoothing;
  73     m_workBuffer = NULL;
  74 }
  75
  76 void PeakPicking::process( double* src, unsigned int len, vector<int> &onsets )
  77 {
  78     if (len < 4) return;
  79
  80     vector <double> m_maxima;
  81
  82     // Signal conditioning
  83     m_DFSmoothing->process( src, m_workBuffer );
  84
  85     for( unsigned int u = 0; u < len; u++)
  86     {
  87         m_maxima.push_back( m_workBuffer[ u ] );
  88     }
  89
  90     quadEval( m_maxima, onsets );
  91
  92     for( int b = 0; b <  (int)m_maxima.size(); b++)
  93     {
  94         src[ b ] = m_maxima[ b ];
  95     }
  96 }
  97
  98 int PeakPicking::quadEval( vector<double> &src, vector<int> &idx )
  99 {
 100     unsigned int maxLength;
 101
 102     vector <int> m_maxIndex;
 103     vector <int> m_onsetPosition;
 104
 105     vector <double> m_maxFit;
 106     vector <double> m_poly;
 107     vector <double> m_err;
 108
 109     m_poly.push_back(0);
 110     m_poly.push_back(0);
 111     m_poly.push_back(0);
 112
 113     for(  int t = -2; t < 3; t++)
 114     {
 115         m_err.push_back( (double)t );
 116     }
 117     for( unsigned int i = 2; i < src.size() - 2; i++)
 118     {
 119         if( (src[i] > src[i-1]) && (src[i] > src[i+1]) && (src[i] > 0) )
 120         {
 121 //          m_maxIndex.push_back(  i + 1 );
 122             m_maxIndex.push_back(i);
 123         }
 124     }
 125
 126     maxLength = m_maxIndex.size();
 127
 128     double selMax = 0;
 129
 130     for( unsigned int j = 0; j < maxLength ; j++)
 131     {
 132         for (int k = -2; k <= 2; ++k)
 133         {
 134             selMax = src[ m_maxIndex[j] + k ] ;
 135             m_maxFit.push_back(selMax);
 136         }
 137
 138         TPolyFit::PolyFit2(m_err, m_maxFit, m_poly);
 139
 140         double f = m_poly[0];
 141         double h = m_poly[2];
 142
 143         if (h < -Qfilta || f > Qfiltc)
 144         {
 145             idx.push_back(m_maxIndex[j]);
 146         }
 147
 148         m_maxFit.clear();
 149     }
 150
 151     return 1;
 152 }