CLAM-Development
1.1
|
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