If the lines of text are too long you can fix the problem with these instructions.
The filter technique presented here offers continuous processing, not batch processing like fft filtering. If your application is only to replace unfiltered stored data with filtered stored data, 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. In the method presented here, during the same sampling period that sample n+1 is read in, completely filtered sample n is written out.
If a digital filter is used to filter a baseband signal with no carrier frequency, the signal is represented by a sequence of real numbers. If the signal represented has a carrier frequency, the carrier frequency will not be represented in the sequence of numbers, but the amplitude and phase of the modulation of the carrier frequency will be represented. This would require a sequence of complex numbers. The method presented here can handle real or complex numbers.
Some people complain that the method presented here has unwanted noise due to roundoff error. That is because they did not know how to implement this method properly to avoid roundoff noise. Proper implementation is explained here.
If the sequence of numbers you wish to filter is complex, the use of this technique is very simple. If the sequence is real, 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.
This document shows example filter curves produced using the technique described in the article at this site "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.
If software filtering is to be used, whether in a hardware application or a purely software application, the technique presented in this article will probably be easier to understand and to program than the conventional method usually found in textbooks. The technique presented here uses complex arithmetic internal to the filter whether the sequence of numbers is a baseband real sequence or a complex sequence representing amplitude and phase of a carrier signal.
The conventional method is rigorously derived. My method is based on a hunch supported weakly by a reasonableness argument and strongly by impressive accuracy in practice. Many books describe the conventional method. "digital filters" by R. W. Hamming, "digital filters" by Andreas Antoniou, "digital signal processing" by Steven W. Smith, and "digital filters for everyone" by Rusty Allred. None of these authors wrote a program to display a large variety of frequency domain and time domain graphic examples the way I have in this article. Perhaps the conventional method presents some difficulty that would make such a program impractical.
See also the derivation of a pole here and the derivation of the zero from the pole here.
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.
The equations are normalized in a way chosen for convenience. The pole has unity gain only at the frequency of the pole. The zero approaches unity gain at frequencies far from the frequency of the zero.
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 is listed at the end of this article.
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:
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.
The labeling of our plots of the complex plane is very unconventional, perhaps even dishonest. If a pole is near the imaginary axis, that will result in a peak at some frequency f. But the imaginary coordinate of the pole on the complex plane will be near omega=2*pi*f. But we are only interested in the final result, the frequency response that peaks at f. So the tic marks on the plots of the complex plane do not represent omega, but omega/2*pi, or f. The tic mark at 2*pi*f is given the false label of f. Thus if we put a pole near a tic mark our response will peak at the corresponding f. This is confusing, but I find it convenient.
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.
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:
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:
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".
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:
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:
The log plot is a lowpass filter which rolls off at 1E3 and goes down at 12dB per octave:
Next we add a zero at the origin:
The result is a bandpass filter centered at 1E3 with a 6dB per octave slope on either side:
Finally, we add a second zero at the origin, indicated by the number "2" nearby:
The log plot shows a highpass filter at 1E3 with a 12dB per octave slope at low frequencies:
Now we show a bandpass filter with the real component of the pole locations only 0.03 of the imaginary part:
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:
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.
We show an elliptical filter. The normalized pole-zero plot shows 5 poles and 4 zeroes:
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:
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 bandwidth measured from zero to positive frequencies. Similarly, if a bandpass filter has a steady sine wave 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:
The corresponding time domain step response:
Next, poles at 12 degrees from the imaginary axis:
The corresponding time domain step response:
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 elliptic 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.
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:
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.
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.
One kind of multi pole 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.
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:
We look at the same four pole butterworth filter with half the original simulation bandwidth:
We end our bandwidth investigation by showing the same filter at one fourth the original simulation bandwidth:
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?
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:
Linear frequency response:
Log frequency response:
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:
Linear frequency response:
Log frequency response:
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:
Linear frequency response:
Log frequency response:
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.
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. In nearly every mathematical or geometrical application every angle should be expressed as positive to avoid complexity or wrong answers. The correct form of square root is given in the source code at the end of this article. It is the procedure sqrt in the module cxarith.
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.
Linear frequency response:
Log frequency response:
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:
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 elliptical 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:
The pole-zero plot is not clear because of many poles and zeroes almost lumped together near the origin:
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.
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:
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:
The pole locations for bessel filters with various numbers of poles can be found on the internet. But be careful, sometimes the poles listed are normalized to group delay rolloff frequency and sometimes to amplitude rolloff frequency. The poles shown here are normalized to amplitude rolloff frequency. The normalized 5 pole bessel filter used here is shown below. The normalized poles are (-1.5023,0), (-1.308,0.7179), (-0.9576,1.4711). Multiply by 1E3 and 2*pi for a 1kHz 3dB point.
The impulse response of the bessel filter may not seem ideal because it is not quite symmetrical and has a very slight undershoot after the pulse. This could probably be improved by addition of well placed allpass elements. Allpass elements are frequency dependent delay elements with flat amplitude frequency response. A series of them up the imaginary axis can simulate a delay line. An allpass pole-zero pattern is a complex pair of poles with negative real coordinate paired with a pair of zeroes with positive real coordinate of the same magnitude and the same imaginary coordinate. The bessel filter has flat group delay within its passband, but not beyond its passband. Allpass elements could extend the flat group delay beyond the passband.
A very simple smoothing filter that has an impulse response with negligible overshoot or undershoot is poles at (-1,0), (-1,1),(-1,2). If this is scaled by 1E3 and multiplied by 2*pi, the 3dB point is approximately 1000Hz.
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 received. This scheme implements a delay line many samples long. This is called a FIFO, or first in, first out scheme.
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.
Out next example is certainly NOT a smoothing filter. Since we are on the subject of impulse responses, it may be of passing interest to see a pole-zero pattern that produces a positive impulse followed by a similar, but not identical, negative impulse. It serves as a rough approximation to a single cycle sine wave impulse response. Each tic mark is 1000. The pole pattern is a 5 pole filter with a zero at the origin. The normalized poles are at (-0.6,0), (-0.5,1), (-0.4,2). Multiply by 1E3 and 2*pi.
As of Jan 2022 the part of the program that plots isolated poles and zeroes has been modified so that the simulation, the transfer function and the gain formula can be plotted separately without the necessity of recompiling the program.
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:
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 :
Next, a single pole at -5E1, 4E4, :
Now a zero at -5E1 and 4E4:
In the case of isolated zeroes the simulated zero matches the gain function of the zero at the frequency of the zero. But the computed transfer function of the zero does not. This is because the zero is normalized to unity gain at a point on the opposite side of the frequency circle from the frequency of the zero. Since the shape of the transfer function of the zero is different from the shape of the simulated zero, the depth of the transfer function will be lower than the depth of the simulation at the frequency of the zero. This will not produce such a large effect when poles and zeroes are combined to produce a filter. In that case each of the poles and zeroes will be normalized to unity gain at the frequency where the filter will be normalized to unity gain.
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:
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:
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:
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:
Now this ability of the pole and zero to cancel each other out seems astonishing since the pole is an infinite impulse response filter and the zero is a finite impulse response filter.
First we explain phasors. In the following the symbol "*" represents multiplication. If the mathematical constant "e" is raised to an exponent "x" we will represent this as "exp(x)". A carrier signal cos(omega*t) is equal to 0.5*(exp(j*omega*t) + exp(-j*omega*t)). To greatly reduce the simulation bandwidth and therefore the sample rate necessary we only use the term with the positive exponent to represent the signal. This term exp(j*omega*t) equals cos(omega*t) + j*sin(omega*t). This complex number computed with cos and sin functions is what we mean by a phasor.
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:
Next we make the I.F. frequency 2E5, ten times the filter bandwidth:
Finally, we make the I.F. frequency 2E6, 100 times the filter bandwidth:
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:
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 immediately 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.
One place where this filter technique is particularly useful is in block diagram simulation of a super heterodyne system with radio frequency signals, or RF signals with a carrier frequency, intermediate frequency signals, or IF signals with a carrier, and baseband signals with no carrier. Complex samples will be used at RF and IF. Real signals will be used at baseband. If the signal is below the nominal carrier frequency the samples will be at negative frequencies. If the signal is above the nominal carrier frequecy the samples will be at positive frequency. If the negative frequency representation of the local oscillator (LO) signal is multiplied by the positive frequency representation of the RF signal, the product will be an IF signal that would result from single sideband mixing. Simple to do in software, complicated to do in hardware. exp(a)*exp(b) = exp(a+b). exp(j*(omega*t+phi))*exp(-j*omegalo*t) = exp(j*((omega-omegalo)*t+phi)). Thus phase shift is preserved through mixing.
Most of the block diagram components will be trivial to implement. One that might not be would be an IF phase detector intended to provide a signal proportional to a small phase difference between two IF signals of the same frequency. The phase detector should produce a positive or negative real output proportional to positive or negative phase difference of one signal relative to the other signal.
Perhaps the simplest way to implement the phase detector is based on the cross product of two three dimensional vectors. The tails of the vectors are at the origin of the three dimensional coordinate system and the points of the vectors are at two different points away from the origin. But two points on the complex plane are not two three dimensional vectors. We will cheat and use the complex plane to represent two of three dimensions. Then the cross product of two points on the complex plane would represent the magnitude of a third vector perpendicular to the complex plane. But instead of interpreting it as the magnitude of a third vector, we will interpret it as a phase difference. If the two points on the complex plane are in the same direction from the origin, the magnitude of the third vector would be zero. If one of the points is the reference point and the other differs by some angle plus or minus from the reference point, the cross product of the two points will produce a positive or negative result. This is the kind of behavior we want for our phase detector
Let our two IF signals be s1 and s2, with real part s1.r and imaginary part s1.x. Then the output of our phase detector is y:=(s1.r*s2.x) - (s1.x*s2.r). But y will be a linear function of phase difference only for small phase differences. As the phase difference moves from zero to 360 degrees y traces a single cycle of a sine wave.
Let us think about what we have accomplished. Two signals at the IF carrier frequency will be represented as being at zero frequency. But we will still be able to measure small phase differences between them because the samples of the signals are complex numbers. This would not be possible with samples that were real numbers. This shows the value of the kind of filters presented in this article instead of conventional filters that only work with real numbers.
The software that produced the curves shown in this article is listed below. Listed first is the instructions for compiling and operating the software, then the software source code itself.
To compile and execute this software you must have ubuntu linux operating system, and you must install the oo2c ooberon compiler as instructed in the article how to program a computer elsewhere at this website. Be sure to follow the directions to create the scripts ooc and oocm. To use this program you must save this article to disk. But some browsers to strange things to html files when they save them to disk. If you followed the instructions for installing the compiler you will have a directory ~/wkspc/scratch. In that directory "wget waltzballs.org/other/engr/pzex.html". That will get you a clean copy of the article. Then in that directory "unhtml pzex.html > pzex.text". In pzex.text the source code of the program will be in compileable form. Using the vi editor save each module separately in the directory ~/wkspc/prog/src. Thus module cxarith would be saved as cxarith.Mod, fft as fft.Mod, cmpt as cmpt.Mod and cmptst as cmptst.Mod. Then in the directory above, you could follow the directions in the programming article for separate compilation of the modules. Then in the directory above the source code ~/wkspc/prog, issue "bin/cmptst" to run the program.
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 = rm.pi; 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/. *)
(*To understand this program read the article
"complex digital filters using isolated poles and zeroes".*)
IMPORT cx:=cxarith,rm:=RealMath,Out;
CONST maxpoles=20; pi=rm.pi;
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. f1 is the frequency where the filter
will be normalized to unity gain.*)
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);
(*w is the frequency where the filter will be normalized
to unity gain. p.omega is the frequency of the signal
appled to the pole*)
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;
(*zero with gain applied before subtract makes no
numerical difference)
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);
z.xkm1:=z.xk;
cx.rmult(1.0/z.gain,z.xk,z.xk);
cx.rmult(1.0/z.gain,t1,t1);
cx.sub(z.xk,t1,z.yk);
cx.mult(z.norm,z.yk,z.yk);
(*note that gain is not applied to xkm1*)
END zero;*)
(*zero with 64bit subtract makes no difference
PROCEDURE zero*(t,w:REAL;VAR z:zerorec;VAR fr:fltrrec);
VAR rotate,t1,onemeps:cx.complex;
dzxr,dzxx,dt1r,dt1x:LONGREAL;
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);
dzxr:=LONG(z.xk.r);dzxx:=(z.xk.x);
dt1r:=LONG(t1.r);dt1x:=LONG(t1.x);
z.yk.r:=SHORT(dzxr-dt1r);z.yk.x:=SHORT(dzxx-dt1x);
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)*pi)/(2*n));
omega:=
cosh((1.0/n)*arcsinh(1.0/eps))*rm.cos(((2*k-1)*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
pzfltr(t,pzxk,pzyk,fr);
END onep;
PROCEDURE onez*(t:REAL;pzxk:cx.complex;
VAR pzyk:cx.complex;VAR fr:fltrrec);
BEGIN
(*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;
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=rm.pi;strlngth=7;
q=1;er=2;graf1=3;save6=4;plot6=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;mpls=34;add=35;plot1=36;last=37;
(*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;str1:str;
(*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;
outfile.WriteString('%%Page: 1 1');outfile.WriteLn;
invar.Close;
outvar.Close;
Out.String('saved in temp5.eps');Out.Ln;
END save1proc;
PROCEDURE plot1proc;
BEGIN
irts:=ProcessManagement.system("lpr temp2");
Out.String('sent to printer');Out.Ln;
END plot1proc;
PROCEDURE plot6proc;
(*writes plot buffer to the laser printer.
works only for output of save6proc.*)
VAR resr,resw:Msg.Msg;invar,outvar:Files.File;
infile:TextRider.Reader;outfile:TextRider.Writer;
str:ARRAY 40 OF CHAR;
BEGIN
Out.String('works only for output of save6proc');
Out.Ln;
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;
outfile.WriteString('%%Page: 1 1');
outfile.WriteLn;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 printer and saved in temp4');Out.Ln;
first:=TRUE;
END plot6proc;
(*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('/Latin-Modern-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);
(*compute impulse and fft and plot*)
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;
(*the following is a procedure variable that can
become any of several different procedures*)
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);
(*compute theoretical transfer function and plot*)
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;
VAR i:LONGINT;f,w:REAL; first:BOOLEAN;
BEGIN
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; fr.nz:=0;first:=TRUE;fr.firstb:=TRUE;
fr.p[1].first:=TRUE;
fr.f1:=2*pi*fr.fc;fr.fc:=2*pi*fr.fc;
fr.p[1].sigma:=2*pi*fr.bw;
fr.p[1].omega:=fr.fc;
Out.String('plot fft output? y/n');Out.Ln;
In.Identifier(str1);
IF ((str1='y')OR(str1='Y'))THEN
pltfltr(fr,cmpt.onep);END(*IF*);
Out.String('plot transfer funcion? y/n');Out.Ln;
In.Identifier(str1);
IF ((str1='y')OR(str1='Y'))THEN
pltxfr(fr);END(*IF*);
Out.String('plot gain formula? y/n');Out.Ln;
In.Identifier(str1);
IF ((str1='y')OR(str1='Y'))THEN
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(*IF*);
END pgainproc;
PROCEDURE zgainproc;
VAR i:LONGINT;f,w:REAL; first:BOOLEAN;
BEGIN
fr.nzo:=0;
Out.String('enter frequencies, not radian frequencies');Out.Ln;
Out.String('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.nz:=1;fr.n:=0;first:=TRUE;fr.firstb:=TRUE;
fr.z[1].first:=TRUE;
IF fr.fc < 0.0 THEN fr.f1:=fr.fc+maxfreq;
ELSE fr.f1:=fr.fc-maxfreq;END;
fr.f1:=2*pi*fr.f1;fr.fc:=2*pi*fr.fc;
fr.z[1].sigma:=2*pi*fr.bw;
fr.z[1].omega:=fr.fc;
Out.String('plot fft output? y/n');Out.Ln;
In.Identifier(str1);
IF ((str1='y')OR(str1='Y'))THEN
pltfltr(fr,cmpt.onez);END(*IF*);
Out.String('plot transfer funcion? y/n');Out.Ln;
In.Identifier(str1);
IF ((str1='y')OR(str1='Y'))THEN
pltxfr(fr);END(*IF*);
Out.String('plot gain formula? y/n');Out.Ln;
In.Identifier(str1);
IF ((str1='y')OR(str1='Y'))THEN
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(*IF*);
END zgainproc;
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 addpz(VAR fr:cmpt.fltrrec);
PROCEDURE getdata;
VAR i,j:LONGINT;
BEGIN
i:=fr.n;j:=fr.nz;
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 addpz;
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;
Out.String('probably from 100 to 1000');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;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 add 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);
|add:addpz(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[plot1 ]:='plot1';
cmdara[plot6 ]:='plot6';
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[add ]:='add';
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';
(*insert additional commands here*)
cmdara[last ]:='last';
pen(up);first:=TRUE;
END initvar;
PROCEDURE determine(VAR cmdvar:LONGINT);
CONST bell=7;
VAR cmi:LONGINT;strvar:str;
BEGIN
Out.String('enter command');Out.Ln;
In.Identifier(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, plot1 or plot6');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*)
|plot1:plot1proc(*plots to laser printer*)
|plot6:plot6proc(*plots to laser printer*)
|q:;
END(*case*);
UNTIL command=q;
END cmptst.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.