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 
PeakPicking(PPickParams Config)35 PeakPicking::PeakPicking( PPickParams Config )
36 {
37     m_workBuffer = NULL;
38     initialise( Config );
39 }
40 
~PeakPicking()41 PeakPicking::~PeakPicking()
42 {
43     deInitialise();
44 }
45 
initialise(PPickParams Config)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 
deInitialise()69 void PeakPicking::deInitialise()
70 {
71     delete [] m_workBuffer;
72     delete m_DFSmoothing;
73     m_workBuffer = NULL;
74 }
75 
process(double * src,int len,vector<int> & onsets)76 void PeakPicking::process( double* src, 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 (int i = 0; i < len; i++) {
86         m_maxima.push_back( m_workBuffer[ i ] );
87     }
88 
89     quadEval( m_maxima, onsets );
90 
91     for( int b = 0; b < (int)m_maxima.size(); b++) {
92         src[ b ] = m_maxima[ b ];
93     }
94 }
95 
quadEval(vector<double> & src,vector<int> & idx)96 int PeakPicking::quadEval( vector<double> &src, vector<int> &idx )
97 {
98     int maxLength;
99 
100     vector <int> m_maxIndex;
101     vector <int> m_onsetPosition;
102 
103     vector <double> m_maxFit;
104     vector <double> m_poly;
105     vector <double> m_err;
106 
107     m_poly.push_back(0);
108     m_poly.push_back(0);
109     m_poly.push_back(0);
110 
111     for (int t = -2; t < 3; t++) {
112         m_err.push_back( (double)t );
113     }
114 
115     for (int i = 2; i < int(src.size()) - 2; i++) {
116         if ((src[i] > src[i-1]) && (src[i] > src[i+1]) && (src[i] > 0) ) {
117             m_maxIndex.push_back(i);
118         }
119     }
120 
121     maxLength = int(m_maxIndex.size());
122 
123     double selMax = 0;
124 
125     for (int j = 0; j < maxLength ; j++) {
126         for (int k = -2; k <= 2; ++k) {
127             selMax = src[ m_maxIndex[j] + k ] ;
128             m_maxFit.push_back(selMax);
129         }
130 
131         TPolyFit::PolyFit2(m_err, m_maxFit, m_poly);
132 
133         double f = m_poly[0];
134         double h = m_poly[2];
135 
136         if (h < -Qfilta || f > Qfiltc) {
137             idx.push_back(m_maxIndex[j]);
138         }
139 
140         m_maxFit.clear();
141     }
142 
143     return 1;
144 }
145