CLAM-Development  1.4.0
SpectralPeakDetect.cxx
Go to the documentation of this file.
1 /*
2  * Copyright (c) 2001-2004 MUSIC TECHNOLOGY GROUP (MTG)
3  * UNIVERSITAT POMPEU FABRA
4  *
5  *
6  * This program is free software; you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation; either version 2 of the License, or
9  * (at your option) any later version.
10  *
11  * This program is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with this program; if not, write to the Free Software
18  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19  *
20  */
21 
22 #include "Complex.hxx"
23 #include "Spectrum.hxx"
24 #include "SpectralPeakArray.hxx"
25 #include "SpectralPeakDetect.hxx"
26 
27 
28 
29 namespace CLAM {
30 
31  /* Processing object Method implementations */
32 
34  : mInput( "Input spectrum", this ),
35  mOutput( "Output spectral peak array", this )
36  {
38  }
39 
41  : mInput( "Input spectrum", this ),
42  mOutput( "Output spectral peak array", this )
43  {
44  Configure(c);
45  }
46 
48  {}
49 
50 
51  /* Configure the Processing Object according to the Config object */
52 
54  {
55 
57  return true;
58  }
59 
60  /* Setting Prototypes for faster processing */
61 
63  {
64  return true;
65  }
66 
68  {
69  return true;
70  }
71 
73  {
74  return true;
75  }
76 
78  {
79  bool result = false;
80 
81  mOutput.GetData().SetScale( EScale::eLog );
82 
83  if (mInput.GetData().GetScale() != EScale::eLog)
84  {
85  mTmpLinearInSpectrum = mInput.GetData();
86  mTmpLinearInSpectrum.ToDB();
87  result = Do( mTmpLinearInSpectrum, mOutput.GetData() );
88  }
89  else
90  {
91  result = Do( mInput.GetData(), mOutput.GetData() );
92  }
93 
94  mInput.Consume();
95  mOutput.Produce();
96  return result;
97  }
98 
99  /* The unsupervised Do() function */
100 
102  {
103  CLAM_ASSERT(CheckInputType(input), "SpectralPeakDetect::Do() - Type of input data doesn't match expected type.");
104  CLAM_ASSERT(CheckOutputType(out), "SpectralPeakDetect::Do() - Type of output data doesn't match expected type.");
105 
106  TData middleMag;
107  TData leftMag;
108  TData rightMag;
109  const TData samplingRate = input.GetSpectralRange() * TData(2.0);
110  const TData magThreshold = mConfig.GetMagThreshold();
111  const TSize nBins = input.GetSize();
112  const TData maxFreq= mConfig.GetMaxFreq();
113 
114  const DataArray& inMagBuffer=input.GetMagBuffer();
115  const DataArray& inPhaseBuffer=input.GetPhaseBuffer();
116 
117  TSize maxPeaks=mConfig.GetMaxPeaks();
118  if (out.GetnMaxPeaks() != maxPeaks)
119  {
120  out.SetnMaxPeaks(maxPeaks);
121  }
122 
123  out.SetnPeaks(0);
124 
125  DataArray& outMagBuffer=out.GetMagBuffer();
126  DataArray& outFreqBuffer=out.GetFreqBuffer();
127  DataArray& outPhaseBuffer=out.GetPhaseBuffer();
128  DataArray& outBinPosBuffer=out.GetBinPosBuffer();
129  DataArray& outBinWidthBuffer=out.GetBinWidthBuffer();
130 
131  // detection loop
132  TSize nSpectralPeaks = 0;
133  TSize binWidth = 0; // BinWidth is in NumBins
134 
135  int i;
136  for (i = 0; (i < nBins-2) && (nSpectralPeaks < maxPeaks); ++i)
137  {
138  leftMag = inMagBuffer[i];
139  middleMag = inMagBuffer[i+1];
140  rightMag = inMagBuffer[i+2];
141 
142  // local constant detected
143  if (middleMag == leftMag && leftMag == rightMag) {
144 
145  if ((nSpectralPeaks > 0) && (outBinWidthBuffer[nSpectralPeaks-1] == 0)) {
146 
147  outBinWidthBuffer[nSpectralPeaks-1]=TData(binWidth); // store last SpectralPeakBinWidth
148  }
149  binWidth = 1; // Reset Binwidth
150  continue;
151  }
152 
153  // local Minimum detected
154  if ((middleMag <= leftMag) && (middleMag<= rightMag)) {
155  if ((nSpectralPeaks > 0) && (outBinWidthBuffer[nSpectralPeaks-1] == 0)) {
156  outBinWidthBuffer[nSpectralPeaks-1]=TData(binWidth); // store last SpectralPeakBinWidth
157  }
158  binWidth = 0; // Reset Binwidth
159  }
160 
161  TData interpolatedBin;
162  TData spectralPeakFreq=0;
163  TData spectralPeakPhase;
164  TData spectralPeakMag;
165  // local maximum Detected !
166  if ((middleMag >= leftMag) && (middleMag >= rightMag) && (middleMag > magThreshold)) {
167 
168  TSize SpectralPeakPosition = i+1; // middleMag has index i+1
169 
170  // update last BinWidth
171  if ((nSpectralPeaks > 0) && (outBinWidthBuffer[nSpectralPeaks-1] == 0)) {
172 
173  TSize lastSpectralPeakBin = (TSize) (outFreqBuffer[nSpectralPeaks-1]*2* (nBins-1)/samplingRate);
174  TSize tempVal = binWidth - (TSize)((SpectralPeakPosition-lastSpectralPeakBin)/2.0);
175  outBinWidthBuffer[nSpectralPeaks-1]=TData(tempVal);
176  binWidth = (TSize) ((SpectralPeakPosition-lastSpectralPeakBin)/2.0);
177  }
178 
179  // if we get to the end of a constant area then ...
180  if ((middleMag == leftMag) && (middleMag > rightMag) && (nSpectralPeaks > 0)){
181 
182  // update last SpectralPeak BinPosition, it will be located in the middle of the constant area
183  interpolatedBin = (double) outBinPosBuffer[nSpectralPeaks-1] + (double) (i+1-outBinPosBuffer[nSpectralPeaks-1])/2; // center BinPos
184  spectralPeakFreq = interpolatedBin * samplingRate / 2 / (nBins-1);
185 
186  outFreqBuffer[nSpectralPeaks-1]=spectralPeakFreq;
187  outMagBuffer[nSpectralPeaks-1]=input.GetMag(spectralPeakFreq);
188  outPhaseBuffer[nSpectralPeaks-1]=input.GetPhase(spectralPeakFreq);
189  outBinWidthBuffer[nSpectralPeaks-1]=int(SpectralPeakPosition-outBinPosBuffer[nSpectralPeaks-1]);
190  outBinPosBuffer[nSpectralPeaks-1]= interpolatedBin; // interpolated BinPos is stored
191  }
192 
193  else { // add SpectralPeak... BinWidth will be updated in the next turn
194 
195  // Estimating the ``true'' maximum peak (frequency and magnitude) of the detected local maximum
196  // using a parabolic cure-fitting. The idea is that the main-lobe of spectrum of most analysis
197  // windows on a dB scale looks like a parabola and therefore the maximum of a parabola fitted
198  // through a local maxima bin and it's two neighboring bins will give a good approximation
199  // of the actual frequency and magnitude of a sinusoid in the input signal.
200  //
201  // The parabola f(x) = a(x-n)^2 + b(x-n) + c can be completely described using 3 points;
202  // f(n-1) = A1, f(n) = A2 and f(n+1) = A3, where
203  // A1 = 20log10(|X(n-1)|), A2 = 20log10(|X(n)|), A3 = 20log10(|X(n+1)|).
204  //
205  // Solving these equation yields: a = 1/2*A1 - A2 + 1/2*A3, b = 1/2*A3 - 1/2*A1 and
206  // c = A2.
207  //
208  // As the 3 bins are known to be a maxima, solving d/dx f(x) = 0, yields (fractional) bin
209  // position x of the estimated peak. Substituting delta_x for (x-n) in this equation yields
210  // the fractional offset in bins from n where the peak's maximum is.
211  //
212  // Solving this equation yields: delta_x = 1/2 * (A1 - A3)/(A1 - 2*A2 + A3).
213  //
214  // Computing f(n+delta_x) will estimate the peak's magnitude (in dB's):
215  // f(n+delta_x) = A2 - 1/4*(A1-A3)*delta_x.
216 
217  const TData diffFromMax = TData(0.5) * ((leftMag-rightMag) / (leftMag- 2*middleMag + rightMag));
218  interpolatedBin = SpectralPeakPosition+diffFromMax;
219 
220  // store SpectralPeak data
221  spectralPeakFreq = interpolatedBin * samplingRate / 2 / (nBins-1); // interpolate Frequency
222  spectralPeakMag = middleMag-TData(.25)*(leftMag-rightMag)*diffFromMax;
223 
224  TData leftPhase,rightPhase;
225  if (diffFromMax>=0)
226  {
227  leftPhase = inPhaseBuffer[i+1];
228  rightPhase = inPhaseBuffer[i+2];
229  }
230  else
231  {
232  leftPhase = inPhaseBuffer[i];
233  rightPhase = inPhaseBuffer[i+1];
234  }
235 
236  if (fabs(rightPhase-leftPhase)>PI)
237  {
238  if (rightPhase>0)
239  leftPhase+=TData(TWO_PI);
240  else
241  rightPhase+=TData(TWO_PI);
242  }
243 
244  if (diffFromMax >= 0)
245  spectralPeakPhase = leftPhase + diffFromMax*(rightPhase-leftPhase);
246  else
247  spectralPeakPhase = leftPhase + (1+diffFromMax)*(rightPhase-leftPhase);
248 
249  if (spectralPeakFreq>maxFreq)
250  break;
251 
252  outFreqBuffer.AddElem(spectralPeakFreq);
253  outMagBuffer.AddElem(spectralPeakMag);
254  outPhaseBuffer.AddElem(spectralPeakPhase);
255  outBinPosBuffer.AddElem(interpolatedBin);
256  outBinWidthBuffer.AddElem(0); // BinWidth will be set later
257 
258  nSpectralPeaks++;
259 
260  }
261  }
262  binWidth++;
263  }
264 
265  // update the very last binwidth value if it's not set yet
266  if ((nSpectralPeaks > 0) && (outBinWidthBuffer[nSpectralPeaks-1] == 0)){
267 
268  TSize lastSpectralPeakBin = (TSize) (outFreqBuffer[nSpectralPeaks-1] * 2 * nBins / samplingRate);
269  TSize tempVal = binWidth - (TSize)((i-lastSpectralPeakBin)/2.0);
270  outBinWidthBuffer[nSpectralPeaks-1]=TData(tempVal);
271  binWidth = (TSize) ((i-lastSpectralPeakBin)/2.0);
272  }
273 
274 
275  //TODO: I don't know if we could also change nMaxPeaks if we have found less
276  //but then we would be resizing array in every call to the Do and that is not
277  //very nice either.
278 
279  //All this is not necessary, it is not doing anything
280 /*
281  if(nSpectralPeaks>maxPeaks)
282  out.SetnMaxPeaks(nSpectralPeaks);
283  out.SetnPeaks(nSpectralPeaks);
284 */
285  return true;
286  }
287 
288 
289  bool SpectralPeakDetect::CheckInputType(const Spectrum &in)
290  {
291  if (!in.HasScale())
292  return false;
293 
294  if (in.GetScale() != EScale::eLog)
295  return false;
296 
297  if (!in.HasSpectralRange())
298  return false;
299 
300  if (!in.HasMagBuffer())
301  return false;
302 
303  if (!in.HasPhaseBuffer())
304  return false;
305 
306  return true;
307  }
308 
309  bool SpectralPeakDetect::CheckOutputType(const SpectralPeakArray &out)
310  {
311  if (!out.HasScale())
312  return false;
313 
314  if (out.GetScale() != EScale::eLog)
315  return false;
316 
317  if (!out.HasBinWidthBuffer())
318  return false;
319 
320  if (!out.HasFreqBuffer())
321  return false;
322 
323  if (!out.HasBinPosBuffer())
324  return false;
325 
326  if (!out.HasMagBuffer())
327  return false;
328 
329  if (!out.HasPhaseBuffer())
330  return false;
331 
332  return true;
333  }
334 
335 }
336