1 /*******************************************************************************
2
3 "A Collection of Useful C++ Classes for Digital Signal Processing"
4 By Vinnie Falco
5
6 Official project location:
7 https://github.com/vinniefalco/DSPFilters
8
9 See Documentation.cpp for contact information, notes, and bibliography.
10
11 --------------------------------------------------------------------------------
12
13 License: MIT License (http://www.opensource.org/licenses/mit-license.php)
14 Copyright (c) 2009 by Vinnie Falco
15
16 Permission is hereby granted, free of charge, to any person obtaining a copy
17 of this software and associated documentation files (the "Software"), to deal
18 in the Software without restriction, including without limitation the rights
19 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
20 copies of the Software, and to permit persons to whom the Software is
21 furnished to do so, subject to the following conditions:
22
23 The above copyright notice and this permission notice shall be included in
24 all copies or substantial portions of the Software.
25
26 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
27 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
28 FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
29 AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
30 LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
31 OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
32 THE SOFTWARE.
33
34 *******************************************************************************/
35
36 #include "DspFilters/Common.h"
37 #include "DspFilters/ChebyshevII.h"
38
39 namespace Dsp {
40
41 namespace ChebyshevII {
42
43 // "Chebyshev Filter Properties"
44 // http://cnx.org/content/m16906/latest/
45
AnalogLowPass()46 AnalogLowPass::AnalogLowPass ()
47 : m_numPoles (-1)
48 {
49 setNormal (0, 1);
50 }
51
design(int numPoles,double stopBandDb)52 void AnalogLowPass::design (int numPoles,
53 double stopBandDb)
54 {
55 if (m_numPoles != numPoles ||
56 m_stopBandDb != stopBandDb)
57 {
58 m_numPoles = numPoles;
59 m_stopBandDb = stopBandDb;
60
61 reset ();
62
63 const double eps = std::sqrt (1. / (std::exp (stopBandDb * 0.1 * doubleLn10) - 1));
64 const double v0 = asinh (1 / eps) / numPoles;
65 const double sinh_v0 = -sinh (v0);
66 const double cosh_v0 = cosh (v0);
67 const double fn = doublePi / (2 * numPoles);
68
69 int k = 1;
70 for (int i = numPoles / 2; --i >= 0; k+=2)
71 {
72 const double a = sinh_v0 * cos ((k - numPoles) * fn);
73 const double b = cosh_v0 * sin ((k - numPoles) * fn);
74 const double d2 = a * a + b * b;
75 const double im = 1 / cos (k * fn);
76 addPoleZeroConjugatePairs (complex_t (a / d2, b / d2),
77 complex_t (0, im));
78 }
79
80 if (numPoles & 1)
81 {
82 add (1 / sinh_v0, infinity());
83 }
84 }
85 }
86
87 //------------------------------------------------------------------------------
88
89 //
90 // Chebyshev Type I low pass shelf prototype
91 // From "High-Order Digital Parametric Equalizer Design"
92 // Sophocles J. Orfanidis
93 // http://www.ece.rutgers.edu/~orfanidi/ece521/hpeq.pdf
94 //
95
AnalogLowShelf()96 AnalogLowShelf::AnalogLowShelf ()
97 : m_numPoles (-1)
98 {
99 setNormal (doublePi, 1);
100 }
101
design(int numPoles,double gainDb,double stopBandDb)102 void AnalogLowShelf::design (int numPoles,
103 double gainDb,
104 double stopBandDb)
105 {
106 if (m_numPoles != numPoles ||
107 m_stopBandDb != stopBandDb ||
108 m_gainDb != gainDb)
109 {
110 m_numPoles = numPoles;
111 m_stopBandDb = stopBandDb;
112 m_gainDb = gainDb;
113
114 reset ();
115
116 gainDb = -gainDb;
117
118 if (stopBandDb >= fabs(gainDb))
119 stopBandDb = fabs (gainDb);
120 if (gainDb<0)
121 stopBandDb = -stopBandDb;
122
123 const double G = std::pow (10., gainDb / 20.0 );
124 const double Gb = std::pow (10., (gainDb - stopBandDb) / 20.0);
125 const double G0 = 1;
126 const double g0 = pow (G0 , 1. / numPoles);
127
128 double eps;
129 if (Gb != G0 )
130 eps = sqrt((G*G-Gb*Gb)/(Gb*Gb-G0*G0));
131 else
132 eps = G-1; // This is surely wrong
133
134 const double b = pow (G/eps+Gb*sqrt(1+1/(eps*eps)), 1./numPoles);
135 const double u = log (b / g0);
136 const double v = log (pow (1. / eps+sqrt(1+1/(eps*eps)),1./numPoles));
137
138 const double sinh_u = sinh (u);
139 const double sinh_v = sinh (v);
140 const double cosh_u = cosh (u);
141 const double cosh_v = cosh (v);
142 const double n2 = 2 * numPoles;
143 const int pairs = numPoles / 2;
144 for (int i = 1; i <= pairs; ++i)
145 {
146 const double a = doublePi * (2 * i - 1) / n2;
147 const double sn = sin (a);
148 const double cs = cos (a);
149 addPoleZeroConjugatePairs (complex_t (-sn * sinh_u, cs * cosh_u),
150 complex_t (-sn * sinh_v, cs * cosh_v));
151 }
152
153 if (numPoles & 1)
154 add (-sinh_u, -sinh_v);
155 }
156 }
157
158 //------------------------------------------------------------------------------
159
setup(int order,double sampleRate,double cutoffFrequency,double stopBandDb)160 void LowPassBase::setup (int order,
161 double sampleRate,
162 double cutoffFrequency,
163 double stopBandDb)
164 {
165 m_analogProto.design (order, stopBandDb);
166
167 LowPassTransform (cutoffFrequency / sampleRate,
168 m_digitalProto,
169 m_analogProto);
170
171 Cascade::setLayout (m_digitalProto);
172 }
173
setup(int order,double sampleRate,double cutoffFrequency,double stopBandDb)174 void HighPassBase::setup (int order,
175 double sampleRate,
176 double cutoffFrequency,
177 double stopBandDb)
178 {
179 m_analogProto.design (order, stopBandDb);
180
181 HighPassTransform (cutoffFrequency / sampleRate,
182 m_digitalProto,
183 m_analogProto);
184
185 Cascade::setLayout (m_digitalProto);
186 }
187
setup(int order,double sampleRate,double centerFrequency,double widthFrequency,double stopBandDb)188 void BandPassBase::setup (int order,
189 double sampleRate,
190 double centerFrequency,
191 double widthFrequency,
192 double stopBandDb)
193 {
194 m_analogProto.design (order, stopBandDb);
195
196 BandPassTransform (centerFrequency / sampleRate,
197 widthFrequency / sampleRate,
198 m_digitalProto,
199 m_analogProto);
200
201 Cascade::setLayout (m_digitalProto);
202 }
203
setup(int order,double sampleRate,double centerFrequency,double widthFrequency,double stopBandDb)204 void BandStopBase::setup (int order,
205 double sampleRate,
206 double centerFrequency,
207 double widthFrequency,
208 double stopBandDb)
209 {
210 m_analogProto.design (order, stopBandDb);
211
212 BandStopTransform (centerFrequency / sampleRate,
213 widthFrequency / sampleRate,
214 m_digitalProto,
215 m_analogProto);
216
217 Cascade::setLayout (m_digitalProto);
218 }
219
setup(int order,double sampleRate,double cutoffFrequency,double gainDb,double stopBandDb)220 void LowShelfBase::setup (int order,
221 double sampleRate,
222 double cutoffFrequency,
223 double gainDb,
224 double stopBandDb)
225 {
226 m_analogProto.design (order, gainDb, stopBandDb);
227
228 LowPassTransform (cutoffFrequency / sampleRate,
229 m_digitalProto,
230 m_analogProto);
231
232 Cascade::setLayout (m_digitalProto);
233 }
234
setup(int order,double sampleRate,double cutoffFrequency,double gainDb,double stopBandDb)235 void HighShelfBase::setup (int order,
236 double sampleRate,
237 double cutoffFrequency,
238 double gainDb,
239 double stopBandDb)
240 {
241 m_analogProto.design (order, gainDb, stopBandDb);
242
243 HighPassTransform (cutoffFrequency / sampleRate,
244 m_digitalProto,
245 m_analogProto);
246
247 Cascade::setLayout (m_digitalProto);
248 }
249
setup(int order,double sampleRate,double centerFrequency,double widthFrequency,double gainDb,double stopBandDb)250 void BandShelfBase::setup (int order,
251 double sampleRate,
252 double centerFrequency,
253 double widthFrequency,
254 double gainDb,
255 double stopBandDb)
256 {
257 m_analogProto.design (order, gainDb, stopBandDb);
258
259 BandPassTransform (centerFrequency / sampleRate,
260 widthFrequency / sampleRate,
261 m_digitalProto,
262 m_analogProto);
263
264 m_digitalProto.setNormal (((centerFrequency/sampleRate) < 0.25) ? doublePi : 0, 1);
265
266 Cascade::setLayout (m_digitalProto);
267 }
268
269 }
270
271 }
272