RTXI  2.1
burg.cpp
Go to the documentation of this file.
1 //
2 // File = burg.cpp
3 //
4
5 #include <fstream>
6 #include <stdlib.h>
7 #include <iostream>
8
9 #include "burg.h"
10 #include "complex.h"
11 #include "matrix_t.h"
13 #include "yulewalk.h"
14
15 using namespace std;
16
17 //==============================================
18
19 template <class T>
20 BurgMethod<T>::BurgMethod(int est_order, T* x, int sig_len)
21 {
22  int i, j, k, n;
23  matrix<T> kappa(1, est_order, 1, est_order);
24  T numer, denom;
25  Ar_Order = est_order;
26  Noise_Seed = 31415927; // arbitrary default
27
28  A_Coeffs = new T[est_order + 1];
29  double* rho = new double[est_order + 1];
30  T* err_bak = new T[sig_len];
31  T* err_fwd = new T[sig_len];
32  T* new_bak = new T[sig_len];
33  T* new_fwd = new T[sig_len];
34
35  rho[0] = 0.0;
36  for (n = 0; n < sig_len; n++) {
37  rho[0] += mag_sqrd(x[n]);
38  }
39  rho[0] /= sig_len;
41  std::cout << "rho[0] = " << rho[0] << std::endl;
42  std::cout << " enter amount to be subtracted from this value" << std::endl;
45  for (n = 0; n < sig_len - 1; n++) {
46  err_bak[n] = x[n];
47  err_fwd[n + 1] = x[n + 1];
48  }
49  for (k = 1; k <= est_order; k++) {
50  numer = 0.0;
51  denom = 0.0;
52  for (n = k; n < sig_len; n++) {
53  numer += err_fwd[n] * conj(err_bak[n - 1]);
54  denom += mag_sqrd(err_fwd[n]) + mag_sqrd(err_bak[n - 1]);
55  }
56  kappa[k][k] = -2.0 * numer / denom;
57  rho[k] = (1.0 - mag_sqrd(kappa[k][k])) * rho[k - 1];
58  if (est_order == 1)
59  break;
60  if (k > 1) {
61  // update coeffs for prediction error filter
62  for (j = 1; j < k; j++)
63  kappa[j][k] = kappa[j][k - 1] + kappa[k][k] * conj(kappa[k - j][k - 1]);
64  }
65  //
66  // update prediction errors
67  for (i = k; i <= n - 2; i++) {
68  new_bak[i] = err_bak[i - 1] + conj(kappa[k][k]) * err_fwd[i];
69  new_fwd[i + 1] = err_fwd[i + 1] + kappa[k][k] * err_bak[i];
70  }
71  for (i = k; i <= n - 2; i++) {
72  err_bak[i] = new_bak[i];
73  err_fwd[i + 1] = new_fwd[i + 1];
74  }
75  } // end of loop over k
76  Drv_Noise_Var = rho[est_order];
77  for (i = 1; i <= est_order; i++) {
78  A_Coeffs[i] = kappa[i][est_order];
79  std::cout << "A_Coeff[" << i << "] = " << A_Coeffs[i] << std::endl;
80  }
81  delete[] rho;
82  delete[] err_bak;
83  delete[] err_fwd;
84  delete[] new_bak;
85  delete[] new_fwd;
86  return;
87 };
88
89 template <class T>
91 //==================================================
92 // Function to dump AR parameters to output
93 // stream indicated by uout
94
95 template <class T>
96 void
98 {
99  uout << "Drv_Noise_Var = " << Drv_Noise_Var << std::endl;
100  for (int indx = 0; indx <= Ar_Order; indx++) {
101  uout << "a[" << indx << "] = " << A_Coeffs[indx] << std::endl;
102  }
103  return;
104 };
105 //==================================================
106 // Function to copy AR parameters to array
107 // pointed to by a_coeff
108
109 template <class T>
110 void
111 BurgMethod<T>::GetParameters(int* ar_order, T* a_coeff)
112 {
113  *ar_order = Ar_Order;
114  for (int indx = 0; indx <= Ar_Order; indx++) {
115  a_coeff[indx] = A_Coeffs[indx];
116  }
117  return;
118 };
119 //---------------------------------------------------
120 template <class T>
121 double
123 {
124  return (Drv_Noise_Var);
125 }
126
127 template class BurgMethod<complex>;
128 template class BurgMethod<double>;
complex conj(const complex _z)
Definition: complex.h:108
void GetParameters(int *ar_order, T *a_coeff)
Definition: burg.cpp:111
~BurgMethod(void)
Definition: burg.cpp:90
double GetDrivingVariance(void)
Definition: burg.cpp:122
BurgMethod(int est_ar_order, T *sig_seq, int seq_len)
Definition: burg.cpp:20
void DumpParameters(std::ostream &uout)
Definition: burg.cpp:97
double mag_sqrd(const complex _z)
Definition: complex.h:128