qm-dsp  1.8
DetectionFunction.cpp
Go to the documentation of this file.
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 "DetectionFunction.h"
17 #include <cstring>
18 
20 // Construction/Destruction
22 
24  m_window(0)
25 {
29  m_magPeaks = NULL;
30 
31  initialise( Config );
32 }
33 
35 {
36  deInitialise();
37 }
38 
39 
41 {
42  m_dataLength = Config.frameLength;
43  m_halfLength = m_dataLength/2 + 1;
44 
45  m_DFType = Config.DFType;
46  m_stepSize = Config.stepSize;
47 
48  m_whiten = Config.adaptiveWhitening;
51  if (m_whitenRelaxCoeff < 0) m_whitenRelaxCoeff = 0.9997;
52  if (m_whitenFloor < 0) m_whitenFloor = 0.01;
53 
54  m_magHistory = new double[ m_halfLength ];
55  memset(m_magHistory,0, m_halfLength*sizeof(double));
56 
57  m_phaseHistory = new double[ m_halfLength ];
58  memset(m_phaseHistory,0, m_halfLength*sizeof(double));
59 
60  m_phaseHistoryOld = new double[ m_halfLength ];
61  memset(m_phaseHistoryOld,0, m_halfLength*sizeof(double));
62 
63  m_magPeaks = new double[ m_halfLength ];
64  memset(m_magPeaks,0, m_halfLength*sizeof(double));
65 
67 
68  m_magnitude = new double[ m_halfLength ];
69  m_thetaAngle = new double[ m_halfLength ];
70  m_unwrapped = new double[ m_halfLength ];
71 
73  m_windowed = new double[ m_dataLength ];
74 }
75 
77 {
78  delete [] m_magHistory ;
79  delete [] m_phaseHistory ;
80  delete [] m_phaseHistoryOld ;
81  delete [] m_magPeaks ;
82 
83  delete m_phaseVoc;
84 
85  delete [] m_magnitude;
86  delete [] m_thetaAngle;
87  delete [] m_windowed;
88  delete [] m_unwrapped;
89 
90  delete m_window;
91 }
92 
93 double DetectionFunction::processTimeDomain(const double *samples)
94 {
95  m_window->cut(samples, m_windowed);
96 
99 
100  if (m_whiten) whiten();
101 
102  return runDF();
103 }
104 
105 double DetectionFunction::processFrequencyDomain(const double *reals,
106  const double *imags)
107 {
108  m_phaseVoc->processFrequencyDomain(reals, imags,
110 
111  if (m_whiten) whiten();
112 
113  return runDF();
114 }
115 
117 {
118  for (unsigned int i = 0; i < m_halfLength; ++i) {
119  double m = m_magnitude[i];
120  if (m < m_magPeaks[i]) {
121  m = m + (m_magPeaks[i] - m) * m_whitenRelaxCoeff;
122  }
123  if (m < m_whitenFloor) m = m_whitenFloor;
124  m_magPeaks[i] = m;
125  m_magnitude[i] /= m;
126  }
127 }
128 
130 {
131  double retVal = 0;
132 
133  switch( m_DFType )
134  {
135  case DF_HFC:
136  retVal = HFC( m_halfLength, m_magnitude);
137  break;
138 
139  case DF_SPECDIFF:
140  retVal = specDiff( m_halfLength, m_magnitude);
141  break;
142 
143  case DF_PHASEDEV:
144  // Using the instantaneous phases here actually provides the
145  // same results (for these calculations) as if we had used
146  // unwrapped phases, but without the possible accumulation of
147  // phase error over time
148  retVal = phaseDev( m_halfLength, m_thetaAngle);
149  break;
150 
151  case DF_COMPLEXSD:
153  break;
154 
155  case DF_BROADBAND:
156  retVal = broadband( m_halfLength, m_magnitude);
157  break;
158  }
159 
160  return retVal;
161 }
162 
163 double DetectionFunction::HFC(unsigned int length, double *src)
164 {
165  unsigned int i;
166  double val = 0;
167 
168  for( i = 0; i < length; i++)
169  {
170  val += src[ i ] * ( i + 1);
171  }
172  return val;
173 }
174 
175 double DetectionFunction::specDiff(unsigned int length, double *src)
176 {
177  unsigned int i;
178  double val = 0.0;
179  double temp = 0.0;
180  double diff = 0.0;
181 
182  for( i = 0; i < length; i++)
183  {
184  temp = fabs( (src[ i ] * src[ i ]) - (m_magHistory[ i ] * m_magHistory[ i ]) );
185 
186  diff= sqrt(temp);
187 
188  // (See note in phaseDev below.)
189 
190  val += diff;
191 
192  m_magHistory[ i ] = src[ i ];
193  }
194 
195  return val;
196 }
197 
198 
199 double DetectionFunction::phaseDev(unsigned int length, double *srcPhase)
200 {
201  unsigned int i;
202  double tmpPhase = 0;
203  double tmpVal = 0;
204  double val = 0;
205 
206  double dev = 0;
207 
208  for( i = 0; i < length; i++)
209  {
210  tmpPhase = (srcPhase[ i ]- 2*m_phaseHistory[ i ]+m_phaseHistoryOld[ i ]);
211  dev = MathUtilities::princarg( tmpPhase );
212 
213  // A previous version of this code only counted the value here
214  // if the magnitude exceeded 0.1. My impression is that
215  // doesn't greatly improve the results for "loud" music (so
216  // long as the peak picker is reasonably sophisticated), but
217  // does significantly damage its ability to work with quieter
218  // music, so I'm removing it and counting the result always.
219  // Same goes for the spectral difference measure above.
220 
221  tmpVal = fabs(dev);
222  val += tmpVal ;
223 
224  m_phaseHistoryOld[ i ] = m_phaseHistory[ i ] ;
225  m_phaseHistory[ i ] = srcPhase[ i ];
226  }
227 
228  return val;
229 }
230 
231 
232 double DetectionFunction::complexSD(unsigned int length, double *srcMagnitude, double *srcPhase)
233 {
234  unsigned int i;
235  double val = 0;
236  double tmpPhase = 0;
237  double tmpReal = 0;
238  double tmpImag = 0;
239 
240  double dev = 0;
241  ComplexData meas = ComplexData( 0, 0 );
242  ComplexData j = ComplexData( 0, 1 );
243 
244  for( i = 0; i < length; i++)
245  {
246  tmpPhase = (srcPhase[ i ]- 2*m_phaseHistory[ i ]+m_phaseHistoryOld[ i ]);
247  dev= MathUtilities::princarg( tmpPhase );
248 
249  meas = m_magHistory[i] - ( srcMagnitude[ i ] * exp( j * dev) );
250 
251  tmpReal = real( meas );
252  tmpImag = imag( meas );
253 
254  val += sqrt( (tmpReal * tmpReal) + (tmpImag * tmpImag) );
255 
256  m_phaseHistoryOld[ i ] = m_phaseHistory[ i ] ;
257  m_phaseHistory[ i ] = srcPhase[ i ];
258  m_magHistory[ i ] = srcMagnitude[ i ];
259  }
260 
261  return val;
262 }
263 
264 double DetectionFunction::broadband(unsigned int length, double *src)
265 {
266  double val = 0;
267  for (unsigned int i = 0; i < length; ++i) {
268  double sqrmag = src[i] * src[i];
269  if (m_magHistory[i] > 0.0) {
270  double diff = 10.0 * log10(sqrmag / m_magHistory[i]);
271  if (diff > m_dbRise) val = val + 1;
272  }
273  m_magHistory[i] = sqrmag;
274  }
275  return val;
276 }
277 
279 {
280  return m_magnitude;
281 }
282 
double phaseDev(unsigned int length, double *srcPhase)
double processFrequencyDomain(const double *reals, const double *imags)
Process a single frequency-domain frame, provided as frameLength/2+1 real and imaginary component val...
double whiteningFloor
Window< double > * m_window
double specDiff(unsigned int length, double *src)
unsigned int m_stepSize
unsigned int m_halfLength
double processTimeDomain(const double *samples)
Process a single time-domain frame of audio, provided as frameLength samples.
bool adaptiveWhitening
void processTimeDomain(const double *src, double *mag, double *phase, double *unwrapped)
Given one frame of time-domain samples, FFT and return the magnitudes, instantaneous phases...
void initialise(DFConfig Config)
#define DF_HFC
PhaseVocoder * m_phaseVoc
unsigned int frameLength
unsigned int m_dataLength
double complexSD(unsigned int length, double *srcMagnitude, double *srcPhase)
unsigned int stepSize
#define NULL
Definition: Filter.h:20
static double princarg(double ang)
The principle argument function.
double whiteningRelaxCoeff
#define DF_PHASEDEV
DetectionFunction(DFConfig Config)
#define DF_COMPLEXSD
#define DF_BROADBAND
complex< double > ComplexData
Definition: MathAliases.h:23
void processFrequencyDomain(const double *reals, const double *imags, double *mag, double *phase, double *unwrapped)
Given one frame of frequency-domain samples, return the magnitudes, instantaneous phases...
double * getSpectrumMagnitude()
double HFC(unsigned int length, double *src)
double broadband(unsigned int length, double *srcMagnitude)
#define DF_SPECDIFF
void cut(T *src) const
Definition: Window.h:61