SpectralPeakDetect.cxx

Go to the documentation of this file.
00001 /*
00002  * Copyright (c) 2001-2004 MUSIC TECHNOLOGY GROUP (MTG)
00003  *                         UNIVERSITAT POMPEU FABRA
00004  *
00005  *
00006  * This program is free software; you can redistribute it and/or modify
00007  * it under the terms of the GNU General Public License as published by
00008  * the Free Software Foundation; either version 2 of the License, or
00009  * (at your option) any later version.
00010  *
00011  * This program is distributed in the hope that it will be useful,
00012  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00013  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00014  * GNU General Public License for more details.
00015  *
00016  * You should have received a copy of the GNU General Public License
00017  * along with this program; if not, write to the Free Software
00018  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00019  *
00020  */
00021 
00022 #include "Complex.hxx"
00023 #include "Spectrum.hxx"
00024 #include "SpectralPeakArray.hxx"
00025 #include "SpectralPeakDetect.hxx"
00026 
00027 
00028 
00029 namespace CLAM {
00030 
00031         /* Processing  object Method  implementations */
00032 
00033         SpectralPeakDetect::SpectralPeakDetect()
00034                 : mInput( "Input spectrum", this ),
00035                   mOutput( "Output spectral peak array", this )
00036         {
00037                 Configure(SpectralPeakDetectConfig());
00038         }
00039 
00040         SpectralPeakDetect::SpectralPeakDetect(const SpectralPeakDetectConfig &c = SpectralPeakDetectConfig())
00041                 : mInput( "Input spectrum", this ),
00042                   mOutput( "Output spectral peak array", this )
00043         {
00044                 Configure(c);
00045         }
00046 
00047         SpectralPeakDetect::~SpectralPeakDetect()
00048         {}
00049 
00050 
00051         /* Configure the Processing Object according to the Config object */
00052 
00053         bool SpectralPeakDetect::ConcreteConfigure(const ProcessingConfig& c)
00054         {
00055 
00056                 CopyAsConcreteConfig(mConfig, c);
00057                 return true;
00058         }
00059 
00060         /* Setting Prototypes for faster processing */
00061 
00062         bool SpectralPeakDetect::SetPrototypes(Spectrum& inputs,const SpectralPeakArray& out)
00063         {
00064                 return true;
00065         }
00066 
00067         bool SpectralPeakDetect::SetPrototypes()
00068         {
00069                 return true;
00070         }
00071         
00072         bool SpectralPeakDetect::UnsetPrototypes()
00073         {
00074                 return true;
00075         }
00076 
00077         bool  SpectralPeakDetect::Do(void)
00078         {
00079                 bool result = false;
00080 
00081                 mOutput.GetData().SetScale( EScale::eLog );
00082 
00083                 if (mInput.GetData().GetScale() != EScale::eLog)
00084                 {
00085                         mTmpLinearInSpectrum = mInput.GetData();
00086                         mTmpLinearInSpectrum.ToDB();
00087                         result = Do( mTmpLinearInSpectrum, mOutput.GetData() );
00088                 }
00089                 else
00090                 {
00091                         result = Do( mInput.GetData(), mOutput.GetData() );
00092                 }
00093 
00094                 mInput.Consume();
00095                 mOutput.Produce();
00096                 return result;
00097         }
00098 
00099         /* The  unsupervised Do() function */
00100 
00101         bool  SpectralPeakDetect::Do(const Spectrum& input, SpectralPeakArray& out)
00102         {
00103                 CLAM_ASSERT(CheckInputType(input), "SpectralPeakDetect::Do() - Type of input data doesn't match expected type.");
00104                 CLAM_ASSERT(CheckOutputType(out), "SpectralPeakDetect::Do() - Type of output data doesn't match expected type.");
00105 
00106                 TData middleMag;
00107                 TData leftMag;
00108                 TData rightMag;
00109                 const TData samplingRate = input.GetSpectralRange() * TData(2.0);
00110                 const TData magThreshold = mConfig.GetMagThreshold();
00111                 const TSize nBins = input.GetSize();
00112                 const TData maxFreq= mConfig.GetMaxFreq();
00113 
00114                 const DataArray& inMagBuffer=input.GetMagBuffer();
00115                 const DataArray& inPhaseBuffer=input.GetPhaseBuffer();
00116 
00117                 TSize maxPeaks=mConfig.GetMaxPeaks();
00118                 if (out.GetnMaxPeaks() != maxPeaks)
00119                 {
00120                         out.SetnMaxPeaks(maxPeaks);
00121                 }
00122 
00123                 out.SetnPeaks(0);
00124                 
00125                 DataArray& outMagBuffer=out.GetMagBuffer();
00126                 DataArray& outFreqBuffer=out.GetFreqBuffer();
00127                 DataArray& outPhaseBuffer=out.GetPhaseBuffer();
00128                 DataArray& outBinPosBuffer=out.GetBinPosBuffer();
00129                 DataArray& outBinWidthBuffer=out.GetBinWidthBuffer();
00130                 
00131                 // detection loop 
00132                 TSize nSpectralPeaks = 0;
00133                 TSize binWidth = 0;      // BinWidth is in NumBins
00134 
00135                 int i;
00136                 for (i = 0; (i < nBins-2) && (nSpectralPeaks < maxPeaks); ++i)
00137                 {
00138                         leftMag         = inMagBuffer[i];
00139                         middleMag       = inMagBuffer[i+1];
00140                         rightMag        = inMagBuffer[i+2];
00141 
00142                         // local constant detected 
00143                         if (middleMag == leftMag && leftMag == rightMag) {
00144 
00145                                 if ((nSpectralPeaks > 0) && (outBinWidthBuffer[nSpectralPeaks-1] == 0)) {
00146                         
00147                                         outBinWidthBuffer[nSpectralPeaks-1]=TData(binWidth); // store last SpectralPeakBinWidth
00148                                 }
00149                                 binWidth = 1; // Reset Binwidth
00150                                 continue;
00151                         }
00152                 
00153                         // local Minimum detected 
00154                         if ((middleMag <= leftMag) && (middleMag<= rightMag)) {
00155                                 if ((nSpectralPeaks > 0) && (outBinWidthBuffer[nSpectralPeaks-1] == 0)) {
00156                                         outBinWidthBuffer[nSpectralPeaks-1]=TData(binWidth); // store last SpectralPeakBinWidth
00157                                 }
00158                                 binWidth = 0; // Reset Binwidth
00159                         }
00160                         
00161                         TData interpolatedBin;
00162                         TData spectralPeakFreq=0;
00163                         TData spectralPeakPhase;
00164                         TData spectralPeakMag;
00165                         // local maximum Detected ! 
00166                         if ((middleMag >= leftMag) && (middleMag >= rightMag) && (middleMag > magThreshold)) {  
00167 
00168                                 TSize SpectralPeakPosition = i+1;       // middleMag has index i+1
00169 
00170                                 // update last BinWidth 
00171                                 if ((nSpectralPeaks > 0) && (outBinWidthBuffer[nSpectralPeaks-1] == 0)) {                       
00172                         
00173                                         TSize lastSpectralPeakBin = (TSize) (outFreqBuffer[nSpectralPeaks-1]*2* (nBins-1)/samplingRate);
00174                                         TSize tempVal = binWidth - (TSize)((SpectralPeakPosition-lastSpectralPeakBin)/2.0);
00175                                         outBinWidthBuffer[nSpectralPeaks-1]=TData(tempVal);
00176                                         binWidth = (TSize) ((SpectralPeakPosition-lastSpectralPeakBin)/2.0);
00177                                 }       
00178 
00179                                 // if we get to the end of a constant area then ...
00180                                 if ((middleMag == leftMag) && (middleMag > rightMag) && (nSpectralPeaks > 0)){
00181                         
00182                                         // update last SpectralPeak BinPosition, it will be located in the middle of the constant area 
00183                                         interpolatedBin = (double) outBinPosBuffer[nSpectralPeaks-1] + (double) (i+1-outBinPosBuffer[nSpectralPeaks-1])/2; // center BinPos 
00184                                         spectralPeakFreq = interpolatedBin * samplingRate / 2 / (nBins-1);
00185                         
00186                                         outFreqBuffer[nSpectralPeaks-1]=spectralPeakFreq;
00187                                         outMagBuffer[nSpectralPeaks-1]=input.GetMag(spectralPeakFreq);
00188                                         outPhaseBuffer[nSpectralPeaks-1]=input.GetPhase(spectralPeakFreq);
00189                                         outBinWidthBuffer[nSpectralPeaks-1]=int(SpectralPeakPosition-outBinPosBuffer[nSpectralPeaks-1]);
00190                                         outBinPosBuffer[nSpectralPeaks-1]= interpolatedBin; // interpolated BinPos is stored    
00191                                 }
00192 
00193                                 else { //  add SpectralPeak... BinWidth will be updated in the next turn
00194                         
00195                                         // Estimating the ``true'' maximum peak (frequency and magnitude) of the detected local maximum 
00196                                         // using a parabolic cure-fitting. The idea is that the main-lobe of spectrum of most analysis 
00197                                         // windows on a dB scale looks like a parabola and therefore the maximum of a parabola fitted 
00198                                         // through a local maxima bin and it's two neighboring bins will give a good approximation 
00199                                         // of the actual frequency and magnitude of a sinusoid in the input signal.
00200                                         //
00201                                         // The parabola f(x) = a(x-n)^2 + b(x-n) + c can be completely described using 3 points; 
00202                                         // f(n-1) = A1, f(n) = A2 and f(n+1) = A3, where 
00203                                         // A1 = 20log10(|X(n-1)|), A2 = 20log10(|X(n)|), A3 = 20log10(|X(n+1)|).
00204                                         //
00205                                         // Solving these equation yields: a = 1/2*A1 - A2 + 1/2*A3, b = 1/2*A3 - 1/2*A1 and 
00206                                         // c = A2.
00207                                         //
00208                                         // As the 3 bins are known to be a maxima, solving d/dx f(x) = 0, yields (fractional) bin 
00209                                         // position x of the estimated peak. Substituting delta_x for (x-n) in this equation yields 
00210                                         // the fractional offset in bins from n where the peak's maximum is.
00211                                         //
00212                                         // Solving this equation yields: delta_x = 1/2 * (A1 - A3)/(A1 - 2*A2 + A3).
00213                                         // 
00214                                         // Computing f(n+delta_x) will estimate the peak's magnitude (in dB's):
00215                                         // f(n+delta_x) = A2 - 1/4*(A1-A3)*delta_x.
00216 
00217                                         const TData diffFromMax = TData(0.5) * ((leftMag-rightMag) / (leftMag- 2*middleMag + rightMag));
00218                                         interpolatedBin = SpectralPeakPosition+diffFromMax;
00219                         
00220                                         // store SpectralPeak data 
00221                                         spectralPeakFreq  = interpolatedBin * samplingRate / 2 / (nBins-1); // interpolate Frequency 
00222                                         spectralPeakMag   = middleMag-TData(.25)*(leftMag-rightMag)*diffFromMax;
00223 
00224                                         TData leftPhase,rightPhase;
00225                                         if (diffFromMax>=0)
00226                                         {
00227                                                 leftPhase = inPhaseBuffer[i+1];
00228                                                 rightPhase = inPhaseBuffer[i+2];
00229                                         }
00230                                         else
00231                                         {
00232                                                 leftPhase = inPhaseBuffer[i];
00233                                                 rightPhase = inPhaseBuffer[i+1];
00234                                         }
00235 
00236                                         if (fabs(rightPhase-leftPhase)>PI)
00237                                         {
00238                                                 if (rightPhase>0)
00239                                                         leftPhase+=TData(TWO_PI);
00240                                                 else
00241                                                         rightPhase+=TData(TWO_PI);
00242                                         }
00243 
00244                                         if (diffFromMax >= 0)
00245                                                 spectralPeakPhase = leftPhase + diffFromMax*(rightPhase-leftPhase);
00246                                         else
00247                                                 spectralPeakPhase = leftPhase + (1+diffFromMax)*(rightPhase-leftPhase);
00248 
00249                                         if (spectralPeakFreq>maxFreq)
00250                                                 break;
00251 
00252                                         outFreqBuffer.AddElem(spectralPeakFreq);
00253                                         outMagBuffer.AddElem(spectralPeakMag);
00254                                         outPhaseBuffer.AddElem(spectralPeakPhase);
00255                                         outBinPosBuffer.AddElem(interpolatedBin);
00256                                         outBinWidthBuffer.AddElem(0); // BinWidth will be set later
00257 
00258                                         nSpectralPeaks++;
00259                                 
00260                                 }
00261                         }
00262                         binWidth++;
00263                 }
00264                 
00265                 // update the very last binwidth value if it's not set yet 
00266                 if ((nSpectralPeaks > 0) && (outBinWidthBuffer[nSpectralPeaks-1] == 0)){                        
00267                 
00268                         TSize lastSpectralPeakBin = (TSize) (outFreqBuffer[nSpectralPeaks-1] * 2 * nBins / samplingRate);
00269                         TSize tempVal = binWidth - (TSize)((i-lastSpectralPeakBin)/2.0);
00270                         outBinWidthBuffer[nSpectralPeaks-1]=TData(tempVal);
00271                         binWidth = (TSize) ((i-lastSpectralPeakBin)/2.0);
00272                 }
00273 
00274 
00275                 //TODO: I don't know if we could also change nMaxPeaks if we have found less
00276                 //but then we would be resizing array in every call to the Do and that is not
00277                 //very nice either.
00278                 
00279                 //All this is not necessary, it is not doing anything
00280 /*
00281                 if(nSpectralPeaks>maxPeaks)
00282                         out.SetnMaxPeaks(nSpectralPeaks);
00283                 out.SetnPeaks(nSpectralPeaks);
00284 */
00285                 return true;
00286         }
00287 
00288 
00289         bool SpectralPeakDetect::CheckInputType(const Spectrum &in)
00290         {
00291                 if (!in.HasScale())
00292                         return false;
00293 
00294                 if (in.GetScale() != EScale::eLog)
00295                         return false;
00296 
00297                 if (!in.HasSpectralRange())
00298                         return false;
00299 
00300                 if (!in.HasMagBuffer())
00301                         return false;
00302 
00303                 if (!in.HasPhaseBuffer())
00304                         return false;
00305 
00306                 return true;
00307         }
00308 
00309         bool SpectralPeakDetect::CheckOutputType(const SpectralPeakArray &out) 
00310         {
00311                 if (!out.HasScale())
00312                         return false;
00313 
00314                 if (out.GetScale() != EScale::eLog)
00315                         return false;
00316 
00317                 if (!out.HasBinWidthBuffer())
00318                         return false;
00319 
00320                 if (!out.HasFreqBuffer())
00321                         return false;
00322 
00323                 if (!out.HasBinPosBuffer())
00324                         return false;
00325 
00326                 if (!out.HasMagBuffer())
00327                         return false;
00328 
00329                 if (!out.HasPhaseBuffer())
00330                         return false;
00331 
00332                 return true;
00333         }
00334 
00335 }
00336 
Generated by  doxygen 1.6.3