Open Ephys GUI
 All Classes Functions Variables
State.h
1 /*******************************************************************************
2 
3 "A Collection of Useful C++ Classes for Digital Signal Processing"
4  By Vincent Falco
5 
6 Official project location:
7 http://code.google.com/p/dspfilterscpp/
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 Vincent 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 #ifndef DSPFILTERS_STATE_H
37 #define DSPFILTERS_STATE_H
38 
39 #include "Common.h"
40 #include "Biquad.h"
41 
42 #include <stdexcept>
43 
44 namespace Dsp {
45 
46 /*
47  * Various forms of state information required to
48  * process channels of actual sample data.
49  *
50  */
51 
52 //------------------------------------------------------------------------------
53 
54 /*
55  * State for applying a second order section to a sample using Direct Form I
56  *
57  * Difference equation:
58  *
59  * y[n] = (b0/a0)*x[n] + (b1/a0)*x[n-1] + (b2/a0)*x[n-2]
60  * - (a1/a0)*y[n-1] - (a2/a0)*y[n-2]
61  */
63 {
64 public:
65  DirectFormI ()
66  {
67  reset();
68  }
69 
70  void reset ()
71  {
72  m_x1 = 0;
73  m_x2 = 0;
74  m_y1 = 0;
75  m_y2 = 0;
76  }
77 
78  template <typename Sample>
79  inline Sample process1 (const Sample in,
80  const BiquadBase& s,
81  const double vsa) // very small amount
82  {
83  double out = s.m_b0*in + s.m_b1*m_x1 + s.m_b2*m_x2
84  - s.m_a1*m_y1 - s.m_a2*m_y2
85  + vsa;
86  m_x2 = m_x1;
87  m_y2 = m_y1;
88  m_x1 = in;
89  m_y1 = out;
90 
91  return static_cast<Sample> (out);
92  }
93 
94 protected:
95  double m_x2; // x[n-2]
96  double m_y2; // y[n-2]
97  double m_x1; // x[n-1]
98  double m_y1; // y[n-1]
99 };
100 
101 //------------------------------------------------------------------------------
102 
103 /*
104  * State for applying a second order section to a sample using Direct Form II
105  *
106  * Difference equation:
107  *
108  * v[n] = x[n] - (a1/a0)*v[n-1] - (a2/a0)*v[n-2]
109  * y(n) = (b0/a0)*v[n] + (b1/a0)*v[n-1] + (b2/a0)*v[n-2]
110  *
111  */
113 {
114 public:
115  DirectFormII ()
116  {
117  reset ();
118  }
119 
120  void reset ()
121  {
122  m_v1 = 0;
123  m_v2 = 0;
124  }
125 
126  template <typename Sample>
127  Sample process1 (const Sample in,
128  const BiquadBase& s,
129  const double vsa)
130  {
131  double w = in - s.m_a1*m_v1 - s.m_a2*m_v2 + vsa;
132  double out = s.m_b0*w + s.m_b1*m_v1 + s.m_b2*m_v2;
133 
134  m_v2 = m_v1;
135  m_v1 = w;
136 
137  return static_cast<Sample> (out);
138  }
139 
140 private:
141  double m_v1; // v[-1]
142  double m_v2; // v[-2]
143 };
144 
145 //------------------------------------------------------------------------------
146 
147 /*
148  * Transposed Direct Form I and II
149  * by lubomir i. ivanov (neolit123 [at] gmail)
150  *
151  * Reference:
152  * http://www.kvraudio.com/forum/viewtopic.php?p=4430351
153  *
154  */
155 
156 // I think this one is broken
158 {
159 public:
161  {
162  reset ();
163  }
164 
165  void reset ()
166  {
167  m_v = 0;
168  m_s1 = 0;
169  m_s1_1 = 0;
170  m_s2 = 0;
171  m_s2_1 = 0;
172  m_s3 = 0;
173  m_s3_1 = 0;
174  m_s4 = 0;
175  m_s4_1 = 0;
176  }
177 
178  template <typename Sample>
179  inline Sample process1 (const Sample in,
180  const BiquadBase& s,
181  const double vsa)
182  {
183  double out;
184 
185  // can be: in += m_s1_1;
186  m_v = in + m_s1_1;
187  out = s.m_b0*m_v + m_s3_1;
188  m_s1 = m_s2_1 - s.m_a1*m_v;
189  m_s2 = -s.m_a2*m_v;
190  m_s3 = s.m_b1*m_v + m_s4_1;
191  m_s4 = s.m_b2*m_v;
192 
193  m_s4_1 = m_s4;
194  m_s3_1 = m_s3;
195  m_s2_1 = m_s2;
196  m_s1_1 = m_s1;
197 
198  return static_cast<Sample> (out);
199  }
200 
201 private:
202  double m_v;
203  double m_s1;
204  double m_s1_1;
205  double m_s2;
206  double m_s2_1;
207  double m_s3;
208  double m_s3_1;
209  double m_s4;
210  double m_s4_1;
211 };
212 
213 //------------------------------------------------------------------------------
214 
216 {
217 public:
219  {
220  reset ();
221  }
222 
223  void reset ()
224  {
225  m_s1 = 0;
226  m_s1_1 = 0;
227  m_s2 = 0;
228  m_s2_1 = 0;
229  }
230 
231  template <typename Sample>
232  inline Sample process1 (const Sample in,
233  const BiquadBase& s,
234  const double vsa)
235  {
236  double out;
237 
238  out = m_s1_1 + s.m_b0*in + vsa;
239  m_s1 = m_s2_1 + s.m_b1*in - s.m_a1*out;
240  m_s2 = s.m_b2*in - s.m_a2*out;
241  m_s1_1 = m_s1;
242  m_s2_1 = m_s2;
243 
244  return static_cast<Sample> (out);
245  }
246 
247 private:
248  double m_s1;
249  double m_s1_1;
250  double m_s2;
251  double m_s2_1;
252 };
253 
254 //------------------------------------------------------------------------------
255 
256 // Holds an array of states suitable for multi-channel processing
257 template <int Channels, class StateType>
259 {
260 public:
261  ChannelsState ()
262  {
263  }
264 
265  const int getNumChannels() const
266  {
267  return Channels;
268  }
269 
270  void reset ()
271  {
272  for (int i = 0; i < Channels; ++i)
273  m_state[i].reset();
274  }
275 
276  StateType& operator[] (int index)
277  {
278  assert (index >= 0 && index < Channels);
279  return m_state[index];
280  }
281 
282  template <class Filter, typename Sample>
283  void process (int numSamples,
284  Sample* const* arrayOfChannels,
285  Filter& filter)
286  {
287  for (int i = 0; i < Channels; ++i)
288  filter.process (numSamples, arrayOfChannels[i], m_state[i]);
289  }
290 
291 private:
292  StateType m_state[Channels];
293 };
294 
295 // Empty state, can't process anything
296 template <class StateType>
297 class ChannelsState <0, StateType>
298 {
299 public:
300  const int getNumChannels() const
301  {
302  return 0;
303  }
304 
305  void reset ()
306  {
307  throw std::logic_error ("attempt to reset empty ChannelState");
308  }
309 
310  template <class FilterDesign, typename Sample>
311  void process (int numSamples,
312  Sample* const* arrayOfChannels,
313  FilterDesign& filter)
314  {
315  throw std::logic_error ("attempt to process empty ChannelState");
316  }
317 };
318 
319 //------------------------------------------------------------------------------
320 
321 }
322 
323 #endif