RTXI  2.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
impresp.cpp
Go to the documentation of this file.
1 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 //
3 // File = impresp.cpp
4 //
5 //
6 
7 #include <math.h>
8 #include <stdlib.h>
9 #include <iostream>
10 
11 #include "impresp.h"
12 #include "complex.h"
13 
14 #ifdef _DEBUG
15 extern std::ofstream DebugFile;
16 #endif
17 
18 //===========================================================
19 // constructors
20 
21 ImpulseResponse::ImpulseResponse(FilterTransFunc* trans_func, int num_resp_pts,
22  double delta_time)
23 {
24  complex k_sub_r_denom, k_sub_r_numer, s_sub_r;
25  int r, n;
26 
27  Delta_Time = delta_time;
28  Num_Resp_Pts = num_resp_pts;
29  Trans_Func = trans_func;
30 
34 
35 #ifdef _DEBUG
36  DebugFile << "Num_Poles = " << Num_Poles << std::endl;
37  DebugFile << "Num_Zeros = " << Num_Zeros << std::endl;
38 #endif
39 
40  R_Max = (Num_Poles + 1) >> 1;
41  Order_Is_Odd = Num_Poles % 2;
42 
43 #ifdef _DEBUG
44  DebugFile << "R_Max = " << R_Max << std::endl;
45 #endif
46 
47  K_Sub_R = new complex[Num_Poles + 1];
48  Sigma = new double[Num_Poles + 1];
49  Omega = new double[Num_Poles + 1];
50 
51  //-----------------------------------
52  // compute Kr coefficients
53 
54  for (r = 1; r <= Num_Poles; r++) {
55 #ifdef _DEBUG
56  DebugFile << "r = " << r << std::endl;
57 #endif
58  k_sub_r_denom = complex(1.0, 0.0);
59  k_sub_r_numer = complex(1.0, 0.0);
60  s_sub_r = Trans_Func->GetPole(r);
61 #ifdef _DEBUG
62  DebugFile << "s_sub_r = " << s_sub_r << std::endl;
63 #endif
64  Sigma[r] = real(s_sub_r);
65  Omega[r] = imag(s_sub_r);
66  for (n = 1; n <= Num_Poles; n++) {
67  if (n == r)
68  continue;
69 #ifdef _DEBUG
70  DebugFile << "pole[" << n << "] = " << (Trans_Func->GetPole(n))
71  << std::endl;
72  DebugFile << "difference term = " << (s_sub_r - (Trans_Func->GetPole(n)))
73  << std::endl;
74 #endif
75  k_sub_r_denom *= (s_sub_r - (Trans_Func->GetPole(n)));
76 #ifdef _DEBUG
77  DebugFile << "k_sub_r_denom = " << k_sub_r_denom << std::endl;
78 #endif
79  }
80  for (n = 1; n <= Num_Zeros; n++) {
81 #ifdef _DEBUG
82  DebugFile << "zero[" << n << "] = " << (Trans_Func->GetZero(n))
83  << std::endl;
84  DebugFile << "difference term = " << (s_sub_r - (Trans_Func->GetZero(n)))
85  << std::endl;
86 #endif
87  k_sub_r_numer *= (s_sub_r - (Trans_Func->GetZero(n)));
88 #ifdef _DEBUG
89  DebugFile << "k_sub_r_numer = " << k_sub_r_numer << std::endl;
90 #endif
91  }
92  K_Sub_R[r] = k_sub_r_numer / k_sub_r_denom;
93 #ifdef _DEBUG
94  DebugFile << "K_Sub_R[" << r << "] = " << K_Sub_R[r] << std::endl;
95 #endif
96  }
97 
98  return;
99 };
100 
101 //=========================================================
102 void
104 {
105  int resp_indx;
106  double h_of_t, time, delta_t;
107 
108  Response_File = new std::ofstream("imp_anal.txt", std::ios::out);
109 
110  //-----------------------------------------------
111  // compute samples of impulse response
112 
113  delta_t = Delta_Time;
114 
115  for (resp_indx = 0; resp_indx < Num_Resp_Pts; resp_indx++) {
116  time = delta_t * resp_indx;
117  h_of_t = ComputeSample(time);
118  (*Response_File) << time << ", " << h_of_t << std::endl;
119  }
120  Response_File->close();
121  return;
122 }
123 
124 //=========================================================
125 double
127 {
128  int r;
129  double cos_part, sin_part;
130  double h_of_t;
131 
132  h_of_t = 0.0;
133 
134  for (r = 1; r <= (Num_Poles >> 1); r++) {
135  cos_part =
136  2 * real(K_Sub_R[r]) * exp(Sigma[r] * time) * cos(Omega[r] * time);
137  sin_part =
138  2 * imag(K_Sub_R[r]) * exp(Sigma[r] * time) * sin(Omega[r] * time);
139  h_of_t += (cos_part - sin_part);
140  }
141  //---------------------------------------------
142  // add the real exponential component
143  // present in odd-order responses
144 
145  if (Order_Is_Odd == 1) {
146  h_of_t += (real(K_Sub_R[R_Max]) * exp(Sigma[R_Max] * time));
147  }
148  h_of_t *= H_Sub_Zero;
149 
150  return (h_of_t);
151 }
double Delta_Time
Definition: impresp.h:25
int Num_Resp_Pts
Definition: impresp.h:26
std::ofstream * Response_File
Definition: impresp.h:27
double * Omega
Definition: impresp.h:30
double H_Sub_Zero
Definition: impresp.h:28
complex GetPole(int pole_indx)
Definition: filtfunc.cpp:71
void GenerateResponse(void)
Definition: impresp.cpp:103
int Order_Is_Odd
Definition: impresp.h:32
int GetNumPoles(void)
Definition: filtfunc.cpp:206
std::ofstream DebugFile
int GetNumZeros(void)
Definition: filtfunc.cpp:217
complex * K_Sub_R
Definition: impresp.h:29
FilterTransFunc * Trans_Func
Definition: impresp.h:24
ImpulseResponse(FilterTransFunc *trans_func, int num_resp_pts, double delta_t)
Definition: impresp.cpp:21
float GetHSubZero(void)
Definition: filtfunc.cpp:253
double * Sigma
Definition: impresp.h:30
complex GetZero(int zero_indx)
Definition: filtfunc.cpp:90
double imag(const complex &_z)
Definition: complex.h:103
double real(const complex &_z)
Definition: complex.h:98
double ComputeSample(double time)
Definition: impresp.cpp:126