DSPRelated.com
Free Books
   Free Books    Spectral Audio Signal Processing  

Hilbert Transform Design Example

We will now use the window method to design a complex bandpass filter which passes positive frequencies and rejects negative frequencies.

Since every real signal $ x(n)$ possesses a Hermitian spectrum $ X(\omega)$ , i.e., $ X(-\omega) = \overline{X(\omega)}$ , it follows that, if we filter out the negative frequencies, we will destroy this spectral symmetry, and the output signal will be complex for every nonzero real input signal (excluding dc and half the sampling rate). In other terms, we want a filter which produces a ``single sideband'' (SSB) output signal in response to any real input signal. The Hermitian spectrum of a real signal can be viewed as two sidebands about dc (with one sideband being the ``conjugate flip'' of the other). See §2.3.3 for a review of Fourier symmetry-properties for real signals.

An ``analytic signal'' in signal processing is defined as any signal $ x(n)$ having only positive or only negative frequencies, but not both (typically only positive frequencies). In principle, the imaginary part of an analytic signal is computed from its real part by the Hilbert transform (defined and discussed below). In other words, one can ``filter out'' the negative-frequency components of a signal $ x(n)$ by taking its Hilbert transform $ y(n) = {\cal H}_n\{x\}$ and forming the analytic signal $ x_a(n) = x(n) + j y(n)$ . Thus, an alternative problem specification is to ask for a (real) filter which approximates the Hilbert transform as closely as possible for a given filter order.

Primer on Hilbert Transform Theory

We need a Hilbert-transform filter $ {\cal H}$ to compute the imaginary part $ y(t)$ of the analytic signal $ x_a(t)$ given its real part $ x(t)$ . That is,

$\displaystyle x_a(t) \eqsp x(t) + j\,y(t)$ (5.14)

where $ y = {\cal H}\{x\}$ . In the frequency domain, we have

$\displaystyle X_a(\omega) \eqsp X(\omega) + j\,Y(\omega) \eqsp \left[1+j\,H(\omega)\right]\,X(\omega)$ (5.15)

where $ H(\omega)$ denotes the frequency response of the Hilbert transform $ {\cal H}$ . Since by definition we have $ X_a(\omega)=0$ for $ \omega<0$ , we must have $ j\,H(\omega)=-1$ for $ \omega<0$ , so that $ H(\omega)=j$ for negative frequencies (an allpass response with phase-shift $ +90$ degrees). To pass the positive-frequency components unchanged, we would most naturally define $ H(\omega)=0$ for $ \omega>0$ . However, conventionally, the positive-frequency Hilbert-transform frequency response is defined more symmetrically as $ H(\omega)=-j$ for $ \omega>0$ , which gives $ j\,H(\omega)=+1$ and $\mbox{re\ensuremath{\left\{X_a(\omega)\right\}}}=2X(\omega)$ , i.e., the positive-frequency components of $ X$ are multiplied by $ 2$ .

In view of the foregoing, the frequency response of the ideal Hilbert-transform filter may be defined as follows:

$\displaystyle H(\omega) \isdefs \left\{\begin{array}{ll} \quad\! j, & \omega<0 \\ [5pt] \quad\!0, & \omega=0 \\ [5pt] -j, & \omega>0 \\ \end{array} \right. \protect$ (5.16)

Note that the point at $ \omega = 0$ can be defined arbitrarily since the inverse-Fourier transform integral is not affected by a single finite point (being a ``set of measure zero'').

The ideal filter impulse response $ h(n)$ is obtained by finding the inverse Fourier transform of (4.16). For discrete time, we may take the inverse DTFT of (4.16) to obtain the ideal discrete-time Hilbert-transform impulse response, as pursued in Problem 10. We will work with the usual continuous-time limit $ h(t)=1/(\pi t)$ in the next section.

Hilbert Transform

The Hilbert transform $ y(t)$ of a real, continuous-time signal $ x(t)$ may be expressed as the convolution of $ x$ with the Hilbert transform kernel:

$\displaystyle h(t) \isdefs \frac{1}{\pi t},\qquad t\in(-\infty,\infty) \protect$ (5.17)

That is, the Hilbert transform of $ x$ is given by

$\displaystyle y = h \ast x. % \qquad \hbox{(Hilbert transform of $x$)}. $ (5.18)

Thus, the Hilbert transform is a non-causal linear time-invariant filter.

The complex analytic signal $ x_a(t)$ corresponding to the real signal $ x(t)$ is then given by

$\displaystyle x_a(t)$ $\displaystyle \isdef$ $\displaystyle x(t) + j y(t)$  
$\displaystyle %\qquad \hbox{(Analytic signal corresponding to $x$)}\\ [10pt] $ $\displaystyle =$ $\displaystyle \frac{1}{\pi}\displaystyle\int_0^{\infty} X(\omega)e^{j\omega t}\, d\omega$ (5.19)

To show this last equality (note the lower limit of 0 instead of the usual $ -\infty$ ), it is easiest to apply (4.16) in the frequency domain:

$\displaystyle X_a(\omega)$ $\displaystyle \isdef$ $\displaystyle X(\omega) + j Y(\omega)$  
  $\displaystyle \isdef$ $\displaystyle (X_++X_-) + j (-j X_+ + j X_-)$ (5.20)
  $\displaystyle =$ $\displaystyle 2X_+(\omega)$ (5.21)

Thus, the negative-frequency components of $ X$ are canceled, while the positive-frequency components are doubled. This occurs because, as discussed above, the Hilbert transform is an allpass filter that provides a $ 90$ degree phase shift at all negative frequencies, and a $ -90$ degree phase shift at all positive frequencies, as indicated in (4.16). The use of the Hilbert transform to create an analytic signal from a real signal is one of its main applications. However, as the preceding sections make clear, a Hilbert transform in practice is far from ideal because it must be made finite-duration in some way.

Filtering and Windowing the Ideal
Hilbert-Transform Impulse Response

Let $ h_i(t)$ denote the convolution kernel of the continuous-time Hilbert transform from (4.17) above:

$\displaystyle h_i(t) = \frac{1}{\pi t} \protect$ (5.22)

Convolving a real signal $ x(t)$ with this kernel produces the imaginary part of the corresponding analytic signal. The way the ``window method'' for digital filter design is classically done is to simply sample the ideal impulse response to obtain $ h_i(nT)$ and then window it to give $ \hat{h}_i(nT)=w(n)h_i(nT)$ . However, we know from above (e.g., §4.5.2) that we need to provide transition bands in order to obtain a reasonable design. A single-sideband filter needs a transition band between dc and $ \Omega_w$ , or higher, where $ \Omega_w$ denotes the main-lobe width (in rad/sample) of the window $ w$ we choose, and a second transition band is needed between $ \pi-\Omega_w$ and $ \pi$ .

Note that we cannot allow a time-domain sample at time 0 in (4.22) because it would be infinity. Instead, time 0 should be taken to lie between two samples, thereby introducing a small non-integer advance or delay. We'll choose a half-sample delay. As a result, we'll need to delay the real-part filter by half a sample as well when we make a complete single-sideband filter.

The matlab below illustrates the design of an FIR Hilbert-transform filter by the window method using a Kaiser window. For a more practical illustration, the sampling-rate assumed is set to $ f_s=22,050$ Hz instead of being normalized to 1 as usual. The Kaiser-window $ \beta $ parameter is set to $ 8$ , which normally gives ``pretty good'' audio performance (cf. Fig.3.28). From Fig.3.28, we see that we can expect a stop-band attenuation better than $ 50$ dB. The choice of $ \beta $ , in setting the time-bandwidth product of the Kaiser window, determines both the stop-band rejection and the transition bandwidths required by our FIR frequency response. 

  M = 257;    % window length = FIR filter length (Window Method)
  fs = 22050; % sampling rate assumed (Hz)
  f1 = 530;   % lower pass-band limit = transition bandwidth (Hz)
  beta = 8;   % beta for Kaiser window for decent side-lobe rejection
Recall that, for a rectangular window, our minimum transition bandwidth would be $ 2f_s/M \approx 172$ Hz, and for a Hamming window, $ 4f_s/M \approx 343$ Hz. In this example, using a Kaiser window with $ \beta=8$ ( $ \alpha=\beta/\pi\approx 2.5$ ), the main-lobe width is on the order of $ 5f_s/M \approx 430$ Hz, so we expect transition bandwidths of this width. The choice $ f_1=530$ above should therefore be sufficient, but not ``tight''.5.8 For each doubling of the filter length (or each halving of the sampling rate), we may cut $ f_1$ in half.

Matlab, Continued

Given the above design parameters, we compute some derived parameters as follows: 

  fn = fs/2;             % Nyquist limit (Hz)
  f2 = fn - f1;          % upper pass-band limit
  N = 2^(nextpow2(8*M)); % large FFT for interpolated display
  k1 = round(N*f1/fs);   % lower band edge in bins
  if k1<2, k1=2; end;    % cannot have dc or fn response
  kn = N/2 + 1;          % bin index at Nyquist limit (1-based)
  k2 = kn-k1+1;          % high-frequency band edge
  f1 = k1*fs/N           % quantized band-edge frequencies
  f2 = k2*fs/N
Setting the upper transition band the same as the low-frequency band ( $ f_2 = f_s/2 - f_1$ ) provides an additional benefit: the symmetry of the desired response about $ f_s/4$ cuts the computational expense of the filter in half, because it forces every other sample in the impulse response to be zero [224, p. 172].5.9

Kaiser Window

With the filter length $ M=257$ and Kaiser window $ \beta=8$ as given above, we may compute the Kaiser window itself in matlab via 

  w = kaiser(M,beta)'; % Kaiser window in "linear phase form"
The spectrum of this window (zero-padded by more than a factor of 8) is shown in Fig.4.9 (full magnitude spectrum) and Fig.4.10 (zoom-in on the main lobe).

Figure 4.9: The Kaiser window magnitude spectrum.
\includegraphics[width=0.8\twidth]{eps/KaiserFR}

Figure 4.10: Kaiser window magnitude spectrum near dc.
\includegraphics[width=0.8\twidth]{eps/KaiserZoomFR}

Windowing a Desired Impulse Response Computed by the
Frequency Sampling Method

The next step is to apply our Kaiser window to the ``desired'' impulse response, where ``desired'' means a time-shifted (by 1/2 sample) and bandlimited (to introduce transition bands) version of the ``ideal'' impulse response in (4.22). In principle, we are using the frequency-sampling method (§4.4) to prepare a desired FIR filter of length $ N=4096$ as the inverse FFT of a desired frequency response prepared by direct Fourier intuition. This long FIR filter is then ``windowed'' down to length $ M=257$ to give us our final FIR filter designed by the window method.

If the smallest transition bandwidth is $ f_1$ Hz, then the FFT size $ N$ should satisfy $ N\gg f_s/f_1$ . Otherwise, there may be too much time aliasing in the desired impulse response.5.10 The only non-obvious part in the matlab below is ``.^8'' which smooths the taper to zero and looks better on a log magnitude scale. It would also make sense to do a linear taper on a dB scale which corresponds to an exponential taper to zero. 

H = [ ([0:k1-2]/(k1-1)).^8,ones(1,k2-k1+1),...
      ([k1-2:-1:0]/(k1-1)).^8, zeros(1,N/2-1)];
Figure 4.11 shows our desired amplitude response so constructed.

Figure 4.11: Desired single-sideband-filter frequency response.
\includegraphics[width=0.8\twidth]{eps/hilbertFRIdeal}

Now we inverse-FFT the desired frequency response to obtain the desired impulse response: 

h = ifft(H); % desired impulse response
hodd = imag(h(1:2:N));     % This should be zero
ierr = norm(hodd)/norm(h); % Look at numerical round-off error
% Typical value: ierr = 4.1958e-15
% Also look at time aliasing:
aerr = norm(h(N/2-N/32:N/2+N/32))/norm(h);
% Typical value: 4.8300e-04
The real part of the desired impulse response is shown in Fig.4.12, and the imaginary part in Fig.4.13.

Figure 4.12: Real part of the desired single-sideband-filter impulse response.
\includegraphics[width=0.8\twidth]{eps/hilbertIRRealIdeal}

Figure 4.13: Desired Hilbert-transform-filter impulse response.
\includegraphics[width=0.8\twidth]{eps/hilbertIRImagIdeal}

Now use the Kaiser window to time-limit the desired impulse response: 

% put window in zero-phase form:
wzp = [w((M+1)/2:M), zeros(1,N-M), w(1:(M-1)/2)];
hw = wzp .* h; % single-sideband FIR filter, zero-centered
Hw = fft(hw);  % for results display: plot(db(Hw));
hh = [hw(N-(M-1)/2+1:N),hw(1:(M+1)/2)]; % caual FIR
% plot(db(fft([hh,zeros(1,N-M)]))); % freq resp plot

Figure 4.14 and Fig.4.15 show the normalized dB magnitude frequency response of our final FIR filter consisting of the $ M$ nonzero samples of hw.

Figure 4.14: Frequency response of the Kaiser-windowed impulse response.
\includegraphics[width=0.8\twidth]{eps/KaiserHilbertFR}

Figure 4.15: Zoomed frequency response of the Kaiser-windowed impulse response.
\includegraphics[width=0.8\twidth]{eps/KaiserHilbertZoomedFR}

More General FIR Filter Design

We have looked at more than just FIR filter design by the window method and frequency-sampling technique. The general steps were

  1. Prepare a desired frequency-response that ``seems achievable''
  2. Inverse FFT
  3. Window the result (time-limit it)
  4. FFT that to see how it looks
In general, the end result will be a somewhat smoothed version of what we originally asked for in the frequency domain. However, this smoothing will be minimal if we asked for a truly ``doable'' desired frequency response. Because the above steps are fast and non-iterative, they can be used effectively in response to an interactive user interface (such as a set of audio graphic-equalizer sliders that are interpolated to form the desired frequency response), or even in a real-time adaptive system.


Comparison to Optimal Chebyshev FIR Filter

Let's now compare the window-method design using the Kaiser window to the optimal equiripple FIR filter design given by the Remez multiple exchange algorithm.

Note, by the way, that many filter-design software functions, such as firpm have special modes for designing Hilbert-transform filters [224].

It turns out that the Remez exchange algorithm has convergence problems for filters larger than a few hundred taps. Therefore, the FIR filter length $ M=257$ was chosen above to be small enough to work out in this comparison. However, keep in mind that for very large filter orders, the Remez exchange method may not be an option. There are more recently developed methods for optimal Chebyshev FIR filter design, using ``convex optimization'' techniques, that may continue to work at very high orders [218,22,153]. The fast nonparametric methods discussed above (frequency sampling, window method) will work fine at extremely high orders.

The following Matlab command will try to design the FIR Hilbert-transform filter of the desired length with the desired transition bands: 

 hri = firpm(M-1, [f1,f2]/fn, [1,1], [1], 'Hilbert');
Instead, however, we will use a more robust method [228] which uses the Remez exchange algorithm to design a lowpass filter, followed by modulation of the lowpass impulse-response by a complex sinusoid at frequency $ f_s/4$ in order to frequency-shift the lowpass to the single-sideband filter we seek: 
tic; % remember the current time
hrm = firpm(M-1, [0,(f2-fs/4)/fn,0.5,1], [1,1,0,0], [1,10]);
dt = toc; % design time dt can be minutes
hr = hrm .* j .^ [0:M-1]; % modulate lowpass to single-sideband
The weighting [1,10] in the call to firpm above says ``make the pass-band ripple $ 10$ times that of the stop-band.'' For steady-state audio spectra, pass-band ripple can be as high as $ 0.1$ dB or more without audible consequences.5.11 The result is shown in Fig.4.16 (full amplitude response) and Fig.4.17 (zoom-in on the dc transition band). By symmetry the high-frequency transition region is identical (but flipped):

Figure 4.16: Frequency response of the optimal Chebyshev FIR filter designed by the Remez exchange algorithm.
\includegraphics[width=0.8\twidth]{eps/OptimalHilbertFR}

Figure 4.17: Transition region of the optimal Chebyshev frequency response.
\includegraphics[width=0.8\twidth]{eps/OptimalHilbertZoomedFR}

In this case we did not normalize the peak amplitude response to 0 dB because it has a ripple peak of about 1 dB in the pass-band. Figure 4.18 shows a zoom-in on the pass-band ripple.

Figure 4.18: Pass-Band ripple for optimal Chebyshev frequency response.
\includegraphics[width=0.8\twidth]{eps/OptimalHilbertZoomedPB}

Conclusions

We can note the following points regarding our single-sideband FIR filter design by means of direct Fourier intuition, frequency-sampling, and the window-method:

  • The pass-band ripple is much smaller than 0.1 dB, which is ``over designed'' and therefore wasting of taps.

  • The stop-band response ``droops'' which ``wastes'' filter taps when stop-band attenuation is the only stop-band specification. In other words, the first stop-band ripple drives the spec ($ -80$ dB), while all higher-frequency ripples are over-designed. On the other hand, a high-frequency ``roll-off'' of this nature is quite natural in the frequency domain, and it corresponds to a ``smoother pulse'' in the time domain. Sometimes making the stop-band attenuation uniform will cause small impulses at the beginning and end of the impulse response in the time domain. (The pass-band and stop-band ripple can ``add up'' under the inverse Fourier transform integral.) Recall this impulsive endpoint phenomenon for the Chebyshev window shown in Fig.3.33.

  • The pass-band is degraded by early roll-off. The pass-band edge is not exactly in the desired place.

  • The filter length can be thousands of taps long without running into numerical failure. Filters this long are actually needed for sampling rate conversion [270,218].

We can also note some observations regarding the optimal Chebyshev version designed by the Remez multiple exchange algorithm:

  • The stop-band is ideal, equiripple.

  • The transition bandwidth is close to half that of the window method. (We already knew our chosen transition bandwidth was not ``tight'', but our rule-of-thumb based on the Kaiser-window main-lobe width predicted only about $ 23$ % excess width.)

  • The pass-band is ideal, though over-designed for static audio spectra.

  • The computational design time is orders of magnitude larger than that for window method.

  • The design fails to converge for filters much longer than 256 taps. (Need to increase working precision or use a different method to get longer optimal Chebyshev FIR filters.)

Practical application of Hilbert-transform filters to phase-vocoder analysis for additive synthesis is discussed in §G.10.

Next Section:
Generalized Window Method
Previous Section:
Window Method for FIR Filter Design