Digital FIR Hilbert Transformers: Fundamentals and Efficient Design Methods

1. Introduction

IIR Hilbert transformers perform a phase approximation. This means that the phase response of the system is approximated to the desired values in a given range of frequencies. The magnitude response allows passing all the frequencies, with the magnitude obtained around the desired value within a given tolerance [6], [7]. On the other hand, FIR Hilbert transformers perform a magnitude approximation. In this case the system magnitude response is approximated to the desired values in a given range of frequencies. The advantage is that their phase response is always maintained in the desired value over the complete range of frequencies [8].

Whereas IIR Hilbert transformers can present instability and they are sensitive to the rounding in their coefficients, FIR filters can have exact linear phase and their stability is guaranteed. Moreover, FIR filters are less sensitive to the coefficients rounding and their phase response is not affected by this rounding. Because of this, FIR Hilbert transformers are often preferred [8]-[15]. Nevertheless, the main drawback of FIR filters is a higher complexity compared with the corresponding IIR filters. Multipliers, the most costly elements in DSP implementations, are required in an amount linearly related with the length of the filter. A linear phase FIR Hilbert transformer, which has an anti-symmetrical impulse response, can be designed with either an odd length (Type III symmetry) or an even length (Type IV symmetry). The number of multipliers m is given in terms of the filter length L as$$m\approx C\cdot L$$, where C = 0.25 for a filter with Type III symmetry or C = 0.5 for a filter with Type IV symmetry.

The design of optimum equiripple FIR Hilbert transformers is usually performed by Parks-McClellan algorithm. Using the MATLAB Signal Processing Toolbox, this becomes a straightforward procedure through the function firpm. However, for small transition bandwidth and small ripples the resulting filter requires a very high length. This complexity increases with more stringent specifications, i.e., narrower transition bandwidths and also smaller pass-band ripples. Therefore, different techniques have been developed in the last 2 decades for efficient design of Hilbert transformers, where the highly stringent specifications are met with an as low as possible required complexity. The most representative methods are [9]-[15], which are based in very efficient schemes to reduce complexity in FIR filters.

Methods [9] and [10] are based on the Frequency Response Masking (FRM) technique proposed in [16]. In [9], the design is based on reducing the complexity of a half-band filter. Then, the Hilbert transformer is derived from this half-band filter. In [10], a frequency response corrector subfilter is introduced, and all subfilters are designed simultaneously under the same framework. The method [11] is based on wide bandwidth and linear phase FIR filters with Piecewise Polynomial-Sinusoidal (PPS) impulse response. These methods offer a very high reduction in the required number of multiplier coefficients compared to the direct design based on Parks-McClellan algorithm. An important characteristic is that they are fully parallel approaches, which have the disadvantage of being area consuming since they do not directly take advantage of hardware multiplexing.

The Frequency Transformation (FT) method, proposed first in [17] and extended in [18], was modified to design FIR Hilbert transformers in [12] based on a tapped cascaded interconnection of repeated simple basic building blocks constituted by two identical subfilters. Taking advantage of the repetitive use of identical subfilters, the recent proposal [13] gives a simple and efficient method to design multiplierless Hilbert transformers, where a combination of the FT method with the Pipelining-Interleaving (PI) technique of [19] allows getting a time-multiplexed architecture which only requires three subfilters. In [14], an optimized design was developed to minimize the overall number of filter coefficients in a modified FT-PI-based structure derived from the one of [13], where only two subfilters are needed. Based on methods [13] and [14], a different architecture which just requires one subfilter was developed in [15].

In this chapter, fundamentals on digital FIR Hilbert transformers will be covered by reviewing the characteristics of analytic signals. The main connection existing between Hilbert transformers and half-band filters will be highlighted but, at the same time, the complete introductory explanation will be kept as simple as possible. The methods to design low-complexity FIR filters, namely FRM [16], FT [17] and PPS [11], as well as the PI architecture [19], which are the cornerstone of the efficient techniques to design Hilbert transformers presented in [9]-[15], will be introduced in a simplified and concise way. With such background we will provide an extensive revision of the methods [9]-[15] to design low-complexity efficient FIR Hilbert transformers, including MATLAB routines for these methods.

2. Complex signals, analytic signals and Hilbert transformers

A real signal is a one-dimensional variation of real values over time. A complex signal is a two-dimensional signal whose value at some instant in time can be specified by a single complex number. The variation of the two parts of the complex numbers, namely the real part and the imaginary part, is the reason for referring to it as two-dimensional signal [20]. A real signal can be represented in a two-dimensional plot by presenting its variations against time. Similarly, a complex signal can be represented in a three-dimensional plot by considering time as a third dimension.

Real signals always have positive and negative frequency spectral components, and these components are generally real and imaginary. For any real signal, the positive and negative parts of its real spectral component always have even symmetry around the zero-frequency point, i.e., they are mirror images of each other. Conversely, the positive and negative parts of its imaginary spectral component are always anti-symmetric, i.e., they are always negatives of each other [1]. This conjugate symmetry is the invariant nature of real signals.

Complex signals, on the other hand, are not restricted to these spectral conjugate symmetry conditions. The special case of complex signals which do not have a negative part neither in their real nor in their imaginary spectral components are known as analytic signalsor also as quadrature signals [2]. An example of analytic signal is the complex exponential signalx_c(t), presented in Figure 1, and described by

$$x_c(t)=e^{j{\omega }_{0}t}=x_r(t)+j{x}_i(t)=\cos ({\omega }_{0}t)+j\sin ({\omega }_{0}t).$$E1

The real part and the imaginary part of the analytic signal are related trough the Hilbert transform. In simple words, given an analytic signal, its imaginary part is the Hilbert transform of its real part. Figure 1 shows the complex signal x_c(t), its real part x_r(t) and its imaginary part, x_i(t). Figure 2 presents the frequency spectral components of these signals. It can be seen that the real part x_r(t) and the imaginary part x_i(t), both real signals, preserve the spectral conjugate symmetry. The complex signal x_c(t) does not have negative parts neither in its real spectral component nor in its imaginary spectral component. For this reason, analytic signals are also referred as one-side spectrum signals. Finally, Figure 3 shows the Hilbert transform relation between the real and imaginary parts of x_c(t).

Figure 1.

The Hilbert transform and the analytic signal of x_r(t) = cos(ω₀t), ω₀= 2π.

Figure 2.

From left to right, frequency spectrum of x_r(t), x_i(t) and x_c(t).

Figure 3.

Hilbert transform relations between x_r(t) and x_i(t) to generate x_c(t).

The motivation for creating analytic signals, or in other words, for eliminating the negative parts of the real and imaginary spectral components of real signals, is that these negative parts have in essence the same information than the positive parts due to the conjugate symmetry previously mentioned. The elimination of these negative parts reduces the required bandwidth for the processing. For the case of DSP applications, it is possible to form a complex sequence x_c(n) given as follows,

$$x_c(n)=x_r(n)+j{x}_i(n),$$E2

with the special property that its frequency spectrum X_c(e^jω) is equal to that of a given real sequence x(n) for the positive Nyquist interval and zero for the negative Nyquist interval, i.e.,

$$X_c(e^{j\omega })=\{\begin{array}{c}X(e^{j\omega })\text{ for }0\le \omega <\pi ,\\ \text{ 0 for }-\pi \le \omega <0.\end{array}$$E3

Although analyticity has no formal meaning for sequences [2], the same terminology, i.e., analytic sequence, will be applied for complex sequences whose frequency spectrum is one-sided, like in (3).

If X_r(e^jω) and X_i(e^jω) respectively denote the frequency spectrums of x_r(n) and x_i(n), then

$$X_c(e^{j\omega })=X_r(e^{j\omega })+j{X}_i(e^{j\omega }).$$E4

The spectrums of x_r(n) and x_i(n) can be readily deduced as

$$X_r(e^{j\omega })=\frac{1}{2}[X_c(e^{j\omega })+X_c^{*}(e^{-j\omega })],$$E5

$$j{X}_i(e^{j\omega })=\frac{1}{2}[X_c(e^{j\omega })-X_c^{*}(e^{-j\omega })],$$E6

where X_c^*(e^jω) is the complex conjugate of X_c(e^jω). Note that (6) gives an expression for jX_i(e^jω), which is the frequency spectrum of the imaginary signal jx_i(n). Also, note that X_r(e^jω) and X_i(e^jω) are both complex-valued functions in general. However, X_r(e^jω) is conjugate symmetric, i.e., X_r(e^jω) = X_r^*(e^–jω). Similarly, jX_i(e^jω) is conjugate anti-symmetric, i.e., jX_i(e^jω) = –jX_i^*(e^–jω). These relations are illustrated in Figure 4.

From (5) and (6) we obtain

$$X_c(e^{j\omega })=2X_r(e^{j\omega })-X_c^{*}(e^{-j\omega }),$$E7

$$X_c(e^{j\omega })=2j{X}_i(e^{j\omega })+X_c^{*}(e^{-j\omega }),$$E8

and since X_c^*(e^–jω) = 0 for 0 <ω<π (see Figure 4b), eqs. (3), (7) and (8) give

Figure 4.

Decomposition of an unilateral spectrum. Solid and dashed lines are, respectively, the real and imaginary parts.

$$X_c(e^{j\omega })=\{\begin{array}{c}2X_r(e^{j\omega })\text{ for }0\le \omega <\pi ,\\ \text{ 0 for }-\pi \le \omega <0.\end{array}$$E9

$$X_c(e^{j\omega })=\{\begin{array}{c}2j{X}_i(e^{j\omega })\text{ for }0\le \omega <\pi ,\\ \text{ 0 for }-\pi \le \omega <0.\end{array}$$E10

Thus

$$X_i(e^{j\omega })=-j{X}_r(e^{j\omega })\text{ for }0\le \omega <\pi .$$E11

On the other hand, from (4), and since X_c(e^jω) = 0 for –π$$\le$$ω< 0, we have

$$X_i(e^{j\omega })=j{X}_r(e^{j\omega })\text{ for }-\pi \le \omega <0.$$E12

Therefore, (11) and (12) can be expressed as

$$X_i(e^{j\omega })=\{\begin{array}{c}-j{X}_r(e^{j\omega })\text{ for }0\le \omega <\pi ,\\ j{X}_r(e^{j\omega })\text{ for }-\pi \le \omega <0,\end{array}$$E13

or

$$X_i(e^{j\omega })=H(e^{j\omega })X_r(e^{j\omega }),$$E14

where

$$H(e^{j\omega })=\{\begin{array}{c}-j\text{ for }0\le \omega <\pi ,\\ j\text{ for }-\pi \le \omega <0.\end{array}$$E15

According to (14), x_i(n) can be obtained by processing x_r(n) with a linear time-invariant discrete-time system whose frequency response H(e^jω) is given in (15). This frequency response has unity magnitude, a phase angle of –π/2 radians for 0 <ω<π, and a phase angle of π/2 radians for –π<ω< 0. A system of this type is commonly referred to as Hilbert transformer or sometimes as 90-degree phase shifter.

The impulse response h(n) of a Hilbert transformer is [2]

$$h(n)=\{\begin{array}{c}\frac{2}{n\pi }{\sin }^2\left(\frac{n\pi }{2}\right)\text{ for }n\ne 0,\\ \text{ 0 for }n=0.\end{array}$$E16

This impulse response is not absolutely summable and thus the frequency response of (15) is ideal. However, approximations to the ideal Hilbert transformer can be obtained with IIR or FIR systems. Thus, Hilbert transformers are considered a special class of filter.

IIR Hilbert transformers have phase error as well as magnitude error in approximating the ideal frequency response. Basically, these filters can be designed by using two all-pass systems whose phase responses differ by approximately π/2 over some well-defined portion of the band 0 < |ω|<π. By taking the outputs of the two all-pass filters as the real and imaginary parts of a complex signal it can be found that the spectrum of such signal nearly vanishes over much of the negative frequency interval. As such, the outputs of the two all-pass filters are quite nearly a Hilbert transformer.

FIR Hilbert transformers with constant group delay can be easily designed. The π/2 phase shift is realized exactly, with an additional linear phase component required for a causal FIR system. By evaluating (16) over some positive and negative values of n, it can be seen that the impulse response is anti-symmetric. Therefore, FIR Hilbert transformers are based on either Type III (i.e., anti-symmetric impulse response with odd length L) or Type IV (i.e., anti-symmetric impulse response with even length L) symmetry. Filters with Type III symmetry have amplitude equal to zero in ω = 0 and ω = π and filters with Type IV symmetry have amplitude equal to zero only in ω = 0. Thus, the FIR approximation is acceptable over a given range of frequencies (a pass-band region) which does not include these extremes.

The exactness of the phase of Type III and Type IV FIR systems is a compelling motivation for their use in approximating Hilbert transformers. Additionally, whereas IIR Hilbert transformers can present instability and they are sensitive to rounding error in their coefficients, FIR filters have guaranteed stability, are less sensitive to the coefficients rounding and their phase response is not affected by this rounding. Because of this, FIR Hilbert transformers are often preferred [8]-[15]. The rest of this chapter will be focused on the design of FIR Hilbert transformers.

2.1. Basic design of FIR Hilbert transformers with MATLAB

Since the phase requirement in FIR Hilbert transformers is accomplished, the design of a FIR Hilbert transformer consists on finding the impulse response h(n), for n = 0 to L–1, which satisfies the following magnitude response specification,

$$(1-\delta )\le \left|H(e^{j\omega })\right|\le (1-\delta )\text{ for }{\omega }_L\le \omega \le {\omega }_H,$$E17

where δ is the allowed pass-band ripple, ω_L is the lower pass-band edge and ω_H is given as ω_H =π– ω_L if the desired Hilbert transformer is a Type III filter or ω_H =π if it is Type IV. The values ω_L and ω_H can be made to approach 0 and π, respectively, as closely as desired by increasing the length L of the filter. For Hilbert transformers, the value ω_L/2π is considered the transition band.

The design of optimum equiripple Type III and Type IV FIR linear phase Hilbert transformers is usually performed by Parks-McClellan algorithm. With the MATLAB Signal Processing Toolbox this becomes a straightforward procedure through the function firpm. The order of the filter, L–1, must be estimated in advance. A useful formula to estimate L, presented in [10], is

$$\begin{array}{l}L\approx \phi (\delta ,{\omega }_L)=[0.002655{\left({\log }_{10}(\delta )\right)}^3+0.031843{\left({\log }_{10}(\delta )\right)}^2-\mathrm{...}\\ 0.554993{\log }_{10}(\delta )-0.049788]/(\frac{{\omega }_L}{2\pi })+1.\end{array}$$E18

Example 1. The following code example illustrates the design of a Type III Hilbert transformer with δ= 0.01, ω_L= 0.1π and ω_H= π– ω_L= 0.9π using the MATLAB Signal Processing Toolbox. From (18), L is estimated as L$$\approx$$24.3266. For convenience, we use L = 4k+3 with k integer and the closest value for L, higher than the estimated value, is chosen. In this case we use L = 4*6+3 = 27.

d = 0.01; w_L = 0.1*pi; w_H = 0.9*pi; L = 27; h = firpm(L-1,[w_L/pi w_H/pi],[1 1],'hilbert'); [H w] = freqz(h,1); figure; stem(0:length(h)-1, h,'fill'); figure; plot(w/pi, abs(H))

Figure 5a shows the impulse response and Figure 5b shows the magnitude response of the obtained Hilbert transformer.

As we mentioned earlier, the most expensive elements in digital filters are multipliers. For a Type III Hilbert transformer, the number of multipliers, m, is (L+1)/4 if L = 4k+3, or (L–1)/4 if L = 4k+1, with k integer. In this last case the impulse response values h(0) and h(L–1) are zero. For a Type IV Hilbert transformer, the number of multipliers is L/2. This number can be simplified as

$$m\approx C\cdot L,$$E19

where C = 0.25 for a Type III Hilbert transformer or C = 0.5 for a Type IV Hilbert transformer.

Figure 5.

a) Impulse response and (b) Magnitude response of a digital FIR Hilbert transformer.

It is worth highlighting the following point. A Type III Hilbert transformer can be derived from a Type IV Hilbert transformer by adding one zero-valued impulse response sample between each two impulse response samples of the Type IV Hilbert transformer. This is equivalent to replacing each z^–1 in the transfer function of this filter by z^–2. The opposite procedure can be straightforwardly followed to derive a Type IV Hilbert transformer from a Type III Hilbert transformer.

From (18), it can be seen that the length L increases inversely proportional with the transition band ω_L/2π. Hence, for cases with narrow transition band and small ripple the implementation cost becomes very high. As a simple example, the estimated length for a Hilbert transformer withδ= 0.001 and ω_L= 0.001π is L$$\approx$$3661.2. It implies approximately 915.3 multipliers, which is impractical. Therefore, the design of FIR Hilbert transformers with strict specifications requires specialized techniques to reduce the number of multipliers and the computational complexity of the filter. Before of the revision of these techniques, in the next section we will highlight the relation between Hilbert transformers and half-band filters because this relation is crucial for developing the specialized techniques to design low-complexity FIR Hilbert transformers.

3. The Hilbert transformer and its relation with the half-band filter

Half-band filters have their transfer function given by

$$H_{Hb}(z)={\displaystyle \sum _{n=0}^{2M}{h}_{Hb}(n)z^{-n}},$$E20

where h_Hb(n) is the impulse response, M is an odd integer and 2Misthe filter order [21]. The coefficients are symmetric with respect to the central coefficient h_Hb(M), namely, they accomplish the following relation,

$$h_{Hb}(2M-n)=h_{Hb}(n)\text{ for }n=0,1,\mathrm{...},2M.$$E21

The length of the filter, L, is an odd number given as L = 2M + 1 with M = 1, 3, 5… etc. In a lineal phase half-band filter, almost a half of the coefficients are zero. Figure 6 shows the procedure to design a half-band filter. It starts with the transfer function of a Type II lineal phase filter, i.e., a filter with symmetric impulse response whose length is even. For this filter we have,

$$Q(z)={\displaystyle \sum _{n=0}^{M}q(n)z^{-n}\text{ for }q(M-n)=q(n)}.$$E22

First, samples with value zero are introduced between the q(n) impulse samples (see Figure 6a and 6c). This generates a transfer function with Type I symmetry, whose order is 2M, given by

Figure 6.

Low-pass half-band filter design. (a) Impulse response of a Type II filter with order M, q(n). (b) Zero-phase frequency response of the Type II filter with order M, Q(ω). (c) Impulse response of a Type I filter with order 2M, f(n). (d) Zero-phase frequency response of the Type I filter with order 2M, F(ω). (e) Impulse response of a half-band filter, h_Hb(n). (f) Zero-phase frequency response of the half-band filter, H_Hb(ω).

$$F(z)={\displaystyle \sum _{n=0}^{2M}f(n)z^{-n}=Q(z^2)={\displaystyle \sum _{n=0}^{M}q(n)z^{-2n}}}.$$E23

Then, the zero-valued sample with index n = M is replaced by the value 1/2 (See Figure 6e). Thus we obtain

$$H_{Hb}(z)=\frac{1}{2}{z}^{-M}+F(z)=\frac{1}{2}{z}^{-M}+{\displaystyle \sum _{n=0}^{M}q(n)z^{-2n}}.$$E24

The coefficients of the half-band filter are obtained from (20) and (24) as

$$h_{Hb}(M)=\frac{1}{2},h_{Hb}(n)=q\left(\frac{n}{2}\right)\text{ for }n\text{ even,}$$E25

$$h_{Hb}(n)=0,\text{ for }n\text{ odd and }n\ne M,$$E26

and the zero-phase frequency response of H_Hb(z) is

$$H_{Hb}(\omega )=\frac{1}{2}+Q(2\omega )=\frac{1}{2}+F(\omega ).$$E27

Based on the previous relations, the design of a half-band low-pass filter with pass-band frequency ω_p and pass-band ripple δ can be carried out by designing the Q(z) filter, such that its zero-phase frequency response, Q(ω), oscillates within 1/2 ± δ over the range of frequencies [0, 2ω_p] (See Figure 6b). Since Q(z) has a Type II transfer function, it presents a fixed zero in z = –1 (ω = π). Note from Figure 6b that Q(ω) oscillates within –1/2 ± δ over the range of frequencies [2π – 2ω_p, 2π]. The corresponding zero-phase frequency response of the filter F(z), given as F(ω) = Q(2ω), remains within 1/2 ± δ over the range of frequencies [0, ω_p] and within –1/2 ± δ over the range of frequencies [π – ω_p, π] (See Figure 6d). Finally, H_Hb(ω) oscillates around 1 over the range of frequencies [0, ω_p] with tolerance δ and around 0 over the range of frequencies [π – ω_p, π] with the same tolerance δ (See Figure 6f). Note that, as a low-pass filter, the half-band filter has the relations ω_s = π–ω_p and δ_p = δ_s.

Example 2. The following code example illustrates the design of a half-band filter with δ_p=δ_s= 0.005, ω_p= 0.4π and ω_s= π– ω_p= 0.5π using the MATLAB Filter Design Toolbox. Eq. (18) can be applied to estimate the filter length L, by substituting ω_L = (π/2) – ω_p and δ= 2δ_p. The length L is estimated as L$$\approx$$24.3266. Since L must be represented as L = 2M+1 with M odd to avoid zero-valued impulse response samples in the left and right extremes of the impulse response, we use L =2*13 + 1 = 27.

3.1. Hilbert transformer derived from a half-band filter

A Hilbert transformer filter can be designed from a half-band filter. First,the sample with value 1/2, located in the index n = M, is replaced by the value 0 (See Figures 7a and 7d). Thus we obtain

$$\tilde{H}(z)=H_{Hb}(z)-\frac{1}{2}{z}^{-M}$$E28

.

The subtraction of this coefficient, drawn in Figure 7d, causes the zero-phase frequency response of the half-band filter to be shifted downwardly, as illustrated in Figure 7c. The filter $$\tilde{H}(\omega )$$must be shifted by π/2 on the horizontal direction in the frequency domain. This is accomplished by multiplying the remaining coefficients, $$\tilde{h}(n)$$, by (j)^–n. The effect of this action produces a Hilbert transformer with odd length and a pass-band gain equal to 0.5. Therefore, a Hilbert transformer with unitary gain will be obtained by scaling all the coefficients by 2. This is illustrated in Figure 8. The transfer function of the Hilbert transformer given in terms of the transfer function of a half-band filter is

$$H(z)=2\left[H_{Hb}(jz)-\frac{1}{2}{(jz)}^{-M}\right]=2{\displaystyle \sum _{\begin{array}{l}n=0\\ n\ne M\end{array}}^{2M}{h}_{Hb}(n){(jz)}^{-n}}$$E29

.

The impulse response of the Hilbert transformer is related with the impulse response of the half-band filter through the following expression,

$$h(n)=\{\begin{array}{c}0;n=2k+1,\\ 2{(-1)}^{k-1}{h}_{Hb}(n);n=2k\text{ with }k=0,1,2,\mathrm{...},M.\end{array}$$E30

It was mentioned earlier that the useful bandwidth in a Hilbert transformer is restricted to some range given as 0 < ω_L ≤ ω ≤ ω_H < π, where ω_H= π – ω_L. The relation of the low-pass edge frequency ω_L with the band-edge frequencies of the half-band filter, ω_p and ω_s, is given by

$${\omega }_L=(\pi /2)-{\omega }_p$$E31

$${\omega }_L={\omega }_s-(\pi /2)$$E32

where ω_s = π – ω_p.

Figure 7.

a) Ideal zero-phase frequency response of the half-band filter H_Hb(ω). (b) Impulse response of the half-band filter h_Hb(n). (c) Ideal zero-phase frequency response resulting of the subtraction of the central coefficient, located at n = M,$$\tilde{H}(\omega )=H_{Hb}(\omega )-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$$. (d) Impulse response,$$\tilde{h}(n)$$.

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 obtained from H_Hb(ω). (b) Impulse response h(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δ_pwith regard to the specifications in Example 1.

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

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(z^M). In general, a filter G(z) becomes expanded-by-M by replacing every of its elements z^–1 by z^–M 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²) 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 zero-phase 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(z^M) and its complementary filter, G_c(z^M), 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

$$G_c(z)=z^{-(L_G-1)/2}-G(z),$$E33

where L_G is the length of the filter G(z). Since the frequency response of these filters is periodical, two non-periodic masking filters, H_Ma(z) and H_Mc(z), are respectively cascaded with G(z^K) and G_c(z^K) 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_{Ma}(z)+[z^{-M(L_G-1)/2}-G(z^M)]H_{Mc}(z).$$E34

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(x)={\displaystyle \sum _{k=0}^{M}q(k)x^k}.$$E35

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$$\in$$X_p$$\cup$$X_s, where X_p is the range of values [x_p,l, x_p,u] and X_s is the range of values [x_s,l, x_s,u], 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(x)-{\delta }_p\le Q(x)\le D(x)+{\delta }_p\text{,}\text{for}{x}_{p,l}\le x\le x_{p,u},$$E36

$$D(x)-{\delta }_s\le Q(x)\le D(x)+{\delta }_s,\text{for}{x}_{s,l}\le x\le x_{s,u},$$E37

$$x_{p,l}\le G(\omega )\le x_{p,u},\text{for}\omega \in {\Omega }_p,$$E38

$$x_{s,l}\le G(\omega )\le x_{s,u},\text{for}\omega \in {\Omega }_s$$E39

Usually, D(x) = 1 for x$$\in$$X_p and D(x) = 0 for x$$\in$$X_s. 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 L_G (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)={\displaystyle \sum _{k=0}^{M}q(k)z^{-(M-k)(L_G-1)/2}{\left[G(z)\right]}^k}.$$E40

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 2N+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,

$$H(z)={\displaystyle \sum _{m=1}^M{z}^{-N_m}{H}_m(z)+z^{-\widehat{N}}\widehat{H}(z)},$$E41

where,

$$H_m(z)=h_m(N-N_m)z^{-(N-N_m)}+{\displaystyle \sum _{n=0}^{(N-N_m)-1}{h}_m(n)[z^{-n}+z^{-(2(N-N_m)-n)}]}.$$E42

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

$$h_m(n)={\displaystyle \sum _{r=0}^L{a}_m^{(L)}(r)n^r\times \sin [{\omega }_c(n-(N-N_m))]}.$$E43

In addition, $$z^{-\widehat{N}}\widehat{H}(z)$$is a conventional direct-form transfer function with non-zero impulse response coefficients$$\widehat{h}(n)$$, with n = N–c+1, N–c+2, …, N–c+T, where c =$$\lceil T/2\rceil$$ and T is the number of additional coefficients at the center of the filter. Given the values N_m 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₁(n) and x₂(n), that are filtered by two identical filters H(z). Thus, two corresponding sequences of independent outputs, y₁(n) and y₂(n), are obtained. An alternative form for this purpose is the multirate implementation using H(z²) 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₁(n) to the second input x₂(n). In this way, H(z²) is used to implement H²(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(z^K). 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(z^K) 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 half-band filter H_Hb(z) [26]. Consider a half-band filter H_a(z) as a band-edge shaping filter whose transfer function is given by

$$H_a(z)=\frac{1}{2}{z}^{-2K+1}+A(z),$$E44

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

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

$$H_{Hb}(z)=H_a(z^M)H_{Ma}(z)+[z^{-M(2K-1)}-H_a(z^M)]H_{Mc}(z).$$E46

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

$$H_{Hb}(z)=\left[\frac{1}{2}{z}^{-M(2K-1)}+A(z^M)\right]H_{Ma}(z)+\left[\frac{1}{2}{z}^{-M(2K-1)}-A(z^M)\right]H_{Mc}(z).$$E47

If M is odd, then either [$$\frac{1}{2}{z}^{-M(2K-1)}+A(z^M)$$] or[$$\frac{1}{2}{z}^{-M(2K-1)}-A(z^M)$$] has a transition band centered at π/2, just as desired for a half-band filter.

If M is given by the form

$$M=4k+1,\text{ with }k=\{0,1,2,\mathrm{...}\},$$E48

then the pass-band of H_Ma(z) is greater than the pass-band of H_Mc(z). With ω_p and ω_s denoting the band-edge frequencies of H_Hb(z), these values can be express as [26]

$${\omega }_p=\frac{2\pi m+\theta }{M},{\omega }_s=\frac{2\pi m+\varphi }{M},$$E49

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

$$m=\lfloor \frac{{\omega }_{p}M}{2\pi }\rfloor ,\theta ={\omega }_{p}M-2\pi m,\varphi ={\omega }_{s}M-2\pi m$$E50

where $$\lfloor x\rfloor$$represents the integer part of x, whereas θ and ϕ are the pass-band and stop-band edges of H_a(z). The pass-band and stop-band edge frequencies of the masking filter H_Ma(z), θ_Ma and ϕ_Ma, as well as the pass-band and stop-band edge frequencies of the masking filter H_Mc(z), θ_Mc and ϕ_Mc, are given by

$${\theta }_{Ma}={\omega }_p=\frac{2\pi m+\theta }{M},{\varphi }_{Ma}=\frac{2\pi (m+1)-\varphi }{M},{\theta }_{Mc}=\frac{2\pi m-\theta }{M},{\varphi }_{Mc}={\omega }_s=\frac{2\pi m+\varphi }{M}$$E51

If M is given by the form

$$M=4k+3,\text{ con }k=\{0,1,2,\mathrm{...}\}.$$E52

then the passband of H_Mc(z) is greater than the pass-band of H_Ma(z). In this case the frequencies ω_p and ω_s are given by

$${\omega }_p=\frac{2\pi m-\varphi }{M},{\omega }_s=\frac{2\pi m-\theta }{M},$$E53

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

$$m=\lceil \frac{{\omega }_{s}M}{2\pi }\rceil ,\theta =2\pi m-{\omega }_{s}M,\varphi =2\pi m-{\omega }_{p}M,$$E54

where $$\lceil x\rceil$$represents the rounding operation to the closest integer greater than x. The values, θ_Ma, ϕ_Ma, θ_Mc and ϕ_Mc, are given by

$${\theta }_{Ma}=\frac{2\pi (m-1)+\varphi }{M},{\varphi }_{Ma}={\omega }_s=\frac{2\pi m-\theta }{M},{\theta }_{Mc}={\omega }_p=\frac{2\pi m-\varphi }{M},{\varphi }_{Mc}=\frac{2\pi m+\theta }{M}.$$E55

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

$$H_{Ma}(z)=h_{Ma}\left((L_{Ma}-1)/2\right)+{\displaystyle \sum _{k=1}^{(L_{Ma}-1)/2}{h}_{Ma}\left(k+(L_{Ma}-1)/2\right)\left[z^{k-(L_{Ma}-1)/2}+z^{-k-(L_{Ma}-1)/2}\right]}.$$E56

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

$$\begin{array}{l}B(z)=h_{Ma}\left(1+(L_{Ma}-1)/2\right)\left[z^{1-(L_{Ma}-1)/2}+z^{-1-(L_{Ma}-1)/2}\right]+\\ h_{Ma}\left(3+(L_{Ma}-1)/2\right)\left[z^{3-(L_{Ma}-1)/2}+z^{-3-(L_{Ma}-1)/2}\right]+\mathrm{...}\end{array}$$E57

$$\begin{array}{l}C(z)=h_{Ma}\left((L_{Ma}-1)/2\right)+\\ h_{Ma}\left(2+(L_{Ma}-1)/2\right)\left[z^{2-(L_{Ma}-1)/2}+z^{-2-(L_{Ma}-1)/2}\right]+\\ h_{Ma}\left(4+(L_{Ma}-1)/2\right)\left[z^{4-(L_{Ma}-1)/2}+z^{-4-(L_{Ma}-1)/2}\right]+\mathrm{...}\end{array}$$E58

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

$$H_{Ma}(z)=B(z)+C(z).$$E59

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

$$H_{Mc}(z)=z^{-(L_{Ma}-1)/2}-H_{Ma}(-z).$$E60

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).$$E61

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_{Hb}(z)=\frac{1}{2}{z}^{-[M(2K-1)+(L_{Ma}-1)/2]}+z^{-M(2K-1)}B(z)+A(z^M)[2C(z)-z^{-(L_{Ma}-1)/2}].$$E62

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 }{{\omega }_s-{\omega }_p}}.$$E63

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

1. Find if M can be expressed either as in (48) or (52). Then obtain the pass-band and stop-band edge frequencies of the band-edge shaping filter H_a(z), θ and ϕ, as well as the ones of the masking filter H_Ma(z), θ_Ma and ϕ_Ma, using (50) and (51) if M is expressed as in (48), or (54) and (55) if M is expressed as in (52).Design these filters with a ripple approx. 15% less than the desired ripple.

2. Obtain A(z) from H_a(z) using (44) and (45). Then obtain B(z) and C(z) from H_Ma(z) using (56), (57) and (58).

3. Synthesize the overall structure of (62) in terms of A(z^M), B(z) and C(z).

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^–jπn/2 (see section 3.1). The transfer function of the Hilbert transformer H(z) is given by [9]

$$H(z)=2{(jz)}^{-M(2K-1)}B(jz)+2A\left({(jz)}^M\right)[2C(jz)-{(jz)}^{-(L_{Ma}-1)/2}].$$E64

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(j^Mz^M) has only real coefficients since the imaginary unit generated by (jz)^–n with n odd is eliminated by zero-valued coefficients in these indexes n. Second, note that all the coefficients in [2C(jz) – (jz)^{–(LMa–1)/2}] are real and all the coefficients in B(jz) are imaginary when L_Ma is expressed as 4k+1, with k integer. Third, the term (jz)^–M(2K–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)^–M(2K–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 H_b(z) = 2(jz)^–M(2K–1)B(jz) and in the second branch we have the cascade of H₁(z^M) and H_M(z), where H₁(z^M) = 2A(j^Mz^M) and H_M(z) = [2C(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 H_b(z) can be seen as a low-order Hilbert transformer, the filter H₁(z^M) as a band-edge shaping filter and H_M(z) as a masking filter. The basic filter H_b(z) provides a low order approximation (with wide transition bandwidth) to the desired specification. The cascaded connection of H₁(z^M) and H_M(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(z)=H_1(z^M)H_M(z)+H_b(z).$$E65

Let the lengths of H_b(z), H₁(z) and H_M(z) be L_b, L₁ and L_M, respectively. The length of H₁(z^M)H_M(z) is ML₁ + L_M – M. The delay introduced by H_b(z) and the delay introduced by H₁(z^M)H_M(z) must be the same; otherwise, pure delays must be introduced into the shorter-delay 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 L_b and of (ML₁ + L_M – M) must be the same. Furthermore, H₁(z^M)H_M(z)must have anti-symmetrical impulse response.

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

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

Since the length of H_b(z) is even, the length of H₁(z^M) H_M(z), i.e.,MN₁ + N_M – M, must also be even. If M is odd, N₁andN_Mmust have different parities. By considering the gain of |H_M(e^jω)| in the vicinity of ω = 0, it is clear that H_M(z) has symmetrical impulse response. Thus, H₁(z) must have anti-symmetrical impulse response to satisfy the condition that H₁(z^M) H_M(z) must have anti-symmetrical impulse response.

The band-edges of H_b(z) and H_M(z) are the same. Let the band-edge of H_b(z) be θ_b and let the band-edge of H₁(z) be θ₁. It can be seen from Figure 12 that the value θ_b satisfies θ_b ≤ (2π – θ₁)/M. For an arbitrary value θ₁, 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₁(z^M), H_b(z) and H_M(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 12.

Magnitude responses of the subfilters for even length H_b(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.

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^Mz^M), 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 L_P = 2N and anti-symmetric impulse response of the form p(2N – 1 – n) = –p(n). Its frequency response is expressed as

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

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

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

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

Using the equivalence cos(Ωn) = T_n{cos(Ω)} [17], where T_n{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\left\{x\right\}=1,T_1\left\{x\right\}=x\text{and}{T}_n\left\{x\right\}=2x{T}_{n-1}\left\{x\right\}-T_{n-2}\left\{x\right\},$$E68

the zero-phase term can be rewritten as

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

where α(n) are obtained from $$\tilde{d}(n)$$using the coefficients of the Chebyshev polynomials. Based on the equivalence given as

$$\cos \left(2x\right)=1-2{\sin }^2\left(x\right)=1+2{\left(j\cdot \sin x\right)}^2,$$E70

the zero-phase term can be expressed by

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

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

$$G(e^{j\omega })=e^{-j\left(2M\omega /2\right)}G(\omega ),$$E72

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

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

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 {\displaystyle \sum _{n=1}^{M}c(n)\sin (\omega n)},$$E74

resulting in

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

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[{\displaystyle \sum _{n=1}^{M}c(n)\sin (\omega n)}\right].$$E76

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){\displaystyle \sum _{n=0}^{N-1}{z}^{-2M(N-1-n)}\alpha (n){\left[H_1(z)\right]}^n},H_1(z)=z^{-2M}+2G^2(z)$$E77

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

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

$$(1-\delta )\le \left|P(\Omega )\right|\le (1+\delta ),\text{ for }{\Omega }_L\le \Omega \le \pi ,$$E78

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

$$v_d-{\delta }_G\le \left|G_0(\omega )\right|\le 1,\text{ for }{\omega }_L\le \omega \le \pi -{\omega }_L,$$E79

$$v_d=\frac{1}{2}+\frac{1}{2}\sin \left(\frac{{\Omega }_L}{2}\right),{\delta }_G=\frac{1}{2}-\frac{1}{2}\sin \left(\frac{{\Omega }_L}{2}\right).$$E80

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₁(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.

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_G\approx {\text{L}}_G=\varphi ({\delta }_G/v_d,{\omega }_L),$$E81

$$L_P\approx {\text{L}}_P=\varphi (\delta ,{\Omega }_L),$$E82

Figure 16.

PI-based structure with one subfilter [15].

where L_G and L_P are the respective approximations to the lengths of the subfilter, L_G, and the prototype filter, L_P, whereas v_d and δ_G are given in (80). Clearly, L_G and L_P 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_G,m_P),$$E83

$$m_G\approx C\cdot L_G,$$E84

$$m_P=L_P/2,$$E85

where m_P and m_G are taken from (19).

Substituting v_d and δ_G from (80) in (81), and using the approximations (81) and (82) respectively in (84) and (85) we have

$$m(\delta ,{\omega }_L,{\Omega }_L)\approx f\left(C\cdot \varphi (\frac{1-\sin ({\Omega }_L/2)}{1+\sin ({\Omega }_L/2)},{\omega }_L),\frac{\text{1}}{\text{2}}\cdot \varphi (\delta ,{\Omega }_L)\right).$$E86

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 Ω^*_L for Ω_L, with 0 < Ω_L<π, such that m(δ, ω_L, Ω_L) is minimized. This optimization problem is given as

$$\begin{array}{l}\min \underset{{\Omega }_L\in ℝ}{m(\delta ,{\omega }_L,{\Omega }_L)}\\ \text{such that 0}<{\Omega }_L<\pi ,\end{array}$$E87

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{array}{l}\varphi (\delta ,{\omega }_L)=\frac{1}{2}+\left[\frac{1.101{\left[-{\log }_{10}(\delta )\right]}^{1.1}}{(\frac{{\omega }_L}{2\pi })}+1\right]\cdot [\frac{2}{3\pi }\arctan \{[2.325{\left[0.30103-{\log }_{10}(\delta )\right]}^{-0.445}\cdot \mathrm{...}\\ {(\frac{{\omega }_L}{2\pi })}^{-1.39}]\left[\frac{1}{0.5-(\frac{{\omega }_L}{2\pi })}\right]\}+\frac{1}{6}].\end{array}$$E88

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 Ω^*_L = 0.2237π and the lengths for the prototype filter and the subfilter are, respectively, L_P = 14 and L_G = 31. The value Ω^*_L 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].

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.

5.3. Hilbert transformer design based on Piecewise Polynomial Sinusoidal technique

This method to synthesize Hilbert transformers, first proposed in [25] and then detailed in [11], is based on the previous method [31] and on a modification to the method [32]-[33]. The method stems from a windowing technique for Type III and Type IV FIR filters of order 2N and 2N-1, respectively. The basic windowing technique for Type III filters is expressed as

$$h(n)=w(n)h_0(n),$$E89

where

$$h_0(n)=\{\begin{array}{l}0\text{ for }n\text{ odd}\\ \frac{1-\cos ((n-N)\pi )}{(n-N)\pi }\text{ for n even}\end{array}$$E90

is an ideal impulse response, w(n) is the window function and N is odd. Similarly, an ideal odd-order impulse response can be expressed as

$$h_0(n)=\frac{1}{(n-(2N-1)/2)\pi }.$$E91

The overall transfer function can be expressed in the following way,

$$F(z)=\{\begin{array}{l}{\displaystyle \sum _{n=0}^{2N}W(n)[1-\cos ((n-N)\pi )]z^{-n},}\text{ Type III}\\ {\displaystyle \sum _{n=0}^{2N-1}W(n)z^{-n},\text{ Type IV}}\end{array}$$E92

with

$$W(n)=\{\begin{array}{l}w(n-N)/[\pi (n-N)],\\ w(n-(2N-1/2))/[\pi (n-(2N-1/2))].\end{array}$$E93

W(n) is an impulse response satisfying W(2N − n) = −W(n) for n = 0, 1,...,N − 1, W(N) = 0, and W(2N − 1 − n) = −W(n) for n = 0, 1,...,N, for Types III and IV respectively. Additionally, w(n) is a positive and symmetric window function for all n.

A way to generate the piecewise-polynomial-sinusoidal transfer function F(z) is to consider the following transfer function:

$$G(z)={\displaystyle \sum _{n=0}^{N_s}W(n)z^{-n}(e^{j(n-N)\pi })}={\displaystyle \sum _{n=0}^{N_s}W(n)\left[1-(\cos ((n-N)\pi )+j\sin ((n-N)\pi ))\right]z^{-n}},$$E94

where N_s is 2N or 2N − 1 for Types III and IV, respectively. From (92) it follows that F(z) can be expressed as

$$F(z)=\mathrm{Re}\left\{{\displaystyle \sum _{n=0}^{N_s}W(n)z^{-n}+z^{-n}(e^{j(n-N)\pi })}\right\}=\mathrm{Re}\left\{{\displaystyle \sum _{n=0}^{N_s}W(n)\left[1-z^{-N\pi }{(z^{-1}{e}^{j\pi })}^n\right]}\right\}.$$E95

In [32] a wideband FIR filter is obtained by first generating the envelope filter with W(n) for n =0, 1,..., 2N [2N − 1] as the impulse-response coefficients. This impulse response is designed to become piecewise polynomial. The coefficient values of this envelope filter are modified by multiplying them with [1−cos(n−N)π] for Type III, which gives F(z). From (92) it can be seen that for Type IV the piecewise-polynomial-sinusoidal impulse response coincides with the piecewise-polynomial impulse response, i.e., the constant part of the real part in (95).

Let us consider the simpler case Type IV filter. The overall transfer function, denoted by H(z), is constructed as presented previously in sub-section 4.3. We will repeat the equation here for convenience,

$$H(z)={\displaystyle \sum _{m=1}^M{z}^{-N_m}{H}_m(z)+z^{-\widehat{N}}\widehat{H}(z)}.$$E96

Recall that the integers N_m in the delay terms z^–Nm satisfy N₁ = 0 and N_m+1>N_m for m = 1, 2, …, M – 1, and the order of H_m(z) is 2(N – N_m) – 1. The impulse response is given for n = 0, 1, 2, …,N − N_m–1by Lth order polynomials as follows

$$h_m(n)={\displaystyle \sum _{r=0}^L{a}_m^{(L)}(r)n^r}.$$E97

Additionally, $$z^{-\widehat{N}}\widehat{H}(z)$$is a conventional [2N − 1]th-order Type IV direct-form transfer function with the additional impulse response coefficients at n = N −c,..., N −1, where c =$$\lceil T/2\rceil$$ and T is the number of additional coefficients at the center of the filter. The delay terms in (96) are used to shift the center of the symmetry at the desired location, which occurs at n = (2N − 1)/2.

In order to indicate that the overall filter has a piecewise-polynomial impulse response the time interval [0,N − 1] is divided into the following M subintervals:

$$X_M=[N_m,N_{m+1}-1]\text{ for }m=1,2,\mathrm{...},M-1$$E98

and

$$X_M=[N_M,N].$$E99

First, we have thatX₁ =[0,N₂− 1] because N₁ = 0, Secondly, the overall impulse response can be studied up to n = N − 1 because of the odd symmetry. The impulse response on X_m can be expressed as

$$h(n)={\displaystyle \sum _{m=1}^M{\tilde{h}}_m(n)},$$E100

where

$${\tilde{h}}_M(n)={\displaystyle \sum _{k=1}^m\left[{\displaystyle \sum _{r=0}^L\left[a_k^{(L)}(r){(n-N_M)}^r\right]}\right]}$$E101

for m = 1, 2, …, M – 1 and

$${\tilde{h}}_M(n)={\displaystyle \sum _{k=1}^M\left[{\displaystyle \sum _{r=0}^L\left[a_k^{(L)}(r){(n-N_M)}^r\right]}\right]+h´(n)},$$E102

which equals the overall impulse response and where h’(n)is a conventional direct-form Type IV filter with nonzero coefficients for n = N − c,..., N − 1, in which c = $$\lceil T/2\rceil$$and T is the number of separately generated additional center coefficients. The slices N_ms should be chosen so that |N₂−N₁| ≠|N₃−N₂| ≠ …≠|N_M −N_M−1|,where N₁ = 0 and M is the number of subintervals in the overall impulse response.

Based on the above equations, in each X_m for m =1, 2,..., M, a separate piecewise-polynomial impulse response can be generated. In addition, in the X_M, there are additional center coefficients, which are of great importance for fine-tuning the overall filter to meet the given criteria.

Given the filter criteria as well as the design parameters M,N, L, N_m’s, and the number of center coefficients included in$$\widehat{H}(z)$$, the overall problem is solvable by using linear programming.

6. Conclusion

In this chapter we have studied the Hilbert transform relations existing among the real part and the imaginary part of complex analytic signals. The importance of these signals has been highlighted in terms of spectral efficiency, i.e., the analytic signals do not have spectral components in their negative-frequency side. For discrete-time sequences, this characteristic holds for the negative-frequency side in every Nyquist period.

The Hilbert transformer has been introduced as a special type of FIR filter which is the key processing system to generate analytic signals. The design of such important filter is, of course, straightforward with the aid of an important filter design tool: the MATLAB Signal Processing Toolbox. However, this direct design method, shown as a very simple and convenient MATLAB code, cannot be efficiently applied for more stringent and realistic specifications. We have presented a concise explanation of the relation of Hilbert transformers and half-band filters because this relation, as has been observed from literature, is one of the most important characteristics to overcome this problem.

The efficient methods to design low-complexity FIR Hilbert transformers with strict specifications have been detailed. Three methods have been analyzed, namely, Frequency-Response Masking (FRM), Frequency Transformation (FT) and Piecewise-Polynomial Sinusoidal (PPS). These schemes are based on three different approaches to design efficient FIR filtering. FRM is a periodical subfilter based method, FT is an identical subfilter based method and PPS is a piecewise-polynomial based method. Additionally, it has been observed that FRM and PPS are fully parallel approaches and do not take direct advantage of hardware multiplexing. On the other hand, we have shown that FT allows area-efficient architectures by multiplexing a simple subfilter.

Finally, the FRM and the time-multiplexed FT approach have been illustrated in MATLAB, with the aid of the Signal Processing Toolbox. Even though the underlying theory on the efficient techniques to design FIR Hilbert transformers is specialized, the MATLAB codes have been preserved in a simple and as clear as possible presentation. The presented codes allow a clearer understanding on such specialized techniques and, at the same time, can serve as a basis for more elaborated algorithms and further research on this fertile area.

Acknowledgement

Special thanks to Miriam G. Cruz Jimenez for her valuable assistance in the development of this chapter.