00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
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