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 "LPC_AutoCorrelation.hxx" 00023 #include "Array.hxx" 00024 #include "Audio.hxx" 00025 #include "Assert.hxx" 00026 #include "CLAM_Math.hxx" 00027 #include "ProcessingFactory.hxx" 00028 #include "SpecTypeFlags.hxx" 00029 00030 namespace CLAM 00031 { 00032 00033 namespace Hidden 00034 { 00035 static const char * metadata[] = { 00036 "key", "LPC_AutoCorrelation", 00037 "category", "Analysis", 00038 "description", "LPC_AutoCorrelation", 00039 0 00040 }; 00041 static FactoryRegistrator<ProcessingFactory, LPC_AutoCorrelation> reg = metadata; 00042 } 00043 00044 void LPCConfig::DefaultInit() 00045 { 00046 AddAll(); 00047 UpdateData(); 00048 // MRJ: Seems that eleven is a 'wise' number 00049 SetOrder( 11 ); 00050 } 00051 00052 LPC_AutoCorrelation::LPC_AutoCorrelation() 00053 : mAudioIn("AudioIn",this) 00054 , mLPModelOut("LPModelOut",this) 00055 , mSpectrumOut("SpectrumOut",this) 00056 { 00057 Configure(LPCConfig()); 00058 } 00059 00060 LPC_AutoCorrelation::~LPC_AutoCorrelation() 00061 { 00062 } 00063 00064 bool LPC_AutoCorrelation::Do() 00065 { 00066 const Audio & audio = mAudioIn.GetData(); 00067 LPModel & lpc = mLPModelOut.GetData(); 00068 Spectrum & spectrum = mSpectrumOut.GetData(); 00069 00070 lpc.UpdateModelOrder(mCurrentConfig.GetOrder()); 00071 bool ok = Do(audio, lpc); 00072 CLAM::SpecTypeFlags flags; 00073 flags.bMagPhase=1; 00074 flags.bComplex = 0; 00075 spectrum.SetSize( audio.GetSize()/2+1 ); 00076 spectrum.SetSpectralRange( audio.GetSampleRate()/2 ); 00077 spectrum.SetType(flags); 00078 lpc.ToSpectrum(spectrum); 00079 00080 mAudioIn.Consume(); 00081 mLPModelOut.Produce(); 00082 mSpectrumOut.Produce(); 00083 return ok; 00084 } 00085 00086 bool LPC_AutoCorrelation::Do( const Audio& in, LPModel& out ) 00087 { 00088 bool mustUpdateData = false; 00089 if ( !out.HasFilterCoefficients() ) 00090 { 00091 out.AddFilterCoefficients(); 00092 mustUpdateData =true; 00093 } 00094 if ( !out.HasReflectionCoefficients() ) 00095 { 00096 out.AddReflectionCoefficients(); 00097 mustUpdateData = true; 00098 } 00099 if ( !out.HasAvgSqrFilterError() ) 00100 { 00101 out.AddAvgSqrFilterError(); 00102 mustUpdateData = true; 00103 } 00104 if ( mustUpdateData ) 00105 out.UpdateData(); 00106 00107 return Do( in, out.GetFilterCoefficients(), 00108 out.GetReflectionCoefficients(), 00109 out.GetAvgSqrFilterError() ); 00110 } 00111 00112 /* this is the "clean" version: 00113 bool LPC_AutoCorrelation::Do( const Audio& in, DataArray& A, DataArray& K, TData& E ) 00114 { 00115 00116 if( !AbleToExecute() ) return true; 00117 00118 TData * inBuffer = in.GetBuffer().GetPtr(); 00119 TData * outBuffer = out.GetBuffer().GetPtr(); 00120 TData norm = 1 / outBuffer[0]; 00121 for (int k = 0; k < out.GetSize(); k++ ) 00122 { 00123 for (int n = 0; n < in.GetSize(); n++ ) 00124 { 00125 if( n < k ) // k is out of the segment 00126 outBuffer[ k ] += 0; 00127 else 00128 outBuffer[ k ] += inBuffer[ n ] * inBuffer[ n - k ] ; 00129 } 00130 00131 outBuffer[ k ] *= in.GetSize() ; 00132 outBuffer[ k ] *= norm; 00133 } 00134 } 00135 */ 00136 00137 // The following does the same, but more efficient, by removing the condition 00138 // from the for loop 00139 bool LPC_AutoCorrelation::Do( const Audio& in, DataArray& A, DataArray& K, TData& E ) 00140 { 00141 if ( !AbleToExecute() ) 00142 return false; 00143 00144 DataArray R; // autocorrelation coefficients 00145 00146 R.Resize( mCurrentConfig.GetOrder() ); 00147 R.SetSize( mCurrentConfig.GetOrder() ); 00148 00149 ComputeAutocorrelation( in.GetBuffer(), R ); 00150 00151 if ( fabs( R[0] ) <= 1e-6 ) 00152 { 00154 DataArray R2( R.GetPtr()+1, R.Size()-1 ); 00155 SolveSystemByLevinsonDurbin( R2, A, K, E ); 00156 } 00157 else 00158 SolveSystemByLevinsonDurbin( R, A, K, E ); 00159 00160 return true; 00161 } 00162 00163 bool LPC_AutoCorrelation::ConcreteConfigure( const ProcessingConfig& cfg ) 00164 { 00165 CopyAsConcreteConfig( mCurrentConfig, cfg ); 00166 00167 CLAM_ASSERT( mCurrentConfig.HasOrder(), 00168 "Invalid configuration object: it must have the 'Order' attribute available" ); 00169 00170 return true; 00171 } 00172 00174 static inline CLAM::TData dot_product( const CLAM::TData* in1, 00175 const CLAM::TData* in2, 00176 const CLAM::TData* endIn ) 00177 { 00178 CLAM::TData accum = 0.0; 00179 while ( in1 != endIn ) 00180 accum+= (*in1++)*(*in2++); 00181 return accum; 00182 } 00183 00184 void LPC_AutoCorrelation::ComputeAutocorrelation(const Array<TData>& signal, 00185 Array<TData>& acCoeffs) 00186 { 00187 //unsigned size = pow(2.,Round(log10(2.*signal.GetSize()-1.)/log10(2.))); 00188 int k = 0; 00189 TData N = TData( signal.Size() ); 00190 const TData *inBuffer = signal.GetPtr(); 00191 const TData *endInBuffer = signal.GetPtr() + signal.Size(); 00192 TData *outBuffer = acCoeffs.GetPtr(); 00193 const TData *endOutBuffer = acCoeffs.GetPtr()+acCoeffs.Size(); 00194 00195 const TData *inBuffer2 = NULL; 00196 00197 *outBuffer = dot_product( inBuffer, inBuffer, endInBuffer ); 00198 *outBuffer *= N; 00199 TData norm = 1.0/ *outBuffer; 00200 *outBuffer *= norm; 00201 outBuffer++; 00202 k++; 00203 00204 while( outBuffer != endOutBuffer ) 00205 { 00206 inBuffer2 = inBuffer; 00207 inBuffer += k; 00208 00209 *outBuffer = dot_product( inBuffer, inBuffer2, endInBuffer ); 00210 *outBuffer*= N; 00211 *outBuffer *= norm; 00212 00213 inBuffer = signal.GetPtr(); 00214 outBuffer++; 00215 k++; 00216 } 00217 } 00218 00219 void LPC_AutoCorrelation::SolveSystemByLevinsonDurbin( const Array<TData>& R, 00220 Array<TData>& A, 00221 Array<TData>& K, 00222 TData& E) 00223 { 00224 unsigned order = mCurrentConfig.GetOrder(); 00225 CLAM_ASSERT( A.Size() == order, 00226 "A coefficient array size mismatch!" ); 00227 CLAM_ASSERT( K.Size() == order, 00228 "K coefficient array size mismatch!" ); 00229 00230 std::vector <TData> Ap(order); 00231 E = R[0]; 00232 A[0] = 1; 00233 for( unsigned i = 1 ; i < order; i++ ) 00234 { 00235 K[i] = R[i]; 00236 for(unsigned j=1; j<i; j++ ) 00237 K[i] += A[j] * R[i-j]; 00238 00239 K[i] = - K[i] / E; 00240 00241 for( unsigned j=1; j<i; j++ ) 00242 Ap[j] = A[i-j]; 00243 00244 for( unsigned j=1; j<i; j++ ) 00245 A[j] += K[i] * Ap[j]; 00246 00247 A[i] = K[i]; 00248 E *= (1-K[i]*K[i]); 00249 } 00250 } 00251 } 00252