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