Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members

burg.h

Go to the documentation of this file.
00001 /*
00002  * SPUC - Signal processing using C++ - A DSP library
00003  * 
00004  * This program is free software; you can redistribute it and/or modify
00005  * it under the terms of the GNU General Public License as published by
00006  * the Free Software Foundation; either version 2, or (at your option)
00007  * any later version.
00008  * 
00009  * This program is distributed in the hope that it will be useful,
00010  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00011  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00012  * GNU General Public License for more details.
00013  * 
00014  * You should have received a copy of the GNU General Public License
00015  * along with this program; if not, write to the Free Software
00016  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
00017 */
00018 #include <math.h>
00019 #include <vec.h>
00020 namespace SPUC {
00024 template <class T> Vector<T> burg( Vector<T>& x, int P) {
00025         const double EPS = 1e-30;
00026                 Vector<T> a(P);
00027                 int N = x.size();
00028                 T* ef = new T[N];
00029                 T* ef_prev = new T[N];
00030                 T* eb = new T[N];
00031                 T* eb_prev = new T[N];
00032                 T*   aa = new T[P];
00033                 T* rc = new T[N];
00034                 
00035                 T gamma;
00036                 
00037                 T num, den ;
00038                 int t, p ;
00039                 int i;
00040                 
00041                 for (i=0;i<N;i++) {     ef[i] = eb[i] = x[i];  }
00042                 
00043                 for (p=0;p<P;p++)  {
00044                         num = 0.0;
00045                         den = 0.0;
00046                         
00047                         for (t=p+1;t<N;t++) { 
00048                                 den += ef[t]*ef[t] + eb[t-1]*eb[t-1];
00049                                 num += ef[t] * eb[t-1]; 
00050                         }
00051                         
00052                         if (magsq(den)>EPS) gamma=2.0*num/den;
00053                         else gamma  = 0;
00054                         
00055                         a[p] = gamma;
00056                         rc[p]=-a[p];
00057                         
00058                         for (i=0;i<N;i++) {
00059                                 ef_prev[i] = ef[i];
00060                                 eb_prev[i] = eb[i];
00061                         }
00062                         
00063                         for (t=1;t<N;t++)  {
00064                                 ef[t] +=  -(conj(gamma) * eb_prev[t-1]);
00065                                 eb[t]  = eb_prev[t-1] - (gamma * ef_prev[t]);
00066                         }
00067           
00068                         if (p>0) for(t=0;t<p;t++) a[t]= aa[t]  + rc[p]*aa[p-t-1];
00069                         for (i=0;i<p+1;i++) aa[i] = a[i];
00070                 }
00071         
00072                 delete [] ef;
00073                 delete [] eb;
00074                 delete [] ef_prev;
00075                 delete [] eb_prev;
00076                 delete [] aa;
00077                 delete [] rc;
00078                 return(a);
00079         }
00080 }
00081 

Generated on Fri Sep 16 11:02:25 2005 for spuc by  doxygen 1.4.4