Winflt 4up

FIR Digital Filter Design MUS421/EE367B Lecture 6 FIR Digital Filter Design Julius O. Smith III ([email protected]) Center for Computer Research in Music and Acoustics (CCRMA) Department of Music, Stanford University Stanford, California 94305 March 28, 2005

In a previous lecture, we looked at many windows, and examined their properties. Today, we are going to see how these windows can be used to design Finite Impulse Response (FIR) digital filters. FFT processors implement long FIR filters more efficiently than any other method (using Overlap-Add). We need flexible ways to design all kinds of FIR filters for use in FFT processors.

Outline • Ideal Lowpass Filter • Optimal Least Squares in Time Domain • Window Method • Optimal Methods • Case Study: FIR Hilbert-Transform Design

1

2

Ideal Lowpass Filter

Impulse Response of Ideal LPF

Amplitude Response:

∆

Hideal hideal(n) = DTFT−1 n Z ωc ∆ 1 ejωndω = 2π −ωc

1




Gain

sin(ωcn) πn = 2fcsinc(2fcn), =

0

0

Frequency

fc

fs/2

n∈Z

Problems: • Infinitely long

• Gain = 1 for f < fc

• Non-causal

• Gain = 0 for f > fc

• Cannot be shifted to make it causal • Cannot be implemented in practice Conclusion: • We must accept some compromise(s) in the design of any practical lowpass filter. • Managing such trade-offs is the topic of digital filter design. 3

4

Optimal (but Poor) Least-Squares Impulse Response Design

Length L = 30 Least-Squares LPF Design

Length 30 FIR Amplitude Response

Let

Sum of squared errors: L−1 ∞ 2 2 X X ∆ ˆ ˆ ˆ = h(n) − h(n) = h(n) − h(n) J2(h) + c2 n=−∞ ∆

n=0

P0

P∞

2

2

where c2 = n=−∞ |h(n)| + n=L |h(n)| does not ˆ Note that J2(h) ˆ ≥ c2 . depend on h.

Result: The error is minimized (in the least-squares sense) by simply matching the first L terms in the desired impulse response.

1

0.5

0

0

0.05

0.1

0.15 0.2 0.25 0.3 0.35 Normalized Frequency (cycles/sample)

0.4

0.45

0

0.05

0.1

0.15 0.2 0.25 0.3 0.35 Normalized Frequency (cycles/sample)

0.4

0.45

0.5

0 −10 Magnitude (dB)

• h(n) = ideal filter impulse response. ˆ • h(n) = length L causal FIR filter (to be designed).

Magnitude (Linear)

1.5

−20 −30 −40 −50

• Desired cut-off frequency is fc = 1/4 • Stopband attenuation only around 10 dB

Optimal least-squares FIR filter: ( h(n), 0 ≤ n ≤ L − 1 ∆ ˆ h(n) = 0, otherwise

• Passband gain has a “resonance” near cut-off

p

Also optimal under any L norm with any error weighting: L−1 p X ˆ + cp ≥ cp ˆ = w(n) h(n) − h(n) Jp(h)

• Filter is quite poor for audio use ˆ = H ∗ asincL •H

n=0

5

6

Length L = 71 Least-Squares LPF Design

Length L = 71 Optimal Chebyshev LPF Design

Length 71 FIR Amplitude Response

remez(70,[0 0.5 0.6 1],[1 1 0 0]); 1 Length 71 FIR Amplitude Response 1.5 0.5

0

0

0.05

0.1

0.15 0.2 0.25 0.3 0.35 Normalized Frequency (cycles/sample)

0.4

0.45

0.5

Magnitude (Linear)

Magnitude (Linear)

1.5

1

0.5

0

Magnitude (dB)

−10

0

−20

0

0.05

0.1

0.15 0.2 0.25 0.3 0.35 Normalized Frequency (cycles/sample)

0.4

0.45

0

0.05

0.1

0.15 0.2 0.25 0.3 0.35 Normalized Frequency (cycles/sample)

0.4

0.45

−30 −40

0.5

0

−50

−10

−70 0

0.05

0.1

0.15 0.2 0.25 0.3 0.35 Normalized Frequency (cycles/sample)

0.4

0.45

• Stopband attenuation still only around 10 dB • Corner-frequency resonance is higher, though narrower • “Gibb’s phenomenon” in the frequency domain What we’re doing wrong:

Magnitude (dB)

−60

−20 −30 −40 −50 −60 −70

• Stopband attenuation better than 60 dB • No corner-frequency resonance

• Desired filter specification asks for too much (e.g, an infinite roll-off rate).

• Transition region from passband to stopband is critical

• Time-domain-least-squares is a poor choice of error criterion for audio work.

• Error is “equiripple” in both stopband (visible) and passband (not visible)

7

8

• Impulse response is slightly “impulsive” at the endpoints (not shown).

Lowpass Filter Specifications Lowpass Filter Specifications Amplitude

1 + δp

Ideal lowpass filter response

δp

1 1 − δp Pass-band

Stop-band TW

δs 0

0

Frequency

ωp

Transition band

ωs

• δs : stopband ripple (≤ 0.001 = −60 dB typical) • δp : passband ripple (≤ 0.1 dB typical) • ωs: stopband edge frequency • ωp: passband edge frequency • TW: transition width = ωs − ωp • SBA: stop-band attenuation = −20 log(δs) 9

Example Ripple Calculations in Matlab Let’s first consider the passband ripple spec, ±0.1 dB. Converting that to linear ripple amplitude gives, in Matlab, format long; dp=10^(0.1/20)-1 dp = 0.01157945425990 Let’s check it: >> 20*log10(1+dp) ans = 0.10000000000000 >> 20*log10(1-dp) ans = -0.10116471483635 Ok, close enough. Now let’s set the stopband ripple to 1/10 times the passband ripple and see where we are:

10

So, that’s about 60 dB stop-band rejection, which is not too bad. Now let’s set the stopband ripple to 1/100 times the passband ripple: >> ds=dp/100; >> 20*log10(ds) ans = -78.72623816882052 So now we’re close to 80dB SBA, which is getting into the ”high fidelity” zone. Ideal Lowpass Filter The ideal lowpass filter is defined by the following specifications: ∆

• ωs = ωp = ωc = 2πfc ⇒ TW = 0 • δp = δs = 0 ⇒ SBA = ∞

>> ds=dp/10; >> 20*log10(ds) ans = -58.72623816882052 11

12

The Window Method for FIR Filter Design • In practical FFT processors, we are limited to a finite duration impulse response.

Example A typical windowed ideal lowpass filter response is depicted in the following diagram: Length 101 Lowpass FIR Filter - Hamming Window

• An “obvious” method for FIR filter design is to simply window the ideal impulse response, just like we window ideal sinusoids in spectrum analysis:

1 0.8 0.6

∆

ˆ w (n) = w(n) · h h ideal(n)

Amplitude

0.4

where w is a zero phase window of length M (odd).

0.2 0

• We saw that the rectangular window is optimal in the time-domain least-squares sense, but a poor choice for typical audio filter design. What about other windows?

-0.2 -0.4 -0.6 -60

-40

-20

0 time (samples)

20

40

60

• Windowing is multiplication in the time domain ↔ convolution in the frequency domain. • Thus, in the window method, the ideal frequency response is convolved with the window transform: ˆ w (ω) = (W ∗ H H ideal)(ω) The convolution ‘smears’ the response, introducing deviations in both the pass-band and stop-band. 13

14

• At this point in the design, we have a non-causal (zero phase) filter. • length 25 FIR filter • ωc = π/4 • hideal(n) =

n∈Z

• We will examine several windows:

• This shift in the time domain results in a linear phase term in the frequency domain:

Rectangular

Amplitude

−jω M2−1 ˆ ˆ H(ω) H causal(ω) = e

sin(ωc n) πn

Hanning

0.3

0.3

0.2

0.2

Amplitude

• In order to implement this in a real time situation, we need to shift the filter by an amount equal to half its length.   M −1 ˆh ˆ causal(n) = h n − 2

Examples

0.1 0

0 -0.1 0

10 20 time (samples) Hamming

0.3

0.3

0.2

0.2

Amplitude

Amplitude

-0.1 0

0.1 0 -0.1 0

15

0.1

10 20 time (samples) Kaiser

0.1 0 -0.1 0

10 20 time (samples)

16

10 20 time (samples)

Amplitude Response of Examples

Rectangular

-20 -40 -60

0

0

-20

−10

-40 −20

-60 -2 0 2 Frequency (radians/sample) Hamming

-2 0 2 Frequency (radians/sample) Kaiser

−30

−40

0 Amplitude - dB

0

Amplitude − dB

Amplitude - dB

Amplitude - dB

Hamming Window

Hanning

0

Amplitude - dB

Effect of Changing Filter Length, M = 11, 15, 41

-20 -40 -60

−50

-20 −60

-40 -60

-2 0 2 Frequency (radians/sample)

−70

−3

-2 0 2 Frequency (radians/sample)

−2

−1

0 1 Frequency (radians/sample)

17

Summarizing:

Ideal highpass, bandpass, and bandstop filters can all be considered as modifications of an ideal lowpass filter. A highpass filter can be constructed by shifting a lowpass filter by an amount equal to half the sampling rate. |H|

i

0 |H|

• To design a high-pass filter, first design a finite length, causal lowpass filter • Use a lowpass cutoff frequency equal to π−(highpass cutoff) • Design the lowpass filter as in the previous section • Modulate the impulse reponse by (−1)n to exchange DC and fs/2

ag replacements −π

3

18

Other Types of Filters

−2π

2

π

ω

Alternatively, if the passband of the prototype lowpass ˆ filter H(ω) is very flat about unity gain, a highpass filter can be designed by simply subtracting from 1:

ω

Hhp(ω) = 1 − Hlp(ω) ↔ δ(n) − hlp(n)

2π

∆

−2π

i

−π

π

0

2π

This is equivalent to a modulation by e−jπn in the time domain. Hence, hhp(n) = hlp(n)e−jπn = hlp(n)(−1)n 19

A bandpass can be derived by designing the appropriate lowpass filter (with cutoff frequency equal to half the width of the bandpass), and then modulating it to the desired center frequency ωc. In the frequency domain, this means shifting two copies of the lowpass prototype, one to the left by −ωc, and the other to the right by ωc. The two shifted copies are added together to give a bandpass filter having a real impulse response. 20

|H|

Summary of the Window Method

−ω0c ωc |H|

−π

π

ω

Basic Idea: • Start with Ideal Lowpass-Filter Impulse Response

−π

−ω0

0

ωl ω0 ωu

πω

To design a bandpass filter: • Design a lowpass prototype filter of the desired width (cutoff ωc = ωu − ω0) • Shift to the left and right by ω0 and sum: hbp(n) = 2 cos(nω0)hlp(n) A bandstop filter can be constructed by first designing a bandpass with similar specifications. In the frequency domain, we desire a filter with a transfer function equal to Hbs = 1 − Hbp. Hence, in the time domain we have:

hbs(n) =

(

1 − hbp(0) n = 0 −hbp(n) n = ±1, ±2, . . . , ± M2−1

• Sinc function sin(πt/T )/(πt/T ) • Infinite duration • Non-causal • Window the infinite duration sinc to get a finite duration filter • Window length determines transition width • Window type determines the sidelobe level (SLL) • Rectangular window yields the optimal filter in the unweighted least squares sense • Shift the noncausal (zero phase) filter to get a causal (linear phase) filter • The linear phase slope equals half the filter length • Other types of filters (highpass, bandpass, bandstop) are designed by first designing a low pass, and applying transformations

21

22

Important Points (Window Method)

Optimal FIR Digital Filter Design

• Easy to design (see fir1 and fir2 in Octave and Matlab SPTB) • Can design extremely long FIR filters without numerical problems • Not usually optimum in any sense General Remarks (Window-Method) • Transition width TW determined by window length M • As window gets longer, its main lobe narrows • Narrower main lobe ⇒ narrower transition band • Tighter specs demand longer filters • Length M ∝ 1/TW in general • Pass-band and stop-band ripple are determined by window side lobes • We do not have individual control over passband and stopband ripple (This is a limitation of the window method) • Ripple determined by window used • Ripple normally largest around transition band

• Optimal Chebyshev FIR Design • Lp Norm Minimization • Least squares problems • Min-Max problems • Least Squares Optimization • PseudoInverse • FIR Filter Design • Linear Phase • Complex Filter Design • Other Applications Books including digital filter design: • Boyd and Vandenberghe • Parks and Burrus • Rabiner and Gold • Oppenheim and Schafer • Steiglitz

23

24

• Strum and Kirk

Optimal Chebyshev Filters

See Music 421 References1 • Minimize the maximum of the error • Stopband and passband errors are equiripple • Remez Multiple Exchange (Parks-McClellan algorithm) • Available in Octave and the Matlab SP Tool Box: remez(), cremez() Problems with Remez Exchange • Convergence unlikely for FIR filters longer than a few hundred taps or so • Fundamentally real polynomial approximation on the unit circle • Applies only to linear phase filters • Cannot design individual (non-minimum phase) “sub-filters” in polyphase filter bank form • Really need optimal weighted complex approximation e.g., minimum phase FIR filter design http://ccrma.stanford.edu/˜jos/refs421/

PSfrag replacements 25

replacements ω π 0 |H| Example: Chebyshev window and its transform: 2π −π Chebyshev Window: M = 101, Ripple = -40 dB −2π 1 ω 0.9 0.7 Amplitude

|H| π ωc ω0 ωu ωl −π −ωc −ω0

0.8

0.6 0.5 0.4 0.3 0.2 0.1

0

-100

0

-50

50

26

ω π 0 |H| 2π −π −2π ω |H| π ωc ω0 ωu ωl −π −ωc −ω0

0

DFT of Chebyshev Window: M = 11, Ripple = -40 dB

-10 -20 dB Magnitude

1

-30 -40 -50 -60 -70 -3

-2

0

-1

freq

100

Time (samples)

27

28

1

2

3

• L∞-norm

Lp norms The Lp norm of an N -dimensional vector x is defined as

∆

k x kp =

N −1 X

|xi|p

i=0

! p1

∆

k x k∞ = lim

p→∞

N −1 X i=0

|xi|p

! p1

In the limit as p → ∞, the Lp norm of x is dominated by the maximum element of x.

Special cases • L1 norm ∆

k x k1 =

N −1 X

|xi|

i=0

• Sum of the absolute values of the elements • “City block” distance • L2 norm

v uN −1 uX ∆ t |xi|2 k x k2 = i=0

• “Euclidean” distance • Minimized by “Least Squares” techniques

29

30

Filter Design using Lp Norms

Least-Squares Linear-Phase FIR Filter Design

In terms of Lp norm minimization, the FIR filter design problem becomes: min kW (ωk ) [H(ωk ) − D(ωk )]kp h

Where

Let the FIR filter length be L + 1 samples, with L even, and suppose we’ll initially design it to be centered about the time origin (“zero-phase”). Then the frequency response is given on our frequency grid ωk by L/2 X hne−jωk n H(ωk ) = n=−L/2

• ωk = suitable discrete set of frequencies • H(ωk ) = frequency response of our FIR filter • D(ωk ) = desired (complex) frequency response • W (ωk ) = (optional) weighting function Lets look at some specific cases:

Enforcing even symmetry in the impulse response, i.e., hn = h−n, gives a zero phase FIR filter which we can later right-shift L/2 samples to make a causal, linear phase filter. In this case, the frequency response reduces to a sum of cosines: L/2 X H(ωk ) = h0 + 2 hn cos(ωk n), k = 0, 1, 2, . . . , N − 1 n=1

or in matrix form:  

 h0 H(ω0) 1 2 cos(ω0 ) . . . 2 cos[ω0(L − 1)]      H(ω1)   1 2 cos(ω1 ) . . . 2 cos[ω1(L − 1)]   h1   =  ..   .. ...   .  . hL/2 H(ωN −1) 1 2 cos(ωN −1 ) . . . 2 cos[ωN −1(L − 1)] {z } | {z | A

x

(Note that Remez exchange algorithms are also based on this formulation internally.) 31

32

    

}

Matrix Formulation: Optimal L2 Design, Cont’d

Least Squares Optimization

In matrix notation, our filter design problem can be stated

xˆ = arg min k Ax − b k2 = arg min k Ax − b k22

min kAx − x

bk22

where, for zero-phase filters,   1 2 cos(ω0) . . . 2 cos [ω0(L − 1)] 2 cos(ω1) . . . 2 cos [ω1(L − 1)]  ∆ 1  A= . .  1 2 cos(ωN −1) . . . 2 cos [ωN −1(L − 1)] ∆

x=h and b = [D(ωk )] is the desired frequency response at the specified frequencies. The function firls in the Matlab Signal Processing Tool Box implements frequency-domain weighted-least-squares FIR filter design (not available in Octave as of 2/2003).

33

PSfrag replacements ω Geometrical π Interpretation of Least Squares 0 |H| Typically, the number of frequency constraints is much greater than2π the number of design variables (filter taps). −π In these cases, we have an overdetermined system of −2π equations (more equations than unknowns). Therefore, ω we cannot generally satisfy all the equations, and we are 0 left with minimizing some error criterion to find the |H| “optimal compromise” solution. π In the case ofωcleast-squares approximation, we are minimizing the ω0 Euclidean distance, which suggests the following geometrical interpretation: ωu ωl −π −ωc column-space of A −ω0 Aˆ x

∆

x

x

Hence we can minimize k Ax − b k22 = (Ax − b)T (Ax − b) Expanding this, we have: (Ax − b)T (Ax − b) = (bT − xT AT )(Ax − b) This is quadratic in x, hence it has a global minimum which we can find by taking the derivative, setting it to zero, and solving for x. Doing this yields: AT Ax − AT b = 0 PSfrag replacements These are the ωfamous normal equations whose solution is given by: π 0   |H| xˆ = (AT A)−1AT b 2π The matrix −π ∆ A† = (AT A)−1AT −2π is known as the ω pseudo-inverse of the matrix A. 0 34 |H| π ωc ω0 ωu ωl −π −ωc column-space of A −ω0 Aˆ x

This diagram suggests that the error vector b − Ax is orthogonal to the column space of the matrix A, hence it must be orthogonal to each column in A. AT (b − Aˆ x) = 0 ⇒ AT Aˆ x = AT b This is how the orthogonality principle can be used to derive the fact that the best least squares solution is given by xˆ = (AT A)−1AT b = A†b Note that the pseudo-inverse A† projects the vector b onto the column space of A.

Aˆ x=b+e Minimizexkek2 = kb − Axk2

In Matlab: •x = A \ b • x = pinv(A) * b

35

36

Chebyshev Optimal Linear Phase filter Design Consider the L∞-norm minimization problem: min k Ax − b k∞ x

We said earlier that the ∞ norm of a vector is the maximum element of that vector, hence minimizing the ∞ norm is minimizing the maximum element of a vector: min max |aTk x − bk | x

where

aTk

k

2π −π −2π Hence we are minimizing a linear objective, subject to a ω linear constraints. This is known as a linear set of 0 programming problem (linprog in Matlab). |H| π ωc ω0 ωu ωl −π −ωc −ω0 Linear Constraints

denotes the kth row of the matrix A.

We can then write this as:

Linear Objective

min t s.t. |aTk x − bk | < t   x ∆ If we introduce a new variable x˜ = , then t t = f T x˜ = [ 0 . . . 0 1 ]˜ x, and our optimization problem becomes: min f T x˜ x˜ s.t. [ aTk 0 ] · x˜ − bk < x˜f T x˜ 37

38

More General Real FIR Filters

Complex FIR Filter Design

So far we have looked at the design of linear phase FIR filters. In this case, A, x and b are all real.

In linear-phase filter design, we assumed symmetry of our filter coefficients [h(n) = h(−n)] ⇒

Sometimes we may want to design filters and specify both the magnitude and the phase of the frequency response.

• The filter frequency response became a sum of cosines (“zero phase”)

Examples:

• The matrix A was real

• Minimum phase filters

• The desired magnitude response b was real

• Inverse filters

• The final zero-phase filter xˆ could be right-shifted L/2 samples to get a corresponding causal linear-phase FIR filter

• Fractional delay filters • Interpolation polyphase subfilters

Now we would like to specify a complex frequency response. This means that: • b is complex • A is complex • We still want x (our filter coefficients) to be real If we try to use ’\’ or pinv in Matlab, we will generally get a complex result for xˆ 39

40

Optimal FIR Filters — Arbitrary Magnitude and Phase Specificiations

Summarizing our problem: min k Ax − b k2 x

where, A ∈ CN xM , b ∈ CN x1, and x ∈ RM x1 Hence we have, min k[R(A) + jI(A)] x − [R(b) + jI(b)]k22 x

which can be written as: min k R(A)x − R(b) + j [I(A)x + I(b)] k22 x

or

    2

R(A) R(b)

min x− I(A) I(b) 2 x

which is written in terms of only real variables. Hence, we can use the standard least squares solvers in Matlab and end up with a real solution.

cremez (Matlab Signal Processing Toolbox) performs complex L∞ FIR filter design: • Documented online at The Mathworks (search for cremez) • Convergence theoretically guaranteed for arbitrary magnitude and phase specifications versus frequency. • Reduces to Parks-McClellan algorithm (Remez second algorithm) as a special case. • Written by Karam and McClellan. See “Design of Optimal Digital FIR Filters with Arbitrary Magnitude and Phase Responses,” by Lina J. Karam and James H. McClellan, ISCAS-96. The paper may be downloaded at http://www.eas.asu.edu/~karam/papers/iscas96 km.html

Related paper “Design of Fractional Delay Filters Using Convex Optimization” (Mohonk-97, Music 421 handout): http://ccrma.stanford.edu/~jos/eps/mohonk97.ps.gz 41

42

Second-Order Cone Problems (SOCP)

Hilbert Transform Filter Design Example: A Case Study in FIR Filter Design

• A linear function is minimized over the intersection of an affine set and the product of second-order (quadratic) cones • Nonlinear, convex problem including linear and (convex) quadratic programs as special cases • Solved by efficient primal-dual interior-point methods

Design Example: • Ideal Single-Sideband Filter (Hilbert transformer) • Window Method using Kaiser Window • Optimal Chebyshev using Remez Exchange

• Number of iterations required to solve a problem grows at most as the square root of the problem size • Typical number of iterations ranges between 5 and 50, almost independent of the problem size See “Applications of second-order cone programming,” by M. Lobo, L. Vandenberghe, S. Boyd, and H. Lebret, Linear Algebra and its Applications, 284:193-228, November 1998, Special Issue on Linear Algebra in Control, Signals and Image Processing. Available online: http://www.stanford.edu/~boyd/socp.html

For more on this topic, see also Prof. Boyd’s course on convex optimization: EE 3642 and text book (Boyd and Vandenberghe, Cambridge U. Press, 2004). 2

http://www.stanford.edu/class/ee364/

43

44

Problem Statement

Hilbert Transform ∆

y(t) = (hi ∗ x)(t)

• Design a “negative-frequency filter” which converts a real signal into its complex “analytic signal” counterpart by filtering out negative frequencies.

∆

hi(t) =

• Equivalently, design a “single sideband” (SSB) filter which outputs a single-sideband signal in response to a complex-conjugate pair of sidebands (a Hermitian spectrum).

(Hilbert transform of x)

1 πt

(Hilbert transform “kernel”)

   −j, ω > 0 ∆ j, ω < 0 Hi(ω) =   1, ω = 0

(Hilbert frequency response)

∆

• We could also call it a “positive-frequency-pass” filter.

xa(t) = x(t) + jy(t) (Analytic signal from x) Z ∞ = π1 X(ω)ejωtdω (Note: Lower limit usually −∞)

• The imaginary part of such a filter is called a Hilbert transform filter.

0

Proof:

• The ideal Hilbert transform is an allpass filter that delays positive-frequency components by 90 degrees, and advances negative-frequency components by 90 degrees.

Xa(ω) = X(ω) + jY (ω)

(By linearity of Fourier transform)

= X+ + X− (Apply frequency response) +j[−jX+ + jX−] = 2X+(ω)

46

45

PSfrag replacements

Kaiser Window Kaiser window in causal, linear-phase form:

ω π 0 |H| 2π −π −2π ω

Kaiser Window Transform, FFT size = 2048, Window length = 257

0 −20

−40

w = kaiser(M,beta); % M=length=257, beta=8

wzp = [w((M+1)/2:M),zeros(1,N-M),w(1:(M-1)/2)]; Kaiser window transform: W = fft(wzp);

47

|H| π ωc ω0 ωu ωl −π −ωc −ω0

−60 Gain (dB)

Zero-pad by a factor of 8 and convert to zero-phase form (N = FFT size = 2048):

−80

−100

−120

−140

−160 −0.5

−0.4

−0.3

0

−0.2 −0.1 0.1 0.2 Normalized Frequency (cycles/sample)

48

0.3

0.4

0.5

eplacements ω π Kaiser window transform, close-up on main lobe 0 |H| Close Up on Kaiser LF Response, FFT size = 2048, Window length = 257 0 2π −π −2π −20 ω

Oversimplified Window Method Sample the ideal Hilbert-transform kernel h(t) = 1/πt to get ∆ ˆ i(n) = h

|H| π ωc ω0 ωu ωl −π −ωc −ω0

Gain (dB)

−40

1 πnT

(Sampled Hilbert transform kernel)

∆ ˆ w (n) = ˆ i(n) h w(n)h −60

ˆ w (ω) = (W ∗ H)(n) ˆ H (Smoothed ideal frequency response)

−80

Design Parameters:

−100

−120 −0.05

(Windowed ideal impulse response)

−0.04

−0.03

0

−0.02 −0.01 0.01 0.02 Normalized Frequency (cycles/sample)

0.03

0.04

0.05

fs = 22050; T = 1/fs; M = 257; N = 8*(M-1); beta = 8;

% % % % %

sampling rate (Hz) sampling period (sec) FIR filter length = window length for interpolated spectral displays beta for Kaiser window

Filter Design: n = [-N/2+0.5:N/2-0.5]; % Time axis (avoid t=0) hi = T ./ (pi*n*T); % Sampled Hilbert kernel hr = sin(pi*n) ./ (pi*n); % 1/2 sample delay filter h = (hr + j*hi)/2; % Sampled ideal final filter plot(f,fftshift(max(-100,20*log10(abs(fft(h)))))); grid; ... 50

49

eplacements

Gain (dB)

Gain (dB)

ω π 0 Hwp = [Hw(N/2+2:N), Hw(1:N/2+1)]; % Neg. freqs on left FFT of Truncated(2048) Ideal Impulse Response |H| Hwpn = abs(Hwp); Hwpn = Hwpn/max(Hwpn); plot(f,20*log10(Hwpn)); grid; 2π PSfrag replacements FFT of Truncated Analytically Specified "Ideal" Impulse Response ... −π 20 ω −2π FIR π Frequency Response by the Window Method 0 ω 0 |H| -20 |H| Length 257 Analytic−Derived FIR Frequency Response, FFT size = 2048 0 2π π -40 −π ωc −20 −2π ω0 -60 ω ωu −40 ωl -80 |H| −60 −π π -100 −ωc 0 -0.6 -0.4 -0.2 0.2 0.4 0.6 Normalized Frequency (cycles/sample) ω −80 c −ω0 ω0 ωu −100 h = [h(N/2+1:N),h(1:N/2)]; % zero-phase form (almost) hw = wzp .* h; % Apply window to ideal impulse response ωl −120 Hw = fft(hw); % Frequency response we really get −π −ωc −140 % Compute total stopband attenuation: −0.1 0 0.1 0.2 0.3 0.4 0.5 −ω0 −0.5 −0.4 −0.3 −0.2Normalized Frequency (cycles/sample) ierr = norm(Hw(N/2+2:N))/norm(Hw) = 0.0378 = 3.8 percent

For spectral displays, rotate buffer so negative frequencies are on the left and DC is at k = N/2: 51

52

eplacements ω π Zoom in on low-frequency transition band 0 |H| Analytic−Derived Low−Frequency Transition Band, f1/fs=0.023926 (marked by "+") 2π 20 −π 0 −2π ω −20

Since we are working with sampled data, the ideal single-side-band filter impulse response is Z π ∆ 1 u(ω)ejωndω h(n) = IDTFTn(u) = 2π −π

where u(ω) is the unit step function in the frequency domain: ( 1, ω ≥ 0 ∆ u(ω) = 0, ω < 0

−40 Gain (dB)

|H| π ωc ω0 ωu ωl −π −ωc −ω0

Ideal Bandlimited Impulse Response

−60

−80

We can evaluate the IDTFT directly as follows: π Z π 1 ejπn − 1 1 jωn h(n) = ejωndω = e = 2π 0 2πjn 2πjn 0 ( 0, n even, n 6= 0 (−1)n − 1 = = 1 2πjn j πn , n odd

−100

−120

−140 −0.05

−0.04

−0.03

0

−0.02 −0.01 0.01 0.02 Normalized Frequency (cycles/sample)

0.03

0.04

0.05

• Poorly positioned transition region from positive to negative frequencies.

For n = 0, going back to the original integral gives Z π 1 1 ejω·0dω = . h(0) = 2π 0 2

• Need to first bandpass-filter the ideal impulse response.

Thus, the real part of h(n) is δ(n)/2, and the imaginary part is 1/(πn) for odd n and zero otherwise (as a result of bandlimiting). There is no aliasing. However, we also

• High-frequency transition identical, by symmetry

54

53

PSfrag replacements need a transition band. We can work this out analytically, or we can simply “draw” it in the frequency domain, take a (large) inverse DFT, and window that. Since the latter approach is more general, we’ll do that. Window Method applied to Bandlimited Ideal Impulse Response

ω π 0 |H| 2π −π −2π ω

Desired Frequency Response

Desired Frequency Response, FFT size = 2048

0

Further Design Parameters: f1 = 530; % lower passband limit (Hz) N = 2^(nextpow2(8*fs/f1)) % FFT size from xition width if N<8*M, N = 8*(M-1); end; % spectral interpolation

Bandlimited “Ideal” Hilbert transformer frequency response (FR), including transition bands ∆ matched to FFT size N such that N  fs/∆:

Gain (dB)

−50

|H| π ωc ω0 ωu ωl −π −ωc −ω0

−100

−150 −0.5

−0.4

−0.3

0

−0.2 −0.1 0.1 0.2 Normalized Frequency (cycles/sample)

H = [ ([0:k1-2]/(k1-1)).^8,ones(1,k2-k1+1), ... ([k1-2:-1:0]/(k1-1)).^8,zeros(1,N/2-1)]; Hp = [H(N/2+2:N), H(1:N/2+1)]; f = [-0.5 + 1/N : 1/N : 0.5]; normalized freq axis plot(f,Hp); grid; ...

55

56

0.3

0.4

0.5

PSfrag replacements ω π 0 |H| 2π −π −2π ω

h = ifft(H); % zero-phase form % measure of round-off error: ierr = norm(imag(h(1:2:N)))/norm(h) % zero? ierr =

|H| π ωc ω0 ωu ωl −π −ωc −ω0

4.1958e-15 % rough estimate of time-aliasing error: aerr = norm(h(N/2-N/32:N/2+N/32))/norm(h) aerr = 4.8300e-04

Real Part of Desired Impulse Response

Ideal SSB Filter Impulse Response, Real Part, FFT size = 2048 0.6

0.5

0.4

Amplitude

Desired Impulse Response

0.3

0.2

0.1

0 −0.1 −150

−100

0

−50

50

100

150

Time (samples)

% for plots, prefer negative times on the left: hp = [h(N/2+2:N), h(1:N/2+1)];

58

57

eplacements ω π Imaginary Part of Desired Impulse Response 0 |H| Ideal SSB Filter Impulse Response, Imaginary Part, FFT size = 2048 2π 0.4 −π 0.3 −2π ω 0.2

ω π 0 |H| 2π −π −2π ω

Final FIR Filter Frequency Response

Length 257 FIR filter Frequency Response, FFT size = 2048

0 −20

−40

0.1

−60

0 −0.1

−0.2

−0.3

−0.4 −150

−100

−50

0 Time (samples)

50

100

150

|H| π ωc ω0 ωu ωl −π −ωc −ω0

Gain (dB)

Amplitude

|H| π ωc ω0 ωu ωl −π −ωc −ω0

PSfrag replacements

−80

−100

−120

−140

−160

−180 −0.5

−0.4

−0.3

0

−0.2 −0.1 0.1 0.2 Normalized Frequency (cycles/sample)

Apply window to ideal impulse response: hw = wzp .* h; % Final FIR coefficients Hw = fft(hw); % Frequency response we’ll see

Total stopband energy: ierr = norm(Hw(N/2+2:N))/norm(Hw) disp(sprintf([’Stop-band energy is %g percent of’,... total spectral energy’], 100*ierr)); 59

60

0.3

0.4

0.5

eplacements

Remez Multiple Exchange Method

Zoom-in on Transition Region of Final Frequency ω Response π 0 |H| Low−Frequency Transition Band, f1/fs=0.023926 (marked by "+") 2π 20 −π 0 −2π ω −20

• Original reference: J. H. McClellan, T. W. Parks, and L. R. Rabiner, “A Computer Program for Designing Optimum FIR Linear Phase Digital Filters,” IEEE Trans. Audio Electro., vol AU-21, Dec. 1973, pp. 506-526.

−40 Gain (dB)

|H| π ωc ω0 ωu ωl −π −ωc −ω0

• Remez exchange algorithm for optimal Chebyshev (equiripple) FIR filter design.

• Text reference: L. R. Rabiner and B. Gold, Theory and Application of Digital Signal Processing, Prentice-Hall, 1975.

−60

−80

−100

• Pertinent recent reference: A. Reilly, G. Frazer, and B. Boashash, “Analytic Signal Generation—Tips and Traps,” IEEE Tr. Signal Processing, vol. 42, no. 11, Nov. 1994.

−120

−140 −0.05

−0.04

−0.03

0

−0.02 −0.01 0.01 0.02 Normalized Frequency (cycles/sample)

0.03

0.04

0.05

• Transition band better positioned

Design an optimal lowpass filter of the desired width (takes quite a while):

• ≈ 100 dB stob-band attenuation

hrm = remez(M-1, [0,(f2-fs/4)/fn,0.5,1], ... [1,1,0,0], [1,10]);

The weighting [1,P] says make the passband ripple P times that of stopband. For steady-state audio spectra, 61

62

PSfrag replacements

Modulate it up to where it’s a single-sideband filter. I.e., right-shift by π/2 in the frequency domain. I.e., rotate the frequency response counterclockwise along the unit circle by 90 degrees: hr = hrm .* j .^ [0:M-1];

63

Optimal Chebyshev FIR(257) Frequency Response

ω π 0 |H| 2π −π −2π ω |H| π ωc ω0 ωu ωl −π −ωc −ω0

Length 257 Optimal Chebyshev FIR Single−Sideband−Filter Freq. Response 50

0 −50 Gain (dB)

passband ripple can be 0.1 dB or more. However, consider an FM signal in which a sinusoid is sweeping back and forth in the passband. In that case, the passband ripple generates AM sidebands, so a spec more like that in the stopband may be called for. Here, we allow the passband ripple to be 10 times the stopband ripple, as a compromise.

−100

−150

−200 −0.5

−0.4

−0.3

0

−0.2 −0.1 0.1 0.2 Normalized Frequency (cycles/sample)

64

0.3

0.4

0.5

eplacements

PSfrag replacements

Low−Frequency Transition Band, f1/fs=0.023926 (marked by "+") 50

0 −50 Gain (dB)

|H| π ωc ω0 ωu ωl −π −ωc −ω0

ω π 0 |H| 2π −π −2π ω

Zoom-In on Transition Band

−100

−150

−200 −0.05

−0.04

−0.03

0

−0.02 −0.01 0.01 0.02 Normalized Frequency (cycles/sample)

0.03

0.04

0.05

|H| π ωc ω0 ωu ωl −π −ωc −ω0

Passband Ripple

−4

4

Passband Ripple

x 10

3

2

1 Gain (dB)

ω π 0 |H| 2π −π −2π ω

0 −1

−2

−3

−4 0.18

0.2

0.22

0.24 0.26 0.28 Normalized Frequency (cycles/sample)

65

66

Initial Impulse Response

Summary of Trade-Offs

0.3

0.32

eplacements ω π 0 |H| 2π −π −2π ω

Window-method filter:

First 50 impulse−response coefficients

x 10

• Stopband ”droops” which ”wastes” filter taps if the worst-case stopband attenuation is the only stopband spec. However, such ”roll-off” is natural in the frequency domain and corresponds to a ”smoother pulse” in the time domain.

6

4

2 Amplitude

|H| π ωc ω0 ωu ωl −π −ωc −ω0

−4

8

0

• Passband is degraded by early roll-off. Edge is off. • Filter length can be thousands of taps as needed for full audio band.

−2

−4

Optimal Remez-Exchange filter:

−6

−8

0

5

10

15

20

25 Index

30

35

40

45

50

• Stopband is ideal, equiripple. • Transition region close to half that of the window method. • Pass-band is ideal, equiripple. • Time-domain impulses from passband ripple are small. • Design time orders of magnitude larger than that for window method.

67

68

• Fails to converge for filters much longer than 256 taps. (Need to increase working precision to get longer filters.)

Matlab hilbert() function The Matlab Hilbert function returns an analytic signal given its real part. PSfrag replacements Nh ω = M-2; % This looks best delta = [1,zeros(1,Nh)]; % zero-phase impulse π hh = hilbert(delta); % zeros negative-freq FFT bins 0 = fft(hh); % FIR frequency response Hh

|H| 2π −π −2π ω

−50

−100

Gain (dB)

|H| π ωc ω0 ωu ωl −π −ωc −ω0

Frequency Response of Length 256 hilbert() using Matlab

0

−150

−200

−250

−300

−350 −0.5

−0.4

−0.3

0

−0.2 −0.1 0.1 0.2 Normalized Frequency (cycles/sample)

0.3

0.4

0.5

Looks pretty good, but let’s zero-pad the impulse response to interpolate the frequency response: 70

69

PSfrag replacements

ω π = [hh(Lh/2+1:Lh), zeros(1,N-Lh), hh(1:Lh/2)]; hhzp Hhzp 0 = fft(hhzp); % Frequency response

|H| 2π −π −2π ω

−50

|H| π ωc ω0 ωu ωl −π −ωc −ω0

−100

Gain (dB)

|H| π ωc ω0 ωu ωl −π −ωc −ω0

Interpolated Frequency Response of Length 256 hilbert() using Matlab

0

−150

−200

−250

−300

−350 −0.5

−0.4

−0.3

0

−0.2 −0.1 0.1 0.2 Normalized Frequency (cycles/sample)

0.3

0.4

0.5

Low−Frequency Transition Band, f1/fs=0.023926 (marked by "+") 50

0 −50

−100 Gain (dB)

eplacements

ω π 0 |H| 2π −π −2π ω

−150

−200

−250

−300

−350 −0.05

−0.04

−0.03

0

−0.02 −0.01 0.01 0.02 Normalized Frequency (cycles/sample)

0.03

0.04

0.05

Problems: • Transition bandwidth only 2 bins wide. • Massive time aliasing. • Only justifiable for periodic signals with period exactly equal to filter length Nh. In this case only, time aliasing is equivalent to overlap-add from adjacent periods.

71

72

More generally, any real signal given to hilbert() must be interpreted as precisely one period of a periodic signal.

73