UNIT 4:

DIGITAL SYSTEM MODELS

 

 

 

 

Summary

 

            This unit is concerned with the description of digital systems, it introduces the concepts of a linear time-invariant system, convolution, the general system difference equation and its implementation, and the frequency response characteristics of systems.

 

            When you have worked through this unit you should:

·      be able to define the characteristics of a linear, time-invariant system

·      understand what is represented by an impulse response and a frequency response of a system

·      be able to explain how and why systems may be simulated with discrete convolution

·      understand how a system can be described with the general difference equation

·      be able to describe one direct form method for the implementation of the general difference equation in an algorithm

·      be able to state the z-transform of a sampled sequence

·      appreciate that the z-transform may be used to obtain a polynomial statement of the general difference equation.

·      understand how the response of a system may be stated in terms of poles and zeros

·      understand that the response can be evaluated for any particular frequency to get the frequency response graph of a system

·      know how to use an algorithm for implementing an LTI system and calculating its response

·      be able to relate the parameters of a simple resonator to its implementation

 

Concepts

 

            We are concerned with systems that transform one digital signal into another, in particular with systems that perform a linear processing that does not change with time (is time-invariant).  We seek to describe the behaviour of such systems and to find a general algorithmic description suitable for digital implementation.

 

            Generally we describe systems that modify signals in terms of their Frequency Response (a graph of system response against frequency) and their Phase Response (a graph of phase shift against frequency).  A linear time-invariant (LTI) system may also be described in the time-domain by means of its Impulse response, that is: the graph of the output of the system after an impulse is presented on its input.  This is possible because we can describe any digital signal as a sequence of impulses of a height equal to the signal amplitude occurring at the sampling instants.  The output of a LTI system to any signal can therefore be predicted by a suitable linear combination of impulse responses starting at each sample instant (this is a consequence of the principle of superposition).

 

            Mathematically, this process can be described as taking a signal x[] and an impulse response h[] and calculating each output sample y[n] with the formula:

                                                   

1

            Assuming h[] is zero at negative times, and samples are counted from 1.  The operation is analogous to laying the impulse response backwards alongside the input signal prior to n and then cross multiplying and  adding.  This process is known as convolution. y[n] is said to be generated by the convolution of x[] and h[].

 

            Convolution with the impulse response is the time-domain equivalent of multiplying by the frequency response.  You can see this by considering convolution of h[] with a sinusoid signal given by:

            so that the output of the system is given by:

 

            which can be rewritten as:

 

            That is the output is also a sinusoid, at the same frequency as the input, but with an amplitude scaling and a phase change given by H(f).

 

            The impulse response of a simple resonator is simply an exponentially decaying sinewave of a frequency equal to the natural frequency and of a rate of decay related to the bandwidth.

 

            Systems that have some kind of 'memory' can have impulse responses of infinite duration (called Infinite-Impulse-Response or IIR systems), others have an impulse response of finite duration (called Finite-Impulse-Response or FIR systems).  Our simple resonator theoretically continues vibrating for ever, so is an IIR system.

 

            To describe an LTI system algorithmically, we need a formulation of its behaviour that allows both finite and infinite impulse responses.  We do this by noting that the characteristics of an LTI system must be expressed in a constant linear relationship between previous output samples and previous input samples.  That is, given input samples x[] and output samples y[], some linear combination of output samples is related to some linear combination of input samples:

                                  

2

            This is the general difference equation of an LTI system.  It defines the behaviour of an LTI system in terms of a set of coefficients a[0..p] that operate on current and previous input samples, and a set of coefficients b[1..q] that operate on previous output samples.  b[0] is always taken as 1.  This is perhaps easier to see with this reformulation:

 

 

            This can be readily implemented as a computer algorithm.

 

            The general difference equation leads to a simple classification of LTI systems:

 

            (i)         systems with the coefficients b[1..q] are zero, in which case output samples are generated by a linear combination of previous input samples.  These are called non-recursive or moving average systems.

 

            (ii)        systems in which a[0] is 1, and a[1..p] are zero, in which case the output is formed from a single input and a linear combination of past output samples.  These are called simply recursive or autoregressive systems.

 

            (iii)       systems with arbitrary a[] and b[] coefficients are sometimes called autoregressive/moving average or ARMA systems.  They are of course also recursive.

 

            The direct form implementation of the general difference equation requires two memory arrays: one for past inputs x[] and one for past outputs y[].  It can be shown that the same result can be achieved with only a single memory of an 'internal' waveform s[].  This direct form II is the basis for the LTISystem() implementation below.

 

            From the description of an LTI system in terms of a[] and b[] coefficients we can also determine its frequency response characteristics (note that we leave to Unit 7 the reverse problem of the determination of a[] and b[] from the frequency response).  This is normally performed using a mathematical transformation of the a[] and b[] coefficients akin to the discrete Fourier transform, called the Z transform.

 

            Put simply, the Z transform converts a sequence of digital samples at successive sampling instants to a polynomial of some operator z-1, in which the amplitudes of the samples become the coefficients of the polynomial.  Thus a finite signal x[1..n] has a Z transform:

 

           

            that is

            One way of thinking of z-1 is that it represents a unit delay. In other words multiplying by z-n delays a z-transformed signal by n samples.  This analogy allows to view a z-transformed signal as scaled impulses delayed appropriately and added together.

 

            The conversion of a sequence to a polynomial is useful because it allows us to write the general difference equation as two polynomials:

                                  

3

            Where X(z) is the z-transform of the input x[], and Y(z) the z-transform of the output y[].  The cross multiplication and collection of terms in z-n is equivalent to convolution on the original signals.  This polynomial form of the difference equation allows us to form the response of the system:

                                                

4

            that is, as a ratio of two polynomials in z.  We can then see that this response will have peaks at the roots[1] of the denominator polynomial - these are called the poles of the system, and dips at the roots of the numerator polynomial - these are called the zeros of the system.  Our LTI system can be described completely (with the addition of an overall gain factor) by the location of the poles and zeros of the system.  These will in general be complex numbers and can be plotted on the complex plane on a diagram called a z-plane diagram.

 

            To obtain the frequency response of our system we need to put into the system sinewaves at different frequencies but of unit amplitude and plot the output amplitude.  We do this by substituting z-1 in our polynomial for the complex sine eiW where W is the angular frequency (= 2πf/Fs) at which we require the response.  We can now obtain a numerical value for the response for a range of frequencies between 0 and the sampling rate (2p).  What we are actually calculating is the product of the distances between each system zero and a point on the unit circle specified by eiW divided by the product of the distances between each pole and that point.  Thus the nearer the poles are to the unit circle the higher the peaks in the response and the closer the zeros are to the unit circle the lower the dips.  The magnitude of this ratio at a given frequency is the magnitude of the frequency response, while the argument of this ratio is the phase response.

 

Algorithms

 

            Algorithm 4.1 Convolution

 

// conv.cpp -- implementation of Convolution of 2 waveforms

//

// C++ (c) 1996 Mark Huckvale University College London

 

#include "tools.h"

#include "conv.h"

 

// The function Convolve performs convolution of two

// time domain waveforms

 

Waveform Convolve(

       const Waveform& x,   // input waveform

       const Waveform& h    // impulse response

       )                    // returns convolved sequence

{

       Waveform y(x.count(),x.rate());

 

       // for each output sample

       for (int n=1;n<=y.count();n++) {

              y[n] = 0;

              // sum of product with reverse IR

              for (int k=0;k<h.count();k++)

                     y[n] += x[n-k] * h[k+1];

                     // exploits bad index capability of x[]

       }

 

       return y;

}

 

           
Algorithm 4.2 Linear time-invariant system model

 

// ltisys.h -- definition for Linear system class

//

// C++ (c) 1996 Mark Huckvale University College London

 

// The LTISystem class provides a convenient implementation

// of a linear time-invariant system including estimation

// of frequency response

 

#ifndef _LTISys_H

#define _LTISys_H

 

#include "wavedoub.h"

#include "complex.h"

 

class LTISystem {

public:       // public for demonstration purposes

       WaveDouble    a;     // numerator

       WaveDouble    b;     // denominator

       WaveDouble    s;     // state memory

public:

       // constructors

       LTISystem(int na,int nb);

       LTISystem(const LTISystem& ltis);

       ~LTISystem();

       // assignment

       LTISystem& operator= (const LTISystem& ltis);

       // clear system

       void clear();

       // run system

       double operator() (double ival);

       // frequency response

       Complex response(double freq) const;

};

 

#endif

 

 

// ltisys.cpp -- implementation for Linear system class

//

// C++ (c) 1996 Mark Huckvale University College London

 

// The LTISystem class provides a convenient implementation

// of a linear time-invariant system including estimation

// of frequency response

 

#include "tools.h"

#include "ltisys.h"

 

// constructors

LTISystem::LTISystem(int na,int nb)

: a(na+1,1,0),                    // a[0]..a[na]

  b(nb+1,1,0),                    // b[0]..b[nb] (but b[0]=1.0)

  s((na>nb)?na+1:nb+1,1,0)        // s[0]..s[max(na,nb)]

{

       clear();

}

LTISystem::LTISystem(const LTISystem& ltis)

{

       a = ltis.a;

       b = ltis.b;

       s = ltis.s;

}

LTISystem::~LTISystem()

{

}

 

// assignment

LTISystem& LTISystem::operator= (const LTISystem& ltis)

{

       if (this == &ltis)

              return *this;

       a = ltis.a;

       b = ltis.b;

       s = ltis.s;

       return *this;

}

 

// clear

void LTISystem::clear()

{

       for (int i=0;i<s.count();i++) s[i]=0;

}

 

// run sample through system

double LTISystem::operator() (double ival)

{

       int i;

 

       // shift state memory

       for (i=s.count()-1;i>0;i--) s[i] = s[i-1];

 

       // compute s[0] from y[] coefficients

       s[0] = ival;

       for (i=1;i<b.count();i++)

              s[0] -= b[i] *s[i];

 

       // compute output from x[] coefficients

       double out = 0;

       for (i=0;i<a.count();i++)

              out += a[i] * s[i];

 

       return out;

}

 

            Algorithm 4.3 LTI System response

 

// calculate frequency response

Complex LTISystem::response(double freq) const

{

       int i;

       WaveComplex omega(s.count(),1,0);

       Complex z=exp(Complex(0,2*PI*freq));

 

       // initialise polynomial of complex frequency

       omega[0] = 1.0;

       omega[1] = z;

       for (i=2;i<s.count();i++)

              omega[i] = omega[i-1] * z;

 

       // compute response of numerator

       Complex rnum(0,0);

       for (i=0;i<a.count();i++)

              rnum += a[i] * omega[i];

 

       // compute response of denominator

       Complex rden=omega[0];

       for (i=1;i<b.count();i++)

              rden += b[i] * omega[i];

 

       // compute ratio

       if (mag(rden)==0)

              return 1.0E5;        // i.e. infinity

       else

              return rnum/rden;

}

 

            Algorithm 4.4 Digital resonator

 

A mathematical justification for the formula used to build a digital resonator is also available.

 

// reson.cpp -- implementation of simple digital resonator

//

// C++ (c) 1996 Mark Huckvale University College London

 

#include "tools.h"

#include "ltisys.h"

#include "reson.h"

 

// The Resonator function sets up a linear system consisting

// of a complex pole pair

 

LTISystem Resonator(

       double freq,  // natural frequency

                     // (fraction of sample rate)

       double bwid   // bandwidth

                     // (fraction of sampling rate)

       )             // returns linear system

{

       LTISystem ltis(0,2);

 

       // get parameters as angles

       double wf = 2*PI*freq;

       double wb = 2*PI*bwid;

 

       // estimate pole radius from bandwidth

       //   (rule of thumb)

       double r = 1.0 - wb/2;

 

       // set up numerator

       ltis.a[0] = 1.0;

 

       // set up denominator: quadratic formula

       ltis.b[1] = -2.0 * r * cos(wf);

       ltis.b[2] = r * r;

 

       // adjust numerator for unity gain at DC

       ltis.a[0] /= mag(ltis.response(0));

 

       return ltis;

}

 

 

Bibliography

 

            Meade & Dillon, Signals and Systems, Chapter 5.

 

            Lynn & Fuerst, Introductory Digital Signal Processing, Chapters 2 & 4.

 

            Orfanidis, Introduction to Signal Processing, Sections 5.1, 5.4, 6.1-6.3.

Example Programs

 

Example 4.1 Demonstration of convolution

 

// conv_t.cpp -- demonstration of convolution

 

#include "tools.h"

#include "conv.h"

 

int main()

{

       int i;

 

       Graph gr(2,3,"Demonstration of Convolution");

 

       // create input signal and display

       Waveform x(10,1);

       for (i=1;i<=10;i++) x[i]=0;

       x[3]=1;

       x.plot(gr,1,"Input Waveform 1");

 

       // create impulse response

       Waveform h(4,1);

       h[1]=1; h[2]=2; h[3]=-2; h[4]=-1;

       h.plot(gr,2,"Impulse Response");

 

       // do convolution

       Waveform y = Convolve(x,h);

       y.plot(gr,3,"Convolved Output 1");

 

       // try a different input

       for (i=1;i<=10;i++) x[i]=i%5;

       x.plot(gr,4,"Input Waveform 2");

 

       // same impulse response

       h.plot(gr,5,"Impulse Response");

 

       // do convolution

       y = Convolve(x,h);

       y.plot(gr,6,"Convolved Output 2");

      

       gr.close();

}

 

 

 

 

            Example 4.2 Equivalence of filtering and convolution

 

// reson_t.cpp -- demonstration of resonator design & convolution

 

#include "tools.h"

#include "conv.h"

#include "ltisys.h"

#include "reson.h"

 

const double ILIMIT=1.0E-5;       // length limit for impulse response

 

int main()

{

       int i;

       Graph gr(3,2,"Simple Resonator System");

 

       // create resonator

       LTISystem ltis=Resonator(0.1,0.01);

 

       // calculate frequency response

       Spectrum fresp(500,1000);

       for (i=0;i<500;i++)

              fresp[i] = ltis.response(i/1000.0);

       fresp.plotLogMag(gr,1,"Magnitude Response");

       fresp.plotPhase(gr,3,"Phase Response");

 

       // calculate impulse response

       Waveform iresp(0,1);

       double lval = 0;                  // last output

       double oval = ltis(1.0);          // put in unit pulse

       while ((fabs(oval) > ILIMIT) || (fabs(oval-lval) > ILIMIT)) {

              iresp += oval;             // append sample

              lval = oval;               // remember sample

              oval = ltis(0.0);          // get next sample

       }

       iresp.plot(gr,5,"Impulse Response");

 

       // create a ramp waveform

       Waveform ramp(1000,1);

       for (i=1;i<=1000;i++)

              ramp[i] = i % 100;

       ramp.plot(gr,2,"Ramp Waveform");

 

       // filter a ramp waveform

       Waveform framp(1000,1);

       ltis.clear();

       for (i=1;i<=1000;i++)

              framp[i] = ltis(ramp[i]);

       framp.plot(gr,4,"Filtered Ramp");

 

       // convolve ramp with impulse response

       Waveform cramp=Convolve(ramp,iresp);

       cramp.plot(gr,6,"Convolved Ramp");

 

       gr.close();

}

 

 

 

Exercises

 

4.1       Explain the first three values of 'Convolved Output 2' (they are: 1, 4 & 5).

 

4.2       Modify example program 4.1 so that the impulse response is a square wave.

 

4.3       Modify example program 4.1 so that the impulse response matches the input waveform in each case.

 

4.4       Adapt example program 4.2 to send a stream of pulses through the resonator rather than the ramp wave, displaying the input and output signals.

 

4.5       Adapt example program 4.2 to generate the output from a resonator with a natural frequency of 25Hz and a bandwidth of 7.5Hz when fed with a pulse train at 10Hz lasting 1 second and sampled at 1,000Hz.  Display the result.  (Hint: create a Signal of 1000 samples at 1000 samples/sec, create a Resonator of centre frequency (25/1000) and bandwidth (7.5/1000), then run the resonator into the signal, putting in pulses every 100 samples)

 


    [1]The roots of a polynomial in x are those values of x for which the polynomial is zero.  If the polynomial is factored into a form (x-a)(x-b)(x-c)..(x-p), then the roots are the values a,b,c,..,p.