Digital filters
Lecturer: Dr. Bui Ha Duc     Email: ducbh@hcmute.edu.vn
Department of Mechatronics
1
Signal Filter
       Signal Filter is a frequency selective LTI system
                                                            𝑀−1
                                                       1
                                                  𝑦 𝑛 =  𝑥[𝑛 − 𝑘]
                                                       𝑀
                                                            𝑘=0
Analog filters                             Digital filters:
• Electronic components are cheap.         • Better performance
• Large dynamic range in amplitude and     • Programmable.
  frequency.                               • One hardware design can implement
• Real-time.                                 many different filters.
• Low stability of resistors, capacitors   • The characteristics of digital filters
  and inductors due to temperature.          are predictable.
• Difficult to get the components          • Can be evaluated in simulation before
  accuracy as calculated by the formula      implementing in hardware.
    2
Digital filters
       Digital filters are a very important part of DSP:
           signal separation: noise, interference, other signals
            E.g. measuring the electrical activity of a baby's heart (EKG) while still in
            the womb; extracting EOG
           signal restoration
            E.g. recovering an audio recording made with poor equipment,
            calibration
    3
    Remove noise on raw RSSI data
4
    Remove Eye blinks from EEG signal
5
    Remove noise on raw GPS data
6
Ideal filters
    Output of a digital filter:
                                  𝒀 𝒆𝒋𝝎 = 𝑯 𝒆𝒋𝝎 𝑿(𝒆𝒋𝝎 )
                             Frequency Filter        Input
                                                     frequencies
7
Ideal Filter
       The output of ideal filter is
           distortionless response over one or more frequency bands
           zero response at the remaining frequencies
                                                             −𝒋𝝎𝒏𝒅 , 𝝎 ≤ 𝝎
                                             𝑯   𝒆𝒋𝝎   =൝  𝒆               𝒄
                                                          𝒐,    𝝎𝒄 ≤ 𝝎 ≤ 𝝅
                             Ideal lowpass filter
    8
Ideal Filter
       Invert Fourier transform of an ideal lowpass filter
                          𝒋𝝎
                                           𝒔𝒊𝒏(𝝎𝒄 (𝒏 − 𝒏𝒅 ))
                    𝑯 𝒆        → 𝒉𝒍𝒑   𝒏 =
                                              𝝅(𝒏 − 𝒏𝒅 )
               → Not practical
    9
Practical filter
    Practical filter is usually done by minimizing some
     approximation error between the non-ideal filter and a
     prototype ideal filter
    A good filter should have only a small ripple in the passband,
     high attenuation in the stopband, and very narrow transition
     bands
    10
Distortion of signals passing through LTI
systems
    Distortionless response systems
      𝑦 𝑛 = 𝐺𝑥 𝑛 − 𝑛𝑑 ,       𝐺 and 𝑛𝑑 𝑎𝑟𝑒 𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡𝑠
    Taking Fourier transform on both sides:
                    𝒀 𝒆𝒋𝝎 = 𝑮𝒆−𝒋𝝎𝒏𝒅 𝑿(𝒆𝒋𝝎 )
    Frequency response:
                                                           MATLAB
         𝑯   𝒆𝒋𝝎   =   𝑮𝒆−𝒋𝝎𝒏𝒅   = 𝑯   𝒆𝒋𝝎   𝒆𝒋∠𝑯(𝒆^𝒋𝝎 )   freqz(b,a)
                           𝑯 𝒆𝒋𝝎   =𝑮                      abs()
                        ∠𝑯 𝒆𝒋𝝎 = −𝝎𝒏𝒅                      angle()
    11
Distortion of signals passing through LTI
systems
    A filter introduces magnitude distortion 𝑯 𝒆𝒋𝝎   ≠𝑮
    A filter introduces phase or delay distortion if ∠𝑯 𝒆𝒋𝝎 ≠ −𝝎𝒏𝒅
    12
 Distortion of signals passing through LTI
 systems
Example:
                        1               1
Input x[𝑛] = cos 𝜔𝑛 −     cos   3𝜔𝑛 +     cos(5𝜔𝑛)
                        3               5
Suppose that the output y[n] given by:
   y[𝑛] = 𝑐1 cos 𝜔𝑛 + 𝜙1 − 𝑐2 cos 3𝜔𝑛 + 𝜙2 + 𝑐3 cos(5𝜔𝑛 + 𝜙3 )
  13
Distortion of signals passing through LTI
systems
14
Distortion of signals passing through LTI
systems
15
Example
    Determine the magnitude and phase response of this
     filter
    16
Example
17
Magnitude and phase response functions and input–output signals for the filter
18
Group delay
    Group delay is a convenient way to check the linearity of
     phase response
    Group delay is defined as the negative of the slope of the
     phase:
                                𝑑(∠𝑯 𝒆𝒋𝝎
                       𝜏𝑔𝑑 = −
                                    𝑑𝜔
    Linear-phase response have constant group delay
    19
Group delay
    Analyze this frequency response
    20
Filter design procedure
1.    Specification: Specify the desired frequency response
      function characteristics to address the needs of a specific
      application.
2.    Approximation: Approximate the desired frequency
      response function
3.    Quantization: Quantize the filter coefficients at the
      required fixed-point arithmetic representation.
4.    Verification: Check whether the filter satisfies the
      performance requirements by simulation or testing with real
      data.
5.    Implementation: Implement the system obtained in
      hardware, software.
 21
Filter specifications
     Practical filter:
     • the passband responses are not perfectly flat
     • the stopband responses cannot completely reject (eliminate) bands
       of frequencies
     • the transition between passband and stopband regions takes place
       over a finite transition band.
22
Filter specifications
23
Decibel
    A bel (in honor of Alexander Graham Bell) means that
     the power is changed by a factor of ten
        E.g 3 bel of amplification produces an output signal with
         10×10×10 = 1000 times the power of the input
    dB =1/10 bel, every ten decibels mean that the power has
     changed by a factor of ten
                                        𝑃2
                         𝑑𝐵 = 10𝑙𝑜𝑔10
                                        𝑃1
    Every twenty decibels mean that the amplitude has
     changed by a factor of ten.
                                       𝐴2
                        𝑑𝐵 = 20𝑙𝑜𝑔10
                                       𝐴1
    24
Filter specifications
    Absolute specifications
        In the passband
              1 − 𝛿𝑝 ≤ 𝐻 𝑒 𝑗𝜔     ≤ 1 + 𝛿𝑝 ,        0 ≤ 𝜔 ≤ 𝜔𝑝
        In the stopband
                    𝐻 𝑒 𝑗𝜔      ≤ 𝛿𝑠 ,      𝜔𝑠 ≤ 𝜔 ≤ 𝜋
    Relative specifications
        passband ripple
                                     1+𝛿𝑝
                    𝐴𝑝 = 20𝑙𝑜𝑔10               𝑑𝐵
                                     1−𝛿𝑝
        stopband attenuation
                                1 + 𝛿𝑝
                𝐴𝑠 = 20𝑙𝑜𝑔10           ≈ −20𝑙𝑜𝑔10 (𝛿𝑠 )   dB
                                  𝛿𝑠
    25
Filter specifications
Then
             𝐴𝑝 ≤ 𝐻 𝑒 𝑗𝜔 , 𝑖𝑛 𝑑𝐵 ≤ 0,         0 ≤ 𝜔 ≤ 𝜔𝑝
               𝐻 𝑒 𝑗𝜔 , 𝑖𝑛 𝑑𝐵 ≤ −𝐴𝑠 ,        𝜔𝑠 ≤ 𝜔 ≤ 𝜋
    In a well designed filter, typically 𝐴𝑝 ≈ 0 and 𝐴𝑠 ≫ 1
    The smaller 𝛿𝑝 and 𝛿𝑠 , and the narrower the transition band,
     the more complicated (higher order) is the required filter.
    26
Example:
 A filter specified have a magnitude response in the
  passband within ±0.001, and the minimum stopband
  attenuation is given as 𝛿𝑠 =0.001.
 𝐴𝑝 , 𝐴𝑠 =?
                        1 + 0.001
         𝐴𝑝 = 20𝑙𝑜𝑔10             = 0.0174 𝑑𝐵
                        1 − 0.001
             𝐴𝑠 = −20𝑙𝑜𝑔10 (0.001)=60 dB
    27
    Example:
    A lowpass digital filter is specified by the following
     relative specifications:
          𝜔𝑝 = 0.3π, 𝐴𝑝 = 0.5 dB; 𝜔𝑠 = 0.5π, 𝐴𝑠 = 40 dB.
𝛿𝑝 , 𝛿𝑠 = ?
    28
FINITE IMPULSE RESPONSE FILTERS
    FIR system:
                          𝑀
                   𝑦 𝑛 =  𝑏𝑘 𝑥[𝑛 − 𝑘]
                          𝑘=0
                         𝑏𝑛 , 𝑛 = 0, … , 𝑀
            and    ℎ 𝑛 =ቊ
                          0,   𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒
    31
Characteristics of FIR filters
    FIR filters are always stable
    The filter has finite memory
    The design of linear-phase filters can be guaranteed
    FIR filters can be easily implemented on most DSP
     processors, but may result in high computational cost
    32
FIR filters with linear phase
    Frequency response
                                  ∞
                     𝑯 𝒆𝒋𝝎 =  𝒉[𝒌]𝒆−𝒋𝝎𝒌
                                𝒌=−∞
    An FIR filter has linear phase if its coefficients /impulse
     response satisfy the following symmetric condition
            ℎ 𝑛 = ±ℎ 𝑀 − 𝑎 ,             𝑀 𝑖𝑠 𝑎 𝑖𝑛𝑡𝑒𝑔𝑒𝑟
    33
Moving Average Filters
    Moving Average Filters:
                                  𝑀−1
                            1
                       𝑦 𝑛 =  𝑥[𝑛 − 𝑘]
                            𝑀
                                  𝑘=0
    Advantages:
        Easy to understand and implement
        Less computation time
    34
Moving Average Filters
    Noise Reduction vs. Step Response
Increasing the number of points lead to?
    35
Moving Average Filters
    Frequency Response of moving average filter
                          𝜔𝑀
                      sin(    )
             𝐻 𝑒 𝑗𝜔 =      2 𝑒 −𝑗𝜔(𝑀−1)/2
                            𝜔
                      𝑀𝑠𝑖𝑛( )
                            2
                                  𝜔𝑀
                              sin(    )
                   𝐻 𝑒 𝑗𝜔   =      2
                                    𝜔
                              𝑀𝑠𝑖𝑛( )
                                    2
    36
Moving Average Filters
                                                    Transitional band?
                                                    Stopband attenuation?
                                                   The moving average filter is
                                                   a good smoothing filter but
                                                   a bad low-pass-filter
     Frequency Response of moving average filter
37
 Multiple pass
Multiple-pass
moving average
filters involve passing
the input signal
through a moving
average filter two or
more times
   38
Example of Moving Average filter
The file djw6576.txt contains the weekly opening value
x[n], 0 ≤ n ≤ N−1, of the Dow Jones Industrial Average
for N = 600 weeks starting in January 1965
    Investigate the effect of changing M
    39
Summary of Moving Average filter
    Optimal filter for the following tasks:
        Reducing random noise while retaining the sharp step
         response
        Useful for time domain encoded signals
    Worst filter concerning frequency encoded signals (no
     frequency separation capabilities !)
    40
Windowed-sinc filter
    Windowed-sinc filters are used to separate different
     bands of frequencies.
    Advantages:
      stable
      very good overall performance in the frequency domain.
     (poor performance in time domain.)
      if done by standard convolution: needs much computation
       power.
      less execution time necessary when realized by FFT-programm.
    41
Windowed-sinc filter
    Motivation
         Ideal lowpass filter ⟶ infinite
         impulse response
    42
Windows-sinc filter
If we make impulse response
become finite                       result from the abrupt discontinuity
Need to remove the abrupt end! But how?
43
    Solution: multiplying the modified-sinc with a window to
     reduce the abruptness of the truncated ends and thereby
     improve the frequency response
    44
     X
45
46
Design Windowed-sinc filters
    2 parameters must be determined:
        the cutoff frequency, 𝑓𝑐
        the order of the filter kernel, M.
    After fc and M have been selected, the filter kernel is
     calculated from the relation:
                         𝑠𝑖𝑛(2𝜋𝑓𝑐 (𝑛 − 𝑛𝑑 ))
               ℎ 𝑛 =𝐴                        x w[n]
                             𝜋(𝑛 − 𝑛𝑑 )
    47
     • Cut-off frequency according to the sampling rate
     • Filter length M defines the transition bandwith
                                         4
                             𝑀≈
                                   𝑏𝑎𝑛𝑑𝑤𝑖𝑑𝑡ℎ
48
49
The order of the filter kernel
           𝑏𝑛 , 𝑛 = 0, … , 𝑀
    ℎ 𝑛 =ቊ
            0,   𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒
    M is the order of the filter
    A filter with order M has L = M+1 coefficients
    M should be an even integer
    50
The order of the filter kernel
    Type-I FIR filter:
        has a symmetric impulse response
        even order M
    Consider the case M = 4
    51
The order of the filter kernel
    In general, the frequency response of a type-I FIR filter
     can be expressed as:
    52
The order of the filter kernel
    Type-II FIR filter:
        has a symmetric impulse response
        odd order M
    Consider the case M = 5
    53
The order of the filter kernel
    In general, the frequency response of a type-II FIR filter
     can be expressed as:
    At 𝜔 = 𝜋, 𝐴 𝑒 𝑗𝜔 = 0, independent of h[k]
                  filters with a frequency response that is nonzero
                  at ω = π, for example highpass filter, band stop
                  filter, cannot be implemented
    54
Design of FIR filters using windows
     A window useful for filter design should have a Fourier
     transform that has a narrow mainlobe, a small relative peak
     sidelobe, and good side lobe roll-off
55
FIR filter design using windows
Procedure:
1. Given the design specifications {𝜔𝑝 , 𝜔𝑠 , 𝐴𝑝 , 𝐴𝑠 }, determine the ripples 𝛿𝑝 and 𝛿𝑠 and
set 𝛿 = min{𝛿𝑝 , 𝛿𝑠 }.
2. Determine the cutoff frequency of the ideal lowpass prototype by 𝜔𝑐 = (𝜔𝑝 +
𝜔𝑠 )/2.
3. Determine the design parameters 𝐴𝑝 , 𝐴𝑠 and ∆𝜔 = 𝜔𝑠 − 𝜔𝑝 .
4. Choose the window function that provides the smallest stopband attenuation
greater than A.
5. Determine the required value of M = L−1 by selecting the corresponding value of ω
from the column labeled “exact ω”. If M is odd, we may increase it by one to have a
flexible type-I filter.
6. Determine the impulse response of the ideal lowpass filter by
                                         𝑠𝑖𝑛(𝜔𝑐 (𝑛 − 𝑀/2))
                               ℎ𝑑 𝑛 =
                                              𝜋(𝑛 − 𝑀/2)
7. Compute the impulse response ℎ[𝑛] = ℎ𝑑 [𝑛]𝑤[𝑛] using the chosen window.
8. Check whether the designed filter satisfies the prescribed specifications; if not,
increase the order M and go back to step 6.
 56
Design of FIR filters using windows
     Properties of commonly used windows (L = M + 1)
57
Design of FIR filters by frequency sampling
Example:
Design a lowpass linear-phase FIR filter to satisfy the
following specifications:
𝜔𝑝 = 0.25𝜋, 𝜔𝑠 = 0.35𝜋, 𝐴𝑝 = 0.1 𝑑𝐵, 𝐴𝑠 = 50 𝑑𝐵.
 58
Kaiser window
    Kaiser window is an windows with adjustable parameters
     designed to optimize the trade-off between main-lobe width
     and sidelobe height
                             𝑛−𝛼 2
                 𝐼0 𝛽 1 − 𝛼
            𝑤𝑛 =                   ,   0≤𝑛≤𝑀
                       𝐼0 (𝛽)
                           0,  𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒
Where
                                       ∞   (𝑥/2)𝑚
         𝛼 = 𝑀/2   and   𝐼0 𝑥 = 1   + σ𝑚=1
                                             𝑚!
python function: signal.kaiserord(L,beta)
    59
Kaiser window
60
Kaiser window
61
 Example
Design a lowpass filter to remove the high frequency noise
above 1,600 Hz with the following filter specifications:
passband frequency range: 0- 1,600Hz; passband ripple: 0.02
dB; stopband frequency range: 1,800 - 4,000 Hz; stopband
attenuation: 50 dB.
 Use the designed lowpass filter to filter the noisy speech
  (nspeech.mat)
    62
Summary of Windowed-sinc filters
    incredible performance
    But trade-off in execution time
    63
PARKS–MCCLELLAN METHOD
    One of the most popular optimal design method used in
     the industry due to its efficiency and flexibility
    The algorithms minimize the maximum error between
     the desired frequency response and actual frequency
     response
                  𝐸 = 𝑚𝑎𝑥 𝐻𝑑 𝑒 𝑗𝜔 − 𝐻 𝑒 𝑗𝜔
    64
     Parks–McClellan algorithm
65
    Calculate the weights
    Use Remez algorithm to calculate filter coefficients
         f_cutoff = 500
         order = cycles*fs // f_cutoff
         # Build the filter using the Remez method
         # Pass band: [0, 0.75 * f_cutoff] (begin to transition at 75% of our cutoff), gain: 1
         # Stop band: [f_cutoff, fs/2] (attenuate everything above the cutoff), gain: 0
         hpm = scipy.signal.remez(order,
                          [0, 0.75 * f_cutoff, f_cutoff, fs/2], # Our frequency bands
                          [1, 0], # Our target gain for each band
                          fs=fs)
         # Apply the filter to a signal
         y = scipy.signal.convolve(x, hpm)
    66
INFINITE IMPULSE RESPONSE FILTERS
    IIR system:
                        𝑀                    𝑀
              𝑦 𝑛 =  𝑏𝑘 𝑥[𝑛 − 𝑘] −  𝑎𝑘 𝑦[𝑛 − 𝑘]
                       𝑘=0                  𝑘=1
    ℎ 𝑛 is a infinite sequence
    The objective in IIR filter design is to determine the coefficients
     so that its frequency response H(ejω) approximates an ideal
     desired response Hd(ejω)
    67
IIR filters
    Step 1 Convert the discrete-time design
     specifications into continuous-time specifications.
    Step 2 Design a continuous-time filter, that is, obtain a
     system function Hc(s) that satisfies the continuous-time
     specifications.
    Step 3 Convert Hc(s) to an appropriate system function
     H(z), which meets the specifications, using a continuous-
     time to discrete-time transformation.
    68
Characteristics of IIR filters
    IIR filters cannot have linear phase
    IIR filter not always stable
    IIR filters involves both the magnitude and phase
     responses
    IIR filters require a much lower filter order than FIR
     filters
         IIR filters are suitable for lower cost system
    69
IIR filters
  Consider the following IIR system:
                 1
          𝑦[𝑛] = (𝑥[𝑛] − 𝑥[𝑛 − 2]) + 𝑦[𝑛 − 1]
                 2
1. What is the order of this system?
2. What are b and a for this system?
3. What is its impulse response?
    70
Transformation of continuous-time filters
to discrete-time IIR filters
    Impulse-invariance transformation
                       ℎ 𝑛 = 𝑇𝑑 ℎ𝑐 𝑛𝑇𝑑
𝑇𝑑 is the sampling period
    This transformation is known as impulse-invariance
     because it preserves the shape of the impulse
     response
    Frequency response:
                    𝑗𝜔
                                𝜔
                𝐻 𝑒    = 𝐻𝑐 𝑗       ,     𝜔 <𝜋
                                𝑇𝑑
If
                  𝐻𝑐 𝑗Ω = 0,        Ω ≥ 𝜋/𝑇𝑑
    71
Transformation of continuous-time filters
to discrete-time IIR filters
    Bilinear transformation
        invertible nonlinear mapping between the s-plane and the z-
         plane
                                      2 1 − 𝑧 −1
                          𝑠=𝑓 𝑧 =
                                      𝑇𝑑 1 + 𝑧 −1
    Td has no effect on the design process since the bilinear
     transformation does not involve any sampling
     operation.
    72
Transformation of continuous-time filters
to discrete-time IIR filters
73
Common IIR filters
74
Butterworth filter
    Butterworth lowpas filter (LPF) was proposed by Stephen
     Butterworth in 1930
    any filter with its frequency response H(Ω) that satisfies
     the equation
     is a low pass filter with order N and cutoff frequency Ωc.
     Laplace transform H(s) :
    75
Butterworth filter
    Design IIR low pass filter with specifications as
    76
Butterworth filter
                     N=2
77
Butterworth filter
     Design plots for the eighth-order digital lowpass Butterworth filter
78
Butterworth filter
    We note that Butterworth filters have the maximum
     flatness, close to the ideal response.
    Butterworth filters have a very smooth group-delay
     response
    Python functions for Butterworth filters:
        b, a = scipy.signal.butter(order, fc, fs)
    79
Butterworth filter
    Order of the filter is determined based on pass-band
     ripple, stop-band attenuation and transition band
           fc = 500 # Cutoff at 500 Hz
           fstop = 1000 # End the transition band at 1000 Hz
           Ap = 3 # 3dB ripple
           As = 60 # 60 dB attenuation of the stop-band
           # Get the smallest order for the filter
           order = scipy.signal.buttord(fc, fstop, Ap, As , fs=fs)
           # Now construct the filter
           b, a = scipy.signal.butter(order, fc, fs=fs)
    80
Butterworth filter
 Use scipy.signal.butter to construct a high-pass filter
  (fs=44100) with the following properties:
1. fc = 50
2. a transition width of 250 Hz,
3. at most 1 dB of pass-band loss, and
4. stop-band attenuation of 60 dB.
What is the order of the resulting filter?
    81
Chebyshev filter
    named after the mathematician Pafnuty Chebyshev
    There are two types of Chebyshev filters:
        Type 1: have a steeper transition, higher pass-band ripple than
         Butterworth filters
             # Get the order and discard the second output (natural frequency)
             order, _ = scipy.signal.cheb1ord(fc, fstop, ripple, attenuation, fs=fs)
             # Build the filter
             b, a = scipy.signal.cheby1(order, ripple, fc, fs=fs)
        Type 2: have a steeper transition, lower stop-band attenuation
         than Butterworth filters
             order2, _ = scipy.signal.cheb2ord(fc, fstop, ripple, attenuation, fs=fs)
             b2, a2 = scipy.signal.cheby2(order2, attenuation, fstop, fs=fs)
    82
Chebyshev filter
     Design plots for the fifth-order digital lowpass Chebyshev filter
83
Elliptic filters
    Elliptic filters have extremely steep transitions but allow
     ripple in both the pass- and stop-bands
    Python code to create Elliptic filters
         # Get the order and discard the second output (natural frequency)
         order, _ = scipy.signal.ellipord(fc, fstop, ripple, attenuation, fs=fs)
         # Build the filter
         b, a = scipy.signal.ellip(order, ripple, attenuation, fc, fs=fs)
    84
Exercise
Implement a low-pass filter at fs = 44100 with the following
specification:
         Cut-off frequency Hz
         No more than 1.5 dB of ripple in the pass-band (less is acceptable)
         At least 80 dB of attenuation in the stop-band (more is acceptable)
         A transition band of no more than 500 Hz (less is acceptable)
a.        Which of the IIR filter(s) should be used to implement a
          filter with these constraints?
b.        Determine the order of the filter. If there are multiple filter
          types that will work, which one gives the smallest order?
c.        Implement an FIR filter with the given constraints. What
          order does the filter need to be to satisfy the constraints?
     85