**4. Efficient methods to design FIR filters**

It is known that the complexity of FIR digital filters increases in an inverse proportion with the transition bandwidth. A simple example in Sub-section 2.1 was given for the case of a Hilbert transformer with small ripple and bandwidth. Several efficient techniques have been developed to efficiently design FIR filters with strict specifications, such that the resulting filter accomplishes the desired specification with a lower complexity than the direct design obtained with the Parks-McClellan algorithm. For Hilbert transformers, however, these methods have been specially adapted, since Hilbert transform filtering has special characteristics as we saw in Section 2. Prior to starting the review of the special methods to design FIR Hilbert transformers with low complexity, we will briefly introduce in this section the general techniques which are the origin of these methods. These techniques are Frequency-Response Masking (FRM), Frequency Transformation (FT) and Piecewise Polynomial Sinusoidal (PPS).

#### **4.1. Frequency-Response Masking technique**

458 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 1

π ω

π/2

Hilbert transformers designed from half-band filters have odd length. In these cases, there is a coefficient of value zero between each coefficient of its impulse response. Thus, a Hilbert transformer with even length can be obtained by eliminating these zero-valued coefficients.

**Figure 8.** Hilbert transformer design, (a) Ideal zero-phase frequency response of the Hilbert transformer

(a) (b)

0 0

*h*(*n*)

½

M 2M

n

**Example 3.** The following code example illustrates a simple way to obtain the impulse response of a Hilbert transformer from the impulse response of a half-band filter. If we consider that the code of Example 2 has been previously run, it can be assumed that h\_half and L are already defined. The resulting Hilbert transformer coefficients in h are the same as the ones obtained in Example 1, since the half-band filter generated in Example 2 accomplish the relations *ωL* = (*π*/2) – *ωp* and *δ*= 2*δp*with regard to the specifications in

It is known that the complexity of FIR digital filters increases in an inverse proportion with the transition bandwidth. A simple example in Sub-section 2.1 was given for the case of a Hilbert transformer with small ripple and bandwidth. Several efficient techniques have been developed to efficiently design FIR filters with strict specifications, such that the resulting filter accomplishes the desired specification with a lower complexity than the direct design obtained with the Parks-McClellan algorithm. For Hilbert transformers, however, these methods have been specially adapted, since Hilbert transform filtering has special characteristics as we saw in Section 2. Prior to starting the review of the special methods to design FIR Hilbert transformers with low complexity, we will briefly introduce in this section the general techniques which are the origin of these methods. These techniques are Frequency-Response Masking (FRM), Frequency Transformation (FT) and Piecewise

middle\_sample = [zeros(1,(L-1)/2) 1/2 zeros(1,(L-1)/2)];

obtained from *HHb*(*ω*). (b) Impulse response *h*(*n*).

0


*H*(*ω*)

1



m = (ones(1,L)\*i).^(-(index-1));

Polynomial Sinusoidal (PPS).

h = 2\*(h\_half.\*m - middle\_sample.\*m);

**4. Efficient methods to design FIR filters** 

Example 1.

index = [1:L];

The FRM technique, introduced in [16], uses the so-called expanded-by-*M* filters as basic building blocks, where the transfer functions have the form *G*(*zM*). In general, a filter *G*(*z*) becomes expanded-by-*M* by replacing every of its elements *z*–1 by *z*–*<sup>M</sup>* or, in other words, by inserting *M*–1 zero-valued impulse response samples between two of its original impulse response samples. The periodic frequency response of these filters has *M* periods in the frequency range [0, 2*π*].

Figures 6a to 6d show how a filter *F*(*z*) = *Q*(*z*2) has a compressed-by-two frequency response in comparison with the original filter *Q*(*z*), whose impulse response is depicted in Figure 6a. The zero-phase frequency response of *Q*(*z*), presented in Figure 6b, shows a period that covers the frequency range [0, 2*π*]. On the other hand, the filter *F*(*z*) has a very similar zerophase frequency response, with the only difference that its period covers the frequency range [0, *π*]. The transition bandwidth of this expanded-by-2 filter *F*(*z*) is a half of the transition bandwidth of the filter *Q*(*z*). However, now this expanded filter has a replica of the frequency response of *Q*(*z*) (which covers the range of frequency [2*π*, 4*π*]) over the range [*π*,2*π*]. The number of multipliers of *F*(*z*) is the same as the one of *Q*(*z*). This can be seen in Figure 6c, where several impulse response samples are zero-valued.

The main idea of the FRM technique is using an expanded filter *G*(*zM*) and its complementary filter, *Gc*(*zM*), to form the transition band of a desired filter *H*(*z*). Because of that, these filters are so-called band-edge shaping filters. The complementary filter is given as

$$\mathbf{G}\_{\mathbf{c}}(\mathbf{z}) = \mathbf{z}^{-(L\_{\mathbb{C}} - 1)/2} - \mathbf{G}(\mathbf{z}),\tag{33}$$

where *LG* is the length of the filter *G*(*z*). Since the frequency response of these filters is periodical, two non-periodic masking filters, *HMa*(*z*) and *HMc*(*z*), are respectively cascaded with *G*(*zK*) and *Gc*(*zK*) to eliminate the unwanted periodic replicas of frequency response. The overall filter formed with the FRM technique is given as

$$H(z) = G(z^M)H\_{M\mathfrak{a}}(z) + [z^{-M(L\_{\mathbb{C}}-1)/2} - G(z^M)]H\_{M\mathfrak{a}}(z) \,. \tag{34}$$

Extensive information about the basic FRM method can be found in [16] and [23].

#### **4.2. Frequency Transformation technique**

The FT technique, first studied in [22] and then generalized in [17], is based on the repetitive use of an identical simple subfilter *G*(*z*). Let us consider *G*(*ω*) as the zero-phase frequency response of *G*(*z*) and an amplitude change function *Q*(*x*) given as

$$Q(\mathbf{x}) = \sum\_{k=0}^{M} q(k)\mathbf{x}^{k} \,. \tag{35}$$

The function *Q*(*x*) allows changing the values *x = G*(*ω*)to new values *y = Q*(*x*). Basically, the new amplitude values *y* = *Q*(*x*) must approximate the desired values *d* = *D*(*x*) for *x Xp Xs*, where *Xp* is the range of values [*xp*,*<sup>l</sup>*, *xp*,*<sup>u</sup>*] and *Xs* is the range of values [*xs*,*<sup>l</sup>*, *xs*,*<sup>u</sup>*], such that the zero-phase frequency response of the overall filter *H*(*z*) achieves the desired values *d* with a maximum absolute pass-band deviation *δp* over the pass-band region Ω*p*, as well as a maximum absolute stop-band deviation *δs* over the stop-band region Ω*s*. This characteristic is reached if the following conditions are simultaneously met,

$$D(\mathbf{x}) - \mathcal{S}\_p \le \mathbf{Q}(\mathbf{x}) \le D(\mathbf{x}) + \mathcal{S}\_{p^{\prime}} \text{ for } \mathbf{x}\_{p,l} \le \mathbf{x} \le \mathbf{x}\_{p,u^{\prime}} \tag{36}$$

$$D(\mathbf{x}) - \mathcal{S}\_s \le Q(\mathbf{x}) \le D(\mathbf{x}) + \mathcal{S}\_{s'} \quad \text{for } \mathbf{x}\_{s,l} \le \mathbf{x} \le \mathbf{x}\_{s,u'} \tag{37}$$

$$\mathbf{x}\_{p,l} \le \mathbf{G}(\boldsymbol{\phi}) \le \mathbf{x}\_{p,u'} \quad \text{for } \boldsymbol{\phi} \in \boldsymbol{\Omega}\_{p'} \tag{38}$$

$$\mathbf{x}\_{s,l} \le \mathbf{G}(o\rho) \le \mathbf{x}\_{s,\mu'} \text{ for } o \in \Omega\_s \tag{39}$$

Usually, *D*(*x*) = 1 for *x Xp* and *D*(*x*) = 0 for *x Xs*. Basically, two problems can be solved from this approach:

*Problem 1.* Given *M*, the number of subfilters, find the optimal coefficients of *Q*(*x*) and the optimal coefficients of *G*(*z*) to meet the conditions (36) to (39) with the minimum length *LG* (which must be odd).

*Problem 2.* Given the subfilter *G*(*z*), find the optimal coefficients of *Q*(*x*) to meet the conditions (36) to (39) with the minimum value *M*.

The overall filter formed with the FT technique is given as

$$H(z) = \sum\_{k=0}^{M} q(k) z^{-(M-k)(L\_G - 1)/2} \left[ G(z) \right]^k. \tag{40}$$

Detailed information about the FT method can be found in [16] and [23].

#### **4.3. Piecewise Polynomial Sinusoidal technique**

In the PPS technique, introduced in [24] for wide-band Type I filters, extended in [25] for Hilbert transformers and detailed in [11] for both cases, the idea is to divide the impulse response of a wideband filter into sub-responses and to generate each sub-response with polynomials with a given degree. For wideband linear-phase filters, the impulse response has a narrow main lobe and the side lobes have very rapid change in sign. Therefore, it is taken advantage of sinusoids in such a way that the polynomial pieces follow the polynomial-sinusoidal shapes to decrease the number of polynomial pieces and, as a consequence, to reduce the number of coefficients.

The overall transfer function *H*(*z*) for a desired Type I filter with length 2*N*+1 is constructed with *M* parallel branches connected and delayed with *z*–*Nm* in order to keep the center of symmetry at the same location for all the sub-impulse responses. These sub-responses are modulated with a sinusoidal function and finally an arbitrary number of separately generated filter coefficients is added as follows,

Digital FIR Hilbert Transformers: Fundamentals and Efficient Design Methods 461

$$H(\mathbf{z}) = \sum\_{m=1}^{M} \mathbf{z}^{-N\_m} H\_m(\mathbf{z}) + \mathbf{z}^{-\hat{N}} \hat{H}(\mathbf{z}) \, , \tag{41}$$

where,

460 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 1

is reached if the following conditions are simultaneously met,

conditions (36) to (39) with the minimum value *M*.

The overall filter formed with the FT technique is given as

**4.3. Piecewise Polynomial Sinusoidal technique** 

consequence, to reduce the number of coefficients.

generated filter coefficients is added as follows,

*k*

Detailed information about the FT method can be found in [16] and [23].

*M*

0

from this approach:

(which must be odd).

where *Xp* is the range of values [*xp*,*<sup>l</sup>*, *xp*,*<sup>u</sup>*] and *Xs* is the range of values [*xs*,*<sup>l</sup>*, *xs*,*<sup>u</sup>*], such that the zero-phase frequency response of the overall filter *H*(*z*) achieves the desired values *d* with a maximum absolute pass-band deviation *δp* over the pass-band region Ω*p*, as well as a maximum absolute stop-band deviation *δs* over the stop-band region Ω*s*. This characteristic

, , ( ) ( ) ( ) , for , *Dx Qx Dx x x x <sup>p</sup> p pl pu*

, , ( ) ( ) ( ) , for , *Dx Qx Dx x x x <sup>s</sup> s sl su*

, , ( ) , for , *<sup>p</sup> <sup>l</sup> p u <sup>p</sup> xG x*

, , ( ) , for *s l s u <sup>s</sup> xG x*

Usually, *D*(*x*) = 1 for *x Xp* and *D*(*x*) = 0 for *x Xs*. Basically, two problems can be solved

*Problem 1.* Given *M*, the number of subfilters, find the optimal coefficients of *Q*(*x*) and the optimal coefficients of *G*(*z*) to meet the conditions (36) to (39) with the minimum length *LG*

*Problem 2.* Given the subfilter *G*(*z*), find the optimal coefficients of *Q*(*x*) to meet the

( )( 1)/ 2

() () ( ) *<sup>G</sup>*

In the PPS technique, introduced in [24] for wide-band Type I filters, extended in [25] for Hilbert transformers and detailed in [11] for both cases, the idea is to divide the impulse response of a wideband filter into sub-responses and to generate each sub-response with polynomials with a given degree. For wideband linear-phase filters, the impulse response has a narrow main lobe and the side lobes have very rapid change in sign. Therefore, it is taken advantage of sinusoids in such a way that the polynomial pieces follow the polynomial-sinusoidal shapes to decrease the number of polynomial pieces and, as a

The overall transfer function *H*(*z*) for a desired Type I filter with length 2*N*+1 is constructed with *M* parallel branches connected and delayed with *z*–*Nm* in order to keep the center of symmetry at the same location for all the sub-impulse responses. These sub-responses are modulated with a sinusoidal function and finally an arbitrary number of separately

*MkL k*

*Hz qkz G z* . (40)

 

 

 

(36)

(37)

(38)

(39)

$$H\_m(z) = h\_m(N - N\_m)z^{-(N - N\_m)} + \sum\_{n=0}^{(N - N\_m) - 1} h\_m(n)[z^{-n} + z^{-(2\langle N - N\_m \rangle - n)}].\tag{42}$$

The integers *Nm* in the delay terms *z*–*Nm* satisfy *N*1 = 0 and *Nm*+1>*Nm* for *m* = 1, 2, …, *M* – 1, and the order of *Hm*(*z*) is 2(*N – Nm*). The impulse response is given as

$$\mathcal{H}\_m(\mathfrak{n}) = \sum\_{r=0}^L a\_m^{(L)}(r) n^r \times \sin\left[\mathcal{O}\_c(n - (N - N\_m))\right]. \tag{43}$$

In addition, <sup>ˆ</sup> ˆ ( ) *<sup>N</sup> z Hz* is a conventional direct-form transfer function with non-zero impulse response coefficients ˆ *h n*( ) , with *n* = *N*–*c*+1, *N*–*c*+2, …, *N*–*c*+*T*, where *c* = *T* / 2 and *T* is the number of additional coefficients at the center of the filter. Given the values *Nm* for *m* = 1, 2, …, *M* – 1, *M* and *L*, the objective is finding the polynomial coefficients such that the error with respect to a desired amplitude characteristic is minimized. An extensive explanation on this method can be reviewed in [11].

#### **4.4. Pipelining-Interleaving architecture**

The Pipelining-Interleaving (PI) technique developed in [19] provides efficient structures of FIR digital filters to avoid the repetitive use of an identical filter. Suppose that we have two sequences of independent signals, *x*1(*n*) and *x*2(*n*), that are filtered by two identical filters *H*(*z*). Thus, two corresponding sequences of independent outputs, *y*1(*n*) and *y*2(*n*), are obtained. An alternative form for this purpose is the multirate implementation using *H*(*z*2) as shown in Figure 9. This structure uses a single filter to implement two identical filters. The clock rate for this implementation must be twice the data rate [19]. If only one sequence of input signal is filtered, it is possible to connect the first output sequence *y*1(*n*) to the second input *x*2(*n*). In this way, *H*(*z*2) is used to implement *H*2(*z*).

**Figure 9.** Filtering of two independent sequences using one filter.

The PI structure of the Figure 9 can be extended to implement the filtering of *K* different signals, each one filtered by an identical filter *H*(*z*), with *K* being an arbitrary positive integer. From this, it is possible to implement the filtering of one signal with *K* identical filters in cascade. Figure 10a presents the general structure to filter a signal using one filter *H*(*zK*). Figure 10b shows the equivalent structure, which consists of the filtering of a signal by a cascade of *K* identical filters *H*(*z*) [19].

In the structure shown in Figure 10a, the clock rate of *H*(*zK*) must be *K* times the data rate. Clearly, for high data rate applications, *K* must be chosen as a relatively small integer, otherwise a very high clock rate will be required. More details on this time-multiplexed architecture can be found in [19].

## **5. Efficient Methods to design FIR Hilbert transformers**

It has been mentioned earlier that the design of low-complexity FIR Hilbert transformers with stringent specifications requires special efficient methods. In the following we will review the most representative and useful methods, which are based on the techniques revised in the previous section.

**Figure 10.** Filtering of a sequence with *K* identical cascaded filters, (a) PI architecture with only one filter, (b) equivalent single-rate structure.

#### **5.1. Hilbert transformer design based on Frequency Response Masking**

This method, proposed in [9], relies on the special case of FRM for the synthesis of a halfband filter *HHb*(*z*) [26]. Consider a half-band filter *Ha*(*z*) as a band-edge shaping filter whose transfer function is given by

$$H\_a(z) = \frac{1}{2}z^{-2K+1} + A(z),\tag{44}$$

$$A(z) = \sum\_{k=1}^{K} a(2k - 1 + 2K - 1) \left[ z^{2k - 1 - 2K + 1} + z^{-2k + 1 - 2K + 1} \right] \tag{45}$$

where *La*= 4*K* – 1 is the length of *Ha*(*z*), with *K* being an integer greater than zero, and *a*(*n*) is the impulse response of *Ha*(*z*). Replacing *G*(*z*) by *Ha*(*z*) in (34) we have

$$\mathcal{H}\_{\rm Idb}(\mathbf{z}) = \mathcal{H}\_{\rm a}(\mathbf{z}^{\prime M})\mathcal{H}\_{\rm Ma}(\mathbf{z}) + [\mathbf{z}^{-\mathcal{M}(2\mathcal{K}-1)} - \mathcal{H}\_{\rm a}(\mathbf{z}^{\prime M})]\mathcal{H}\_{\rm Mc}(\mathbf{z}).\tag{46}$$

We can express the transfer function of the overall half-band filter as

$$H\_{Hb}(\mathbf{z}) = \left[\frac{1}{2}\mathbf{z}^{-M(2K-1)} + A(\mathbf{z}^{M})\right]H\_{Ma}(\mathbf{z}) + \left[\frac{1}{2}\mathbf{z}^{-M(2K-1)} - A(\mathbf{z}^{M})\right]H\_{Mc}(\mathbf{z})\,. \tag{47}$$

If *M* is odd, then either [ <sup>1</sup> (2 1) <sup>2</sup> ( ) *<sup>M</sup> K M z Az* ] or[ <sup>1</sup> (2 1) <sup>2</sup> ( ) *<sup>M</sup> K M z Az* ] has a transition band centered at *π*/2, just as desired for a half-band filter.

If *M* is given by the form

462 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 1

**5. Efficient Methods to design FIR Hilbert transformers** 

*K*

*K*

*K*

*K*

architecture can be found in [19].

revised in the previous section.

*K*

<sup>2</sup> *x n*( )

<sup>3</sup> *x n*( )

*K*

*K*

( ) *<sup>K</sup> x n*

filter, (b) equivalent single-rate structure.

transfer function is given by

In the structure shown in Figure 10a, the clock rate of *H*(*zK*) must be *K* times the data rate. Clearly, for high data rate applications, *K* must be chosen as a relatively small integer, otherwise a very high clock rate will be required. More details on this time-multiplexed

It has been mentioned earlier that the design of low-complexity FIR Hilbert transformers with stringent specifications requires special efficient methods. In the following we will review the most representative and useful methods, which are based on the techniques

**Figure 10.** Filtering of a sequence with *K* identical cascaded filters, (a) PI architecture with only one

1( ) *<sup>K</sup> y n*

<sup>2</sup> *y* ( ) *n*

( ) *a* ( ) *b*

( ) *<sup>K</sup> y n*

This method, proposed in [9], relies on the special case of FRM for the synthesis of a halfband filter *HHb*(*z*) [26]. Consider a half-band filter *Ha*(*z*) as a band-edge shaping filter whose

> 1 2 1 <sup>2</sup> ( ) ( ), *<sup>K</sup> H z z Az <sup>a</sup>*

( ) (2 1 2 1) ,

where *La*= 4*K* – 1 is the length of *Ha*(*z*), with *K* being an integer greater than zero, and *a*(*n*) is

(2 1) () ( ) () [ ( )] ( ). *<sup>M</sup> MK M H z Hz H z z Hz H z Hb a Ma a Mc*

*Az a k K z z*

(44)

*H z*( )

<sup>1</sup> *x n*( ) <sup>2</sup> *x n*( )

<sup>3</sup> *x n*( )

( ) *<sup>K</sup> x n*

*H z*( )

<sup>1</sup> *z*

<sup>1</sup> *y* ( ) *n*

<sup>2</sup> *y n*( )

( ) *<sup>K</sup> y n*

<sup>1</sup> *z*

<sup>1</sup> *z*

*H z*( )

*H z*( )

2 12 1 2 12 1

(45)

(46)

*kK kK*

**5.1. Hilbert transformer design based on Frequency Response Masking** 

1

the impulse response of *Ha*(*z*). Replacing *G*(*z*) by *Ha*(*z*) in (34) we have

*k*

*K*

<sup>1</sup> *z* <sup>1</sup> *z*

<sup>1</sup> *z* <sup>1</sup> *z*

<sup>1</sup> *z* <sup>1</sup> *z*

<sup>1</sup> *x n*( ) <sup>1</sup> *<sup>y</sup>* ( ) *<sup>n</sup>*

*K* ( ) *<sup>K</sup> H z*

$$M = 4k + 1, \quad \text{with} \; k = \{0, 1, 2, \ldots\} \; , \tag{48}$$

then the pass-band of *HMa*(*z*) is greater than the pass-band of *HMc*(*z*). With *ωp* and *ωs* denoting the band-edge frequencies of *HHb*(*z*), these values can be express as [26]

$$
\alpha\_p = \frac{2\pi m + \theta}{M}, \qquad \alpha\_s = \frac{2\pi m + \phi}{M}, \tag{49}
$$

where *m* is an integer less than *M*. The values *m*, *θ* and *ϕ* can be calculated as

$$m = \left\lfloor \frac{\alpha\_p M}{2\pi} \right\rfloor, \qquad \theta = \alpha\_p M - 2\pi m, \qquad \phi = \alpha\_s M - 2\pi m \tag{50}$$

where *x* represents the integer part of *x*, whereas *θ* and *ϕ* are the pass-band and stopband edges of *Ha*(*z*). The pass-band and stop-band edge frequencies of the masking filter *HMa*(*z*), *θMa* and *ϕMa*, as well as the pass-band and stop-band edge frequencies of the masking filter *HMc*(*z*), *θMc* and *ϕMc*, are given by

$$
\theta\_{\rm Ma} = \alpha\_p = \frac{2\pi m + \theta}{M}, \phi\_{\rm Ma} = \frac{2\pi(m+1) - \phi}{M}, \theta\_{\rm Mc} = \frac{2\pi m - \theta}{M}, \phi\_{\rm Mc} = \alpha\_s = \frac{2\pi m + \phi}{M} \tag{51}
$$

If *M* is given by the form

$$M = 4k + 3, \quad \text{con } k = \{0, 1, 2, \ldots\} \text{ .} \tag{52}$$

then the passband of *HMc*(*z*) is greater than the pass-band of *HMa*(*z*). In this case the frequencies *ωp* and *ωs* are given by

$$
\alpha\_p = \frac{2\pi m - \phi}{M}, \alpha\_s = \frac{2\pi m - \theta}{M}, \tag{53}
$$

To calculate *m*, *θ* and *ϕ* we use the following relations,

$$
\sigma m = \left[\frac{\alpha\_{\rm s} M}{2\pi}\right]\_{\rm r} \theta = 2\pi m - \alpha\_{\rm s} M \,\,\,\phi = 2\pi m - \alpha\_{\rm p} M \,\,\tag{54}
$$

where *x* represents the rounding operation to the closest integer greater than *x*. The values, *θMa*, *ϕMa*, *θMc* and *ϕMc*, are given by

$$\theta\_{\rm Ma} = \frac{2\pi(m-1) + \phi}{M}, \phi\_{\rm Ma} = \alpha\_{\rm s} = \frac{2\pi m - \theta}{M}, \theta\_{\rm Mc} = \alpha\_{p} = \frac{2\pi m - \phi}{M}, \phi\_{\rm Mc} = \frac{2\pi m + \theta}{M}.\tag{55}$$

If *HMa*(*z*) is a Type I filter with length *LMa* = 4*k* + 1, where *k* is an integer, we can write

$$H\_{\rm Ma}(\mathbf{z}) = h\_{\rm Ma} \left( \left( L\_{\rm Ma} - 1 \right) / 2 \right) + \sum\_{k=1}^{(L\_{\rm Ma} - 1)/2} h\_{\rm Ma} \left( k + \left( L\_{\rm Ma} - 1 \right) / 2 \right) \left[ z^{k - (L\_{\rm Ma} - 1)/2} + z^{-k - (L\_{\rm Ma} - 1)/2} \right]. \tag{56}$$

Now we define the transfer functions *B*(*z*) and *C*(*z*) as

$$\begin{split} B(z) &= h\_{\rm Ma} \left( 1 + \left( L\_{\rm Ma} - 1 \right) / 2 \right) \Bigg[ z^{1 - \left( L\_{\rm Ma} - 1 \right) / 2} + z^{-1 - \left( L\_{\rm Ma} - 1 \right) / 2} \Bigg] + \\ &h\_{\rm Ma} \left( 3 + \left( L\_{\rm Ma} - 1 \right) / 2 \right) \Bigg[ z^{3 - \left( L\_{\rm Ma} - 1 \right) / 2} + z^{-3 - \left( L\_{\rm Ma} - 1 \right) / 2} \Bigg] + \dots \end{split} \tag{57}$$

$$\begin{aligned} \mathbf{C}(\mathbf{z}) &= h\_{\mathrm{Ma}} \left( \left( L\_{\mathrm{Ma}} - \mathbf{1} \right) / 2 \right) + \\ &h\_{\mathrm{Ma}} \left( 2 + \left( L\_{\mathrm{Ma}} - \mathbf{1} \right) / 2 \right) \bigg[ z^{2 - \left( L\_{\mathrm{Ma}} - \mathbf{1} \right) / 2} + z^{-2 - \left( L\_{\mathrm{Ma}} - \mathbf{1} \right) / 2} \bigg] + \\ &h\_{\mathrm{Ma}} \left( 4 + \left( L\_{\mathrm{Ma}} - \mathbf{1} \right) / 2 \right) \bigg[ z^{4 - \left( L\_{\mathrm{Ma}} - \mathbf{1} \right) / 2} + z^{-4 - \left( L\_{\mathrm{Ma}} - \mathbf{1} \right) / 2} \bigg] + ... \end{aligned} \tag{58}$$

with *hMa*(*n*) as coefficients of the filter *HMa*(*z*). Replacing (57) and (58) in (56), we have

$$H\_{Ma}(z) = B(z) + C(z) \, \text{.}\tag{59}$$

In the half-band filter design, the masking filters are related by

$$H\_{Mc}(\mathbf{z}) = \mathbf{z}^{-(L\_{Ma}-1)/2} - H\_{Ma}(-\mathbf{z}) \cdot \tag{60}$$

From (59) and (60), and noting that *B*(–*z*) = –*B*(*z*) and that *C*(–*z*) = *C*(*z*), we have

$$H\_{Mc}(z) = z^{-(L\_{Ma}-1)/2} + B(z) - C(z) \,. \tag{61}$$

Once known the transfer function of the masking filters from (59) and (61), the overall transfer function of the half-band filter can be obtained by substituting (59) and (61) in (47). Finally, we obtain

$$H\_{\rm Hb}(z) = \frac{1}{2}z^{-[M(2K-1)+(L\_{\rm Mh}-1)/2]} + z^{-M(2K-1)}B(z) + A(z^M)[2C(z) - z^{-(L\_{\rm Mh}-1)/2}].\tag{62}$$

The half-band filter with desired deviation *δ* and pass-band edge frequency *ωp* can be designed with the FRM technique by applying the following steps:

1. Get the optimal value of *M*, using the following approximation,

$$M\_{opt} \approx \frac{1}{2} \sqrt{\frac{2\pi}{\alpha\_s - \alpha\_p}} \cdot \tag{63}$$

Note that the obtained value must be rounded to an odd integer.

464 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 1

*x* represents the rounding operation to the closest integer greater than *x*. The

(55)

*L L*

1 ( 1)/2 1 ( 1)/2

*Ma Ma*

*L L*

*L L*

*Ma Ma*

2 ( 1)/2 2 ( 1)/2

*L L*

*Ma Ma*

4 ( 1)/2 4 ( 1)/2

*<sup>M</sup>* () () () *H z Bz Cz <sup>a</sup>* . (59)

( 1)/2 ( ) ( ) *Ma <sup>L</sup> Hzz H z Mc Ma* . (60)

( 1)/2 ( ) () () *Ma <sup>L</sup> H z z Bz Cz Mc* . (61)

*M* . (63)

*Ma Ma*

3 ( 1)/2 3 ( 1)/2

   

(57)

(58)

 

( 1)/2 ( 1)/2

*Ma Ma*

*k L k L*

2 ( 1) 2 2 2 *Ma* , *Ma s* , ,. *Mc p Mc m m mm M M MM*

 

*H zh L h kL z z* . (56)

3 ( 1) / 2 ...

4 ( 1) / 2 ...

Once known the transfer function of the masking filters from (59) and (61), the overall transfer function of the half-band filter can be obtained by substituting (59) and (61) in (47).

The half-band filter with desired deviation *δ* and pass-band edge frequency *ωp* can be

2 *opt*

1 2

 *s p*

 <sup>1</sup> [ (2 1) ( 1)/2] (2 1) ( 1)/2 <sup>2</sup> ( ) ( ) ( )[2 ( ) ] *<sup>M</sup> K LMa Ma <sup>M</sup> K M <sup>L</sup> Hz z Hb z Bz Az Cz z* . (62)

 

If *HMa*(*z*) is a Type I filter with length *LMa* = 4*k* + 1, where *k* is an integer, we can write

where

Finally, we obtain

values, *θMa*, *ϕMa*, *θMc* and *ϕMc*, are given by

*Ma Ma Ma Ma Ma*

Now we define the transfer functions *B*(*z*) and *C*(*z*) as

*k*

1 ( ) ( 1) / 2 ( 1) / 2 *Ma*

*Bz h L z z*

*hL z z*

*hL z z*

*hL z z*

with *hMa*(*n*) as coefficients of the filter *HMa*(*z*). Replacing (57) and (58) in (56), we have

From (59) and (60), and noting that *B*(–*z*) = –*B*(*z*) and that *C*(–*z*) = *C*(*z*), we have

*Ma Ma*

*Ma Ma*

In the half-band filter design, the masking filters are related by

designed with the FRM technique by applying the following steps: 1. Get the optimal value of *M*, using the following approximation,

( 1)/2

*L*

( ) 1 ( 1) / 2

*Ma Ma*

*Ma Ma*

( ) ( 1) / 2 2 ( 1) / 2

*Ma Ma*

*Cz h L*


A Hilbert transformer can be derived from a unity gain half-band filter by subtracting the constant ½ from its transfer function and then modulating the remaining coefficients by *e– <sup>j</sup>πn*/2 (see section 3.1). The transfer function of the Hilbert transformer *H*(*z*) is given by [9]

$$H(\mathbf{z}) = \mathcal{Z}(j\mathbf{z})^{-M(2K-1)}B(j\mathbf{z}) + \mathcal{Z}A\left((j\mathbf{z})^M\right)[2C(j\mathbf{z}) - (j\mathbf{z})^{-(L\_{M\mathbf{z}}-1)/2}].\tag{64}$$

It is worth highlighting that the filter in (64) does not need complex-number arithmetic processing because of the following reasons. First, note that the filter *A*(*jMzM*) has only real coefficients since the imaginary unit generated by (*jz*)*–n* with *n* odd is eliminated by zerovalued coefficients in these indexes *n*. Second, note that all the coefficients in [2*C*(*jz*) – (*jz*)– (*LMa*–1)/2] are real and all the coefficients in *B*(*jz*) are imaginary when *LMa* is expressed as 4*k*+1, with *k* integer. Third, the term (*jz*)–*<sup>M</sup>*(2*<sup>K</sup>*–1) is always imaginary, since its exponent is always odd. From these reasons, we have that if *B*(*jz*) has imaginary coefficients, the term (*jz*)–*<sup>M</sup>*(2*K*–1) makes them real and the overall filter has real coefficients.

The Hilbert transformer from (64) can be seen as a parallel connection of two branches. In the first branch we have *Hb*(*z*) = 2(*jz*)–*<sup>M</sup>*(2*K*–1)*B*(*jz*) and in the second branch we have the cascade of *H*1(*zM*) and *HM*(*z*), where *H*1(*zM*) = 2*A*(*j MzM*) and *HM*(*z*) = [2*C*(*jz*) – (*jz*)–(*LMa*–1)/2]. This structure is presented in Figure 11. Let us review a different point of view of the FRM technique, presented in [10].

The filter *Hb*(*z*) can be seen as a low-order Hilbert transformer, the filter *H*1(*zM*) as a bandedge shaping filter and *HM*(*z*) as a masking filter. The basic filter *Hb*(*z*) provides a low order approximation (with wide transition bandwidth) to the desired specification. The cascaded connection of *H*1(*zM*) and *HM*(*z*) produces a correction term to the transfer function that decreases the transition bandwidth. The transfer function for the overall filter is given by

$$H(\mathbf{z}) = H\_1(\mathbf{z}^M) H\_M(\mathbf{z}) + H\_b(\mathbf{z}) \,. \tag{65}$$

Let the lengths of *Hb*(*z*), *H*1(*z*) and *HM*(*z*) be *Lb*, *L*1 and *LM*, respectively. The length of *H*1(*zM*)*HM*(*z*) is *ML*1 + *LM* – *M*. The delay introduced by *Hb*(*z*) and the delay introduced by *H*1(*zM*)*HM*(*z*) must be the same; otherwise, pure delays must be introduced into the shorterdelay branch to equalize them.

**Figure 11.** Structure for the synthesis of a Hilbert transformer using the FRM technique.

In order to avoid inserting half-sample delay in the implementation, the parities of *Lb* and of (*ML*1 + *LM* – *M*) must be the same. Furthermore, *H*1(*zM*)*HM*(*z*)must have anti-symmetrical impulse response.

Consider the magnitude response of *Hb*(*z*) as |*Hb*(*ejω*)|, as shown in Figure 12a, where *Lb* is even. The computational complexity of *Hb*(*z*) is low since its transition band is wide. Now consider the magnitude response of a transition band correction filter, |*He*(*ejω*)|, as shown in Figure 12b. The Hilbert transformer with sharp transition bandwidth, as shown in Figure 12c, is obtained from the parallel connection of the correction filter with *Hb*(*z*).

The objective is designing a correction filter with very low complexity using FRM technique. Consider the band-edge shaping filter *H*1(*z*) with magnitude response |*H*1(*ejω*)|, as shown in Figure 12d. The complexity of *H*1(*z*) is low because it has a wide transition band. Replacing each delay of *H*1(*z*) by *M* delays, a magnitude response |*H*1(*ejMω*)| is obtained, as shown in Figure 12e. A masking filter *HM*(*z*), with magnitude response |*HM*(*ejω*)| shown in Figure 12f, is used to mask the unwanted pass-band of |*H*1(*ejMω*)|. With this masking, the magnitude response |*He*(*ejω*)|, shown in Figure 12b, is produced. *HM*(*z*) has low complexity because its magnitude response has a wide transition band.

Since the length of *Hb*(*z*) is even, the length of *H*1(*zM*) *HM*(*z*), i.e.,*MN*1 + *NM* – *M*, must also be even. If *M* is odd, *N*1and*NM*must have different parities. By considering the gain of |*HM*(*ejω*)| in the vicinity of *ω* = 0, it is clear that *HM*(*z*) has symmetrical impulse response. Thus, *H*1(*z*) must have anti-symmetrical impulse response to satisfy the condition that *H*1(*zM*) *HM*(*z*) must have anti-symmetrical impulse response.

The band-edges of *Hb*(*z*) and *HM*(*z*) are the same. Let the band-edge of *Hb*(*z*) be *θb* and let the band-edge of *H1*(*z*) be *θ*1. It can be seen from Figure 12 that the value *θb* satisfies *θ<sup>b</sup>* ≤ (2*π* – *θ*1)/*M*. For an arbitrary value *θ*1, it is possible to obtain *θb* if the appropriate value of *M* is known, which is obtained with the objective of minimizing the overall number of coefficients. Finally, the overall filter is designed with a joint simultaneous optimization of *H*1(*zM*), *Hb*(*z*) and *HM*(*z*). For the examples in [10], the algorithm in [27] was used.

In the following we present a simple example to design an efficient FIR Hilbert transformer with stringent specifications based on the FRM technique. The approach presented in this example follows the procedure based on the four steps to design a half-band filter in terms of filter *A*(*z*), *B*(*z*) and *C*(*z*).

The Hilbert transformer is derived using (64). Since this approach does not require a simultaneous optimization for all the filters, it is simple and straightforward. Additionally, the sensitivity to the rounded coefficients is less since every filter is designed separately [28].

**Figure 11.** Structure for the synthesis of a Hilbert transformer using the FRM technique.

12c, is obtained from the parallel connection of the correction filter with *Hb*(*z*).

magnitude response has a wide transition band.

have anti-symmetrical impulse response.

of filter *A*(*z*), *B*(*z*) and *C*(*z*).

impulse response.

In order to avoid inserting half-sample delay in the implementation, the parities of *Lb* and of (*ML*1 + *LM* – *M*) must be the same. Furthermore, *H*1(*zM*)*HM*(*z*)must have anti-symmetrical

( ) *Hb z*

*Y z*( ) *X z*( )

1( ) *<sup>M</sup> H z* ( ) *HM z*

Consider the magnitude response of *Hb*(*z*) as |*Hb*(*ejω*)|, as shown in Figure 12a, where *Lb* is even. The computational complexity of *Hb*(*z*) is low since its transition band is wide. Now consider the magnitude response of a transition band correction filter, |*He*(*ejω*)|, as shown in Figure 12b. The Hilbert transformer with sharp transition bandwidth, as shown in Figure

The objective is designing a correction filter with very low complexity using FRM technique. Consider the band-edge shaping filter *H*1(*z*) with magnitude response |*H*1(*ejω*)|, as shown in Figure 12d. The complexity of *H*1(*z*) is low because it has a wide transition band. Replacing each delay of *H*1(*z*) by *M* delays, a magnitude response |*H*1(*ejMω*)| is obtained, as shown in Figure 12e. A masking filter *HM*(*z*), with magnitude response |*HM*(*ejω*)| shown in Figure 12f, is used to mask the unwanted pass-band of |*H*1(*ejMω*)|. With this masking, the magnitude response |*He*(*ejω*)|, shown in Figure 12b, is produced. *HM*(*z*) has low complexity because its

Since the length of *Hb*(*z*) is even, the length of *H*1(*zM*) *HM*(*z*), i.e.,*MN*1 + *NM* – *M*, must also be even. If *M* is odd, *N*1and*NM*must have different parities. By considering the gain of |*HM*(*ejω*)| in the vicinity of *ω* = 0, it is clear that *HM*(*z*) has symmetrical impulse response. Thus, *H*1(*z*) must have anti-symmetrical impulse response to satisfy the condition that *H*1(*zM*) *HM*(*z*) must

The band-edges of *Hb*(*z*) and *HM*(*z*) are the same. Let the band-edge of *Hb*(*z*) be *θb* and let the band-edge of *H1*(*z*) be *θ*1. It can be seen from Figure 12 that the value *θb* satisfies *θ<sup>b</sup>* ≤ (2*π* – *θ*1)/*M*. For an arbitrary value *θ*1, it is possible to obtain *θb* if the appropriate value of *M* is known, which is obtained with the objective of minimizing the overall number of coefficients. Finally, the overall filter is designed with a joint simultaneous optimization of

In the following we present a simple example to design an efficient FIR Hilbert transformer with stringent specifications based on the FRM technique. The approach presented in this example follows the procedure based on the four steps to design a half-band filter in terms

The Hilbert transformer is derived using (64). Since this approach does not require a simultaneous optimization for all the filters, it is simple and straightforward. Additionally, the sensitivity to the rounded coefficients is less since every filter is designed separately [28].

*H*1(*zM*), *Hb*(*z*) and *HM*(*z*). For the examples in [10], the algorithm in [27] was used.

**Figure 12.** Magnitude responses of the subfilters for even length *Hb*(*z*). Note that *ωL*/2*π* is the desired transition bandwidth.

**Example 4.** The following code example illustrates the design of a Type III Hilbert transformer with *δ*= 0.0001, *ωL*= 0.00125*π* and *ωH*= *π*– *ωL*= 0.99875*π* using the MATLAB Signal Processing Toolbox and the Filter Design Toolbox to generate the half-band filter. The conversion of filters with argument *z* into filters with argument *jz* is performed with the same principle illustrated in the code of example 3.

```
wL=0.00125*pi; wH=0.99875*pi; d=0.0001; % Hilbert transformer specification
wp=pi/2 - wL; ws=pi/2 + wL; %find the half-band band-edge frequencies
dp=d/2; ds=d/2; %find the half-band ripple specification
%--------STEP 1. Optimum M----------
M = (1/2)*round(sqrt(2*pi/(ws-wp))); if mod(M,2)==0; M = M+1; end
%-----------------------------------
%-------STEP 2. Band-edge-Shaping and Masking Filters--------------
if mod((M-1)/4,1)==0; 
 m=floor(wp*M/(2*pi)); theta=(wp*M)-(2*pi*m); phi=(ws*M)-(2*pi*m); 
 theta_Ma=wp; phi_Ma=(2*pi*(m+1)-phi)/M; theta_Mc=((2*pi*m)-theta)/M; 
 phi_Mc=ws; 
else
 m=ceil(ws*M/(2*pi)); theta=(2*pi*m)-(ws*M); phi=(2*pi*m)-(wp*M); 
 theta_Ma=(2*pi*(m-1)+phi)/M; phi_Ma=ws; theta_Mc=wp; 
 phi_Mc=((2*pi*m)+theta)/M;
```

```
end
La=firpmord([theta/pi phi/pi],[1 0], [0.85*dp 0.85*ds]); 
if mod(La,2)==0; La=La+1; end
if mod((La+1)/4,1)~=0; La=La+2; end
K=(La+1)/4; 
[L_Ma,fo,ao,W]=firpmord([theta_Ma/pi phi_Ma/pi],[1 0], [0.85*dp 0.85*ds]); 
if mod(L_Ma,2)==0; L_Ma=L_Ma+1; end
if mod((L_Ma-1)/4,1)~=0; L_Ma=L_Ma+2; end
ha = firhalfband(La-1,theta/pi); 
h_Ma = firpm(L_Ma-1,fo,ao,W); 
%----------------------------------------------------------------
%--------STEP 3. Filters A(z), B(z) and C(z)---------------------
a = ha - [zeros(1,2*K-1) 1/2 zeros(1,2*K-1)]; 
for i=1:L_Ma; if mod(i-1,2)==0; b(i)=0; c(i)=h_Ma(i); else b(i)=h_Ma(i); 
c(i)=0; end
end
delay_b = [zeros(1,M*(2*K-1)) 1 zeros(1,M*(2*K-1))]; 
delay_c = [zeros(1,(L_Ma-1)/2) 1 zeros(1,(L_Ma-1)/2)]; 
%---------------------------------------------------------------
%--------STEP 4. Form the Hilbert transformer (or half-band filter)-------
a_M = upsample(a,M); a_M = a_M(1:end-(M-1)); 
index_a_M = [1:length(a_M)]; 
m_a_M = (ones(1,length(a_M))*j).^(-(index_a_M-1)); 
index_b = [1:length(b)]; 
m_b = (ones(1,length(b))*j).^(-(index_b-1)); 
index_c = [1:length(c)]; 
m_c = (ones(1,length(c))*j).^(-(index_c-1)); 
index_delay_b = [1:length(delay_b)]; 
m_delay_b = (ones(1,length(delay_b))*j).^(-(index_delay_b-1)); 
index_delay_c = [1:length(delay_c)]; 
m_delay_c = (ones(1,length(delay_c))*j).^(-(index_delay_c-1)); 
h = 2*conv(delay_b.*m_delay_b, b.*m_b) +...
 2*conv(a_M.*m_a_M,(2*c.*m_c - delay_c.*m_delay_c)); 
[H w] = freqz(h,1,10000); 
figure; plot(w/pi, abs(H)) 
%----------------------------------------------------------
```
Figure 13 shows the magnitude response of the obtained Hilbert transformer. The overall structure requires 148 coefficients in total, i.e., 69 for *A*(*j MzM*), 39 for *B*(*jz*) and 40 for *C*(*jz*). Clearly, the FRM-based design is a very efficient method comparing to a direct design, like the one presented in example 1, where the estimated length is 4017 and which would require approximately 1005 coefficients.

#### **5.2. Hilbert transformer Design based on Frequency Transformation**

The Frequency Transformation (FT) method, developed in [12] to design FIR Hilbert transformers, allows designing FIR Hilbert transformers using a tapped cascaded interconnection of repeated simple basic building blocks constituted by two identical subfilters. To this end, two simple filters are required, namely, a prototype filter and a subfilter. The number of times that the subfilter is used, as well as the coefficients used between each cascaded subfilter, depends on the prototype filter. Both, the prototype filter and the subfilter are Hilbert transformers. The former is always a Type IV filter whereas the latter can be a Type III or Type IV filter according to the type of the desired Hilbert transformer [12].

The prototype filter must be a Type IV FIR filter, i.e., with even length given as *LP* = 2*N* and anti-symmetric impulse response of the form *p*(2*N* – 1 – *n*) = –*p*(*n*). Its frequency response is expressed as

$$P(e^{j\Omega}) = e^{-j\left((2N-1)\Omega/2 - \pi/2\right)}P(\Omega) \, , \tag{66}$$

where *P*(Ω), the zero-phase term, is given by

468 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 1

%---------------------------------------------------------------- %--------STEP 3. Filters A(z), B(z) and C(z)---------------------

%---------------------------------------------------------------

m\_delay\_b = (ones(1,length(delay\_b))\*j).^(-(index\_delay\_b-1));

m\_delay\_c = (ones(1,length(delay\_c))\*j).^(-(index\_delay\_c-1));

2\*conv(a\_M.\*m\_a\_M,(2\*c.\*m\_c - delay\_c.\*m\_delay\_c));

%----------------------------------------------------------

for i=1:L\_Ma; if mod(i-1,2)==0; b(i)=0; c(i)=h\_Ma(i); else b(i)=h\_Ma(i);

%--------STEP 4. Form the Hilbert transformer (or half-band filter)-------

[L\_Ma,fo,ao,W]=firpmord([theta\_Ma/pi phi\_Ma/pi],[1 0], [0.85\*dp 0.85\*ds]);

La=firpmord([theta/pi phi/pi],[1 0], [0.85\*dp 0.85\*ds]);

end

K=(La+1)/4;

c(i)=0; end

end

if mod(La,2)==0; La=La+1; end

if mod((La+1)/4,1)~=0; La=La+2; end

if mod(L\_Ma,2)==0; L\_Ma=L\_Ma+1; end

ha = firhalfband(La-1,theta/pi); h\_Ma = firpm(L\_Ma-1,fo,ao,W);

if mod((L\_Ma-1)/4,1)~=0; L\_Ma=L\_Ma+2; end

a = ha - [zeros(1,2\*K-1) 1/2 zeros(1,2\*K-1)];

a\_M = upsample(a,M); a\_M = a\_M(1:end-(M-1));

m\_b = (ones(1,length(b))\*j).^(-(index\_b-1));

m\_c = (ones(1,length(c))\*j).^(-(index\_c-1));

h = 2\*conv(delay\_b.\*m\_delay\_b, b.\*m\_b) +...

index\_delay\_b = [1:length(delay\_b)];

index\_delay\_c = [1:length(delay\_c)];

m\_a\_M = (ones(1,length(a\_M))\*j).^(-(index\_a\_M-1));

index\_a\_M = [1:length(a\_M)];

index\_b = [1:length(b)];

index\_c = [1:length(c)];

[H w] = freqz(h,1,10000); figure; plot(w/pi, abs(H))

delay\_b = [zeros(1,M\*(2\*K-1)) 1 zeros(1,M\*(2\*K-1))]; delay\_c = [zeros(1,(L\_Ma-1)/2) 1 zeros(1,(L\_Ma-1)/2)];

$$P(\Omega) = j \cdot \sin\left(\frac{\Omega}{2}\right) \sum\_{n=0}^{N-1} \tilde{d}(n) \cos\left(\Omega n\right),\tag{67}$$

and *Ω* denotes the frequency domain of the prototype filter. The coefficients *d n*( ) can be obtained directly from the impulse response *p*(*n*) [1].

Using the equivalence cos(*Ωn*) = *Tn*{cos(*Ω*)} [17], where *Tn*{*x*}is the nth-degree Chebyshev polynomial defined with the following recursive formulas,

**Figure 13.** Magnitude responses of the FRM-based Hilbert transformer. From left to right: overall magnitude response, transition bandwidth detail and passband ripple detail.

$$T\_0\{\mathbf{x}\} = \mathbf{1}, \ T\_1\{\mathbf{x}\} = \mathbf{x} \quad \text{and} \quad T\_n\{\mathbf{x}\} = 2\mathbf{x}T\_{n-1}\{\mathbf{x}\} - T\_{n-2}\{\mathbf{x}\} \tag{68}$$

the zero-phase term can be rewritten as

$$P(\Omega) = j \cdot \sin\left(\frac{\Omega}{2}\right) \sum\_{n=0}^{N-1} a(n) \left[\cos\left(\Omega\right)\right]^n,\tag{69}$$

where *α*(*n*) are obtained from *d n*( ) using the coefficients of the Chebyshev polynomials. Based on the equivalence given as

$$2\cos\left(2\pi\right) = 1 - 2\sin^2\left(\pi\right) = 1 + 2\left(j \cdot \sin\pi\right)^2,\tag{70}$$

the zero-phase term can be expressed by

$$P(\Omega) = j \cdot \sin\left(\frac{\Omega}{2}\right) \sum\_{n=0}^{N-1} \alpha(n) \left[1 + 2\left(j \cdot \sin\left(\frac{\Omega}{2}\right)\right)^2\right]^n. \tag{71}$$

Consider the case of a Type III subfilter with odd length given as *LG*= 2*M* + 1 and antisymmetric impulse response of the form *g*(2*M* – *n*) = –*g*(*n*). Its frequency response is expressed as

$$\mathcal{G}(e^{j\alpha}) = e^{-j\left(2M\alpha/2\right)}\mathcal{G}(\alpha)\,,\tag{72}$$

where *G*(*ω*) is the zero-phase term, given by

$$G(\omega) = j \cdot \sum\_{n=1}^{M} c(n) \sin(\omega n) \,. \tag{73}$$

The coefficients *c*(*n*) can be obtained directly from *g*(*n*) [1]. Note that the term *G*(*ω*) can be put in (71) by using the following expression,

$$j \cdot \sin\left(\frac{\Omega}{2}\right) = j \cdot \sum\_{n=1}^{M} c(n) \sin(\alpha n) \,\,\,\,\tag{74}$$

resulting in

$$\mathbf{H}(\alpha) = \mathbf{j} \cdot \sum\_{n=1}^{M} \mathbf{c}(n) \sin(\alpha n) \sum\_{n=0}^{N-1} \alpha(n) \left[ 1 + 2 \left( \mathbf{j} \cdot \sum\_{n=1}^{M} \mathbf{c}(n) \sin(\alpha n) \right)^{2} \right]^{n},\tag{75}$$

where *H*(*ω*) is the zero-phase term of the overall filter. Therefore, the frequency transformation is obtained from (74) and is given by

$$\Omega = 2 \sin^{-1} \left[ \sum\_{n=1}^{M} c(n) \sin(\omega n) \right]. \tag{76}$$

Equation (76) implies that the magnitude response of the prototype filter is preserved, but its frequency domain is changed by the subfilter.

The transfer function of the overall Hilbert transformer is given as

$$H(z) = G(z) \sum\_{n=0}^{N-1} z^{-2M\{N-1-n\}} a(n) \left[ H\_1(z) \right]^n, \ H\_1(z) = z^{-2M} + 2G^2(z) \tag{77}$$

with *G*(*z*) being the transfer function of the subfilter.

470 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 1

the zero-phase term can be rewritten as

where *α*(*n*) are obtained from

Based on the equivalence given as

expressed as

resulting in

the zero-phase term can be expressed by

where *G*(*ω*) is the zero-phase term, given by

put in (71) by using the following expression,

transformation is obtained from (74) and is given by

0 1 1 2 { } 1, { } and { } 2 { } { }, *n nn T x T x x T x xT x T x* (68)

*Pj n* , (69)

<sup>2</sup> <sup>2</sup> cos 2 1 2sin 1 2 sin *x x jx* , (70)

*Pj n j* . (71)

*d n*( ) using the coefficients of the Chebyshev polynomials.

*<sup>N</sup> <sup>n</sup>*

 

 

Consider the case of a Type III subfilter with odd length given as *LG*= 2*M* + 1 and antisymmetric impulse response of the form *g*(2*M* – *n*) = –*g*(*n*). Its frequency response is

*n*

1 ( ) ( )sin( ) *M*

The coefficients *c*(*n*) can be obtained directly from *g*(*n*) [1]. Note that the term *G*(*ω*) can be

*n*

 2 1 sin ( )sin( ) *M*

where *H*(*ω*) is the zero-phase term of the overall filter. Therefore, the frequency

*n*

 <sup>1</sup> 1 2sin ( )sin( ) *M*

2 2

1 2

<sup>2</sup> <sup>1</sup>

*H j cn n n j cn n* , (75)

2 /2 ( ) ( ) *<sup>j</sup> j M Ge e G* , (72)

*G j cn n* . (73)

*j j cn n* , (74)

*cn n* . (76)

*n*

*N n*

*n*

2 0 ( ) sin ( ) cos

*n*

*MN M*

 

*nn n*

10 1 ( ) ( )sin( ) ( ) 1 2 ( )sin( )

0 ( ) sin ( ) 1 2 sin

1

For a desired Hilbert transformer specification expressed as in (17), the magnitude response |*P*(Ω)| of the prototype filter must satisfy the following condition,

$$\left| (1 - \delta) \le \left| P(\Omega) \right| \le (1 + \delta), \quad \text{for } \Omega\_\perp \le \Omega \le \pi \right. \tag{78}$$

with *ΩL* being the lower band-edge frequency of the prototype filter. The magnitude response of the subfilter, |*G*(*ω*)|, must fulfill simultaneously

$$\left| \left| \left. v\_{\mathrm{d}} - \delta\_{\mathrm{G}} \right| \leq \left| \mathbb{G}\_{0}(o) \right| \leq 1 \right. \right. \\ \left. \left. \left( \left. o \right| \left. o \right| \leq o \right. \right. \right. \\ \left. \left. o \right| \leq o \leq \left. \pi - o \right|\_{\mathrm{L}} \tag{79}$$

$$
\sigma\_d = \frac{1}{2} + \frac{1}{2}\sin\left(\frac{\Omega\_l}{2}\right), \; \mathcal{S}\_G = \frac{1}{2} - \frac{1}{2}\sin\left(\frac{\Omega\_l}{2}\right). \tag{80}
$$

The design procedure proposed in [12] starts with an arbitrary prototype filter, and then the subfilter is designed accordingly.

Note that the complexity of the subfilter depends almost exclusively on its transition bandwidth, since its ripple specification is considerably relaxed. Similarly, the prototype filter is a low-complexity filter because, even though its ripple specification is strict, its transition bandwidth is relaxed. The relaxed ripple specification of the subfilter makes it suitable to be implemented as a simple, multiplierless system with rounded coefficients [13]. Additionally, it was observed in [13] that the repeated use of identical subfilters can be avoided by taking advantage of the PI technique, which has been introduced in sub-section 4.4. Therefore, a time-multiplexed design with lower area can be obtained.

Figure 14 presents the PI-based architecture proposed in [13]. This structure was straightforwardly derived from [19], where a similar example is given for the sharpening technique of [22]. Additionally, the design of the Hilbert transformer was made multiplierless by applying rounding to the coefficients of the prototype filter and the subfilter.

Instead of choosing an arbitrary prototype filter as in [12], a heuristic search was employed to select the prototype filter and the subfilter, such that the proposed architecture uses a number of coefficients less or equal to 0.25 times the estimated number of multipliers required in a direct design using the Parks-McClellan algorithm.

**Figure 14.** PI-based structure with three subfilters [13].

**Figure 15.** PI-based structure with two subfilters [14].

In [14] the authors observed that the cascaded interconnection of the two subfilters *G*(*z*) required to build *H*1(*z*) can be decoupled and also implemented with the PI technique. Thus, the PI-based architecture shown in Figure 15, which only requires two subfilters, was obtained. The approach of PI-based architectures for FT designs was further developed in [15], and a simple procedure to derive a PI-based structure from a FT-based design with identical subfilters was proposed. With this procedure, the architecture presented in Figure 16 was proposed for Hilbert transformers, where only a simple subfilter is required.

**Figure 16.** PI-based structure with one subfilter [15].

472 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 1

*<sup>K</sup>* + ) *G*(*z*

2

*K*)

+

*z* –1

*z* –1  *K*

 *K*

 *K*

+

*α*(*K*)

*z* –*M*–1

+

. . . . . .

+

+

. . . . . .

+

+

*α*(0)

*z –M–*2

*z –M–*2

*α*(1)

+

+

*z –M–*2

*z* –*M*–1

*z* –*M*–1

*α*(0)

*α*(1)

*z* –1

+

<sup>+</sup> *<sup>α</sup>*(*K*)

*z* –1

*z* –1  *K*

 *K*

 *K*

*z* –1

*M* = *LG* – 1

*z* –1

. . .

*z* –1

. . .

*K* = *LP*/2 – 1

 *K K*

*z* –*KM*

*M* = *LG* – 1

 *K K*

*z –KM–*1

*G*(*z* <sup>2</sup>*<sup>K</sup>* 2 + ) *z* –1

> *z* –1

2

2

+

. . .

*z* –1

 *K*

 *K*

 *K*

*z* –1

+

*G*(*z*) *G*(*z*

*z* –1

*z* –1

**Figure 14.** PI-based structure with three subfilters [13].

*G*(*z*) +

**Figure 15.** PI-based structure with two subfilters [14].

In [14] the authors observed that the cascaded interconnection of the two subfilters *G*(*z*) required to build *H*1(*z*) can be decoupled and also implemented with the PI technique. Thus, the PI-based architecture shown in Figure 15, which only requires two subfilters, was obtained. The approach of PI-based architectures for FT designs was further developed in

*K* = *LP*/2 – 1

2 .

2

*z* –1

. .

 *K*

 *K*

 *K*

*z* –1

+

*z* –1

*z* –1

An important insight proposed in [14] was avoiding the arbitrary selection of the prototype filter as in [12] through the optimized search of the adequate prototype filter, such that a cost metric is minimized. Since multipliers are the most expensive elements in digital filters, reducing the overall number of coefficients is the goal. In general terms, this is equivalent to improve the original heuristic search proposed in [13]. From (78), (79) and (80) it can be observed that the prototype filter and the subfilter can be designed if the frequency Ω*L* is known. The problem consists on finding the optimal frequency Ω*L*.

Consider a function *φ*(*δ*, *ωL*), which can estimate with an acceptable exactitude the length of a HT in terms of its ripple *δ* and its lower passband edge *ωL*, such as the one presented in (18). We have, for the prototype filter and the subfilter,

$$L\_{\mathcal{G}} \approx \mathcal{L}\_{\mathcal{G}} = \phi(\mathcal{S}\_{\mathcal{G}} / \upsilon\_{d \prime} \cdot o\_{\mathcal{L}}) \, \prime \tag{81}$$

 *L* ( , ) *LPP L* , (82)

where *LG* and *LP* are the respective approximations to the lengths of the subfilter, *LG*, and the prototype filter, *LP*, whereas *vd* and *δG* are given in (80). Clearly, *LG* and *LP* are given as functions of the ripple and transition band of the subfilter and of the prototype filter, respectively.

For the previously revised PI-based structures, it is always possible to express the overall number of multipliers *m* in terms of the numbers of multipliers of the prototype filter and the subfilter as,

$$m = f(m\_{\odot'}, m\_P) \, \, \, \, \, \, \tag{83}$$

$$m\_{\mathbb{G}} \approx \mathbb{C} \cdot L\_{\mathbb{G}}\,\prime\,.\tag{84}$$

$$m\_{\rm P} = L\_{\rm P} / \mathcal{D} \,, \tag{85}$$

where *mP* and *mG* are taken from (19).

Substituting *vd* and *δG* from (80) in (81), and using the approximations (81) and (82) respectively in (84) and (85) we have

$$\mathfrak{m}(\delta,\ o\_{\mathbb{L}},\ \Omega\_{\mathbb{L}}) \approx f\left(\mathbb{C} \cdot \phi(\frac{1-\sin(\Omega\_{\mathbb{L}}/2)}{1+\sin(\Omega\_{\mathbb{L}}/2)},\ o\_{\mathbb{L}}),\ \frac{1}{2} \cdot \phi(\delta,\ \Omega\_{\mathbb{L}})\right). \tag{86}$$

Note that, even though *m*(*δ*, *ωL*, Ω*L*) is a function of *δ*, *ωL* and Ω*L*, the values *δ* and *ωL* are known a priori because they are given by the problem at hand (see (17)). Therefore, since the unique unknown is Ω*L*, the approach consists in finding the optimum value Ω\* *<sup>L</sup>* for Ω*L*, with 0 < Ω*L*<*π*, such that *m*(*δ*, *ωL*, Ω*L*) is minimized. This optimization problem is given as

$$\begin{array}{c} \min m(\mathcal{S}, \; \mathcal{O}\_{\mathcal{L}}, \; \Omega\_{\mathcal{L}})\\ \mathcal{\Omega}\_{\mathcal{L}} \; \text{e} \; \mathbb{R} \\ \text{such that } \; 0 < \Omega\_{\mathcal{L}} < \pi\_{\text{\textquotedblleft}} \end{array} \tag{87}$$

where *m*(*δ*, *ωL*, Ω*L*) is given in (86). The result obtained from (86) is an estimation which depends on the exactitude of the function *φ*(*δ*, *ωL*).

Equation (18) was utilized in [17] as the function *φ*(*δ*, *ωL*). However, this function does not give good length estimation for filters with a huge ripple (like the subfilters in the FT method). Using the proposal from [29] as starting point, we have recently derived the following more accurate formula,

$$\begin{split} \phi(\delta,\ o\_{L}) = & \frac{1}{2} + \left[ \frac{1.101 \left[ -\log\_{10}(\delta) \right]^{1.1}}{\left( \frac{o\_{l}}{2\pi} \right)} + 1 \right] \cdot \left[ \frac{2}{3\pi} \arctan \left[ \left2.325 \left[ 0.30103 - \log\_{10}(\delta) \right]^{-0.45} \cdot \dots \right] \right. \\ \left. \left( \frac{o\_{l}}{2\pi} \right)^{-1.39} \right] \left[ \frac{1}{0.5 - \left( \frac{o\_{l}}{2\pi} \right)} \right] + \frac{1}{6} \right]. \end{split} \tag{88}$$

Thus, the FT method consists on finding the optimum value Ω*L* by solving (87). With this value, the prototype filter and the subfilter are designed as given in (78) and (79). Finally, the coefficients *α*(*n*) are found from the prototype filter coefficients by relating (67) and (69) and the overall filter is synthesized by using any of the structures existing in literature [13]- [15] or [17].

**Example 5.** The following code example illustrates the design of a Type III Hilbert transformer with *δ*= 0.004, *ωL*= 0.01*π* and *ωH*= *π*– *ωL*= 0.99*π* using the MATLAB Signal Processing Toolbox. The optimized value for Ω*L* is Ω\* *<sup>L</sup>* = 0.2237*π* and the lengths for the prototype filter and the subfilter are, respectively, *LP* = 14 and *LG* = 31. The value Ω\* *<sup>L</sup>* has been optimized to minimize the number of coefficients in the structure of Figure 16. This code makes use of the MATLAB function ChebyshevPoly.m, which is available online [30].

474 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 1

the subfilter as,

where *mP* and *mG* are taken from (19).

respectively in (84) and (85) we have

 

depends on the exactitude of the function *φ*(*δ*, *ωL*).

2

*L*

 

2

*L*

1.1

following more accurate formula,

2

 

[15] or [17].

*L*

*L*

1.39

1 1 () . 0.5 ( ) <sup>6</sup>

For the previously revised PI-based structures, it is always possible to express the overall number of multipliers *m* in terms of the numbers of multipliers of the prototype filter and

Substituting *vd* and *δG* from (80) in (81), and using the approximations (81) and (82)

Note that, even though *m*(*δ*, *ωL*, Ω*L*) is a function of *δ*, *ωL* and Ω*L*, the values *δ* and *ωL* are known a priori because they are given by the problem at hand (see (17)). Therefore, since the

such that 0 , *L L L*

where *m*(*δ*, *ωL*, Ω*L*) is given in (86). The result obtained from (86) is an estimation which

Equation (18) was utilized in [17] as the function *φ*(*δ*, *ωL*). However, this function does not give good length estimation for filters with a huge ripple (like the subfilters in the FT method). Using the proposal from [29] as starting point, we have recently derived the

1 2 1.101 log ( ) ( , ) 1 arctan 2.325 0.30103 log ( ) ... 2 3 ( )

 

min ( , , )

 

10 0.445

Thus, the FT method consists on finding the optimum value Ω*L* by solving (87). With this value, the prototype filter and the subfilter are designed as given in (78) and (79). Finally, the coefficients *α*(*n*) are found from the prototype filter coefficients by relating (67) and (69) and the overall filter is synthesized by using any of the structures existing in literature [13]-

*L*

 

*m L L <sup>L</sup> <sup>L</sup> f C* . (86)

 1 sin( / 2) <sup>1</sup> 1 sin( / 2) <sup>2</sup> ( , , ) ( , ), ( , ) *<sup>L</sup> L*

 

unique unknown is Ω*L*, the approach consists in finding the optimum value Ω\*

0 < Ω*L*<*π*, such that *m*(*δ*, *ωL*, Ω*L*) is minimized. This optimization problem is given as

*m*

( , ) *m fm m G P* , (83)

*m CL G G* , (84)

/ 2 *m L P P* , (85)

*<sup>L</sup>* for Ω*L*, with

(87)

10

(88)

```
%*********** INITIAL DATA ****************
L_g=31; wl=0.01*pi; L_p=14; Om_L=0.2237*pi; 
%************** SUBFILTER**********************
wH=pi-wl; dG = (1/2) - (1/2)*sin(Om_L/2); 
vd = (( 1 - sin(Om_L/2) )/2) + sin(Om_L/2); 
g = firpm(L_g-1,[wl/pi wH/pi],[vd vd],'hilbert'); 
%*************PROTOTYPE FILTER*****************
[p]=firpm(L_p-1,[Om_L/pi 1],[1 1],'hilbert'); 
%***************BASIC BUILDING BLOCK H1********************
r1=2*conv(g,g); 
delay=[zeros(1,L_g-1) 1 zeros(1,L_g-1)]; 
h1=r1+delay; 
%**************ALPHA COEFFICIENTS FROM CHEBYSHEV POLYNOMIAL **************
N = L_p/2 
for mm=1:N; d(mm)=2*(p((N+1)-mm)); end
D(N)=2*d(N); 
for Mm=fliplr([3:N]); D(Mm-1)=(2*d(Mm-1))+(D(Mm)); D(1)=d(1)+((1/2)*D(2)); 
end
tt=0; 
for nn=fliplr([0:N-1]); tk=ChebyshevPoly(nn); T(nn+1,:)=[zeros(1,tt) tk'] ; 
 tt=tt+1; 
end
ll=sum((D'*[ones(1,N)]).*T); 
alpha=fliplr(ll); 
%**************OVERALL FILTER **************
upper_branch = g; lower_branch = g*alpha(1); h = lower_branch; 
for ii=1:N-1; 
upper_branch = conv(upper_branch,h1); 
lower_branch = conv(h, [zeros(1,L_g-1) 1 zeros(1,L_g-1)]); 
h=lower_branch + alpha(ii+1)*upper_branch;%Overall Hilbert transformer
end
[H w]=freqz(h,1,1000); 
figure 
plot(w/pi,abs(H))
```
Figure 17 shows the magnitude response of the obtained Hilbert transformer. The overall structure requires only 15 coefficients in total, i.e., 7 structural coefficients, *α*(0) to *α*(6), and 8 coefficients for the subfilter *G*(*z*). Compared to a direct design, where the estimated length is 287 and which would require approximately 72 coefficients, The FT-based design achieves almost a 75% of reduction in the number of distinct required coefficients.

**Figure 17.** Magnitude responses of FT-based Hilbert transformer. From left to right: overall magnitude response, transition bandwidth detail and passband ripple detail.
