-
Notifications
You must be signed in to change notification settings - Fork 1
/
Discrete.cpp
161 lines (123 loc) · 3.47 KB
/
Discrete.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
/******************************************************************************
Discrete.cpp
Copyright (c) Roger Burghall 2014..2017
******************************************************************************/
/** \file
'Discrete.cpp' holds the code for implementing discrete, i.e. z-transform, filters.
*/
#include <cmath>
#include "math.h"
#include "Calcs.h"
#include <iostream>
#include <fstream>
#include <complex>
#include <stdlib.h>
#include "Discrete.h"
iir iir1;
/// 'iir::step(...)' implements the Direct Form 1 realisation.
double iir::step(double in)
{
int i;
for(i=2; i>0; i--) x[i] = x[i-1], y[i] = y[i-1];
x[0] = in * k;
y[0] = -a1*y[1] - a2*y[2] + b0*x[0] + b1*x[1] + b2*x[2];
return(y[0]);
}
/** Set the iir filter objects coefficients. */
void iir::setup(double kk, double aa0, double aa1, double aa2, double bb1, double bb2)
{
k = kk;
a1 = aa1, a2 = aa2;
b1 = bb1, b2 = bb2;
}
void iir::write(int stage)
{
char str[256];
std::ofstream out;
if(stage == 0) {
out.open("./out.c", std::ofstream::trunc);
/// Write initial code to 'out.c'
}
else out.open("./out.c", std::ofstream::app);
sprintf(str, "y%1d = -a%1d1*y%1d[1] - a%1d2*y%1d[2] + b%1d0*x%1d[0] + b%1d1*x%1d[1] + b%1d2*x%1d[2];\n",
stage, stage, stage, stage, stage, stage, stage, stage, stage, stage, stage);
out << str;
if(out.is_open( )) { out.close( ); }
}
/** Pre-warp the continuous filter. */
void prewarp(void)
{
/// Not yet designed.
}
/** Transform the continuous filter into a discrete one. */
void iir::transform(filter_class fclass, double Omega_ac, double q)
{
double ro2 = 1.0 / (Omega_ac * Omega_ac);
double roq = 1.0 / (Omega_ac * q);
b0 = 1.0;
switch(fclass) {
case Lowpass:
b1 = 2.0;
b2 = 1.0;
break;
case Highpass:
b0 = ro2;
b1 = -2.0 * ro2;
b2 = ro2;
break;
case Bandpass:
b0 = 1.0;
b1 = 0.0;
b2 = -1.0;
break;
default: cout << "\nUnknown filter type.\n"; return;
}
a0 = ro2 + roq + 1.0;
std::cout << "a0 =" << a0 << "\n";
// k = 1.0 / a0;
a1 = 2.0 - 2.0 * ro2;
a2 = ro2 - roq + 1.0;
a1 /= a0;
a2 /= a0;
b0 /= a0;
b1 /= a0;
b2 /= a0;
std::cout << "(" << b0 << "+(z^-1)*" << b1 << "+(z^-2)*" << b2 << ")/(" << 1.0 << "+(z^-1)*" << a1 << "+(z^-2)*" << a2 << ")\n";
std::cout << "-->" << (b0 + b1 + b2) / (1.0 + a1 + a2) << "\n";
std::cout << "[a0 = 1.0, a1 =" << a1 << ", a2 =" << a2 << "; b0 =" << b0 << ", b1 =" << b1 << ", b2 =" << b2 << "]\n";
/* std::cout << "Direct Form 2 realisation.\n--------------------------\n";
std::cout << "y[n] = b0 * w[n] + b1 * w[n-1] + b2 * w[n-2], w[n] = x[n] - a1 * w[n-1] - a2 * w[n-2]";
*/
}
/** 'normalise(...)' sets 'omega' */
double normalise(double fc, double fs)
{
double omega = 2.0 * M_PI * fc / fs;
return(omega);
}
#if false
int main( )
{
int i;
double fc, fs, omega;
std::cout << "cut-off freq: ";
std::cin >> fc;
std::cout << "sampling freq: ";
std::cin >> fs;
if(fc >= fs/2.0) {
std::cout << "That doesn't make sense! Cut-off must be less than sampling freq / 2.\n";
exit(-1);
}
omega = normalise(fc, fs);
std::cout << "omega =" << omega << "\n";
double omega_ac = tan(omega/2.0);
std::cout << "omega_ac =" << omega_ac << "\n\n";
iir1.transform(omega_ac, 0.7);
// iir1.setup(1.384, 0.1311136, 0.2162924, 0.1311136, -0.829328, 0.307046);
iir1.step(0.0);
std::cout << "\n";
std::cout << "Direct Form 1 realisation.\n--------------------------\n";
for(i=0; i<200; i++) std::cout << iir1.step(1.0) << "\n";
return(0);
}
#endif