Example Digital Filters

Donald Daniel

Jan 2008, revised Sep 2016

up one level

If the lines of text are too long you can fix the problem with these instructions.

introduction

The filter technique presented here offers continuous processing, not batch processing like overlapping fft filtering. If your application is only to replace unfiltered stored data with filtered stored data, overlapping fft methods will be satisfactory. But if your application is the simulation of a complex system of which the filter is only a part, or if your application is a real time application, then you will probably prefer a filter with continuous processing instead of batch processing. Continuous filters require only a minimum number of input samples after a change before the output samples start to change. The older technique of continuous filters using z-transforms is difficult to learn and apply. The newer technique presented here with isolated poles and zeroes is easier.

This document shows example filter curves produced using the technique described in the article "Complex Digital Filters Using Isolated Poles and Zeroes". That article includes an introduction to the concept of digital filters of the simple type that simulate the behavior of analog filters. This kind of digital filter is the most appropriate for computer programs that simulate the behavior of lumped physical systems, whether electrical, mechanical, or acoustical. It is good for data smoothing. It is also the most appropriate for some kinds of audio signal processing. It could be applied to any application requiring digital filtering. If you are new to the subject that article and this one should be read together, they complement each other. To see that article click here. See also the derivation of a pole here and the derivation of the zero from the pole here.

To use the filter technique presented here, you start with a stream of real numbers, add zero imaginary part to convert it to a stream of complex numbers, feed it through the filter. Provided the pole-zero pattern is symmetrical about the real axis as it always is in the analog world, the output stream will have an imaginary part that is zero and can be discarded to produce a stream of real numbers for the output. You do not need a programming language with complex numbers built in. The language Oberon-2 used to program the examples shown here does not have complex numbers. The simple amount of complex arithmetic needed here is easy to program in any programming language.

the equations

In every example in this article, the gain of the filter is normalized to unity gain, zero dB, at some particular frequency. The other article explains how to do this. The pole and zero are given respectively by the following two equations, which are described more fully in the other article. "j" is the square root of minus one. The time interval between samples is "t", a constant.

pz.png

If the complex exponentials had to be computed every time a sample was fed through a pole or zero it would not be very efficient. Fortunately, if the location of the poles and zeroes does not change with time, the exponentials only have to be computed once.

The illustrations presented here were computed using oo2c and filter formulas primarily from the book "Synthesis of Passive Networks" by Ernst A. Guillemin, 1957. The source code of the computer program that computed the curves shown in this article is available as a zip file, if you do not have too much "clickophobia" to unpack a zip file and read the "readme" file. Be forewarned, most public wifi spots will not allow you to download a zip file. The source code is also listed at the end of this article.

simple filters

A pole on the negative real axis of the complex plane is shown in this example where the pole is shown as an "x" at the first tic mark to the left of the origin:

p1pln.png

Frequency response plots start at the origin and move up the vertical axis for positive frequencies, or down for negative frequencies. The frequencies shown in all plots are regular frequencies, they are radian frequencies divided by two pi. Linear plots show both positive and negative frequencies. Log plots show only positive frequencies. In the log plots "1.0E1" is 10 Hz, and "1.0E5" is 100000 Hz. A decade is a factor of 10 change in frequency. The log plots cover 4 decades. Each decade is marked with vertical lines in a 1, 2, 5 sequence. In most of the plots the upper limit of frequency is 1E5, which is half the sample frequency, and is the highest frequency where the digital filter can possibly operate. The log plots are all calculated from products of the pole-zero gain formulas, not measured by FFT's. This is very accurate. Most of the linear plots are FFT measurements of filter performance.

It should be mentioned that because we call our sample rate twice 1E5 that does not mean that calculations are actually done at that rate. In our software the size of the time step variable is just a number, not a time. Our calculations show what would happen if we had a hardware digital filter that actually operated at the chosen rate. But our computer program can compute what would happen at a rate much faster or slower than the software actually does the calculations.

The vertical axis represents the ratio of the output of the filter to its input. The solid curve represents the output of the digital filter. The dashed curve represents the response of a real analog filter that the digital filter simulates. If the simulation is perfect the two curves will overlap and the dashed curve will not be visible. The vertical axis in all plots is in decibels (dB). A decibel is ten times the logarithm to the base ten of a power ratio. Since power is proportional to the square of voltage, pressure or force, a decibel is twenty times the logarithm to the base ten of a ratio of voltage, pressure, or force. When the response is at zero dB, that means the output is the same magnitude as the input.

The next figure is a log plot of the previous example, where the negative pole is scaled to the value -1E3. Even though the pole is on the horizontal axis, it affects response as the frequency moves up the vertical axis. The response is 3dB down at 1E3 on the vertical axis. An octave is a factor of 2 change in frequency. The response is flat at low frequencies and rolls off to a slope of 6dB per octave or 20dB per decade. At the upper limit of frequency, 1E5, the simulation has deviated from the correct answer, the dashed line. Because this filter is flat at low frequencies and down at high frequencies, it is called a lowpass filter.

p1log.png

Next we put two poles in the same location where there was only one pole before. The result is similar, but now it is 6dB down at 1E3 and the slope is 12dB per octave or 40dB per decade:

p2log.png

In audio applications it is sometimes useful to have a filter with a slope of 3dB/octave, or 10dB/decade. Such a filter will convert white noise to pink noise with equal noise power in each octave. Here we give sigma numbers on the real axis that are not divided by two pi. A filter normalized to unity gain at zero frequency, with a pole sigma of -1E2, a zero sigma of -3.16E2, pole -1E3, zero -3.16E3, pole -1E4, zero -3.16E4, pole -1E5 will have the desired slope as shown here:

db3.png

We now put one pole on the negative real axis as before, and add a zero at the origin. Any filter with a zero at the origin models an analog filter that passes only alternating current (AC), blocking direct current (DC). Such behavior is called either "a DC block" or "AC coupled".

p1hppln.png

Since the pole and zero are both on the real axis this filter can be computed with real arithmetic, as shown in the audio programming example at this link. The result is now a highpass filter with a positive 6dB per octave slope at low frequencies and a 3dB point of 1E3:

p1hplog.png

We now present examples to explore the uses of a complex conjugate pair of poles 1E3 from the origin at 45 degrees from the imaginary axis:

p2lppln.png

The log plot is a lowpass filter which rolls off at 1E3 and goes down at 12dB per octave:

p2lplog.png

Next we add a zero at the origin:

p2bppln.png

The result is a bandpass filter centered at 1E3 with a 6dB per octave slope on either side:

p2bplog.png

Finally, we add a second zero at the origin, indicated by the number "2" nearby:

p2bplog.png

The log plot shows a highpass filter at 1E3 with a 12dB per octave slope at low frequencies:

p2bplog.png

Now we show a bandpass filter with the real component of the pole locations only 0.03 of the imaginary part:

p2nbpln.png

This simple combination of a complex pole pair and a zero at the origin is a simulation of any simple resonant system in nature. The log plot shows a very narrow spike, but lower down 6dB per octave slopes on each side. Note that a voltage ratio of 0.03 is about -30dB, and that is where the 6dB per octave slopes start:

p2nblog.png

If you wanted to make a filter with several such narrow peaks at widely separated frequencies, a single filter would not be the best way to go. It would be better to feed the input simultaneously into separate filters and add the outputs.

complex lowpass filter

We show an elliptical filter. The normalized pole-zero plot shows 5 poles and 4 zeroes:

eliptpln.png

We plot the elliptic filter where the filter bandwidth is one tenth the simulation bandwidth. The solid line is the simulated filter using isolated poles and zeroes in series, and the dotted line is the ideal transfer function:

elipt.png

step response

Thus far we have only discussed the frequency domain behavior of filters. Time domain behavior is also important. Suppose a lowpass filter has no input signal, then suddenly a step function signal is applied that rises instantly to a constant value. The time required for the output to rise to the constant value is called the rise time. The rise time is approximately equal to the reciprocal of the bandwidth measured from negative to positive frequencies, or half the reciprocal of the bandwith measured from zero to positive frequencies. Similarly, if a bandpass filter has a steady sinewave suddenly applied in the frequency range of the passband, the risetime of the envelope of the output sinewave will be approximately the reciprocal of the bandwidth.

The next two examples show two pole lowpass filters normalized to unity gain at zero frequency that both have 40dB/decade slopes that pass through -40dB at a frequency of 1E4. The two do not have identical frequency response. The first has a smooth rolloff, the second has a very high sharp peak at the rolloff frequency. Each tic on the complex plane is 1E3. Each pole is a distance of 1E3 from the origin. The input step function has zero imaginary part. The real part of the output is a solid line, the imaginary part is a dashed line. Because the pole/zero plots are symmetrical about the real axis, the output has zero imaginary part. The imaginary part could be discarded to make a stream of real numbers. First, poles at 63 degrees from the imaginary axis:

pln63.png

The corresponding time domain step response:

stp63.png

Next, poles at 12 degrees from the imaginary axis:

pln12.png

The corresponding time domain step response:

stp12.png

how to prevent numerical problems

An all pole filter can be expected to have good numerical performance. A filter with both poles and zeroes may not have good numerical performance if the order that the signal passes through the poles and zeroes is not correct. In this section we will first show how things can go bad with the wrong order, then how good they are with the correct order.

The elliptic filter example above has a bandwidth one tenth of the simulation bandwidth and behaves very well numerically. The signal is passed through all the poles, then all the zeroes, which is most likely to cause numerical difficulties. Later the improvement of a different arrangement will be shown. With 32 bit variables the numerical behavior is good with reducing filter bandwidth until a lower limit when the behavior gets worse rapidly. The following examples show the noisy frequency response and the smooth step response of an eliptic filter the same as above but with the bandwidth reduced to three percent of the simulation bandwidth. Notice the the small amount of noise on the frequency response at high frequencies. The signal is passed through all the poles first, then all the zeroes, which makes the numerical problems most evident.

elipt3.png

stp3.png

The following example shows the frequency response and step response of the same elliptic filter when the filter bandwidth is one percent of the simulation bandwidth.

eliptbd.png

stpbd.png

The numerical problems can be made to disappear by re-arranging the poles and zeroes so that there is always a pole immediately after each zero as the signal is passed through the poles and zeros. Perhaps the older z-transform methods did this automatically. The following plots show the result of this arrangement:

eliptgd.png

stpgd.png

To see better how well the problem has been solved, re-compile the program with the constant dbmax changed from 60 to 120. We now show a ten pole, ten zero chebychev high pass filter with one dB ripple. The ripple at the top is not roundoff noise, but the bad place at 95 dB down is roundoff noise. Presumably with 32 bit variables the technique can be relied on to have good numerical behavior until about 95 dB down.

chhp.png

The code used in my Oberon-2 program to feed the data through the poles and zeroes in series in such a way as to insert a pole after every zero is shown here. fr.n is the number of poles; fr.nz is the number of zeroes. pzxk is the input and pzyk is the output.

PROCEDURE pzfltr*(t:REAL;pzxk:cx.complex;
VAR pzyk:cx.complex; VAR fr:fltrrec);
(*feed data through poles and zeroes in series*)
VAR k:LONGINT;tempxk:cx.complex;
BEGIN
tempxk:=pzxk;
FOR k:=1 TO fr.n DO
IF k < fr.nz+1 THEN
fr.z[k].xk:=tempxk;zero(t,fr.f1,fr.z[k]); 
tempxk:=fr.z[k].yk;END(*IF*);
fr.p[k].xk:=tempxk;pole(t,fr.f1,fr.p[k]);
tempxk:=fr.p[k].yk;
END(*FOR*);
pzyk:=tempxk;END pzfltr;

To see the structure of the records referred to in the above code look for "fltrrec" in the third example when you click here.

effects of simulation bandwidth

One kind of multipole filter is the bessel filter, with maximally flat group delay in the passband, for minimal distortion of the shape of a pulse. Most kinds of filters have pole/zero patterns that must be calculated with complex mathematical formulas. But the most common multipole filter is the butterworth filter which has a pattern defined by a simple geometrical arrangement. The Butterworth lowpass filter has maximally flat amplitude response in the passband. However many poles there are, they are evenly spaced on a semicircle in the left half plane. The point where the semicircle intersects the vertical axis is where the filter is 3dB down. The size of the semicircle is arbitrary and determined by the bandwidth you want. The number of poles in the semicircle is arbitrary and determined by the steepness of the filter slopes you need. In the case of digital filters, more poles take more time to compute for each sample. In the case of analog filters, more poles require more hardware to build. Here is a normalized pole-zero plot of a butterworth filter with 4 poles.

btwpln.png

Next we look at the effect of simulation bandwidth. The simulation bandwidth extends to plus and minus half the sample frequency. Half the sample frequency is also called the Nyquist frequency. We start with a four pole butterworth filter centered at zero Hz with a bandwidth of 2E4, one tenth of the total simulation bandwidth:

btw4.png

We look at the same four pole butterworth filter with half the original simulation bandwidth:

btw4h.png

We end our bandwidth investigation by showing the same filter at one fourth the original simulation bandwidth:

btw4q.png

Some people have suggested ways to make digital filters that do not suffer so much distortion with limited simulation bandwidth. This may miss the point. If your simulation bandwidth is no wider than the bandwidth of the filter, why do you need a filter in the simulation?

transformation to highpass, bandpass, bandstop

Next we show how to transform from a lowpass prototype to highpass, bandpass and stopband filters. The transformations work the same way for lowpass prototypes that are all pole, and for those with both poles and zeroes. We start with a three pole butterworth lowpass filter.

LOW PASS PROTOTYPE

Pole pattern:

lppln.png

Linear frequency response:

lplin.png

Log frequency response:

lplog.png

HIGH PASS FILTER

Next, we apply the lowpass to highpass transformation to the preceding example to produce the corresponding highpass filter. Since our lowpass prototype is butterworth, the pole pattern is unchanged. The butterworth pattern is a special case in this regard because all the poles are equally distant from the origin. Any other type of filter would have pole positions altered by the transformation.

Pole-zero plot:

hppln.png

Linear frequency response:

hplin.png

Log frequency response:

hplog.png

The lowpass to highpass transformation requires a prototype that is symmetrical about the real axis. There will automatically be no poles at the origin since poles must have a negative real part. Choose a reflection frequency, a real number, typically the 3dB point of a lowpass prototype filter. The square of the reflection frequency divided by the complex location of each pole or zero in the lowpass prototype is the new location after the transformation. Guillemin's "Theory of Linear Physical Systems", 1963, p.257 says "If the number of finite zeroes is less than the number of poles, then a number of zeroes, equal to the difference, are at s=infinity." The number of new zeroes to be created at the origin is the number of poles less the number of zeroes in the original lowpass prototype. A zero at the origin would translate to a zero at infinity and disappear.

As another example, if we apply this transformation to the earlier two pole bandpass and the two pole high pass examples, the bandpass remains a bandpass, and the high pass becomes a low pass. The software listed below cannot handle cases with a zero at the origin because that would require dividing by zero.

BAND PASS FILTER

Next we transform our lowpass prototype to a bandpass filter. The lowpass version centered at zero frequency is symmetrical in a linear plot. The bandpass version would not be symmetrical in a linear plot, but is in a log frequency plot, except for errors near the Nyquist frequency, the upper limit of the simulation. If the bandwidth is very small compared with the center frequency, there is very little difference between the shape of a linear plot and a log plot, but if not, the difference can be large. Notice that the original lowpass pole pattern is not merely shifted, but also distorted. Notice that the passband, though narrow, is flat. If the original pole pattern had merely been shifted and not distorted, the passband would not be flat. Therefore the mathematical transformation is necessary to get the desired result.

Pole/zero plot:

bppln.png

Linear frequency response:

bplin.png

Log frequency response:

bplog.png

The lowpass to bandpass transformation is simple. The version of the transformation given here preserves the bandwidth from the negative frequency limit to the positive frequency limit, not from zero to the positive frequency limit. Start with the "lowpass prototype", which does not need to be a lowpass filter, it only needs to be symmetrical about the real axis on the complex plane. The number of poles should be equal to or greater than the number of zeroes. The difference between the number of poles and the number of zeroes will be the number of zeroes at the origin after the transformation. If the location of a pole or zero in the prototype is "a", the corresponding pair of locations after the transformation is "b". The geometrical mean frequency of the bandpass filter is "c", a real number. The geometric mean of two numbers is the square root of the product of the two numbers.

lp2bp.png

Unfortunately, the built-in complex square root routine in your favorite programming language is likely to be the wrong one. There is only one right answer for the square root of a positive real number, but there is more than one right answer for the square root of a complex number. That is to say, there is more than one number that when squared will give the original complex number. Unfortunately, only one of them will work in this application. The form of complex square root needed is the form having the principal root from zero to pi, not the form having the principal root from zero to plus or minus half pi. Thus the square root of 0-j1 on the negative imaginary axis should be -0.707+j0.707, in the upper left hand quadrant.

Having created a pair of poles or a pair of zeroes for each original pole or zero, add the appropriate number of zeroes at the origin as previously noted. When making bandpass filters that will appear narrow in a log frequency plot, "c" will be greater than the bandwidth of the prototype. When making filters that will appear broad in a log frequency plot, the bandwidth of the prototype is greater than "c". The transformation works equally well for all relative sizes of center frequency and bandwidth.

BAND STOP FILTER

The two transformations given above can be combined to make the lowpass to bandstop transformation. One must be careful to move the unity gain normalization frequency away from the stop band. Zero frequency is a good choice. Applying the lowpass to bandpass transformation to the highpass pole-zero pattern shown previously results in the following stop band example. Notice that the width of the 3dB points on the stopband is the same as the width of the 3dB points on the passband of the original lowpass prototype.

First, the pole-zero pattern. There are three overlapping zeroes in both places where one zero is shown. The program used to plot these figures could only indicate the number of overlapping zeroes in the case of zeroes at the origin.

sbpln.png

Linear frequency response:

sblin.png

Log frequency response:

sblog.png

The band stop filter serves as an example of the need for exact mathematical transformations instead of heuristic placement of poles and zeroes. Without the exact mathematical transformation the flat region between the two bandstop regions in the linear frequency plot would not be at the same level as the flat regions outside of the two bandstop regions.

The band stop example is disappointing because the stop band is so narrow at 50dB down. With more poles the situation improves because the filter slopes are steeper, as shown in this example where the only difference is that we started with a ten pole low pass prototype:

sblin2.png

WIDE BAND FILTER

The wide band filter is just another case of the bandpass filter, but with "c", the geometric mean frequency of the filter small compared with the bandwidth of the lowpass prototype. The lowpass to bandpass transformation is used just like before with no difference. In the case of a narrow bandpass filter, the transformation looks so simple that you are not convinced that you need the mathematical formula. The wideband case shows the power of the formula, and proves that you really do need the formula. For this example the lowpass prototype that we will use is not the three pole butterworth used in the previous examples in this section, but instead the elliptial lowpass filter shown in the "complex lowpass filter" section that extends from -1E4 to +1E4. The geometrical mean frequency "c" will be 1E3, which is only one tenth of 1E4. The bandwidth is 2E4, the same as the total bandwidth of the lowpass prototype. The log plot is:

wbelp.png

The pole-zero plot is not clear because of many poles and zeroes almost lumped together near the origin:

wbpln.png

If the plot is scaled up the stuff near the origin will look like a highpass version of the prototype. However, if we had chosen "c" to be 0.7E4 the pole zero plot would not have looked like a bandpass filter or like the combination of a lowpass and a highpass filter, but it would have produced the desired log plot.

smoothing filters

This article is about IIR filters implemented with poles and zeroes. The best known low pass filters of this type are butterworth, elliptical, tschebyscheff (also called chebychev), and bessel. These are usually used for frequency selective filters. Another application for filters is smoothing a noisy stream of data. Smoothing filters perform a weighted running average of noisy data to produce smooth data. The easiest way to get any weighting function you want is with FIR filters, not with IIR filters. The only one of the common IIR filters that is suitable for a smoothing filter is the bessel filter.

The output of a filter is the convolution of the input with the impulse response of the filter. The impulse response is the weighting function for a weighted running average of the input. But for a smoothing filter you want a well behaved impulse response. We are concerned with a low bandwidth signal that has been corrupted with additive wide bandwidth noise. We want to remove as much of the noise as possible without changing the shape of the signal waveform. The impulse response for a 5 pole butterworth filter that is 3dB down at 1kHz is shown here, the input impulse being simulated by one sample of input:

btmpls.png

The butterworth impulse response is not suitable for a smoothing filter because of the ringing after the main impulse. The impulse response for a 5 pole bessel filter that is down 3dB at 1kHz is much better as shown here:

bsmpls.png

The pole locations for bessel filters with various numbers of poles can be found on the internet. The 5 pole bessel filter used here is shown below where each tic is 1E3:

bspln.png

If you want to present the unfiltered data and the filtered data in the same illustration, you will need to delay the unfiltered data an amount equal to the delay of the peak of the impulse response so the two curves will have the same time scale. The unfiltered data will be written into an array equal in size to the number of samples of delay needed. As each new sample is available the active point on the array will be incremented to the next location until it reaches the end of the array and starts over at the beginning of the array. At the active point the oldest sample still in the array will be read out to be plotted in the illustration and the old sample will be replaced in the array by the newest sample that has just been recieved. This scheme implements a delay line many samples long.

A smoothing filter without a separate delay line preserves causality, but not coincidence. A smoothing filter with a separate delay line preserves coincidence, but not causality. With the delay line, the peak of a sudden impulse like the one that produced the impulse response will be coincident in both the unsmoothed and the smoothed data, but in the smoothed data the impulse will start to rise before the event that caused it, thus violating causality. For most data plotting purposes the preservation of coincidence is much more important for meeting the user's heuristic expectations than is the preservation of causality.

properties of poles and zeroes

Before looking at other complex examples, we should look at the properties of the poles and zeroes themselves. Remember that these simulations are of necessity circular functions and must match at the right and left edges of the simulation. First we show a pole at real -5E3 and imaginary 4E4:

p44r.png

The solid curve is the FFT measurement of the simulated pole, and the dashed curve is the transfer function of the actual pole, both normalized to unity gain at the frequency of the pole. The pole is at 4E4. The opposite side of the frequency circle is -6E4, where the simulated pole has flared out, but the actual pole is still monotonically decreasing. Near the frequency of the pole the simulation is quite accurate. On the opposite side of the frequency circle the simulation is not very accurate. The simulation is a good approximation for only about half the frequency circle. If we keep the poles and zeroes and our region of interest near the center of the simulation errors will be negligible. High sample rates widen the outer edges of the simulation and push the opposite side of the frequency circle far from the poles and zeroes.

The transfer function of a pole at a point on the complex plane is the reciprocal of the transfer function of a zero at the same point on the complex plane. Similarly, the gain function of a simulated pole at a point on the complex plane is the reciprocal of the gain function of a simulated zero at the same point on the complex plane, except for a constant vertical shift in dB. Therefore the errors of the simulated zero will be identical to the errors of the simulated pole shown above.

The gain functions of poles and zeroes are used to produce the solid curves in the log frequency plots in this article. Therefore it is of interest to see how accurately the gain functions match the FFT measurements of the simulated poles and zeroes. The simulated poles and zeroes shown are shown here with a solid line and the theoretical gain functions are shown in dashed lines. The dashed lines cannot be seen because the theoretical gain functions exactly overlap the FFT plots of the pole and zero functions. The theoretical gain functions are not to be confused with the ideal transfer functions, which are shown in the other plots in this article. The theoretical gain functions only predict the behavior of the simulated poles and zeros, whereas the transfer functions give the behavior of the actual poles and zeroes in a real filter, not a simulated filter. To start with we show a single pole with a real frequency component of -5E3 and a positive imaginary frequency component of 4E4 :

p44.png

Next, a single pole at -5E1, 4E4, :

p42.png

Now a zero at -5E1 and 4E4:

z42.png

In order for a filter to be physically realizable in analog hardware any complex poles and zeroes must appear in conjugate pairs at equal positive and negative imaginary frequencies. Only that way will a real signal in result in a real signal out, with zero imaginary parts. In the computer, however, we can simulate unrealizable filters with only positive or only negative imaginary frequency components. Now we make an unrealizable notch filter with only positive frequency poles and zeroes by combining the first pole example with the zero example. The notch filter is normalized to unity gain on the opposite side of the frequency circle from the notch. Both the simulated notch and the transfer function are plotted, but the transfer function overlaps the simulation and is invisible:

notch.png

Presumably the close match between the simulation and the actual transfer function is because the errors of the pole and zero are nearly equal and opposite in sign, and nearly cancel out.

The notch example previously shown was unrealizable because it did not have conjugate poles and zeroes. To create a realizable version, we move the same notch to zero frequency, then use a lowpass to bandpass transformation to move it to 1E4. First the normalized pole-zero pattern:

bpnchpln.png

Notice that each pole and the corresponding zero are no longer at the same imaginary frequency in this version, whereas they were in the previous version. Next a log frequency plot:

bpnchlog.png

Since the simulated pole and the simulated zero are both approximations and not exact, it is of interest to see how well they match. To see this the pole at -5E1, 4E4 is plotted with the zero at -5E1, 4E4 as a notch filter normalized as before. They match so well that the plot is a straight line as it should be:

strt.png

use of complex phasors

Now to examine the use of this technique for baseband simulation of narrowband filters centered at high frequencies. Filters such as this are used in radio receivers at the intermediate frequency (I.F.) of the receiver. The technique presented here is used to reduce the simulation bandwidth that would otherwise be required. Complex phasors are used in the data stream. Negative frequencies in the simulation represent positive frequencies in the real world that are below the center frequency of the narrowband filter. The real filter has poles at both positive and negative frequencies but the simulated filter represents only the positive poles in the real filter by means of the lowpass prototype centered at zero. This produces a distortion in the simulated response that we will now investigate. We show a four pole butterworth filter with bandwidth of 2E4 centered at zero in the simulation. We also plot the transfer function of a real filter with various I.F. frequencies. The I.F. frequency will be at zero frequency in the plot, and zero frequency in the real world will be offset to negative frequencies in the plot. The match between overlapping plots would be worse with fewer poles and better with more poles. First, an I.F. frequency of 3E4:

if3.png

Next we make the I.F. frequency 2E5, ten times the filter bandwidth:

if10.png

Finally, we make the I.F. frequency 2E6, 100 times the filter bandwidth:

if100.png

response to carrier signals

The previous section showed that the frequency response of a bandpass filter about its center frequency is similar, but not identical to the response of its lowpass prototype about zero frequency. Similarly the time domain response of a bandpass filter is in a sense similar but not identical to the lowpass prototype.

A step function into a lowpass filter could be said to have a carrier frequency of zero. An analogous signal for a bandpass filter is a sinewave at the center frequency of the bandpass filter that starts abruptly and remains at a constant amplitude. The carrier frequency is the frequency of the sinewave. This signal is the product of a step function with a sinewave. An impulse signal contains all frequencies equally, and could be applied either to a lowpass or a bandpass filter to compare their impulse responses. The impulse response of a lowpass 5 pole butterworth filter has already been shown previously. We now show the impulse response of a bandpass version of the same filter centered at 1E4:

bpmpls.png

Notice that the absolute value of the impulse response of the lowpass filter is similar but not identical to the envelope of the sinewave peaks of the impulse response of the bandpass filter. Similarly, the step response of the lowpass filter would be similar to the step response of a carrier signal at the center frequency of the bandpass filter.

Neither the lowpass or the bandpass filter has a step response that rises imediately to its maximum value. The risetime of the step response is roughly equal to the reciprocal of the bandwidth. The risetime of the step response is often taken to be the time between the 10% point and the 90% point on the step response. It is also roughly the width of the main part of the impulse response. The lowpass filter has a bandwidth of 2E3, from -1E3 to +1E3. The bandpass filter has the same bandwidth. The reciprocal of this is 0.5E3, which is about the risetime of the step response.

Suppose a carrier signal in the passband of the filter is modulated by a rectangular pulse that is zero everywhere except within the pulse, where the amplitude of the pulse is one. For pulses much shorter than the risetime of the step response the output will be similar to the impulse response shown above. If a rectangular pulse much shorter than the risetime of the step response is fed into the filter, the voltage out will be proportional to width of the pulse. In this case the convolution process may be thought of as producing an output that is a weighted running average of the impulse response, the weighting function being the narrow pulse. Since power is proportional to the square of voltage, the signal power out will be proportional to the square of the width of the pulse. Let us assume that the filter has been normalized to zero dB gain in the passband. For pulses longer than the rise time of the step response the output voltage will have time to rise up to the same level as the input voltage.

To interpret this in the frequency domain, notice that the shorter the pulse the wider its frequency spectrum, and the more of the energy falls outside of the passband of the filter, and is filtered out and not passed by the filter. The frequency bandwidth of most of the energy of the pulse spectrum is approximately the reciprocal of the pulse duration. For pulses longer than the rise time of the step response, most of the energy is within the filter bandwidth. For pulses much shorter than the rise time of the step response, much of the energy is outside of the filter bandwidth. The energy in the unfiltered pulse is proportional to the duration of the pulse if the amplitude is held constant. The combination of the two effects described in this paragraph means that the power of the impulse response of the filter to short pulses is proportional to the square of the width of the pulse, the same as before.

software

The software that produced the curves shown in this article is listed below. It appears correct in your browser, but in the raw html file it is not correct because of the process of translating to html. In the introduction section of this article is a link where you can download the software in correct form in a zip file. Listed first is the instructions for compiling and operating the software, then the software source code itself.


This program allows the user to experiment with and
demonstrate digital filters. It is provided in source code
that can be compiled and run on your computer. It is written
in oberon-2 an runs under linux. You can add linux to your
windows computer with a single hard disk. Instructions on
how to do this are provided in the article "how to program a
computer" at http://www.waltzballs.org/other/prog.html. The
article also explains how to install the oberon-2 compiler
oo2c on your linux system.

The program source files must have the suffix ".Mod" added
before they can be compiled in linux, but windows thinks
they are video files if they have that suffix.  The source
code must be compiled in the order cxarith, fft, cmpt then
finally cmptst.

When you run the program you see the main menu. The
choices are:

1. graf1: allows you to specify standard filters and test
them to see how they perform. Linear FFT plots and transfer
function plots are shown. Filter response is plotted as a
solid line, ideal transfer function as a dashed line.

2. graf2: does not test frequency response, but predicts
performance accurately. Linear and log plots are
shown. Solid line is filter response computed as product
of pole-zero gain functions. Ideal transfer function is
dashed line.  This mode does actually test step response
and impulse response.

3. save1: saves a plot to disk as the file temp5.eps.
If you have imagemagick software installed on your computer,
the command "convert temp5.eps temp5.png" will convert it to
a png file suitable for display on the internet. Slightly
better graphics and more convenient use can be attained by
the script "dfcvt" listed along with the source code.

4. save6: saves up to 6 plots to be printed on one sheet
of paper.

5 plot: prints the save6 file.

Within graf1 your first choice is screen width. This
is to show effects of varying data sample rate. Next
you choose between types of filter: A pole, a zero,
a notch, a butterworth lowpass filter, an I.F. filter,
a chebychev lowpass filter, a chebychev highpass filter,
and an elliptical lowpass filter. The number of poles can be
specified on most of the filters, but not on the elliptical
filter, which has 5 poles and 4 zeroes.

Within graf2 your choices are:

1. new: create a custom filter. rp:real pole; rz:real zero;
zo:zero at origin; cp:complex conjugate pair of poles;
cz:complex conjugate pair of zeroes. Normally you will use
small numbers from zero to three so you can see your pole
zero plot on the screen, but you can use large numbers if
you do not want to see the pole-zero plot.

2. scale: apply a scale factor to poles an zeros.  Enter
normalization frequency and scale factor. The normaliztion
frequency is the frequency that will be unity gain in the
plot. This should be the peak of a bandpass, or the shelf
of a highpass, lowpass, or stop band filter.  In addition
to the scale factor, you have a choice of multiplying by 1:
multiplying by one, 2: to multiply by two pi to go from
the complex plane to a frequency plot, or 3: to multiply
by the reciprocal of two pi to go from a frequency response
plot to the complex plane.

3. lp2bp: lowpass to bandpass transformation. Choice
1 does not multiply the center frequency by two pi,
and is apropriate for transforming plots on the complex
plane. Choice 2 does multiply the center frequency by two
pi, and is appropriate for transforming frequency response
plots. This command automatically sets the normalization
frequency to the center frequency, which is what you want
for a bandpass filter, but NOT what you want for a stopband
filter. For a stopband filter you must use setf1 after this
command before plotting.

4. lp2hp: lowpass to highpass transformation. As in the
previous command choice 1 is plain, choice 2 multiplies
the reflection frequency by two pi.

5. setf1: the lp2bp command sets the gain normalization
frequency f1 to the new bandpass frequency. But if the
"lowpass prototype" was a notch filter, not a lowpass filter,
this would not be appropriate, and a new f1 away from the
notch must be picked before plotting. IMPORTANT: If the
program attempts to plot a graph, but the plot "explodes"
and refuses to plot, it is usually because setf1 nees to
be used to move the normalization frequency.

6. pltpln: plot complex plane with poles and zeros.
If greater than 1, the number of zeroes at the origin
is inicate by a number below and to the right of the
origin. This numbering of zeroes only works at the origin.

7. pltlin: linear frequency response plot that is calulated,
not measured.

8. pltlog: log frequency response plot that is calculated,
not measured.

9. prtpz: write a text file "temp6" with pole-zero locations.
The gain term necessary to normalize the gain is shown
as a lumped term, but in the program it is applied as a
collection of distributed terms.

10. step: plot time domain step response.

11. mpls: plot time domain impulse response.

Now we will step you through a tutorial example of how to
use the program. Any part of the example that starts with
"---->" is something that you enter from the keyboard. The
rest is text that the program responds with.

Example 1: Stop band filter. Your linux system will need
the imagemagick package installed for the "display" command
in this example to work.

---->bin/cmptst
enter graf1, graf2, save6, save1 or plot
if graf1 or graf2, enter quit when finished
viewing the graf, enter q to exit program
enter command
---->graf1
enter 1, 2 or 4 do divide screen
--->1
screen fraction=  1.0000
enter graf type:
pgain zgain notch btwth iffltr cheby chhp elipt 
enter command
--->cheby
enter cf bw n db
---->0 2E4 3 2
---->quit
enter graf1, graf2, save6, save1 or plot
if graf1 or graf2, enter quit when finished
viewing the graf, enter q to exit program
enter command
---->graf2
enter q2 to quit this menu
enter operation:
new scale lp2bp lp2hp setf1
pltpln pltlin pltlog prtpz step mpls
enter command
--->scale
enter f1 before scaling, scale factor, choice
choice 1=1, 2=2*pi, 3=1/(2*pi)
--->0 1E-4 3
enter q2 to quit this menu
enter operation:
new scale lp2bp lp2hp setf1
pltpln pltlin pltlog prtpz step mpls
enter command
---->pltpln
---->quit
enter q2 to quit this menu
enter operation:
new scale lp2bp lp2hp setf1
pltpln pltlin pltlog prtpz step mpls
enter command
---->scale
enter f1 before scaling, scale factor, choice
choice 1=1, 2=2*pi, 3=1/(2*pi)
--->0 0.5 1
enter q2 to quit this menu
enter operation:
new scale lp2bp lp2hp setf1
pltpln pltlin pltlog prtpz step mpls
enter command
---->pltpln
---->quit
enter q2 to quit this menu
enter operation:
new scale lp2bp lp2hp setf1
pltpln pltlin pltlog prtpz step mpls
enter command
---->lp2hp
enter reflection frequency of highpass,
f1 normalization frequency
and 1 for plain or 2 for *2pi
all on one line: fm fn choice
---->0.5 9 1
enter q2 to quit this menu
enter operation:
new scale lp2bp lp2hp setf1
pltpln pltlin pltlog prtpz step mpls
enter command
---->pltpln
--->quit
enter q2 to quit this menu
enter operation:
new scale lp2bp lp2hp setf1
pltpln pltlin pltlog prtpz step mpls
enter command
---->lp2bp
enter geometric mean frequency of bandpass
and 1 for plain or 2 for *2pi
--->2 1
enter q2 to quit this menu
enter operation:
new scale lp2bp lp2hp setf1
pltpln pltlin pltlog prtpz step mpls
enter command
--->pltpln
---->quit
enter q2 to quit this menu
enter operation:
new scale lp2bp lp2hp setf1
pltpln pltlin pltlog prtpz step mpls
enter command
---->scale
enter f1 before scaling, scale factor, choice
choice 1=1, 2=2*pi, 3=1/(2*pi)
---->0 1E4 2
enter q2 to quit this menu
enter operation:
new scale lp2bp lp2hp setf1
pltpln pltlin pltlog prtpz step mpls
enter command
---->pltlin
---->quit
enter q2 to quit this menu
enter operation:
new scale lp2bp lp2hp setf1
pltpln pltlin pltlog prtpz step mpls
enter command
--->prtpz
output in temp6
enter q2 to quit this menu
enter operation:
new scale lp2bp lp2hp setf1
pltpln pltlin pltlog prtpz step mpls
enter command
--->q2
enter graf1, graf2, save6, save1 or plot
if graf1 or graf2, enter quit when finished
viewing the graf, enter q to exit program
enter command
--->save1
saved in temp5.eps
enter graf1, graf2, save6, save1 or plot
if graf1 or graf2, enter quit when finished
viewing the graf, enter q to exit program
enter command
--->q
---->display temp5.eps
---->cat temp6
fr.f1=          0.00000000
fr.gain=       3.54497248E+8
n nz nzo =  6  6  0
pole   1      -4.88958398E+3       9.69739688E+4
pole   2      -8.51584688E+4       9.24087500E+4
pole   3      -8.18988232E+3       1.62428062E+5
pole   4      -8.18988525E+3      -1.62428062E+5
pole   5      -8.51584531E+4      -9.24087422E+4
pole   6      -4.88958252E+3      -9.69739688E+4
zero  1          0.00000000       1.25663602E+5
zero  2          0.00000000       1.25663602E+5
zero  3          0.00000000       1.25663602E+5
zero  4          0.00000000      -1.25663602E+5
zero  5          0.00000000      -1.25663602E+5
zero  6          0.00000000      -1.25663602E+5


MODULE cxarith;
(*written 1980 donald daniel.  Originally written in pascal and
translated to Oberon-2 in 2000. This program is free software:
you can redistribute it and/or modify it under the terms of the
GNU General Public License as published by the Free Software
Foundation, version 3 of the License.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General
Public License along with this program.  If not, see
http://www.gnu.org/licenses/.  *)
IMPORT rm:=RealMath;
TYPE complex* = RECORD r*,x*:REAL END;
VAR one-,zero-,jone-,mone-:complex;
PROCEDURE neg*(z1:complex;VAR z2:complex);
BEGIN z2.r:=-z1.r;z2.x:=-z1.x END neg;
PROCEDURE conj*(z1:complex;VAR z2:complex);
BEGIN z2.r:=z1.r; z2.x:=-z1.x END conj;
PROCEDURE add*(z1,z2:complex;VAR z3:complex);
 BEGIN z3.r:=z1.r+z2.r;z3.x:=z1.x+z2.x END add;
PROCEDURE sub*(z1,z2:complex;VAR z3:complex);
 BEGIN z3.r:=z1.r-z2.r;z3.x:=z1.x-z2.x END sub;
PROCEDURE mult*(z1,z2:complex;VAR z3:complex);
 BEGIN z3.r := z1.r*z2.r-z1.x*z2.x;
 z3.x := z1.r*z2.x+z2.r*z1.x END mult;
PROCEDURE rmult*(r:REAL;z2:complex;VAR z3:complex);
BEGIN z3.r:=r*z2.r;z3.x:=r*z2.x;END rmult;
PROCEDURE div*(z1,z2:complex;VAR z3:complex);
(*z3:=z1/z2*)
VAR h:REAL;
BEGIN
IF (ABS(z2.r)>ABS(z2.x)) THEN
h:=z2.x/z2.r;z2.r:=h*z2.x+z2.r;
z3.r:=(z1.r+h*z1.x)/z2.r;
z3.x:=(z1.x-h*z1.r)/z2.r 
ELSE
h:=z2.r/z2.x;z2.x:=h*z2.r+z2.x;
z3.r:=(h*z1.r+z1.x)/z2.x;
z3.x:=(h*z1.x-z1.r)/z2.x END END div;
PROCEDURE cxr*(r:REAL;VAR z:complex);
BEGIN z.r:=r;z.x:=0;END cxr;
PROCEDURE cxj*(r:REAL;VAR z:complex);
BEGIN z.r:=0;z.x:=r;END cxj;
PROCEDURE cpx*(r,x:REAL;VAR z:complex);
BEGIN z.r:=r; z.x:=x  END cpx;
PROCEDURE abs*(z1:complex):REAL;
VAR h:REAL;
 BEGIN
 z1.r:=ABS(z1.r);z1.x:=ABS(z1.x);
 IF z1.x>z1.r THEN h:=z1.r;z1.r:=z1.x;z1.x:=h END;
 IF z1.x=0.0 THEN RETURN z1.r ELSE
 RETURN z1.r*rm.sqrt(1.0+(z1.x/z1.r)*(z1.x/z1.r)) END END abs;
PROCEDURE expj*(x:REAL;VAR z:complex);
(* e^jx *)
 BEGIN z.r:=rm.cos(x); z.x:=rm.sin(x) END expj;
PROCEDURE exp*(z1:complex;VAR z2:complex);
VAR x:REAL;
BEGIN
x:=rm.exp(z1.r);
z2.r:=x*rm.cos(z1.x);z2.x:=x*rm.sin(z1.x);END exp;
PROCEDURE sqrt*(z1:complex;VAR z2:complex);
VAR h:REAL;
BEGIN
h:=rm.sqrt((ABS(z1.r)+abs(z1))/2.0);
IF z1.x#0.0 THEN z1.x:=z1.x/(2.0*h)END;
IF z1.r>=0.0 THEN z1.r:=h
ELSIF z1.x>=0.0 THEN z1.r:=z1.x;z1.x:=h 
ELSE z1.r:=-z1.x;z1.x:=-h END;
(*the else part of the following if statement
adopts the convention that zero to pi, rather
thanzero to plus or minus half pi, is the
principal root*)
IF z1.x>=0.0 THEN z2.r:=z1.r;z2.x:=z1.x 
ELSE z2.r:=-z1.r;z2.x:=-z1.x;END;
END sqrt;

BEGIN
one.r:=1.0;one.x:=0.0;
zero.r:=0.0;zero.x:=0.0;
jone.r:=0.0;jone.x:=1.0;
mone.r:=-1.0;mone.x:=0.0 END cxarith.


MODULE fft;
(*written 1981 donald daniel.  Originally written in pascal
and translated to Oberon-2 in 2000. This program is free software:
you can redistribute it and/or modify it under the terms of the
GNU General Public License as published by the Free Software
Foundation, version 3 of the License.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General
Public License along with this program.  If not, see
http://www.gnu.org/licenses/.  *)

(*this fft operates from 1 to n, not 0 to n-1, so array must
be declared n+1 *)

(*time to freq involves exp(-jwt) so sign is -1.  For impulses
of (1,0) division by n is implicit, otherwize division by
n required.  Zero freq at ar[1].  If desired at ar[n2+1]
then calling prog must follow FFT with: FOR i:=1 TO n2 DO
temp:=ar[i]; ar[i]:=ar[i+n2];ar[i+n2]:=temp;END;*)

(*freq to time involves exp(+jwt) so sign is +1.  No division
by n needed. Calling prog will calculate positive freqs, mirror
complex conjugate to get negative freqs before calling FFT.*)

IMPORT rm:=RealMath,cx:=cxarith;
TYPE cxarray=ARRAY OF cx.complex;

PROCEDURE fft*(n:LONGINT;sign:REAL; VAR ar:cxarray);

PROCEDURE reorder(n:LONGINT; VAR
 ar:cxarray);
VAR nv2,nm1,i,j,k:LONGINT;
 t:cx.complex;
BEGIN
nv2 := n DIV 2;
nm1 := n-1;
j := 1;
FOR i := 1 TO nm1 DO 
 IF i < j THEN 
  t := ar[j];
  ar[j] := ar[i];
  ar[i] := t;
 END;
 k := nv2;
 WHILE k < j DO 
  j := j-k;
  k := k DIV 2;
 END;
 j := j+k
END(*FOR*);END reorder;

PROCEDURE calculate(sign:REAL;
n:LONGINT;VAR ar:cxarray);
CONST pi = 3.1415926; 
VAR 
 u,w,t:cx.complex;
 m,i,j,l,le,le1,ip:LONGINT;ang:REAL;
BEGIN
m := ENTIER((rm.ln(n)/rm.ln(2))+0.5);
FOR l := 1 TO m DO 
 le := ENTIER((rm.power(2.0,l))+0.5);
 le1 := le DIV 2;
 u.r := 1.0; 
 u.x := 0.0;
 ang := sign*pi/le1;
 w.r := rm.cos(ang);
 w.x := rm.sin(ang);
 FOR j := 1 TO le1 DO 
  i := j;
  WHILE i < n DO 
   ip := i+le1;
   cx.mult(ar[ip],u,t);
   cx.sub(ar[i],t,ar[ip]);
   cx.add(ar[i],t,ar[i]);
   i := i+le;
  END(*WHILE*);
  cx.mult(u,w,u)
 END(*FOR*)
END(*FOR*)
END calculate;

BEGIN
IF (sign > 0.0)THEN sign:=1.0 ELSE sign:=-1.0;END;
reorder(n,ar);
calculate(sign,n,ar);
END fft;

BEGIN END fft.


MODULE cmpt;
(*written by donald daniel 2011.  This program is free software:
you can redistribute it and/or modify it under the terms of the
GNU General Public License as published by the Free Software
Foundation, version 3 of the License.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General
Public License along with this program.  If not, see
http://www.gnu.org/licenses/.  *)
IMPORT cx:=cxarith,rm:=RealMath,Out;
CONST maxpoles=20; pi=3.141592654;
TYPE
polerec*=RECORD sigma*, omega*, gain*:REAL;
xk*,yk*,ykm1*,onemeps*,rt1meps*,eps*,gaincx:cx.complex;
first*:BOOLEAN; END;
zerorec*=RECORD sigma*,omega*,gain*:REAL;xk*,yk*,xkm1*,
rt1meps*,norm*,gaincx:cx.complex;first*:BOOLEAN;END;
(*pole and zero first is so fixed terms only calculated once*)
para*=ARRAY maxpoles+1 OF polerec;
zara*=ARRAY maxpoles+1 OF zerorec;
fltrrec*=RECORD fc*,bw*,f1*,db*,gain*:REAL;n*,nz*,nzo*:LONGINT;
        firstb*:BOOLEAN;p*:para;z*:zara END;
(*fltrrec firstb is so location of poles and zeroes
only calculated once*)
VAR s:ARRAY maxpoles+maxpoles OF cx.complex;
(* s  is temporary variable 
includes both poles and zeroes *)

PROCEDURE pgain*(t,w:REAL;VAR p:polerec);
VAR dwt,st,temp,num,denom,quot:cx.complex;
BEGIN
IF ~(p.sigma < 0.0)THEN
Out.String('pole sigma must be negative');Out.Ln;
HALT(1);END;
cx.expj((p.omega-w)*t,dwt);cx.cpx(rm.exp(p.sigma*t),0.0,st);
cx.mult(dwt,st,temp);cx.sub(cx.one,temp,denom);
cx.sub(cx.one,st,num);cx.div(num,denom,quot);
p.gaincx:=quot;
p.gain:=cx.abs(quot);
END pgain;

PROCEDURE zgain*(t,w:REAL;VAR z:zerorec);
VAR dwt,st,temp,num,denom,quot:cx.complex;
BEGIN
cx.expj((z.omega-w)*t,dwt);
cx.cpx(rm.exp(z.sigma*t),0.0,st);
cx.mult(st,dwt,temp);cx.sub(cx.one,temp,num);
cx.add(cx.one,st,denom);cx.div(num,denom,quot);
z.gaincx:=quot;
z.gain:=cx.abs(quot);
END zgain;

PROCEDURE pole*(t,w:REAL;VAR p:polerec;VAR fr:fltrrec);
VAR rotate,t1:cx.complex;
BEGIN
IF p.first THEN
cx.cpx(rm.exp(p.sigma*t),0.0,p.onemeps);
cx.cpx(1.0-rm.exp(p.sigma*t),0.0,p.eps);
cx.expj(p.omega*t,rotate);
cx.mult(p.onemeps,rotate,p.rt1meps);p.ykm1:=cx.zero;
p.first:=FALSE;pgain(t,w,p); END(*if*);
cx.mult(p.ykm1,p.rt1meps,p.yk);cx.mult(p.xk,p.eps,t1);
cx.add(p.yk,t1,p.yk);p.ykm1:=p.yk;
(*note that gain is not applied to ykm1*)
cx.rmult(1.0/p.gain,p.yk,p.yk);END pole;

PROCEDURE zero*(t,w:REAL;VAR z:zerorec;VAR fr:fltrrec);
VAR rotate,t1,onemeps:cx.complex;
BEGIN
IF z.first THEN
cx.cpx(rm.exp(z.sigma*t),0.0,onemeps);
cx.cpx(1.0/(1.0+rm.exp(z.sigma*t)),0.0,z.norm);
cx.expj(z.omega*t,rotate);
cx.mult(onemeps,rotate,z.rt1meps);z.xkm1:=cx.zero;
z.first:=FALSE;zgain(t,w,z); END(*if*);
cx.mult(z.xkm1,z.rt1meps,t1);cx.sub(z.xk,t1,z.yk);
cx.mult(z.norm,z.yk,z.yk);z.xkm1:=z.xk;
(*note that gain is not applied to xkm1*)
cx.rmult(1.0/z.gain,z.yk,z.yk);END zero;

PROCEDURE pzshift(VAR fr:fltrrec);
(*scale and shift poles and zeroes in frequency
no use of s[k] after pzshift*)
VAR k,nt:LONGINT; bwcx:cx.complex;
BEGIN
nt:=fr.n+fr.nz;
fr.f1:=2*pi*fr.f1;fr.fc:=2*pi*fr.fc;fr.bw:=2*pi*fr.bw;
(*At this point we are normalized to low pass w3db=1.0Hz*)
fr.bw:=fr.bw*0.5;(*lowpass to bandpass*)
cx.cpx(fr.bw,0.0,bwcx);
FOR k:=1 TO nt DO cx.mult(s[k],bwcx,s[k]);END;
(*At this point we are centered at 0.0 Hz*)
FOR k:=1 TO nt DO s[k].x:=s[k].x+fr.fc;END;
(*At this point we are centered at fc Hz*)
FOR k:=1 TO fr.n DO
fr.p[k].sigma:=s[k].r;fr.p[k].omega:=s[k].x;
fr.p[k].first:=TRUE;END;
FOR k:=1 TO fr.nz DO
fr.z[k].sigma:=s[fr.n+k].r;fr.z[k].omega:=s[fr.n+k].x;
fr.z[k].first:=TRUE;END;END pzshift;

PROCEDURE lp2hp*(wr:REAL;VAR fr:fltrrec);
VAR i:LONGINT;wrcx,temp:cx.complex;
(*wr squared is divided by the pole-zero pattern.
This transformation is from Ernst Guillemin's books:
"synthesis of passive networks" p.603, or
"theory of linear physical systems" p.273.
Will not work for zeroes at the origin*)
BEGIN
cx.cpx(wr*wr,0.0,wrcx);
FOR i:=1 TO fr.n DO
cx.cpx(fr.p[i].sigma,fr.p[i].omega,temp);
cx.div(wrcx,temp,temp);
fr.p[i].sigma:=temp.r;fr.p[i].omega:=temp.x;
(*some have already been initialized, but new ones 
need to be*)
fr.p[i].first:=TRUE;
END;
FOR i:=1 TO fr.nz DO
cx.cpx(fr.z[i].sigma,fr.z[i].omega,temp);
cx.div(wrcx,temp,temp);
fr.z[i].sigma:=temp.r;fr.z[i].omega:=temp.x;
fr.z[i].first:=TRUE;
END;
FOR i:=1 TO (fr.n-fr.nz) DO
fr.z[i+fr.nz].sigma:=0.0;fr.z[i+fr.nz].omega:=0.0;
fr.z[i+fr.nz].first:=TRUE;
END;
fr.nzo:=fr.n-fr.nz;
fr.nz:=fr.nzo+fr.nz;
END lp2hp;

PROCEDURE lp2bp*(wm:REAL;VAR fr:fltrrec);
VAR w2,diff,s1,s2,temp:cx.complex;i:LONGINT;
nzo:LONGINT;
(*wm is the geometric mean frequency of the bandpass.
This transformation is from Ernst Guillemin's books:
"synthesis of passive networks" p.604 eq.56, or
"theory of linear physical systems" p.274 eq.66, 
but simplified to preserve bandwidth from minus
3db to plus 3db, rather than from 0 to plus 3db*)
BEGIN
nzo:=fr.n-fr.nz;
cx.cpx(wm,0.0,temp); cx.mult(temp,temp,w2);
FOR i:=1 TO fr.n DO
cx.cpx(fr.p[i].sigma,fr.p[i].omega,s1);
cx.mult(s1,s1,temp);cx.sub(temp,w2,diff);
cx.sqrt(diff,temp);
cx.add(s1,temp,s2);
fr.p[i].sigma:=s2.r;fr.p[i].omega:=s2.x;
fr.p[i].first:=TRUE;
cx.sub(s1,temp,s2);
fr.p[i+fr.n].sigma:=s2.r;fr.p[i+fr.n].omega:=s2.x;
(*some have already been initialized, but new ones 
need to be*)
fr.p[i+fr.n].first:=TRUE;
END;
fr.n:=2*fr.n;
FOR i:=1 TO fr.nz DO
cx.cpx(fr.z[i].sigma,fr.z[i].omega,s1);
cx.mult(s1,s1,temp);cx.sub(temp,w2,diff);
cx.sqrt(diff,temp);
cx.add(s1,temp,s2);
fr.z[i].sigma:=s2.r;fr.z[i].omega:=s2.x;
fr.z[i].first:=TRUE;
cx.sub(s1,temp,s2);
fr.z[i+fr.nz].sigma:=s2.r;fr.z[i+fr.nz].omega:=s2.x;
fr.z[i+fr.nz].first:=TRUE;
END;
fr.nz:=2*fr.nz;
FOR i:=1 TO nzo DO
fr.z[i+fr.nz].sigma:=0.0;fr.z[i+fr.nz].omega:=0.0;
fr.z[i+fr.nz].first:=TRUE;
END;
fr.nz:=fr.nz+nzo; fr.nzo:=nzo; fr.f1:=wm;
END lp2bp;

PROCEDURE iffltr*(fm:REAL;VAR fr:fltrrec);
(*transform baseband filter to double sideband
IF filter s2=(s1+-sqrt(s1*s1-4*wc*wc))/2 then
shift so upper sideband is centered on baseband,
for use in plotting transfer functions only*)
VAR i:LONGINT;
BEGIN
fm:=2*pi*fm; lp2bp(fm,fr);
FOR i:=1 TO fr.n DO 
fr.p[i].omega:=fr.p[i].omega-fm;END;
FOR i:=1 TO fr.nz DO
fr.z[i].omega:=fr.z[i].omega-fm;END;
fr.f1:=0.0;
END iffltr;

PROCEDURE pzfltr*(t:REAL;pzxk:cx.complex;
VAR pzyk:cx.complex; VAR fr:fltrrec);
(*feed data through poles and zeroes in series*)
VAR k:LONGINT;tempxk:cx.complex;
BEGIN
tempxk:=pzxk;
FOR k:=1 TO fr.n DO
IF k < fr.nz+1 THEN
fr.z[k].xk:=tempxk;zero(t,fr.f1,fr.z[k],fr); 
tempxk:=fr.z[k].yk;END(*IF*);
fr.p[k].xk:=tempxk;pole(t,fr.f1,fr.p[k],fr);
tempxk:=fr.p[k].yk;
END(*FOR*);
pzyk:=tempxk;END pzfltr;

(*this is numerically bad version that does not put
poles between zeroes
PROCEDURE pzfltr*(t:REAL;pzxk:cx.complex;
VAR pzyk:cx.complex; VAR fr:fltrrec);
(*feed data through poles and zeroes in series*)
VAR k:LONGINT;
BEGIN
fr.p[1].xk:=pzxk;pole(t,fr.f1,fr.p[1],fr);
FOR k:=1 TO fr.n-1 DO
fr.p[k+1].xk:=fr.p[k].yk;pole(t,fr.f1,fr.p[k+1],fr);END;
IF fr.nz>0 THEN
fr.z[1].xk:=fr.p[fr.n].yk;zero(t,fr.f1,fr.z[1],fr);
FOR k:=1 TO fr.nz-1 DO
fr.z[k+1].xk:=fr.z[k].yk;zero(t,fr.f1,fr.z[k+1],fr);END;
pzyk:=fr.z[fr.nz].yk;
ELSE pzyk:=fr.p[fr.n].yk;END (*IF*);END pzfltr;*)

PROCEDURE pzgain*(t,f:REAL;VAR fr:fltrrec;VAR prod:cx.complex);
(*fltr gain at f by pole/zero gain products*)
VAR pn,zn,w1,w:REAL;pg,zg:cx.complex;i:LONGINT;
BEGIN
w1:=fr.f1;w:=2*pi*f;
prod:=cx.one;
FOR i:=1 TO fr.n DO
pgain(t,w1,fr.p[i]);pn:=fr.p[i].gain;
pgain(t,w,fr.p[i]);pg:=fr.p[i].gaincx;
cx.rmult(1/pn,pg,pg);cx.mult(pg,prod,prod);END;
FOR i:=1 TO fr.nz DO
zgain(t,w1,fr.z[i]);zn:=fr.z[i].gain;
zgain(t,w,fr.z[i]);zg:=fr.z[i].gaincx;
cx.rmult(1/zn,zg,zg);cx.mult(zg,prod,prod);END;
END pzgain;

PROCEDURE fltrgain*(t:REAL;VAR fr:fltrrec);
(*lumped gain term in front of poles and zeroes*)
VAR pn,zn,w1:REAL;i:LONGINT;
BEGIN
w1:=fr.f1;
FOR i:=1 TO fr.n DO
pgain(t,w1,fr.p[i]);pn:=fr.p[i].gain;
fr.gain:=fr.gain/pn;
END;
FOR i:=1 TO fr.nz DO
zgain(t,w1,fr.z[i]);zn:=fr.z[i].gain;
fr.gain:=fr.gain/zn;
END;
END fltrgain;

PROCEDURE pzxf*(VAR fr:fltrrec;f:REAL;VAR y:cx.complex);
(*transfer function for filters with poles and zeroes*)
VAR w1,wp,wz,w,quot,prod,vectw,vectwc:cx.complex;k:LONGINT;
BEGIN
cx.cpx(0.0,fr.f1,w1);cx.cpx(0.0,2*pi*f,w);prod:=cx.one;
FOR k:=1 TO fr.n DO
cx.cpx(fr.p[k].sigma,fr.p[k].omega,wp);
cx.sub(w,wp,vectw);cx.sub(w1,wp,vectwc);
cx.div(vectwc,vectw,quot);
cx.mult(quot,prod,prod);END;
FOR k:=1 TO fr.nz DO
cx.cpx(fr.z[k].sigma,fr.z[k].omega,wz);
cx.sub(w,wz,vectw);cx.sub(w1,wz,vectwc);
cx.div(vectw,vectwc,quot);
cx.mult(quot,prod,prod);END;
y:=prod; END pzxf;

PROCEDURE butterworth(n:LONGINT);
(*from spn p.590 eq. 5*)
VAR k:LONGINT;BEGIN
FOR k:=1 TO n DO
cx.expj(pi*(2*k+n-1)/(2*n),s[k]); END;END butterworth;

PROCEDURE chebychev(n:LONGINT;db:REAL);
VAR k:LONGINT;eps,sigma,omega:REAL;
PROCEDURE sinh(x:REAL):REAL;
BEGIN RETURN 0.5*(rm.exp(x)-rm.exp(-x));END sinh;
PROCEDURE cosh(x:REAL):REAL;
BEGIN RETURN 0.5*(rm.exp(x)+rm.exp(-x));END cosh;
PROCEDURE arcsinh(x:REAL):REAL;
BEGIN RETURN rm.ln(x+rm.sqrt(x*x+1.0));END arcsinh;
BEGIN
eps:=rm.sqrt(rm.power(10.0,0.1*db)-1);
FOR k:=1 TO n DO
sigma:=
-sinh((1.0/n)*arcsinh(1.0/eps))*rm.sin(((2*k-1)*rm.pi)/(2*n));
omega:=
cosh((1.0/n)*arcsinh(1.0/eps))*rm.cos(((2*k-1)*rm.pi)/(2*n));
cx.cpx(sigma,omega,s[k]);
END;END chebychev;

PROCEDURE elliptic;
(*special case only*)
BEGIN
cx.cpx(-0.730,0.0,s[1]);
cx.cpx(-0.422,0.761,s[2]);cx.cpx(-0.422,-0.761,s[3]);
cx.cpx(-0.1026,0.974,s[4]);cx.cpx(-0.1026,-0.974,s[5]);
cx.cpx(0,1.27,s[6]);cx.cpx(0,-1.27,s[7]);
cx.cpx(0,1.88,s[8]);cx.cpx(0,-1.88,s[9]);
END elliptic;

PROCEDURE btrwrth*(t:REAL;pzxk:cx.complex;
VAR pzyk:cx.complex;VAR fr:fltrrec);
BEGIN
IF fr.firstb THEN fr.firstb:=FALSE;
fr.nz:=0;
butterworth(fr.n);
pzshift(fr);
END;(*firstb*)
pzfltr(t,pzxk,pzyk,fr);
END btrwrth;

PROCEDURE elipt*(t:REAL;pzxk:cx.complex;
VAR pzyk:cx.complex;VAR fr:fltrrec);
BEGIN
IF fr.firstb THEN fr.firstb:=FALSE;
fr.n:=5;fr.nz:=4;
elliptic;
pzshift(fr);
END;(*firstb*)
pzfltr(t,pzxk,pzyk,fr);
END elipt;

PROCEDURE chebylp*(t:REAL;pzxk:cx.complex;
VAR pzyk:cx.complex;VAR fr:fltrrec);
BEGIN
IF fr.firstb THEN fr.firstb:=FALSE;
fr.nz:=0;
chebychev(fr.n,fr.db);
pzshift(fr);
END;(*firstb*)
pzfltr(t,pzxk,pzyk,fr);
END chebylp;

PROCEDURE chebyhp*(t:REAL;pzxk:cx.complex;
VAR pzyk:cx.complex;VAR fr:fltrrec);
BEGIN
IF fr.firstb THEN fr.firstb:=FALSE;
fr.nz:=0;
chebychev(fr.n,fr.db);
pzshift(fr);
lp2hp(fr.bw,fr);
END;(*firstb*)
pzfltr(t,pzxk,pzyk,fr);
END chebyhp;

PROCEDURE onep*(t:REAL;pzxk:cx.complex;
VAR pzyk:cx.complex;VAR fr:fltrrec);
BEGIN
IF fr.firstb THEN fr.firstb:=FALSE;
fr.n:=1;fr.nz:=0;
fr.f1:=2*pi*fr.f1;fr.fc:=2*pi*fr.fc;
fr.p[1].sigma:=2*pi*fr.p[1].sigma; 
fr.p[1].omega:=fr.fc; fr.p[1].first:=TRUE; 
END;(*firstb*)
pzfltr(t,pzxk,pzyk,fr);
END onep;

PROCEDURE onez*(t:REAL;pzxk:cx.complex;
VAR pzyk:cx.complex;VAR fr:fltrrec);
BEGIN
IF fr.firstb THEN fr.firstb:=FALSE;
fr.n:=0;fr.nz:=1;
fr.f1:=2*pi*fr.f1;fr.fc:=2*pi*fr.fc;
fr.z[1].sigma:=2*pi*fr.z[1].sigma; 
fr.z[1].omega:=fr.fc; fr.z[1].first:=TRUE;
END;(*firstb*)
(*pzfltr will not work here*)
fr.z[1].xk:=pzxk;zero(t,fr.f1,fr.z[1],fr);
pzyk:=fr.z[1].yk;
END onez;

PROCEDURE notch*(t:REAL;pzxk:cx.complex;
VAR pzyk:cx.complex;VAR fr:fltrrec);
BEGIN
IF fr.firstb THEN fr.firstb:=FALSE;
fr.n:=1;fr.nz:=1;
fr.f1:=2*pi*fr.f1;fr.fc:=2*pi*fr.fc;
fr.p[1].sigma:=2*pi*fr.p[1].sigma; 
fr.p[1].omega:=fr.fc;
fr.z[1].sigma:=2*pi*fr.z[1].sigma; 
fr.z[1].omega:=fr.fc;
fr.p[1].first:=TRUE; fr.z[1].first:=TRUE;
END;(*firstb*)
pzfltr(t,pzxk,pzyk,fr);
END notch;

PROCEDURE hpfltr*(t:REAL;pzxk:cx.complex;
VAR pzyk:cx.complex;VAR fr:fltrrec);
(*pzfltr will not work here*)
PROCEDURE highpass(t,fhp,xk:REAL;VAR yk:REAL);
CONST pi2=6.283185;
VAR est,gf:REAL;
BEGIN
IF fr.firstb THEN fr.firstb:=FALSE;
fr.z[0].xkm1.r:=0.0;fr.p[0].ykm1.r:=0.0;
END;
est:=rm.exp(-pi2*fhp*t);
gf:=(1+est)/(1-est);
yk:=(xk-fr.z[0].xkm1.r)/2;fr.z[0].xkm1.r:=xk; 
xk:=yk;yk:=est*fr.p[0].ykm1.r+(1-est)*xk;fr.p[0].ykm1.r:=yk;
yk:=gf*yk;
END highpass;
BEGIN
highpass(t,fr.fc,pzxk.r,pzyk.r);
pzyk.x:=0.0;
END hpfltr;

BEGIN END cmpt.


MODULE cmptst;
(*written 2011 donald daniel.  This program is free software:
you can redistribute it and/or modify it under the terms of the
GNU General Public License as published by the Free Software
Foundation, version 3 of the License.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General
Public License along with this program.  If not, see
http://www.gnu.org/licenses/.  *)

(*section 1. Global declarations*)

IMPORT Msg, In, Out, Files, TextRider, rm:=RealMath, 
OS:ProcessManagement,
 RealStr, IntStr, cx:=cxarith, cmpt, fft;

CONST
pi=3.14159;strlngth=7;
q=1;er=2;graf1=3;save6=4;plot=5;
grid=6;notch=7;btwth=8;elipt=9;
cheby=10;chhp=11;
pgain=12;save1=13;zgain=14;iffltr=15;graf2=16;
rp=17;rz=18;zo=19;cp=20;cz=21;
scale=22;lp2bp=23;pltlin=24;pltlog=25;new=26;
q3=27;q2=28;pltpln=29;prtpz=30;lp2hp=31;setf1=32;
step=33;hpf=34;mpls=35;last=36;
(*insert additional commands before last and
increment the value of last*)
up=1;down=2;

TYPE
str=ARRAY strlngth OF CHAR;
cma=ARRAY last+1 OF str;
pt=PROCEDURE(t:REAL;apxk:cx.complex;
VAR apyk:cx.complex; VAR fr:cmpt.fltrrec);

VAR
t:REAL;first:BOOLEAN;
cmdara:cma;command,penstate,pos:LONGINT;
filevar:Files.File;outfile:TextRider.Writer;
res:Msg.Msg;irts:LONGINT;
fr:cmpt.fltrrec;

(*section 2. Graphics utilities*)

PROCEDURE pento(xr,yr:REAL);
VAR str:ARRAY 40 OF CHAR;
BEGIN
RealStr.RealToFixed(xr,2,str);
outfile.WriteString(str);
outfile.WriteString(' ');
RealStr.RealToFixed(yr,2,str);
outfile.WriteString(str);
IF penstate=up THEN outfile.WriteString(' moveto');
ELSE outfile.WriteString(' lineto');END;
outfile.WriteLn;
END pento;

PROCEDURE width(w:REAL);
VAR str1:ARRAY 40 OF CHAR;
BEGIN
RealStr.RealToFixed(w,2,str1);
outfile.WriteString(str1);
outfile.WriteString(' setlinewidth');
outfile.WriteLn;
outfile.WriteString('stroke');outfile.WriteLn;
END width;

PROCEDURE dash(w:REAL);
(*it is assumed that width will be called
right after dash, because dash has no 'stroke' *)
VAR str:ARRAY 40 OF CHAR;
BEGIN
outfile.WriteString('[ ');
RealStr.RealToFixed(w,2,str);
outfile.WriteString(str);
outfile.WriteString(' ');
RealStr.RealToFixed(w,2,str);
outfile.WriteString(str);
outfile.WriteString(' ] 0 setdash'); outfile.WriteLn;
END dash;

PROCEDURE pen(ps:LONGINT);
BEGIN
penstate:=ps;
END pen;

PROCEDURE initgraf;
BEGIN
filevar:=Files.New('temp1',{Files.write},res);
outfile:=TextRider.ConnectWriter(filevar);
END initgraf;

PROCEDURE showgraf;
VAR resr,resw:Msg.Msg;invar,outvar:Files.File;
infile:TextRider.Reader;outfile:TextRider.Writer;
str:ARRAY 40 OF CHAR;
BEGIN
filevar.Close;
invar:=Files.Old('temp1',{Files.read},resr);
infile:=TextRider.ConnectReader(invar);
outvar:=Files.New('temp2',{Files.write},resw);
outfile:=TextRider.ConnectWriter(outvar);
outfile.WriteString('%!PS');outfile.WriteLn;
LOOP infile.ReadLine(str);
IF infile.res#TextRider.done THEN EXIT END; 
outfile.WriteString(str);
outfile.WriteLn;END;invar.Close;
outfile.WriteString('showpage');outfile.WriteLn;
outvar.Close;irts:=ProcessManagement.system("gs -sDEVICE=x11 temp2");
Out.String('saved in temp2');Out.Ln;
END showgraf;

PROCEDURE save6proc;
(*Stores up to six plots in memory to be plotted
on one sheet of paper*)
VAR r,c:LONGINT;resr,resw:Msg.Msg;
yofst:ARRAY 4 OF REAL;xofst:ARRAY 3 OF REAL;
str:ARRAY 40 OF CHAR;
invar,outvar:Files.File;
infile:TextRider.Reader;outfile:TextRider.Writer;
BEGIN
Out.String('max six plots per page');Out.Ln;
Out.String('rows 1,2,3; cols 1,2');Out.Ln;
Out.String('enter row col for plot');Out.Ln;
In.LongInt(r);In.LongInt(c);
yofst[1]:=530.0;yofst[2]:=320.0;yofst[3]:=110.0;
xofst[1]:=71.0;xofst[2]:=287.0;
invar:=Files.Old('temp1',{Files.read},resr);
infile:=TextRider.ConnectReader(invar);
IF first THEN
outvar:=Files.New('temp3',{Files.write},resw);
outfile:=TextRider.ConnectWriter(outvar);
outfile.WriteString('%!PS');outfile.WriteLn;
ELSE
outvar:=Files.Old('temp3',{Files.write},resw);
outfile:=TextRider.ConnectWriter(outvar);
outfile.SetPos(pos);
END;
outfile.WriteString('gsave');outfile.WriteLn;
RealStr.RealToFixed(xofst[c],2,str);
outfile.WriteString(str);
outfile.WriteString(' ');
RealStr.RealToFixed(yofst[r],2,str);
outfile.WriteString(str);
outfile.WriteString(' translate'); outfile.WriteLn;
outfile.WriteString(' 0.45 0.45 scale');
outfile.WriteLn;
LOOP infile.ReadLine(str);
IF infile.res#TextRider.done THEN EXIT END; 
outfile.WriteString(str);outfile.WriteLn;END;
invar.Close;
outfile.WriteString('grestore');outfile.WriteLn;
pos:=outfile.Pos();outvar.Close;
Out.String('saved in temp3');Out.Ln;first:=FALSE;
END save6proc;

PROCEDURE save1proc;
(*saves one plot to disk in temp5.eps*)
VAR resr,resw:Msg.Msg;
str1,str2:ARRAY 40 OF CHAR;
invar,outvar:Files.File;
infile:TextRider.Reader;outfile:TextRider.Writer;
BEGIN
str1:='';str2:='';
invar:=Files.Old('temp1',{Files.read},resr);
infile:=TextRider.ConnectReader(invar);
outvar:=Files.New('temp5.eps',{Files.write},resw);
outfile:=TextRider.ConnectWriter(outvar);
outfile.WriteString('%!PS');outfile.WriteLn;
outfile.WriteString('%%BoundingBox: 98 68 542 320');
outfile.WriteLn;
LOOP infile.ReadLine(str1);
IF infile.res#TextRider.done THEN EXIT END; 
outfile.WriteString(str1);outfile.WriteLn;END;
outfile.WriteString('showpage');outfile.WriteLn;
outfile.WriteString('%%EOF');outfile.WriteLn;
invar.Close;
outvar.Close;
Out.String('saved in temp5.eps');Out.Ln;
END save1proc;

PROCEDURE plotproc;
(*writes plot buffer to the laser printer*)
VAR resr,resw:Msg.Msg;invar,outvar:Files.File;
infile:TextRider.Reader;outfile:TextRider.Writer;
str:ARRAY 40 OF CHAR;
BEGIN
invar:=Files.Old('temp3',{Files.read},resr);
infile:=TextRider.ConnectReader(invar);
outvar:=Files.New('temp4',{Files.write},resw);
outfile:=TextRider.ConnectWriter(outvar);
LOOP infile.ReadLine(str);
IF infile.res#TextRider.done THEN EXIT END; 
outfile.WriteString(str);
outfile.WriteLn;END;invar.Close;
outfile.WriteString('showpage');outfile.WriteLn;
outvar.Close;
irts:=ProcessManagement.system("lpr temp4");
(*lpr is assumed to send file through Ghostscript
to laser printer*)
Out.String('sent to plotter and saved in temp4');Out.Ln;
first:=TRUE;
END plotproc;

(*section 3. graphics application*)

PROCEDURE^ determine(VAR cmdvar:LONGINT);
(*forward declaration to permit command
interpretation in grafproc*)

PROCEDURE grafproc(choice:LONGINT);
CONST maxgrid=361.0; fmax=100000.0;dbmax=60.0;
fftmax=4096; 
VAR fftsize,n0,n2:LONGINT;
scrfrac,maxfreq,maxdb,mindb:REAL;
ar:ARRAY fftmax+1 OF cx.complex;
(*fftmax+1 because fft operates 1-4096, not 0-4095 *)

PROCEDURE xpos(f:REAL):REAL;(*f in, x out*)
CONST xmargin=145.0;
BEGIN RETURN ((maxgrid/2.0)*(1.0+f/fmax)+xmargin);
END xpos;

PROCEDURE ypos(dbv:REAL):REAL;(*db in, y out*)
CONST ymargin=117.0;
BEGIN
RETURN ( ymargin+(maxgrid/2.0)*((dbv-mindb)/60.0));
END ypos;

PROCEDURE db(v,vref:REAL):REAL;
VAR dbv:REAL;BEGIN
IF((v/vref)<(1.0E-7))THEN dbv:=-140.0
ELSIF((v/vref)>(1.0E1))THEN dbv:=+20.0
ELSE dbv:=20.0*0.43429*rm.ln(v/vref)END;
RETURN dbv;END db;

PROCEDURE initfont;
BEGIN
outfile.WriteString('/Palatino-Roman 16 selectfont');
outfile.WriteLn;
END initfont;

PROCEDURE initsfont;
BEGIN
outfile.WriteString('/Symbol 16 selectfont');
outfile.WriteLn;
END initsfont;

PROCEDURE lbldb(db:REAL);
VAR str:ARRAY 40 OF CHAR;
BEGIN
pen(up);
pento(xpos(-maxfreq)-0.07*maxgrid,ypos(db-2));
outfile.WriteString('(');
RealStr.RealToFixed(db,-1,str);
outfile.WriteString(str);
outfile.WriteString(') show'); outfile.WriteLn;
END lbldb;

PROCEDURE lblfreq;
VAR str:ARRAY 40 OF CHAR;

PROCEDURE writefreq(f,offset:REAL);
BEGIN
pen(up);pento(xpos(f-2.3E3),ypos(mindb-4));
outfile.WriteString('(|) show');outfile.WriteLn;
pen(up);pento(xpos(f-offset),ypos(mindb-8));
outfile.WriteString('(');
RealStr.RealToFloat(f,2,str);
outfile.WriteString(str);
outfile.WriteString(') show');
outfile.WriteLn;END writefreq;

BEGIN
writefreq(-maxfreq,1.0E4);writefreq(0,4.0E3);
writefreq(maxfreq,1.0E4);
pen(up);pento(xpos(-1.8E4),ypos(mindb-13));
outfile.WriteString('(frequency Hz) show');
outfile.WriteLn;
END lblfreq;

PROCEDURE lblgain;
VAR dbh:REAL;
BEGIN
dbh:=(mindb)/2;
pen(up);pento(xpos(-maxfreq)-28.0,ypos(dbh));
outfile.WriteString('90 rotate');outfile.WriteLn;
outfile.WriteString('(gain dB) show');outfile.WriteLn;
outfile.WriteString('-90 rotate');outfile.WriteLn;
END lblgain;

PROCEDURE gridproc(choice:LONGINT);
VAR f,db:REAL;i:LONGINT;
BEGIN
initgraf;
initfont;
IF choice=1 THEN
Out.String('enter 1, 2 or 4 do divide screen');Out.Ln;
In.LongInt(i);CASE i OF 1: scrfrac:=1.0;fftsize:=fftmax;
|2:scrfrac:=0.5;fftsize:=fftmax DIV 2;
|4:scrfrac:=0.25;fftsize:=fftmax DIV 4;END;
ELSE scrfrac:=1.0;fftsize:=fftmax;END;
n2:=fftsize DIV 2;n0:=n2+1;
Out.String('screen fraction=');Out.Real(scrfrac,8,5);Out.Ln;
maxfreq:=fmax*scrfrac;t:=1.0/(2.0*maxfreq);
FOR i:=0 TO (ENTIER(dbmax/10.0))DO
db:=maxdb-i*10.0;
lbldb(db);
pento(xpos(-maxfreq),ypos(db));pen(down);
pento(xpos(maxfreq),ypos(db));pen(up);END;
pento(xpos(-maxfreq),ypos(maxdb));pen(down);
pento(xpos(-maxfreq),ypos(mindb));pen(up);
pento(xpos(maxfreq),ypos(maxdb));pen(down);
pento(xpos(maxfreq),ypos(mindb));pen(up);
lblfreq;
FOR i:=0 TO ENTIER(maxfreq/10000.0)DO
f:=i*10000.0;
pento(xpos(f),ypos(maxdb));pen(down);
pento(xpos(f),ypos(mindb));pen(up);
pento(xpos(-f),ypos(maxdb));pen(down);
pento(xpos(-f),ypos(mindb));pen(up);END;
width(1.0);
lblgain;
END gridproc;

PROCEDURE clipgrid;
VAR str:ARRAY 60 OF CHAR;
width,height:REAL;
BEGIN
RealStr.RealToFixed(xpos(-maxfreq),-1,str);
outfile.WriteString(str);
outfile.WriteString(' ');
RealStr.RealToFixed(ypos(mindb),-1,str);
outfile.WriteString(str);
outfile.WriteString(' ');
width:=xpos(maxfreq)-xpos(-maxfreq);
height:=ypos(10)-ypos(mindb);
RealStr.RealToFixed(width,-1,str);
outfile.WriteString(str);
outfile.WriteString(' ');
RealStr.RealToFixed(height,-1,str);
outfile.WriteString(str);
outfile.WriteString(' rectclip'); outfile.WriteLn;
END clipgrid;

PROCEDURE pltfltr(VAR fr:cmpt.fltrrec;proc:pt);
VAR i:LONGINT;f:REAL;
first:BOOLEAN;xk,yk,temp:cx.complex;
BEGIN
clipgrid;
first:=TRUE;fr.firstb:=TRUE;
FOR i:=1 TO fftsize DO
IF first THEN first:=FALSE;xk:=cx.one;
ELSE xk:=cx.zero;END;
proc(t,xk,yk,fr);
ar[i]:=yk;END;
fft.fft(fftsize,-1.0,ar);first:=TRUE;
FOR i:=1 TO n2 DO
temp:=ar[i];ar[i]:=ar[i+n2];ar[i+n2]:=temp;END;
FOR i:=1 TO fftsize DO
f:=((i-n0)/n2)*maxfreq;
pento(xpos(f),ypos(db(cx.abs(ar[i]),1.0)));
IF first THEN first:=FALSE;pen(down);END;
END;pen(up);width(2.8);
END pltfltr;

PROCEDURE pltxfr(VAR fr:cmpt.fltrrec);
VAR i:LONGINT;f:REAL;
first:BOOLEAN;y:cx.complex;
BEGIN
first:=TRUE;
FOR i:=1 TO fftsize DO
f:=((i-n0)/n2)*maxfreq;
cmpt.pzxf(fr,f,y);
pento(xpos(f),ypos(db(cx.abs(y),1.0)));
IF first THEN first:=FALSE;pen(down);END;
END;pen(up);dash(2.8);width(2.8);
END pltxfr;

PROCEDURE btwthproc;
BEGIN
fr.n:=0;fr.nz:=0;fr.nzo:=0;
Out.String('enter cf bw n');Out.Ln;
In.Real(fr.fc);In.Real(fr.bw);In.LongInt(fr.n);
fr.f1:=fr.fc;
pltfltr(fr,cmpt.btrwrth);
pltxfr(fr);
END btwthproc;

PROCEDURE ifproc;
VAR fm:REAL;
BEGIN
fr.n:=0;fr.nz:=0; fr.nzo:=0;
Out.String('enter I.F. bw n');Out.Ln;
In.Real(fm);In.Real(fr.bw);In.LongInt(fr.n);
fr.fc:=0.0; fr.f1:=fr.fc;
pltfltr(fr,cmpt.btrwrth);
cmpt.iffltr(fm,fr);
pltxfr(fr);
END ifproc;

PROCEDURE eliptproc;
BEGIN
fr.n:=0;fr.nz:=0; fr.nzo:=0; 
Out.String('enter cf bw');Out.Ln;
In.Real(fr.fc);In.Real(fr.bw);
fr.f1:=fr.fc;
pltfltr(fr,cmpt.elipt);
pltxfr(fr);
END eliptproc;

PROCEDURE chebyproc;
BEGIN
fr.n:=0;fr.nz:=0; fr.nzo:=0;
Out.String('enter cf bw n db');Out.Ln;
In.Real(fr.fc);In.Real(fr.bw);In.LongInt(fr.n);
In.Real(fr.db);
fr.f1:=fr.fc;
pltfltr(fr,cmpt.chebylp);
pltxfr(fr);
END chebyproc;

PROCEDURE chhpproc;
BEGIN
fr.n:=0;fr.nz:=0; fr.nzo:=0;fr.fc:=0.0;
Out.String('enter bw n db');Out.Ln;
In.Real(fr.bw);In.LongInt(fr.n);
In.Real(fr.db);
IF fr.fc < 0.0 THEN fr.f1:=fr.fc+maxfreq; 
ELSE fr.f1:=fr.fc-maxfreq;END;
pltfltr(fr,cmpt.chebyhp);
pltxfr(fr);
END chhpproc;

PROCEDURE pgainproc;
(*this procedure does not show accuracy of pole,
just equivalence between simulated pole and 
computed pole gain formula*)
VAR i:LONGINT;f,w:REAL; first:BOOLEAN; 
BEGIN
fr.n:=0;fr.nz:=0; fr.nzo:=0;
Out.String('enter frequencies, not radian frequencies');Out.Ln;
Out.String('negative sigma/2pi and omega/2pi');Out.Ln;
Out.String('enter real and imaginary frequencies');Out.Ln;
In.Real(fr.bw); In.Real(fr.fc);
fr.n:=1; first:=TRUE;fr.firstb:=TRUE;
fr.f1:=fr.fc;
fr.p[1].sigma:=fr.bw;
pltfltr(fr,cmpt.onep);
(*we do not want pltxfr, we want pgain here*)
first:=TRUE;
FOR i:=1 TO fftsize DO
f:=((i-n0)/n2)*maxfreq;
w:=2.0*pi*f;
cmpt.pgain(t,w,fr.p[1]);
pento(xpos(f),ypos(db(fr.p[1].gain,1.0)));
IF first THEN first:=FALSE;pen(down);END;
END;pen(up);dash(2.8);width(2.8);
END pgainproc;

PROCEDURE zgainproc;
(*this procedure does not show accuracy of zero,
just equivalence between simulated zero and 
computed zero gain formula*)
VAR i:LONGINT;f,w:REAL; first:BOOLEAN; 
BEGIN
fr.n:=0;fr.nz:=0; fr.nzo:=0;
Out.String('enter frequencies, not radian frequencies');Out.Ln;
Out.String('negative sigma/2pi and omega/2pi');Out.Ln;
Out.String('enter real and imaginary frequencies');Out.Ln;
In.Real(fr.bw); In.Real(fr.fc);
IF fr.fc < 0.0 THEN fr.f1:=fr.fc+maxfreq; 
ELSE fr.f1:=fr.fc-maxfreq;END;
fr.z[1].sigma:=fr.bw;
pltfltr(fr,cmpt.onez);
(*we do not want pltxfr, we want zgain here*)
first:=TRUE;
FOR i:=1 TO fftsize DO
f:=((i-n0)/n2)*maxfreq;
w:=2.0*pi*f;
cmpt.zgain(t,w,fr.z[1]);
pento(xpos(f),ypos(db(fr.z[1].gain,1.0)));
IF first THEN first:=FALSE;pen(down);END;
END;pen(up);dash(2.8);width(2.8);
END zgainproc;

(*PROCEDURE hpfproc;
BEGIN
Out.String('enter 3db freq');Out.Ln;
In.Real(fr.fc);
pltfltr(fr,cmpt.hpfltr);
END hpfproc;*)

PROCEDURE notchproc;
VAR bwpole,bwzero:REAL; 
BEGIN
fr.n:=0;fr.nz:=0; fr.nzo:=0;
Out.String('enter frequencies, not radian frequencies');Out.Ln;
Out.String('negative sigma/2pi and omega/2pi');Out.Ln;
Out.String('enter sigma-pole sigma-zero fc');Out.Ln;
In.Real(bwpole); In.Real(bwzero); In.Real(fr.fc); 
IF fr.fc < 0.0 THEN fr.f1:=fr.fc+maxfreq; 
ELSE fr.f1:=fr.fc-maxfreq;END;
fr.p[1].sigma:=bwpole;fr.z[1].sigma:=bwzero;
pltfltr(fr,cmpt.notch);
pltxfr(fr);
END notchproc;

PROCEDURE pltpz(fr:cmpt.fltrrec);
VAR i:LONGINT;

PROCEDURE xpz(f:REAL):REAL; (*f in x out*)
BEGIN
RETURN xpos(f*1.6666E4);END xpz;

PROCEDURE ypz(f:REAL):REAL; (*f in y out*)
BEGIN
RETURN ypos(f*10-20.0);END ypz;

PROCEDURE cxplain;
BEGIN
initgraf; 
pen(up);pento(xpz(-3),ypz(0));pen(down);
pento(xpz(3),ypz(0));pen(up);
pento(xpz(2.8),ypz(-0.1));pen(down);
pento(xpz(3),ypz(0));pen(up);
pento(xpz(2.8),ypz(0.1));pen(down);
pento(xpz(3),ypz(0));pen(up);
initsfont;
outfile.WriteString('(s/2p) show');outfile.WriteLn;
pento(xpz(0),ypz(-3));pen(down);
pento(xpz(0),ypz(3));pen(up);
pento(xpz(-0.1),ypz(2.8));pen(down);
pento(xpz(0),ypz(3));pen(up);
pento(xpz(0.1),ypz(2.8));pen(down);
pento(xpz(0),ypz(3));pen(up);
pento(xpz(-0.2),ypz(3.2));
initfont;
outfile.WriteString('(j) show');outfile.WriteLn;
initsfont;
outfile.WriteString('(w/2p) show');outfile.WriteLn;
pento(xpz(-2),ypz(-0.1));pen(down);
pento(xpz(-2),ypz(0.1));pen(up);
pento(xpz(-1),ypz(-0.1));pen(down);
pento(xpz(-1),ypz(0.1));pen(up);
pento(xpz(1),ypz(-0.1));pen(down);
pento(xpz(1),ypz(0.1));pen(up);
pento(xpz(2),ypz(-0.1));pen(down);
pento(xpz(2),ypz(0.1));pen(up);
pento(xpz(-0.1),ypz(-2));pen(down);
pento(xpz(0.1),ypz(-2));pen(up);
pento(xpz(-0.1),ypz(-1));pen(down);
pento(xpz(0.1),ypz(-1));pen(up);
pento(xpz(-0.1),ypz(1));pen(down);
pento(xpz(0.1),ypz(1));pen(up);
pento(xpz(-0.1),ypz(2));pen(down);
pento(xpz(0.1),ypz(2));pen(up);
END cxplain;

PROCEDURE pltx(x,y:REAL);
CONST l=0.1;
BEGIN
pen(up);
pento(xpz(x-l),ypz(y-l));pen(down);
pento(xpz(x+l),ypz(y+l));pen(up);
pento(xpz(x-l),ypz(y+l));pen(down);
pento(xpz(x+l),ypz(y-l));pen(up);
pen(up);
END pltx;

PROCEDURE pltz(x,y:REAL);
CONST r=0.1;
VAR str:ARRAY 40 OF CHAR;
BEGIN
outfile.WriteString('stroke');outfile.WriteLn;
RealStr.RealToFixed(xpz(x),2,str);
outfile.WriteString(str);
outfile.WriteString(' ');
RealStr.RealToFixed(ypz(y),2,str);
outfile.WriteString(str);
outfile.WriteString(' ');
RealStr.RealToFixed(ypz(r)-ypz(0),2,str);
outfile.WriteString(str);
outfile.WriteString(' ');
outfile.WriteString(' 0 360 arc');outfile.WriteLn;
pen(up);
END pltz;

PROCEDURE pltzcnt;
VAR str:ARRAY 40 OF CHAR;
BEGIN
pen(up);pento(xpz(0.2),ypz(-0.5));
initfont;
outfile.WriteString('(');
IntStr.IntToStr(fr.nzo,str);
outfile.WriteString(str);
outfile.WriteString(' ) show');
outfile.WriteLn;
END pltzcnt;

BEGIN
cxplain;
FOR i:=1 TO fr.n DO
pltx(fr.p[i].sigma,fr.p[i].omega);END;
FOR i:=1 TO fr.nz DO
pltz(fr.z[i].sigma,fr.z[i].omega);END;
IF fr.nzo > 1 THEN pltzcnt;END;
width(1.0);
END pltpz;

PROCEDURE newpz(VAR fr:cmpt.fltrrec);

PROCEDURE getdata;
VAR i,j:LONGINT;
BEGIN
i:=0;j:=0;fr.nzo:=0;
REPEAT
Out.String('enter q3 to quit this menu');Out.Ln;
Out.String('enter choice:');Out.Ln;
Out.String('rp rz zo cp cz');Out.Ln;
determine(command);
CASE command OF
er:;
|rp:Out.String('enter negative location');Out.Ln;
INC(i); In.Real(fr.p[i].sigma);fr.p[i].omega:=0.0;
|rz:Out.String('enter location');Out.Ln;
INC(j); In.Real(fr.z[j].sigma);fr.z[j].omega:=0.0;
|zo:Out.String('zero at origin created');Out.Ln;
INC(fr.nzo);
INC(j); fr.z[j].sigma:=0.0;fr.z[j].omega:=0.0;
|cp:Out.String('enter real then imaginary parts');Out.Ln;
INC(i); In.Real(fr.p[i].sigma);In.Real(fr.p[i].omega);
INC(i); fr.p[i].sigma:=fr.p[i-1].sigma;
fr.p[i].omega:=-fr.p[i-1].omega;
|cz:Out.String('enter real then imaginary parts');Out.Ln;
INC(j); In.Real(fr.z[j].sigma);In.Real(fr.z[j].omega);
INC(j); fr.z[j].sigma:=fr.z[j-1].sigma;
fr.z[j].omega:=-fr.z[j-1].omega;
|q3:;
END(*case*);
UNTIL command=q3;fr.n:=i;fr.nz:=j;
END getdata;

BEGIN
getdata;
END newpz;

PROCEDURE scaleproc(VAR fr:cmpt.fltrrec);
VAR i,choice:LONGINT;scale,factor:REAL;
BEGIN
Out.String('enter f1 before scaling, scale factor, choice');
Out.Ln;
Out.String('choice 1=1, 2=2*pi, 3=1/(2*pi)');Out.Ln;
In.Real(fr.f1);In.Real(scale);In.LongInt(choice);
CASE choice OF
1:factor:=1|
2:factor:=2*pi|
3:factor:=1.0/(2*pi);END;
fr.f1:=factor*fr.f1*scale;
FOR i:=1 TO fr.n DO 
fr.p[i].sigma:=factor*scale*fr.p[i].sigma;
fr.p[i].omega:=factor*scale*fr.p[i].omega;END;
FOR i:=1 TO fr.nz DO 
fr.z[i].sigma:=factor*scale*fr.z[i].sigma;
fr.z[i].omega:=factor*scale*fr.z[i].omega;END;
END scaleproc;

PROCEDURE setf1proc(VAR fr:cmpt.fltrrec);
BEGIN
Out.String('enter f1 normalization frequency');Out.Ln;
In.Real(fr.f1);
fr.f1:=2*pi*fr.f1;
END setf1proc;

PROCEDURE pltlinproc(fr:cmpt.fltrrec);
(*linear frequency plot by gain products*)
PROCEDURE pltgflt(VAR fr:cmpt.fltrrec);
VAR i:LONGINT;f:REAL;prod:cx.complex;
first:BOOLEAN;
BEGIN
clipgrid; first:=TRUE;
FOR i:=1 TO fftsize DO
f:=((i-n0)/n2)*maxfreq;
cmpt.pzgain(t,f,fr,prod);
pento(xpos(f),ypos(db(cx.abs(prod),1.0)));
IF first THEN first:=FALSE;pen(down);END;END;
pen(up);width(2.8);
END pltgflt;

BEGIN
gridproc(2); pltgflt(fr); pltxfr(fr); 
END pltlinproc;

PROCEDURE pltlogproc(fr:cmpt.fltrrec);
(*log frequency plot by gain products*)
CONST fl=10.0; npts=200;
VAR left,right,dist:REAL;

PROCEDURE xlpos(f:REAL):REAL;
BEGIN
RETURN left + dist*(rm.ln(f/fl)/rm.ln(fmax/fl));
END xlpos;

PROCEDURE loggrid;
VAR db:REAL;i:LONGINT;f:ARRAY 20 OF REAL;

PROCEDURE writefreq(f:REAL);
VAR str:ARRAY 40 OF CHAR;
BEGIN
pen(up);pento(xlpos(f)-4.5,ypos(mindb-4));
outfile.WriteString('(|) show');outfile.WriteLn;
pen(up);pento(xlpos(f)-20,ypos(mindb-8));
outfile.WriteString('(');
RealStr.RealToFloat(f,2,str);
outfile.WriteString(str);
outfile.WriteString(') show');
outfile.WriteLn; END writefreq;

BEGIN
initgraf; initfont; scrfrac:=1.0;
maxfreq:=fmax*scrfrac;t:=1.0/(2.0*maxfreq);
left:=xpos(-maxfreq);right:=xpos(maxfreq);
dist:=right-left;
FOR i:=0 TO (ENTIER(dbmax/10.0))DO
db:=maxdb-i*10.0;
lbldb(db);
pento(xpos(-maxfreq),ypos(db));pen(down);
pento(xpos(maxfreq),ypos(db));pen(up);END;
pento(xpos(-maxfreq),ypos(maxdb));pen(down);
pento(xpos(-maxfreq),ypos(mindb));pen(up);
pento(xpos(maxfreq),ypos(maxdb));pen(down);
pento(xpos(maxfreq),ypos(mindb));pen(up);
f[1]:=10;f[2]:=20;f[3]:=50;f[4]:=1.0E2;
f[5]:=2.0E2;f[6]:=5.0E2;f[7]:=1.0E3;f[8]:=2.0E3;
f[9]:=5.0E3;f[10]:=1.0E4;f[11]:=2.0E4;f[12]:=5.0E4;
f[13]:=1.0E5;
FOR i:=1 TO 13 DO 
pento(xlpos(f[i]),ypos(maxdb));pen(down);
pento(xlpos(f[i]),ypos(mindb));pen(up);
IF ((i-1) MOD 3)=0 THEN
writefreq(f[i]);END; END; 
pen(up);pento(xpos(0-1.8E4),ypos(mindb-13));
outfile.WriteString('(frequency Hz) show');
outfile.WriteLn;
width(1.0); lblgain;
END loggrid;

PROCEDURE pltlogf(VAR fr:cmpt.fltrrec);
VAR i:LONGINT;f:REAL;prod:cx.complex;
first:BOOLEAN;
BEGIN
clipgrid; first:=TRUE;
FOR i:=1 TO npts DO
f:=fl*rm.exp(i*rm.ln(fmax/fl)/npts);
cmpt.pzgain(t,f,fr,prod);
pento(xlpos(f),ypos(db(cx.abs(prod),1.0)));
IF first THEN first:=FALSE;pen(down);END;END;
pen(up);width(2.8);
END pltlogf;

PROCEDURE pltlogx(VAR fr:cmpt.fltrrec);
VAR i:LONGINT;f:REAL;
first:BOOLEAN;y:cx.complex;
BEGIN
first:=TRUE;
FOR i:=1 TO npts DO
f:=fl*rm.exp(i*rm.ln(fmax/fl)/npts);
cmpt.pzxf(fr,f,y);
pento(xlpos(f),ypos(db(cx.abs(y),1.0)));
IF first THEN first:=FALSE;pen(down);END;
END;pen(up);dash(2.8);width(2.8);
END pltlogx;

BEGIN
loggrid; pltlogf(fr); pltlogx(fr); 
END pltlogproc;

PROCEDURE plttimeproc(fr:cmpt.fltrrec;choice:LONGINT);
VAR nsec,samples,left,right,dist:REAL;

PROCEDURE xcoord(i:LONGINT):REAL;
BEGIN
RETURN left+dist*(i/samples);
END xcoord;

PROCEDURE ycoord(y:REAL):REAL;
BEGIN
RETURN ypos(20*y-30);
END ycoord;

PROCEDURE timegrid;
CONST xofst=15.0;yofst=0.1;
VAR tic,itic:REAL;i,ntic,ix:LONGINT;
str:ARRAY 40 OF CHAR;
BEGIN
IF 
nsec<1.1E-4 THEN tic:=1.0E-5 
ELSIF 
nsec<1.1E-3 THEN tic:=1.0E-4 
ELSIF 
nsec<1.1E-2 THEN tic:=1.0E-3
ELSIF
nsec<1.1E-1 THEN tic:=1.0E-2
ELSE Out.String('duration too long');Out.Ln;HALT(1); END;
ntic:=ENTIER(nsec/tic);
itic:=2*fmax*tic;
pento(left,ycoord(-1));pen(down);
pento(right,ycoord(-1));
pento(right,ycoord(2));
pento(left,ycoord(2));
pento(left,ycoord(-1));pen(up);
pento(left,ycoord(1));pen(down);
pento(right,ycoord(1));pen(up);
pento(left,ycoord(0));pen(down);
pento(right,ycoord(0));pen(up);
pento(left-xofst,ycoord(2-yofst));
outfile.WriteString('(2) show');outfile.WriteLn;
pento(left-xofst,ycoord(1-yofst));
outfile.WriteString('(1) show');outfile.WriteLn;
pento(left-xofst,ycoord(0-yofst));
outfile.WriteString('(0) show');outfile.WriteLn;
pento(left-xofst,ycoord(-1-yofst));
outfile.WriteString('(-1) show');outfile.WriteLn;
FOR i:=1 TO ntic DO ix:=ENTIER(i*itic);
pento(xcoord(ix),ycoord(-1));pen(down);
pento(xcoord(ix),ycoord(2));pen(up);END;
pento(xpos(-4.0E4),ypos(mindb-13));
outfile.WriteString('(each tic=) show');
outfile.WriteLn;
outfile.WriteString('(');
RealStr.RealToFloat(tic,2,str);
outfile.WriteString(str);
outfile.WriteString(') show');outfile.WriteLn;
outfile.WriteString('( seconds) show');
outfile.WriteLn;
width(1);
END timegrid;

PROCEDURE pltstep;
VAR i,k:LONGINT;xk,yk:cx.complex;
BEGIN
FOR k:=1 TO fr.n DO fr.p[k].first:=TRUE;
fr.p[k].ykm1:=cx.zero;END;
FOR k:=1 TO fr.nz DO fr.z[k].first:=TRUE;
fr.z[k].xkm1:=cx.zero;END;
clipgrid;
FOR i:=1 TO ENTIER(samples) DO
xk:=cx.one;
cmpt.pzfltr(t,xk,yk,fr);
ar[i]:=yk;END;
pento(xcoord(0),ycoord(0));pen(down);
FOR i:=1 TO ENTIER(samples) DO
pento(xcoord(i),ycoord(ar[i].r));END;
pen(up);width(2.8);
pento(xcoord(0),ycoord(0));pen(down);
FOR i:=1 TO ENTIER(samples) DO
pento(xcoord(i),ycoord(ar[i].x));END;
pen(up);dash(2.8);width(2.8); 
END pltstep;

PROCEDURE pltmpls;
VAR i,k:LONGINT;xk,yk,mplscx:cx.complex;
mplsx:REAL;
BEGIN
Out.String('enter impulse scale factor');Out.Ln;
In.Real(mplsx);cx.rmult(mplsx,cx.one,mplscx);
FOR k:=1 TO fr.n DO fr.p[k].first:=TRUE;
fr.p[k].ykm1:=cx.zero;END;
FOR k:=1 TO fr.nz DO fr.z[k].first:=TRUE;
fr.z[k].xkm1:=cx.zero;END;
clipgrid;
FOR i:=1 TO ENTIER(samples) DO
IF i=1 THEN xk:=mplscx ELSE xk:=cx.zero END;
cmpt.pzfltr(t,xk,yk,fr);
ar[i]:=yk;END;
pento(xcoord(0),ycoord(0));pen(down);
FOR i:=1 TO ENTIER(samples) DO
pento(xcoord(i),ycoord(ar[i].r));END;
pen(up);width(2.8);
pento(xcoord(0),ycoord(0));pen(down);
FOR i:=1 TO ENTIER(samples) DO
pento(xcoord(i),ycoord(ar[i].x));END;
pen(up);dash(2.8);width(2.8); 
END pltmpls;

BEGIN
initgraf; initfont; scrfrac:=1.0;
maxfreq:=fmax*scrfrac;t:=1.0/(2.0*maxfreq);
left:=xpos(-maxfreq);right:=xpos(maxfreq);
dist:=right-left;
Out.String('enter plot duration seconds');Out.Ln;
Out.String('maximum duration 2E-2');Out.Ln;
In.Real(nsec);samples:=nsec*2*fmax;
IF choice=step THEN timegrid;pltstep;
ELSIF choice=mpls THEN timegrid;pltmpls;END;
END plttimeproc;

PROCEDURE printpz(fr:cmpt.fltrrec);
VAR resw:Msg.Msg;outvar:Files.File;
outfile:TextRider.Writer;
i:LONGINT;
BEGIN
outvar:=Files.New('temp6',{Files.write},resw);
outfile:=TextRider.ConnectWriter(outvar);
outfile.WriteString('fr.f1=');
outfile.WriteReal(fr.f1/(2*pi),20,0);outfile.WriteLn;
fr.gain:=1.0;cmpt.fltrgain(t,fr);
outfile.WriteString('fr.gain=');
outfile.WriteReal(fr.gain,20,0);outfile.WriteLn;
outfile.WriteString('n nz nzo =');
outfile.WriteLInt(fr.n,3);
outfile.WriteLInt(fr.nz,3);
outfile.WriteLInt(fr.nzo,3);outfile.WriteLn;
FOR i:=1 TO fr.n DO 
outfile.WriteString('pole ');
outfile.WriteLInt(i,3);
outfile.WriteReal(fr.p[i].sigma,20,0);
outfile.WriteReal(fr.p[i].omega,20,0);
outfile.WriteLn;END;
FOR i:=1 TO fr.nz DO 
outfile.WriteString('zero');
outfile.WriteLInt(i,3);
outfile.WriteReal(fr.z[i].sigma,20,0);
outfile.WriteReal(fr.z[i].omega,20,0);
outfile.WriteLn; END;
outvar.Close;
Out.String('output in temp6');Out.Ln;
END printpz;

PROCEDURE lp2bpproc(VAR fr:cmpt.fltrrec);
VAR fm:REAL;choice:LONGINT;
BEGIN
Out.String('enter geometric mean frequency of bandpass');
Out.Ln;
Out.String('and 1 for plain or 2 for *2pi');
Out.Ln;
In.Real(fm);In.LongInt(choice);
IF choice=2 THEN fm:=fm*2*pi;END;
cmpt.lp2bp(fm,fr);
END lp2bpproc;

PROCEDURE lp2hpproc(VAR fr:cmpt.fltrrec);
VAR fm,fn:REAL;choice:LONGINT;
BEGIN
Out.String('enter reflection frequency of highpass,');
Out.Ln;
Out.String('f1 normalization frequency');
Out.Ln;
Out.String('and 1 for plain or 2 for *2pi');
Out.Ln;
Out.String('all on one line: fm fn choice');
Out.Ln;
In.Real(fm);In.Real(fn);In.LongInt(choice);
IF choice=2 THEN fm:=fm*2*pi;fn:=fn*2*pi;END;
fr.f1:=fn;
cmpt.lp2hp(fm,fr);
END lp2hpproc;

PROCEDURE menu1;
BEGIN
gridproc(1);
Out.String('enter graf type:');Out.Ln;
Out.String('pgain zgain notch btwth iffltr');
Out.String(' cheby chhp elipt ');Out.Ln;
determine(command); 
CASE command OF 
er:
|btwth:btwthproc|elipt:eliptproc
|cheby:chebyproc|chhp:chhpproc
|pgain:pgainproc|iffltr:ifproc
|zgain:zgainproc|notch:notchproc
(*|hpf:hpfproc*) END;
showgraf;
END menu1;

PROCEDURE menu2;
BEGIN
REPEAT
Out.String('enter q2 to quit this menu');Out.Ln;
Out.String('enter operation:');Out.Ln;
Out.String('new scale lp2bp lp2hp setf1'); Out.Ln;
Out.String('pltpln pltlin pltlog prtpz step mpls'); Out.Ln;
determine(command);
CASE command OF
er:;
|new:newpz(fr);
|scale:scaleproc(fr);
|lp2bp:lp2bpproc(fr)
|lp2hp:lp2hpproc(fr);
|setf1:setf1proc(fr);
|pltpln:pltpz(fr);showgraf;
|pltlin:pltlinproc(fr);showgraf;
|pltlog:pltlogproc(fr);showgraf;
|prtpz:printpz(fr)
|step:plttimeproc(fr,step);showgraf;
|mpls:plttimeproc(fr,mpls);showgraf;
|q2:;
END;
UNTIL command=q2;END menu2;

BEGIN (*grafproc*)
maxdb:=+10.0;mindb:=maxdb-dbmax;
IF choice=1 THEN menu1 ELSE menu2 END;
END grafproc;

(*section 4. command interpreter*)

PROCEDURE initvar;
BEGIN
cmdara[q     ]:='q';
cmdara[er    ]:='er';
cmdara[graf1  ]:='graf1';
cmdara[graf2  ]:='graf2';
cmdara[save1 ]:='save1';
cmdara[save6 ]:='save6';
cmdara[plot  ]:='plot';
cmdara[grid  ]:='grid';
cmdara[btwth ]:='btwth';
cmdara[notch ]:='notch';
cmdara[elipt ]:='elipt';
cmdara[cheby ]:='cheby';
cmdara[chhp  ]:='chhp';
cmdara[pgain ]:='pgain';
cmdara[zgain ]:='zgain';
cmdara[iffltr]:='iffltr';
cmdara[rp    ]:='rp';
cmdara[rz    ]:='rz';
cmdara[cp    ]:='cp';
cmdara[cz    ]:='cz';
cmdara[zo    ]:='zo';
cmdara[scale ]:='scale';
cmdara[pltlin]:='pltlin';
cmdara[new   ]:='new';
cmdara[pltlog]:='pltlog';
cmdara[lp2bp ]:='lp2bp';
cmdara[q3    ]:='q3';
cmdara[q2    ]:='q2';
cmdara[pltpln]:='pltpln';
cmdara[prtpz ]:='prtpz';
cmdara[lp2hp ]:='lp2hp';
cmdara[setf1 ]:='setf1';
cmdara[step  ]:='step';
cmdara[mpls  ]:='mpls';
cmdara[hpf   ]:='hpf';
(*insert additional commands here*)
cmdara[last  ]:='last';
pen(up);first:=TRUE;
END initvar;

PROCEDURE readwd(VAR c:str);
BEGIN
In.Identifier(c);
IF In.Done() & (c # "") THEN
ELSE HALT(1) END;
END readwd;

PROCEDURE determine(VAR cmdvar:LONGINT);
CONST bell=7;
VAR cmi:LONGINT;strvar:str;
BEGIN
Out.String('enter command');Out.Ln;
cmdvar:=er;
readwd(strvar);
FOR cmi:=q TO last DO
IF (cmdara[cmi]=strvar)THEN cmdvar:=cmi;END;END;
IF cmdvar=er THEN
Out.Char(CHR(bell));Out.String('<--<< error');Out.Ln;END;
END determine;

BEGIN (*main*)
initvar; penstate:=up;
REPEAT
Out.String('enter graf1, graf2, save6, save1 or plot');Out.Ln;
Out.String('if graf1 or graf2, enter quit when finished');Out.Ln;
Out.String('viewing the graf, enter q to exit program');Out.Ln;
determine(command);
CASE command OF
er:
|graf1:grafproc(1)(*computes fft graph,displays on screen.
                type quit to terminate ghostscript*)
|graf2:grafproc(2)(*manual pz entry, gain plots*)
|save6:save6proc(*accumulates 6 graphs*)
|save1:save1proc(*saves one graph to disk*)
|plot:plotproc(*plots to laser printer*)
|q:;
END(*case*);
UNTIL command=q;
END cmptst.


up one level