# IIR filter analysis using VHDL.Allpass, multiple delay, and masking filters

**A.Sergiyenko **

National Technical University of Ukraine “КРІ”,

Email: aser@comsys.kpi.ua

*z-transform and transfer function *

z-transform is widely used in DSP due to the fact, that the signal *x(n) = z ^{n}* is the

*of any linear system which is invariant to the signal shift. The eigenfunction means that when such a function is processed by the system then the output signal has the same form but is delayed, and its magnitude is different. We can show that*

**eigenfunction***z*is the eigenfunction of the system which is described by some impulse response

^{n}*h(k)*. Consider the signal

*z*is inputted to the system. Then due to the properties of the linear shift invariant system the output signal is

^{n}(1)

We see that the function *z ^{n}* passes through the system invariably. Due to the formula (1), the

*of the system is*

**transfer function**(2)

The function *H(z)* is the complex function with the magnitude |*H(z)*| and phase *arg(H(z))*. The function *H(z)* has two sets of special points. The first set contains the points * q _{1},…,q_{M}*, for that |

*H(q*| = 0. Therefore, these points are called

_{і})*of the function*

**zeros***H(z)*. The second set contains the points

*r*, for that |

_{1},…,r_{N}*H(r*| = ∞, and therefore, they are called

_{і})*of the function*

**poles**,*H(z)*.

To find the poles and zeros of

*H(z)*, the sum (2) is transferred to the fraction with the positive powers of

*z*. Such a fraction is the following

(3)

The numerator and denominator of the fraction (3) are polynomials with the positive powers of *z*, and * a _{0} = 1*. The multiplier

*z*represents the delay to

^{-L}*L*cycles, and can be taken off. The roots of these polynomials are equal to

*q*and

_{1},…,q_{M}*r*, and are called as

_{1},…,r_{N},*and*

**zeros***, respectively. The function*

**poles***H(z)*can be represented by the sets of zeros and poles as the following fractions

(4)

where *K* is the transfer ratio or the amplification factor of the system.

*Poles and zeros of the transfer function *

Poles and zeros are usually represented in the complex plane as circles and crosses respectively. Fig.1 illustrates an example of pole and zero allocation of some low pass filter. A root point is equal to a couple of real and imaginary numbers: *r = r _{R} + jr_{I} *. It is convenient to represent it in the polar coordinates by the magnitude |

*r*|, and phase

*ω*:

*r =*|

*r*|

*e*. When we consider the discrete signal computing the phase

^{-jω}*ω*represents some frequency of the transfer function. The pole magnitude is not higher than a one for real filters, i.e. |

*r*|

*. The number of poles with the nonzero imaginary parts is always even.*

__<__1*Magnitude-frequency response*

The magnitude of the transfer function |*H(z)*| can be represented in the three dimensional space as the endless elastic diaphragm, which is pulled up the complex plane at the mark one. The diaphragm points, which represent zeros, are pulled down to the complex plane, and points, which represent poles, are pulled up to infinity. If we draw the concentric circles in the diaphragm, then after its stretching we derive the figure like in the Fig.2. The bold line at the surface |*H(z)*| (Fig. 2) is mapped to the complex plane as a circle with the radius of one. This line is derived by

(5)

When we input the sequence * e ^{-jω} *in the formula

*H(z)*, then we derive the discrete Fourier transform, and when we input it into formulas (2) or (3) we get the transfer function

|*H(ω)*| = |*H(e ^{-jω})*|, (6)

This transfer function means the transfer ratio dependence on the input signal frequency *ω*. The function |*H(ω)*| is called the * magnitude frequency response* (MFR) of the system. We can imagine that the bold curve in the Fig.2 is depicted in the cylinder, which is rotated on the axis o|

*H(z)*|. When we untwist this curve and map it in the plane we derive the diagram, which is formed by the unlimited periodical curve with the period of

*2π*, like the curve in the Fig. 3.

The frequency axis of the MFR is normed by dividing to *2π* for the usability, so as a one in this axis represents the sampling frequency.

Similarly the * phase frequency responce* (PFR) of the system is derived:

(7)

It is adopted for the digital filters, that they have the level of MFR higher than √0,5 ≈ 0,7 in the * pass band*. This means that in the pass band the power of the filtered signal is no less than a half of the nominal power. The respective frequency bound is called as the

*(*

**cutoff frequency***ω*in Fig. 3). The maximum level of the filtered signal in the

_{c}*is called the*

**stop band***of the filter (*

**suppression level***H*in Fig.3).

_{s}The DSP systems are usually the chains of sequentially connected units.

Therefore, it is useful to represent MFR with the logarithmical scale, i.e.

*H(ω) = 20*lg|H(e ^{-jωn})|*

Then the resulting MFR of the system is equal to the sum of MFRs of its parts. The level of such MFR is measured in decibels. For example, the level * H(ω)* = 1 represents 0 db, level 0,7 means -3 db., level 0,5 means -6 db, and so on.

*Allpass filters*

Further for the instance we will investigate the transfer function of the form

(8)

and other functions based on it. Function * H _{1}(z)* represents the

*. The poles and zeros of the function*

**allpass filter***H*are illustrated by Fig. 4. The poles and zeros of the function (8) have the equal angles

_{1}(z)*ω*. If length of the pole vector is equal to

*r*, then the length of the zero vector is equal to 1/

*r*. As a result, the zeros fully compensate the poles, and MFR is equal to 1 in the whole frequency band. For this property such a filter is called as an allpass one.

The factors *a* and *b* determine the zero and pole coordinates due to the formulas

*a = – cos(ω); b = r ^{2}.* (9)

The allpass filter has not effect on the signal magnitude, but it distorts cardinally its phase. Moreover, when the factor *b* approaches to a one, then the phase distortion at the frequency *ω* is sharper. If we connect two allpass filters with the characteristics *H _{1}(z), H_{2}(z)* in parallel, i.e.

*H _{3}(z) =(H_{1}(z) + H_{2}(z))/2,* (10)

then the poles and zeros in the resulting characteristic * H _{3}(z)* stop to compensate each other. Both cannels

*H*,

_{1}(z)*H*pass the input signal without its magnitude exchange. But on the frequency, where the phase difference is equal to

_{2}(z)*π, 3π,…,*, the output signals compensate each other and the resulting signal is zero. On the contrary, on the frequency, where the phase difference is equal to

*0,2π, 4π,…,*the resulting signal is doubled after addition. This causes the division to the factor 2 in the formula (10). As a result, the pass band and stop band are present in MFR |

*H(ω)*|.

*Multiple delay filters and masking filters*

One can increase all the delays in * k * times in the linear DSP system. Then instead of the transfer function *H(z)* the function with multiple delays *H(z ^{k})* is considered. For this feature such a filter is named as a multiple delay filter. The effect of delay multiplication consists in scaling MFR, and PFR diagrams in

*k*times along the axis of frequencies.

Let the coefficient

*k*is equal to 2, then in the diagram in the Fig. 3 the marks

*π, 2π*and so on are exchanged to marks

*π/2, π*etc. As a result, the number of both pass bands and stop bands is doubled.

Consider *H _{3}(z)* is a low pass filter, and we construct the multiple delay filter

*H*. Then its MFR |

_{4}(z) = H_{3}(z^{k})*H*| has

_{4}(ω)*k*pass bands. In this situation we can select a single working pass band using, so called,

*with the transfer function*

**masking filter***H*. Such a filter is connected to the multiple delay filter in sequence, i.e. the resulting transfer function is

_{M}(z)*H(z)*=

*H*(11)

_{4}(z)*H_{M}(z)MFR of the masking filter has the stop bands in that frequencies where the pass bands in MFR |*H _{4}(ω)*| are present, which we want to suppress.

*Some useful VHDL functions*

In this article we consider the computation of the transfer functions using VHDL. For this purpose we have to use the mathematical packages IEEE.MATH_REAL, and IEEE.MATH_COMPLEX. In the last package the constants, types, and functions are defined, which operate with the complex data.

Below we consider a set of developed functions, which are useful to calculate the transfer function. The calculation of the sequence * z ^{n}* is implemented using the following functions, which overload each other.

--complex vector Z with magnitude MAG, and phase j*fi function Z(mag,fi:real) return COMPLEX_POLAR is begin return mag*exp(COMPLEX_TO_POLAR(MATH_CBASE_J)*(fi)); end Z; --complex number Z with magnitude 1, and phase j*fi function Z(fi:real) return COMPLEX_POLAR is begin return exp(COMPLEX_TO_POLAR(MATH_CBASE_J)*(fi)); end Z;

The following function calculates * H _{1}(z) * due to the formula (8), but for positive degrees of

*z*, i.e. considering the equation(3). It is useful for the scanning of the space |

*H*| when we search for zeros and poles. By this process we sequentially set different values of the magnitude

_{1}(z)*mag*, аnd phase

*fi*of the complex vector

*z*. It is used for computations of MFR and PFR as well.

calculating (bZ^2+а(b+1)Z+1)/(Z^2+c(b+1)Z+b) function Allpass2(b,a: real; --coefficients mag, fi: real) -- magnitude, angle for z return COMPLEX_POLAR is variable tn,td: COMPLEX_POLAR; begin td:=COMPLEX_TO_POLAR(COMPLEX'(b,0.0)) +a*(b+1.0)*Z(mag,fi) + Z(mag**2,2.0*fi); tn:=COMPLEX_TO_POLAR(COMPLEX'(1.0,0.0)) +a*(b+1.0)*Z(mag,fi) + b*Z(mag**2,2.0*fi); return tn/td ; end Allpass2;

The following function calculates the formula * H _{1}(z^{2})* (8) for the negative degrees of

*z*with a single magnitude, i.e. it calculates the function

*H*, for

_{1}(ω)*ω = fi*. Therefore, it is used only for computing MFR, and PFR.

function Allpass2x2(b,c:real; fi:real) return COMPLEX_POLAR is variable tn,td: COMPLEX_POLAR; begin tn:=COMPLEX_TO_POLAR(COMPLEX'(b,0.0)) +c*(b+1.0)*Z(-2.0*fi) + Z(-4.0*fi); --numerator td:=COMPLEX_TO_POLAR(COMPLEX'(1.0,0.0)) +c*(b+1.0)*Z(-2.0*fi) + b*Z(-4.0*fi); --denominator return tn/td ; end Allpass2x2;

To scan the three dimensional space |*H _{1}(z)*| the following part of the VHDL program is used.

signal clk:bit; signal m,ph:real:=0.0; signal a,b,logm,phase,mag:real:=0.0; begin clk<=not clk after 5ns; --clock generator process(CLK) variable p,phas:real:=0.0; variable Hz: COMPLEX_POLAR; begin a<=-0.5; b<=0.64; -- parameters of H(z) if clk='1' and clk'event then phas:= phas+0.001; -- phase counter (frequency) p:=phas*MATH_PI*2.0; -- normed phase m<= trunc(phas)*0.1+0.1; -- magnitude of Z ph<=phas; -- phase signal end if; Hz:=Allpass2(b,a,m,p); --resulting H(z) mag<=abs(Hz); -- magnitude of H(z) phase<=Hz.ARG; -- phase of H(z) logm<=20.0*log10(abs(Hz)); -- H(z) in decibels end process;

*Investigation of the allpass filters using VHDL*

Let us investigate the function * H _{1}(z)* using the VHDL program, which outputs the diagrams |

*H*| for different values |

_{1}(z)*z*| to estimate the coordinates of poles and zeros of the function (7). Consider

*a = -0.5, b = 0.64*. Then the transfer function is

(12)

After running the described above VHDL program in the system Aldec ActiveHDL the diagram |*H _{1}(z)*| (Fig. 5) was outputted. In this diagram the curves |

*H*| are placed in a set of intervals, marked as

_{1}(m,ph)*m*, where

*m*= 0.1, 0.2,…, is magnitude,

*ph*is phase (frequency) of the vector

*z*. The vector

*z*makes a single turnover when

*ph*is exchanged from

*i*to

*i+1*, where

*i*is the number of the interval. The first diagram |

*H*| on the Fig. 5 is built for

_{1}(m,ph)*m*= 0.1,...0.9, and the second one is for

*m*= 1.0,...1.9. We can see that each interval represents an exclusive curve like ones in the Fig. 2. In the Fig.5 the diagram

*arg(H*for

_{1}(z))*m*= 0.6,...1.1 is shown as well.

The analysis of the diagrams in the Fig. 5 shows the following. We see, that the magnitude of * H _{1}(z)* at

*m*= 1 is really equal to a one. The poles of the function

*H*have the measured phases

_{1}(z)*± 2π*ph = ± 2π*0.1645 = ± 1.0336*, and they have the polar coordinates

*r*and

_{1}= 0.8 ∠1.0336*r*. The zeros of the function

_{2}= 0.8 ∠-1.0336*H*have the same phase and are equal to

_{1}(z)*q*. One can see that

_{1}= 1.3 ∠1.0336, q_{2}= 1.3∠-1.0336*-cos(1.0336) = -0.512 ≈ a, 1/0.8 = = 1.25 ≈ 1.3*, and

*0.8*, i.e. the formulas (9) are true.

^{2}= 0.64 = b

Fig. 5. Involute of the function |H_{1}(z)| (a) and of the function arg(H_{1}(z)) (b)

The analysis of the diagram * arg(H _{1}(z))* shows that for the angles

*ph*= 0 the phase is equal to 0, and for the gradual increase of

*ph*the resulting phase is monotonously decreased. But when

*ph*and norm of

*z*approach to poles the phase sharply approaches to

*±π*.

Let us analyze the function * H _{3}(z)* by

*H*. Here

_{2}(z) = z^{-1}*H*is the transfer function of another allpass filter, because |

_{2}(z) = z^{-1}*H*| = 1 for all

_{3}(z)*z*.

To compute the function |* H _{3}(ω)*| in the previous VHDL program one variable assignment operator is substituted to the following operator, which represents the formula (10) .

*Hz:=(Allpass2(b,a,1.0,p)+Z(-1.0*p))/2.0*(13)

Some correction of the factor

*b*is made, i.e. = 0.625. The resulting function is

(14)

One can show that the multiplication to the multiplier 0.625 can be implemented as a single addition of the multiplicand, which is shifted right to one and three bits respectively. Similarly the multiplication to a factor 0.8125 can be made. This feature simplifies substantially the hardware of the filter, which is represented by the formula (14).

After the running of this program the diagrams |*H _{3}(ω)*| and

*arg(H*are derived for the frequencies

_{3}(ω))*ω=(0, 2π)*or for the normed frequencies

*ph*= (0, 1), which areshown in the Fig.6. We can see that after addition of the functions

*H*and

_{1}(z)*H*the function of the low pass filter is obtained. Due to the diagram this filter has the cut-off frequency

_{2}(z)*f*, and the suppress level of 8.5 db.

_{c}= 0,156 ≈ arg(r_{1}) =1/(2π) 1,0336 = 0,164
Now we can investigate the function * H _{4}(z)*.

(15)

The diagram |*H _{4}(ω)*| (Fig.7,а) is built by the substitution of an operator to the following one as well as in the previous example.

*Hz:=(Allpass2x2(b,a,p)+Z(-2.0*p))/2.0;*

Fig. 7. Diagrams of the magnitude frequency responses |H_{4}(ω)| (a), |H_{M}(ω)| (b), |H(ω)| in the linear and logarithmical scales

We can see that the diagram in Fig.7,a represents the twofold repeated diagram in Fig 6, which is typical to the filters with multiple delays.

Тhe function of the masking filter

(16)

we can prove by the diagram |*H _{M}(ω)*| (Fig.7,b), which is computed using the following operator.

Hz:=(COMPLEX_TO_POLAR(COMPLEX'(1.0,0.0))+Z(-p)) *(COMPLEX_TO_POLAR(COMPLEX'(1.0,0.0))-Z(-7.0*p))/ (COMPLEX_TO_POLAR(COMPLEX'(1.0,0.0))-Z(-p))/14.0;

Finally we can investigate the function (11), which is the composition of functions * H _{4}(z)* and

*H*:

_{M}(z)(17)

The diagram of resulting MFR |*H(ω)*| (Fig.7,c) is generated using the following operator.

Hz:=(Allpass2x2(b,a,p)+Z(-2.0*p)) *COMPLEX_TO_POLAR(COMPLEX'(1.0,0.0))+Z(-p) *COMPLEX_TO_POLAR(COMPLEX'(1.0,0.0))-Z(-7.0*p))/ (COMPLEX_TO_POLAR(COMPLEX'(1.0,0.0))-Z(-p))/28.0;

As we can see in Fig.7, the function * H _{М}(z)* masks the function

*H*. As a result, we derive the transfer function

_{4}(z)*H(z)*of the low pass filter which cutoff frequency is equal to

*f*= 0,06 and suppression level is higher than 23 db.

_{c}Concluding the article we can see that the simple transfer function (12), which can be implemented only with additions, has rather high filtering properties. The FIR filter with the analogous characteristics has at least the order 60. The properties of this function were investigated by using only VHDL, not to consider Matlab or other popular DSP modeling tools.

Further data about allpass filters, masking and multiple delay filter one can find in the references below.

References

1. Chung J. G., Parhi K. K. Pipelined wave digital filter design for narrow-band sharp-transition digital filters // Proc. IEEE Workshop VLSI Signal Processing, La Jolla, CA. –1994. –Р. 501-510.

2. Regalia P., Mitra S.K., Vaidyanathan P.P. The Digital All-Pass Filter: A Versatile Signal Processing Building Block // Proc. IEEE. –1988. –V.76. –№1. – p.19-37.