The impulse invariance method of IIR filter design is based upon the notion that we can design a discrete filter whose timedomain impulse response is a sampled version of the impulse response of a continuous analog filter. If that analog filter (often called the prototype filter) has some desired frequency response, then our IIR filter will yield a discrete approximation of that desired response. The impulse response equivalence of this design method is depicted in Figure 623, where we use the conventional notation of d to represent an impulse function and hc(t) is the analog filter's impulse response. We use the subscript "c" in Figure 623(a) to emphasize the continuous nature of the analog filter. Figure 623(b) illustrates the definition of the discrete filter's impulse response: the filter's timedomain output sequence when the input is a single unityvalued sample (impulse) preceded and followed by all zerovalued samples. Our goal is to design a digital filter whose impulse response is a sampled version of the analog filter's continuous impulse response. Implied in the correspondence of the continuous and discrete impulse responses is the property that we can map each pole on the splane for the analog filter's Hc(s) transfer function to a pole on the zplane for the discrete IIR filter's H(z) transfer function. What designers have found is that the impulse invariance method does yield useful IIR filters, as long as the sampling rate is high relative to the bandwidth of the signal to be filtered. In other words, IIR filters designed using the impulse invariance method are susceptible to aliasing problems because practical analog filters cannot be perfectly bandlimited. Aliasing will occur in an IIR filter's frequency response as shown in Figure 624.
Figure 623. Impulse invariance design equivalence of (a) analog filter continuous impulse response; (b) digital filter discrete impulse response.
Figure 624. Aliasing in the impulse invariance design method: (a) prototype analog filter magnitude response; (b) replicated magnitude responses where HIIR(w) is the discrete Fourier transform of h(n) = hc(nts); (c) potential resultant IIR filter magnitude response with aliasing effects.
From what we've learned in Chapter 2 about the spectral replicating effects of sampling, if Figure 64(a) is the spectrum of the continuous hc(t) impulse response, then the spectrum of the discrete hc(nts) sample sequence is the replicated spectra in Figure 624(b).
In Figure 624(c) we show the possible effect of aliasing where the dashed curve is a desired HIIR(w) frequency magnitude response. However, the actual frequency magnitude response, indicated by the solid curve, can occur when we use the impulse invariance design method. For this reason, we prefer to make the sample frequency fs as large as possible to minimize the overlap between the primary frequency response curve and its replicated images spaced at multiples of ±fs Hz. To see how aliasing can affect IIR filters designed with this method, let's list the necessary impulse invariance design steps and then go through a filter design example.
There are two different methods for designing IIR filters using impulse invariance. The first method, which we'll call Method 1, requires that an inverse Laplace transform as well as a ztransform be performed[12,13]. The second impulse invariance design technique, Method 2, uses a direct substitution process to avoid the inverse Laplace and ztransformations at the expense of needing partial fraction expansion algebra necessary to handle polynomials[14–17]. Both of these methods seem complicated when described in words, but they're really not as difficult as they sound. Let's compare the two methods by listing the steps required for each of them. The impulse invariance design Method 1 goes like this:
Method 1, Step 1: 
Design (or have someone design for you) a prototype analog filter with the desired frequency response.[] The result of this step is a continuous Laplace transfer function Hc(s) expressed as the ratio of two polynomials, such as 


Which is the general form of Eq. (610) with N < M, and a(k) and b(k) are constants. 

Method 1, Step 2: 
Determine the analog filter's continuous timedomain impulse response hc(t) from the Hc(s) Laplace transfer function. I hope this can be done using Laplace tables as opposed to actually evaluating an inverse Laplace transform equation. 
Method 1, Step 3: 
Determine the digital filter's sampling frequency fs, and calculate the sample period as ts = 1/fs. The fs sampling rate is chosen based on the absolute frequency, in Hz, of the prototype analog filter. Because of the aliasing problems associated with this impulse invariance design method, later, we'll see why fs should be made as large as practical. 
Method 1, Step 4: 
Find the ztransform of the continuous hc(t) to obtain the IIR filter's zdomain transfer function H(z) in the form of a ratio of polynomials in z. 
Method 1, Step 5: 
Substitute the value (not the variable) ts for the continuous variable t in the H(z) transfer function obtained in Step 4. In performing this step, we are ensuring, as in Figure 623, that the IIR filter's discrete h(n) impulse response is a sampled version of the continuous filter's hc(t) impulse response so that h(n) = hc(nts), for 0 n . 
Method 1, Step 6: 
Our H(z) from Step 5 will now be of the general form 
Equation 644 

Because the process of sampling the continuous impulse response results in a digital filter frequency response that's scaled by a factor of 1/ts, many filter designers find it appropriate to include the ts factor in Eq. (644). So we can rewrite Eq. (644) as 

Equation 645 

Incorporating the value of ts in Eq. (645), then, makes the IIR filter timeresponse scaling independent of the sampling rate, and the discrete filter will have the same gain as the prototype analog filter.[] 

Method 1, Step 7: 
Because Eq. (644) is in the form of Eq. (625), by inspection, we can express the filter's timedomain difference equation in the general form of Eq. (621) as 
Equation 646 

Choosing to incorporate ts, as in Eq. (645), to make the digital filter's gain equal to the prototype analog filter's gain by multiplying the b(k) coefficients by the sample period ts leads to an IIR filter timedomain expression of the form 

Equation 647 

Notice how the signs changed for the a(k) coefficients from Eqs. (644) and (645) to Eqs. (646) and (647). These sign changes always seem to cause problems for beginners, so watch out. Also, keep in mind that the timedomain expressions in Eq. (646) and Eq. (647) apply only to the filter structure in Figure 618. The a(k) and b(k), or ts · b(k), coefficients, however, can be applied to the improved IIR structure shown in Figure 622 to complete our design. 
[] In a lowpass filter design, for example, the filter type (Chebyshev, Butterworth, elliptic), filter order (number of poles), and the cutoff frequency are parameters to be defined in this step.
[] Some authors have chosen to include the ts factor in the discrete h(n) impulse response in the above Step 4, that is, make h(n) = tshc(nts) [14, 18]. The final result of this, of course, is the same as that obtained by including ts as described in Step 5.
Before we go through an actual example of this design process, let's discuss the other impulse invariance design method.
The impulse invariance Design Method 2, also called the standard ztransform method, takes a different approach. It mathematically partitions the prototype analog filter into multiple singlepole continuous filters and then approximates each one of those by a singlepole digital filter. The set of M singlepole digital filters is then algebraically combined to form an Mpole, Mthordered IIR filter. This process of breaking the analog filter to discrete filter approximation into manageable pieces is shown in Figure 625. The steps necessary to perform an impulse invariance Method 2 design are:
Figure 625. Mathematical flow of the impulse invariance design Method 2.
Method 2, Step 1: 
Obtain the Laplace transfer function Hc(s) for the prototype analog filter in the form of Eq. (643). (Same as Method 1, Step 1.) 
Method 2, Step 2: 
Select an appropriate sampling frequency fs and calculate the sample period as ts = 1/fs. (Same as Method 1, Step 3.) 
Method 2, Step 3: 
Express the analog filter's Laplace transfer function Hc(s) as the sum of singlepole filters. This requires us to use partial fraction expansion methods to express the ratio of polynomials in Eq. (643) in the form of 
Equation 648 

where the individual Ak factors are constants and the kth pole is located at –pk on the splane. We'll denote the kth singlepole analog filter as Hk(s), or 

Equation 649 

Method 2, Step 4: 
Substitute for s + pk in Eq. (648). This mapping of each Hk(s) pole, located at s = –pk on the splane, to the location on the zplane is how we approximate the impulse response of each singlepole analog filter by a singlepole digital filter. (The reader can find the derivation of this substitution, illustrated in our Figure 625, in references [14] through [16].) So, the kth analog singlepole filter Hk(s) is approximated by a singlepole digital filter whose zdomain transfer function is 
Equation 650 

The final combined discrete filter transfer function H(z) is the sum of the singlepoled discrete filters, or 

Equation 651 

Keep in mind that the above H(z) is not a function of time. The ts factor in Eq. (651) is a constant equal to the discretetime sample period. 

Method 2, Step 5: 
Calculate the zdomain transfer function of the sum of the M singlepole digital filters in the form of a ratio of two polynomials in z. Because the H(z) in Eq. (651) will be a series of fractions, we'll have to combine those fractions over a common denominator to get a single ratio of polynomials in the familiar form of 
Equation 652 

Method 2, Step 6: 
Just as in Method 1 Step 6, by inspection, we can express the filter's timedomain equation in the general form of 
Equation 653 

Again, notice the a(k) coefficient sign changes from Eq. (652) to Eq. (653). As described in Method 1 Steps 6 and 7, if we choose to make the digital filter's gain equal to the prototype analog filter's gain by multiplying the b(k) coefficients by the sample period ts, then the IIR filter's timedomain expression will be in the form 

Equation 654 

yielding a final H(z) zdomain transfer function of 

Equation 654' 

Finally, we can implement the improved IIR structure shown in Figure 622 using the a(k) and b(k) coefficients from Eq. (653) or the a(k) and ts.b(k) coefficients from Eq. (654). 
To provide a more meaningful comparison between the two impulse invariance design methods, let's dive in and go through an IIR filter design example using both methods.
6.4.1 Impulse Invariance Design Method 1 Example
Assume that we need to design an IIR filter that approximates a secondorder Chebyshev prototype analog lowpass filter whose passband ripple is 1 dB. Our fs sampling rate is 100 Hz (ts = 0.01), and the filter's 1 dB cutoff frequency is 20 Hz. Our prototype analog filter will have a frequency magnitude response like that shown in Figure 626.
Figure 626. Frequency magnitude response of the example prototype analog filter.
Given the above filter requirements, assume that the analog prototype filter design effort results in the Hc(s) Laplace transfer function of
Equation 655
It's the transfer function in Eq. (655) that we intend to approximate with our discrete IIR filter. To find the analog filter's impulse response, we'd like to get Hc(s) into a form that allows us to use Laplace transform tables to find hc(t).
Searching through systems analysis textbooks we find the following Laplace transform pair:
Equation 656
Our intent, then, is to modify Eq. (655) to get it into the form on the left side of Eq. (656). We do this by realizing that the Laplace transform expression in Eq. (656) can be rewritten as
Equation 657
If we set Eq. (655) equal to the right side of Eq. (657), we can solve for A, a, and w. Doing that,
Equation 658
Solving Eq. (658) for A, a, and w, we first find
Equation 659
Equation 660
so
Equation 661
and
Equation 662
OK, we can now express Hc(s) in the desired form of the left side of Eq. (657) as
Equation 663
Using the Laplace transform pair in Eq. (656), the timedomain impulse response of the prototype analog filter becomes
Equation 664
OK, we're ready to perform Method 1, Step 4, to determine the discrete IIR filter's zdomain transfer function H(z) by performing the ztransform of hc(t). Again, scanning through digital signal processing textbooks or a good math reference book, we find the following ztransform pair where the timedomain expression is in the same form as Eq. (664)'s hc(t) impulse response:
Equation 665
Remember now, the a and w in Eq. (665) are generic and are not related to the a and w values in Eq. (659) and Eq. (661). Substituting the constants from Eq. (664) into the right side of Eq. (665), we get the ztransform of the IIR filter as
Equation 666
Performing Method 1, Step 5, we substitute the ts value of 0.01 for the continuous variable t in Eq. (666), yielding the final H(z) transfer function of
Equation 667
OK, hang in there; we're almost finished. Here are the final steps of Method 1. Because of the transfer function H(z) = Y(z)/X(z), we can crossmultiply the denominators to rewrite the bottom line of Eq. (667) as
or
Equation 668
By inspection of Eq. (668), we can now get the timedomain expression for our IIR filter. Performing Method 1, Steps 6 and 7, we multiply the x(n–1) coefficient by the sample period value of ts = 0.01 to allow for proper scaling as
Equation 669
and there we (finally) are. The coefficients from Eq. (669) are what we use in implementing the improved IIR structure shown in Figure 622 to approximate the original secondorder Chebyshev analog lowpass filter.
Let's see if we get the same result if we use the impulse invariance design Method 2 to approximate the example prototype analog filter.
6.4.2 Impulse Invariance Design Method 2 Example
Given the original prototype filter's Laplace transfer function as
Equation 670
and the value of ts = 0.01 for the sample period, we're ready to proceed with Method 2's Step 3. To express Hc(s) as the sum of singlepole filters, we'll have to factor the denominator of Eq. (670) and use partial fraction expansion methods. For convenience, let's start by replacing the constants in Eq. (670) with variables in the form of
Equation 671
where b = 137.94536, and c = 17410.145. Next, using Eq. (615) with a = 1, we can factor the quadratic denominator of Eq. (671) into
Equation 672
If we substitute the values for b and c in Eq. (672), we'll find that the quantity under the radical sign is negative. This means that the factors in the denominator of Eq. (672) are complex. Because we have lots of algebra ahead of us, let's replace the radicals in Eq. (672) with the imaginary term jR, where j = and R = (b2–4c)/4, such that
Equation 673
OK, partial fraction expansion methods allow us to partition Eq. (673) into two separate fractions of the form
Equation 674
where the K1 constant can be found to be equal to jc/2R and constant K2 is the complex conjugate of K1, or K2 = –jc/2R. (To learn the details of partial fraction expansion methods, the interested reader should investigate standard college algebra or engineering mathematics textbooks.) Thus, Hc(s) can be of the form in Eq. (648) or
Equation 675
We can see from Eq. (675) that our secondorder prototype filter has two poles, one located at p1 = –b/2 – jR and the other at p2 = –b/2 + jR. We're now ready to map those two poles from the splane to the zplane as called out in Method 2, Step 4. Making our substitution for the s + pk terms in Eq. (675), we have the following expression for the zdomain singlepole digital filters,
Equation 676
Our objective in Method 2, Step 5 is to massage Eq. (676) into the form of Eq. (652), so that we can determine the IIR filter's feed forward and feedback coefficients. Putting both fractions in Eq. (676) over a common denominator gives us
Equation 677
Collecting like terms in the numerator and multiplying out the denominator gives us
Equation 678
Factoring the exponentials and collecting like terms of powers of z in Eq. (678),
Equation 679
Continuing to simplify our H(z) expression by factoring out the real part of the exponentials,
Equation 680
We now have H(z) in a form with all the like powers of z combined into single terms, and Eq. (680) looks something like the desired form of Eq. (652). Knowing that the final coefficients of our IIR filter must be real numbers, the question is "What do we do with those imaginary j terms in Eq. (680)?" Once again, Euler to the rescue.[] Using Euler's equations for sinusoids, we can eliminate the imaginary exponentials and Eq. (680) becomes
[] From Euler, we know that sin(ø) = (ejø – e–jø)/2j, and cos(ø) = (ejø + e–jø)/2.
Equation 681
If we plug the values c = 17410.145, b = 137.94536, R = 112.48517, and ts = 0.01 into Eq. (681), we get the following IIR filter transfer function:
Equation 682
Because the transfer function H(z) = Y(z)/X(z), we can again crossmultiply the denominators to rewrite Eq. (682) as
or
Equation 683
Now we take the inverse ztransform of Eq. (683), by inspection, to get the timedomain expression for our IIR filter as
Equation 684
One final step remains. To force the IIR filter gain equal to the prototype analog filter's gain, we multiply the x(n–1) coefficient by the sample period ts as suggested in Method 2, Step 6. In this case, there's only one x(n) coefficient, giving us
Equation 685
that compares well with the Method 1 result in Eq. (669). (Isn't it comforting to work a problem two different ways and get the same result?)
Figure 627 shows, in graphical form, the result of our IIR design example. The splane pole locations of the prototype filter and the zplane poles of the IIR filter are shown in Figure 627(a). Because the splane poles are to the left of the origin and the zplane poles are inside the unit circle, both the prototype analog and the discrete IIR filters are stable. We find the prototype filter's splane pole locations by evaluating Hc(s) in Eq. (675). When s = –b/2 – jR, the denominator of the first term in Eq. (675) becomes zero and Hc(s) is infinitely large. That s = –b/2 – jR value is the location of the lower splane pole in Figure 627(a). When s = –b/2 + jR, the denominator of the second term in Eq. (675) becomes zero and s = –b/2 + jR is the location of the second splane pole.
Figure 627. Impulse invariance design example filter characteristics: (a) splane pole locations of prototype analog filter and zplane pole locations of discrete IIR filter; (b) frequency magnitude response of the discrete IIR filter.
The IIR filter's zplane pole locations are found from Eq. (676). If we multiply the numerators and denominators of Eq. (676) by z,
Equation 686
In Eq. (686), when z is set equal to the denominator of the first term in Eq. (686) becomes zero and H(z) becomes infinitely large. The value of z of
Equation 687
defines the location of the lower zplane pole in Figure 627(a). Specifically, this lower pole is located at a distance of = 0.5017 from the origin, at an angle of q = –Rts radians, or –64.45°. Being conjugate poles, the upper zplane pole is located the same distance from the origin at an angle of q = Rts radians, or +64.45°. Figure 627(b) illustrates the frequency magnitude response of the IIR filter in Hz.
Two different implementations of our IIR filter are shown in Figure 628. Figure 628(a) is an implementation of our secondorder IIR filter based on the general IIR structure given in Figure 622, and Figure 628(b) shows the secondorder IIR filter implementation based on the alternate structure from Figure 621(b). Knowing that the b(0) coefficient on the left side of Figure 628(b) is zero, we arrive at the simplified structure on the right side of Figure 628(b). Looking carefully at Figure 628(a) and the right side of Figure 628(b), we can see that they are equivalent.
Figure 628. Implementations of the impulse invariance design example filter.
Although both impulse invariance design methods are covered in the literature, we might ask, "Which one is preferred?" There's no definite answer to that question because it depends on the Hc(s) of the prototype analog filter. Although our Method 2 example above required more algebra than Method 1, if the prototype filter's sdomain poles were located only on the real axis, Method 2 would have been much simpler because there would be no complex variables to manipulate. In general, Method 2 is more popular for two reasons: (1) the inverse Laplace and ztransformations, although straightforward in our Method 1 example, can be very difficult for higher order filters, and (2) unlike Method 1, Method 2 can be coded in a software routine or a computer spreadsheet.
Upon examining the frequency magnitude response in Figure 627(b), we can see that this secondorder IIR filter's rolloff is not particularly steep. This is, admittedly, a simple loworder filter, but its attenuation slope is so gradual that it doesn't appear to be of much use as a lowpass filter.[] We can also see that the filter's passband ripple is greater than the desired value of 1 dB in Figure 626. What we'll find is that it's not the low order of the filter that contributes to its poor performance, but the sampling rate used. That secondorder IIR filter response is repeated as the shaded curve in Figure 629. If we increased the sampling rate to 200 Hz, we'd get the frequency response shown by the dashed curve in Figure 629. Increasing the sampling rate to 400 Hz results in the much improved frequency response indicated by the solid line in the figure. Sampling rate changes do not affect our filter order or implementation structure. Remember, if we change the sampling rate, only the sample period ts changes in our design equations, resulting in a different set of filter coefficients for each new sampling rate. So we can see that the smaller we make ts (larger fs) the better the resulting filter when either impulse invariance design method is used because the replicated spectral overlap indicated in Figure 624(b) is reduced due to the larger fs sampling rate. The bottom line here is that impulse invariance IIR filter design techniques are most appropriate for narrowband filters; that is, lowpass filters whose cutoff frequencies are much smaller than the sampling rate.
[] A piece of advice: whenever you encounter any frequency representation (be it a digital filter magnitude response or a signal spectrum) that has nonzero values at +fs/2, be suspicious—be very suspicious—that aliasing is taking place.
Figure 629. IIR filter frequency magnitude response, on a linear scale, at three separate sampling rates. Notice how the filter's absolute cutoff frequency of 20 Hz shifts relative to the different fs sampling rates.
The second analytical technique for analog filter approximation, the bilinear transform method, alleviates the impulse invariance method's aliasing problems at the expense of what's called frequency warping. Specifically, there's a nonlinear distortion between the prototype analog filter's frequency scale and the frequency scale of the approximating IIR filter designed using the bilinear transform. Let's see why.
URL http://proquest.safaribooksonline.com/0131089897/ch06lev1sec4
Amazon  