Realizing a discrete dynamic system

The biquad cascade

Although we could implement the difference eq. 1 as shown, the sensitivity of the output to the coefficients leads to numerical inaccuracies as the order of the system \(N\) becomes large, so we will represent the \(N\)th-order system as a sequence of second-order systems as described in section 6.4.

Again, this technique is called a biquad cascade, and it is illustrated in figure 6.18.

Notice that the output of each second-order section (biquad) is the input to the subsequent section. Each biquad implements the same second-order difference equation, but with different coefficients, inputs, outputs, and internal states. For example, the current output \(y_i(n)\) from the \(i\)th section would be \[ y_i(n) = \left[b_{0_i}x_i(n) +b_{1_i}x_i(n-1)+b_{2_i}x_i(n-2) -a_{1_i}y_i(n-1)-a_{2_i}y_i(n-2)\right]/a_{0_i}. \qquad{(3)}\]

Of course, a first- or second-order transfer function would require only one biquad. Depending on the value of \(N\), some of the coefficients of at least one biquad may be zero. We will implement a function to handle any value of \(N\).

The MATLAB Control System Toolbox contains convenient functions for the conversion of continuous models to discrete models:

1. The function c2d(), “continuous to discrete,” computes an equivalent discrete system from a continuous system. For example,

   SYSD = c2d(SYS, T, 'tustin'); % SYS is a continuous system model

2. The function tf2sos(), “transfer function to second-order sections,” breaks a transfer function into second-order biquad sections. For example,

   Cd = tf(SYSD);                % SYSD is a discrete system model
   [b, a] = tfdata(Cd, 'v');     % tf coefficients
   sos = tf2sos(b, a);           % biquad sections

Main function

The only tasks of main() will be the following:

1. Set up and enable the Timer IRQ interrupt.
2. Enter a loop, ending only when a keyboard_backspace is received from the keypad. Use getkey(), that we wrote in lab 3.
3. Signal the timer thread to terminate using the irqThreadRdy flag and wait for the thread to terminate.

The cascade() function

The cascade() function implements the complete dynamic system by passing the measured input through the string of biquad sections. The ISR must pass to cascade() the value of the current input \(x(n)\) measured by the ADC, the array of biquad structures containing the coefficients, the number of biquad sections \(n_s\), and history variables (\(x\)s and \(y\)s) for all sections. Its prototype should be

double cascade(double xin,         // input
               struct biquad *fa,  // biquad array
               int    ns,          // no. segments
               double ymin,        // min output  
               double ymax);       // max output

Here, xin is the current system input, fa is the name of an array of biquad structures, ns is the corresponding number of biquad sections, and ymin and ymax are the saturation limits.

In this example, myFilter would be passed through fa. The value returned by cascade() is the current value of the system output \(y(n)\).

Programming cascade()

An efficient way to code cascade() is to allocate a pointer f in cascade() that will be used to point to elements of the array of biquad structures. Begin the function by setting the pointer to the first element in the array (i.e., the first biquad): f = fa. Variables inside the biquad structure are accessed by using the pointer name, (e.g., f->a0, f->b0, f->x0, f->y1, etc.). (The -> operator is equivalent to dereferencing and then accessing a member—say, (*f).a0—and is typed as a minus sign immediately followed by >.)

Then, loop \(n_s\) times, to cycle through each of the biquad sections in the array. At the beginning of each loop, the output value y0 of the previous biquad must be passed to the input value f->x0 of the current biquad.

Within the loop, coding the output value y0 might look like

y0 = (f->b0 * f->x0 + f->b1 * f->x1 + /* etc */ ) / f->a0;

See eq. 3.

Each time through the loop, after the output value has been computed, the previous values of \(x\) and \(y\) must be updated so they will be correct at the next time step. For example,

f->x2 = f->x1; f->x1 = f->x0; // etc.

At the end of the loop, the pointer f is incremented to advance to the next biquad in the array.

One more point: If the DAC is given a value beyond its range \([-10, +10]\) V, it will saturate its output value appropriately. However, our difference eq. 3 depends on previous values of the output, but it doesn’t saturate. To correct this disparity, cascade() should saturate the output y0 of the final biquad before it is saved for the next iteration.

For example, define the macro as follows:

#define SATURATE(x,lo,hi) ((x) < (lo) ? (lo) : (x) > (hi) ? (hi) : (x))

Pass appropriate values of the xmin and xmax parameters to cascade(). Then, for the last biquad, immediately after y0 is computed, saturate its value as follows

y0 = SATURATE(y0, ymin, ymax);