CLAM-Development  1.4.0
LPC_AutoCorrelation.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 "LPC_AutoCorrelation.hxx"
23 #include "Array.hxx"
24 #include "Audio.hxx"
25 #include "Assert.hxx"
26 #include "CLAM_Math.hxx"
27 #include "ProcessingFactory.hxx"
28 #include "SpecTypeFlags.hxx"
29 
30 namespace CLAM
31 {
32 
33 namespace Hidden
34 {
35  static const char * metadata[] = {
36  "key", "LPC_AutoCorrelation",
37  "category", "Analysis",
38  "description", "LPC_AutoCorrelation",
39  0
40  };
41  static FactoryRegistrator<ProcessingFactory, LPC_AutoCorrelation> reg = metadata;
42 }
43 
45 {
46  AddAll();
47  UpdateData();
48  // MRJ: Seems that eleven is a 'wise' number
49  SetOrder( 11 );
50 }
51 
53  : mAudioIn("AudioIn",this)
54  , mLPModelOut("LPModelOut",this)
55  , mSpectrumOut("SpectrumOut",this)
56 {
58 }
59 
61 {
62 }
63 
65 {
66  const Audio & audio = mAudioIn.GetData();
67  LPModel & lpc = mLPModelOut.GetData();
68  Spectrum & spectrum = mSpectrumOut.GetData();
69 
70  lpc.UpdateModelOrder(mCurrentConfig.GetOrder());
71  bool ok = Do(audio, lpc);
72  CLAM::SpecTypeFlags flags;
73  flags.bMagPhase=1;
74  flags.bComplex = 0;
75  spectrum.SetSize( audio.GetSize()/2+1 );
76  spectrum.SetSpectralRange( audio.GetSampleRate()/2 );
77  spectrum.SetType(flags);
78  lpc.ToSpectrum(spectrum);
79 
80  mAudioIn.Consume();
83  return ok;
84 }
85 
86 bool LPC_AutoCorrelation::Do( const Audio& in, LPModel& out )
87 {
88  bool mustUpdateData = false;
89  if ( !out.HasFilterCoefficients() )
90  {
91  out.AddFilterCoefficients();
92  mustUpdateData =true;
93  }
94  if ( !out.HasReflectionCoefficients() )
95  {
96  out.AddReflectionCoefficients();
97  mustUpdateData = true;
98  }
99  if ( !out.HasAvgSqrFilterError() )
100  {
101  out.AddAvgSqrFilterError();
102  mustUpdateData = true;
103  }
104  if ( mustUpdateData )
105  out.UpdateData();
106 
107  return Do( in, out.GetFilterCoefficients(),
108  out.GetReflectionCoefficients(),
109  out.GetAvgSqrFilterError() );
110 }
111 
112 /* this is the "clean" version:
113 bool LPC_AutoCorrelation::Do( const Audio& in, DataArray& A, DataArray& K, TData& E )
114 {
115 
116  if( !AbleToExecute() ) return true;
117 
118  TData * inBuffer = in.GetBuffer().GetPtr();
119  TData * outBuffer = out.GetBuffer().GetPtr();
120  TData norm = 1 / outBuffer[0];
121  for (int k = 0; k < out.GetSize(); k++ )
122  {
123  for (int n = 0; n < in.GetSize(); n++ )
124  {
125  if( n < k ) // k is out of the segment
126  outBuffer[ k ] += 0;
127  else
128  outBuffer[ k ] += inBuffer[ n ] * inBuffer[ n - k ] ;
129  }
130 
131  outBuffer[ k ] *= in.GetSize() ;
132  outBuffer[ k ] *= norm;
133  }
134 }
135 */
136 
137 // The following does the same, but more efficient, by removing the condition
138 // from the for loop
140 {
141  if ( !AbleToExecute() )
142  return false;
143 
144  DataArray R; // autocorrelation coefficients
145 
146  R.Resize( mCurrentConfig.GetOrder() );
147  R.SetSize( mCurrentConfig.GetOrder() );
148 
149  ComputeAutocorrelation( in.GetBuffer(), R );
150 
151  if ( fabs( R[0] ) <= 1e-6 )
152  {
154  DataArray R2( R.GetPtr()+1, R.Size()-1 );
155  SolveSystemByLevinsonDurbin( R2, A, K, E );
156  }
157  else
158  SolveSystemByLevinsonDurbin( R, A, K, E );
159 
160  return true;
161 }
162 
164 {
166 
167  CLAM_ASSERT( mCurrentConfig.HasOrder(),
168  "Invalid configuration object: it must have the 'Order' attribute available" );
169 
170  return true;
171 }
172 
174 static inline CLAM::TData dot_product( const CLAM::TData* in1,
175  const CLAM::TData* in2,
176  const CLAM::TData* endIn )
177 {
178  CLAM::TData accum = 0.0;
179  while ( in1 != endIn )
180  accum+= (*in1++)*(*in2++);
181  return accum;
182 }
183 
185  Array<TData>& acCoeffs)
186 {
187  //unsigned size = pow(2.,Round(log10(2.*signal.GetSize()-1.)/log10(2.)));
188  int k = 0;
189  TData N = TData( signal.Size() );
190  const TData *inBuffer = signal.GetPtr();
191  const TData *endInBuffer = signal.GetPtr() + signal.Size();
192  TData *outBuffer = acCoeffs.GetPtr();
193  const TData *endOutBuffer = acCoeffs.GetPtr()+acCoeffs.Size();
194 
195  const TData *inBuffer2 = NULL;
196 
197  *outBuffer = dot_product( inBuffer, inBuffer, endInBuffer );
198  *outBuffer *= N;
199  TData norm = 1.0/ *outBuffer;
200  *outBuffer *= norm;
201  outBuffer++;
202  k++;
203 
204  while( outBuffer != endOutBuffer )
205  {
206  inBuffer2 = inBuffer;
207  inBuffer += k;
208 
209  *outBuffer = dot_product( inBuffer, inBuffer2, endInBuffer );
210  *outBuffer*= N;
211  *outBuffer *= norm;
212 
213  inBuffer = signal.GetPtr();
214  outBuffer++;
215  k++;
216  }
217 }
218 
220  Array<TData>& A,
221  Array<TData>& K,
222  TData& E)
223 {
224  unsigned order = mCurrentConfig.GetOrder();
225  CLAM_ASSERT( A.Size() == order,
226  "A coefficient array size mismatch!" );
227  CLAM_ASSERT( K.Size() == order,
228  "K coefficient array size mismatch!" );
229 
230  std::vector <TData> Ap(order);
231  E = R[0];
232  A[0] = 1;
233  for( unsigned i = 1 ; i < order; i++ )
234  {
235  K[i] = R[i];
236  for(unsigned j=1; j<i; j++ )
237  K[i] += A[j] * R[i-j];
238 
239  K[i] = - K[i] / E;
240 
241  for( unsigned j=1; j<i; j++ )
242  Ap[j] = A[i-j];
243 
244  for( unsigned j=1; j<i; j++ )
245  A[j] += K[i] * Ap[j];
246 
247  A[i] = K[i];
248  E *= (1-K[i]*K[i]);
249  }
250 }
251 }
252