LPC_AutoCorrelation.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 "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 
Generated by  doxygen 1.6.3