IIR_null_phase
Weimich

“Non-Causal“ Digital Chebyshev Filter with Linear Phase Response. C Implementation and Octave Script

1. Abbreviation

DSP – Digital Signal Processing
IIR – Infinite Impulse Response
FIR – Finite Impulse Response
TR – Time Reverse
MAC – Multiplications and Additions per Input Sample
SOS – Second Order Section

 

2. Introduction

Often in Digital Signal Processing applications it is necessary to use filters with a linear phase characteristic. In the case is possible to avoid phase distortion of the signal. It is impossible to construct a stable IIR filter with an ideal linear phase characteristic, since the following condition must be met:

Chebyshev_formula1

The (1) means that all the filter poles must be inverted:

Chebyshev_formula2_3

In this case some of the poles are outside the unit circle or all of them must lie on the unit circle of the z plane. It means that the filter is an unstable.
Note: the proof of the formula (1) you can see in [1].
In the article “Bessel Filter in Digital Form?” I am considering the possibility of constructing an IIR filter with a phase characteristic close to linear but only an approximation is obtained.
[1] is showed a non-causal an IIR filter with a “zero” delay which uses reverse time. Sozanski in [2] describes the possibility of implementing such an approach in digital technology. In the article I want to design a Digital IIR Chebyshev Filter with a linear phase using [2].
Figure 2-1 shows a non-causal filter with a zero phase (linear phase).
IIR_null_phase

Figure 2-1: Non-Causal IIR Filter with Null Phase/Delay Response

 

We describe the signal spectra at the output of each block of the non-causal filter using the signal property of the y(n) and y (- n) which have conjugate spectrum
Chebyshev_formula3_1
  • Spectrum of the g(n):
    Chebyshev_formula4
  • Spectrum of the g(-n):
    Chebyshev_formula5
  • Spectrum of the y(-n):
    Chebyshev_formula6
  • Spectrum of the output signal y(n):
    Chebyshev_formula7
(7) matches the filter shown in Figure 2-2.

IIR_null_phase_equal

Figure 2-2: Equivalent filter for Figure 2-1

 

Notes:
  • The filter (7) has a zero phase (linear phase). It is not possible to implement a zero-latency filter. The implementation [2] has a sufficiently large time delay (see below), which occurs when implementing (approximating) the reverse time
  • The filter (7) differs from the base filter H(z) not only in the phase response but also in the magnitude response.
    Base filter H (z):
    Chebyshev_formula8
    Non-causal filter:
    Chebyshev_formula9
    Ripples of the non-causal filter is verdoppelt, and also reduces the bandwidth compared to H(z). Usually the filter bandwidth is defined at -3dB. The non-causal filter has a passband of the original filter -6dB (verdoppelt).

 

3. Time Reverse (TR) Implementation from [2]

See the Figure 2-1. First the input signal x(n) is passed through the first filter H(z) after which the g(n) must be reversed in time and then g(-n) is passed through the second filter H(z) again. We divide g(n) into blocks of length N and process each block separately. See the Figure 3-1.

blocks_definition

Figure 3-1: Block Definitions

 

Next each block is done a time reversion and sent these samples to the second H(z) filter input. But unfortunately it doesn’t work so simple. The second H(z) filter must prepare correct internal state (filter history). Therefore extra Noverlap samples from the subsequent block are added to each block which prepare only the filter state. The second filter output is discarded for these Noverlap samples.
Note: Before starting each subsequent time reverse block we need to clear the history of the second H(z) filter.
Blocks in reverse time you can see in the Figure 3-2.

g_reverse_blocks

Figure 3-2: Reverse Time for g(-n) Blocks

 

At the output second H(z) the y (- n) is got . See the Figure 3-3.

y_reverse_blocks

Figure 3-3: Reverse Time for y(-n) Blocks

 

The subsequent time reversion for y(- n) blocks gets the output signal y(n). See the Figure 3-4.

y_blocks

Figure 3-4: Output y(n) Blocks

 

Conclusions and notes:
  • Processing of the input signal x(n) is possible if the (N + Noverlap) samples is accumulated and then the algorithm shown in Figure 2-1 can be started. Thus the filter has the delay:

    Chebyshev_formula10

    The formula (10) is valid only for the first / starting block and then only the next N samples should be accumulated, since the next block already has Noverlap samples from previous step and the filter time delay is reduced:

    Chebyshev_formula11

    [2] is defined another time delay of the filter:

    Chebyshev_formula12

  • The Noverlap shall be choosed correct. This value should be equal to the length of the impulse response (transient) of the filter H(z), but the filter has an infinite impulse response, which means that the transient takes infinite time. On the other hand the filter stability condition:
    Chebyshev_formula12_1
    It means that h(n) decreases with time and one can find a Noverlap where the values of h(n) can be neglected without a significant error:

    Chebyshev_formula13

    In the next chapter I will analyse at the possibility of determining the transition time for the prototype analog filter which can be used to select Noverlap
    Thus the implementation of the TR is only an approximation since it is necessary to limit the time of the transition process
  • The implementation of the algorithm from Figure 2-1 leads to a sufficiently large time delay of the signal. See the (11) and (12). This imposes a restriction on the use of this method. For example, it is not possible to apply such a filter in control systems where a quick response is required. But in audio and image processing devices it can be effective.
  • It is possible to use a FIR filter with a linear phase. Before choosing a method we need to compare the amount of calculations (multiplications and additions) required for the FIR filter and the filter on the Figure 2-1. Usually the IIR filter is easier to implement but here we pass the signal through the IIR filter twice and need to do additional calculations for Noverlap samples.
  • [2] is calculated the number of the multiplications and additions per input sample (MAC) which required for the method:
    Chebyshev_formula14
  • I think that the method on Figure 2-1 can also be implemented efficiently using Fast Convolution. See the articles: “FIR Filter Implementation Using Octave GNU Tool and C Language“ and “Fast Fourier Transform (FFT) and C Implementation Using the Octave GNU Tool“. Perhaps I will return to this topic in my subsequent articles.

 

4. Estimate the Transient Process of the Analog Filter Prototype

For the time estimation of transients in analog circuits I want to use the method described in [4]. This book is now a bibliographic rarity and I have not seen its English translation.
To begin with consider the simplest RC circuit:

RC_analog_filter

Figure 4-1: RC Analog Low Pass Filter

 

For example: the RC circuit (low pass filter) is described with the first order differential equation:

Chebyshev_formula15

or in the Laplace form:

Chebyshev_formula16

Now let’s feed the Heaviside function 1(t) to the RC circuit input

Chebyshev_formula17

Laplace transform of the 1(t):

Chebyshev_formula18

Then using (16) and (18), we get the filter response at the output:

Chebyshev_formula19

This corresponds to the signal:

Chebyshev_formula20

The transition process takes approximately 5τ:

Chebyshev_formula21

The transition process is almost complete. The exponent values (21) reach about zero.
[4] is defined a generalized time constant as follows:

Chebyshev_formula22

The area under the graph g(t) – g (+∞) is given according to the area of the rectangle [g(t) – g(+∞)] τ.
(22) is correct if
Chebyshev_formula22_1
Let’s check the definition (22) for the RC filter:

Chebyshev_formula23

We will calculate the generalized time constant (22) for the circuit:

Chebyshev_formula24

Let’s use two properties of Laplace Transform:

Chebyshev_formula25_26

Now the step function Heaviside 1(t) will be fed to the input (24). Then the signal at the filter output is:

Chebyshev_formula27

Using (25) and (26) can be written:

Chebyshev_formula28

Note: (28) is true for the case n>m. If n=m then H (+∞)=bn/an

Chebyshev_formula29

Now we’ll calculate:

Chebyshev_formula30

Suppose t → +∞, which corresponds to s → 0 in (30), then (22) can be written using (26):

Chebyshev_formula31

Calculating (31) for (24):

Chebyshev_formula32

Using the 5τ rule for (32):

Chebyshev_formula33

Let’s check the formula (33) for (16):

Chebyshev_formula34

Then:

Chebyshev_formula35

Conclusions and notes
  • Formulas (32) and (33) can be used to select the Noverlap of the digital filter. The (32) and (33) describe a prototype analog filter
  • The proof of the formula (32) is somewhat different from the [4]. Because [4] uses the Carson Transform which is slightly different from the classic Laplace Transform:
    Chebyshev_formula36_37

 

5. Chebyshev Analog Filter Prototype

Now we will develop as an example:
Analog Low Pass Chebyshev Filter with the Fband=1kHz, fourth-order, 0.5dB Ripple in frequency band.
Using the [3] the normalized Laplace transform of the Chebyshev filter can be written:
Chebyshev_formula38_39
Using the (38), (39)
Chebyshev_formula39_1
Then

Chebyshev_formula40

Analog Chebyshev Filter (40) has a magnitude response in Figure 5.1, a phase response in Figure 5-2 and a step response in Figure 5-3.

AnalogChebyshevMagnitudeResponse_4th_order

Figure 5-1: Magnitude Response of the Analog Low Pass Chebyshev Filter, Fband=1kHz, 0.5dB Ripple, 4th order

 

The Magnitude response meets the filter requirements: Fband=1kHz, Ripple = 0.5dB

AnalogChebyshevPhaseResponse_4th_order

Figure 5-2: Phase Response of the Analog Low Pass Chebyshev Filter, Fband=1kHz, 0.5dB Ripple, 4th order

 

The phase response is non-linear.

AnalogChebyshevStepResponse_4th_order

Figure 5-3: Step Response of the Analog Low Pass Chebyshev Filter, Fband=1kHz, 0.5dB Ripple, 4th order

 

The step response has an overshoot about 1.18 (18%).
Now let’s represent (40) as (24) and use the formula (33):

Chebyshev_formula41

Figure 5-3 shows that after 2.4 ms there are two more fluctuations with levels of 0.97 (3%) and 1.016 (1.6%), so it is more reasonable to choose the transition time:

Chebyshev_formula42

 

6. Digital Chebyshev Filter

Using the analog prototype of the Chebyshev filter (40) can be calculated the digital filter. The bilinear transformation shall be used for it:

Chebyshev_formula43

Choose the sampling frequency:

Chebyshev_formula44

To calculate the bilinear transformation, we use the ready-made functions of the Octave GNU Tool:
1. Second Order Section (SOS):
[ZB1, ZA1] = bilinear ([1], [8.6987e-8, 4.1829e-4, 1], 83.3e-6)
ZB1 = [0.016343 0.032686 0.016343]
ZA1 = [1.00000 -1.60636 0.67173]
2. Second Order Section (SOS):
[ZB2, ZA2] = bilinear ([1], [2.9153e-8, 5.8060e-5, 1], 83.3e-6)
ZB2 = [0.052085 0.104169 0.052085]
ZA2 = [1.00000 -1.64645 0.85479]
Then digital Chebyshev Filter has follow z-Transform function with second order sections:
Chebyshev_formula45
Digital Chebyshev Filter (45) has a magnitude response in Figure 6.1, a phase response in Figure 6-2 and a step response in Figure 6-3.

DigitalChebyshevMagnitudeResponse_4th_order

Figure 6-1: Magnitude Response of the Digital Low Pass Chebyshev Filter, Fband=1kHz, 0.5dB Ripple, 4th order, fs=12kHz (Ts=83.3μs)

 

DigitalChebyshevPhaseResponse_4th_order

Figure 6-2: Phase Response of the Digital Low Pass Chebyshev Filter, Fband=1kHz, 0.5dB Ripple, 4th order, fs=12kHz (Ts=83.3μs)

 

DigitalChebyshevStepResponse_4th_order

Figure 6-3: Step Response of the Digital Low Pass Chebyshev Filter, Fband=1kHz, 0.5dB Ripple, 4th order

 

7. “Non-Causal“ Digital Chebyshev Filter with the Linear Phase

Now let’s analyze and implement the filter shown in Figure 2-1. Using (9), we get the Magnitude Response of the filter. See the Figure 7-1.

DigitaNonCausallChebyshevMagnitudeResp_4th_order

Figure 7-1: Magnitude Response of the non-causal Digital Low Pass Chebyshev Filter, Fband=1kHz, 0.5dB Ripple, 4th order, fs=12kHz (Ts=83.3μs). See the Figure 2-1

 

Comparing Figure 7-1 and Figure 6-1, we can see that the bandwidth has decreased: Fband: from 1 kHz to 950 Hz, and the graph has twice steep. For example: f=2kHz: from -38 dB to -76 dB.
Now calculate the Noverlap value using (42), (44):
Chebyshev_formula46
Let’s choose N block size to reduce the MAC value in formula (14) for this: N > > Noverlap. Take for example:
N=3 x Noverlap = 3 x 60 =180 samples
Now the non-causal filter from Figure 2-1 can be implemented using the Time Reverse described above. See the function non_casual_digitalFilter( ) in Octave Script and iir__non_casual_filter( ) in C Code.
Figure 7-2 shows the step response of the non-causal filter. The filter response is very similar to the response of an ideal low-pass filter with a linear phase (sync-filter). The filter response appears before the edge of the step signal. However this digital filter is physically realizable.
DigitaNonCausallChebyshevStepResp_4th_order

Figure 7-2: Step Response of the Non-Causal Low Pass Chebyshev Filter 4th Order. Blue – input, red – output

 

Now we will feed the sum of two sinusoids: 500Hz (inside the bandwidth) and 5kHz (outside the bandwidth) to the filter input. See the Figure 7-3. Only the 500Hz signal remains at the filter output. See the Figure 7-4

Signal500HzPlus5kHz

Figure 7-3: Sum of the two sines 500Hz and 5kHz. Filter Input Signal.

 

DigitaNonCausallChebyshev500HzPlus5kHzResp_4th_order

Figure 7-4: Response on the Sum of the two sines 500Hz and 5kHz. Non-Causal Low Pass Chebyshev Filter 4th Order. Filter Output Signal.

 

Comparing Figure 7-3 and Figure 7-4 we can see that the phase of the input and output signals of 500Hz does not change, which indicates a “zero” signal delay or a linear phase response. In reality the filter has a delay but it is related to the fact that the output signal can be calculated only when N + Noverlap input samples are accumulated. See the formula (10), (11). If we neglect the Tdsp time then for our example:

Chebyshev_formula47

 

8. Octave GNU Script DigitalChebyshevFilter.m

The script supports all plots of the article, generates the help files for the C code, calculates the non-causal digital filter.
Basic Functions:

% Plotting the Magnitude Response for the 4th order analog filter: 1/(a11*s^2+a12*s+a13)(a21*s^2+a22*s+a23)
% Input parameters:
% SA1 – denominator of the second order section
% SA2 – denominator of the second order section
% FhighInHz – Frequency space in Hz
function AnalogMagnitudeResponse(SA1, SA2, FhighInHz)

% Plotting the Phase Response for the 4th order analog filter: 1/(a11*s^2+a12*s+a13)(a21*s^2+a22*s+a23)
% Input parameters:
% SA1 – denominator 1 of the second order section
% SA2 – denominator 2 of the second order section
% FhighInHz – Frequency space in Hz
function AnalogPhaseResponse(SA1, SA2, FhighInHz)

% Plotting the Step Response for the 4th order analog filter: 1/(a11*s^2+a12*s+a13)(a21*s^2+a22*s+a23)
% Input parameters:
% SA1 – denominator 1 of the second order section
% SA2 – denominator 2 of the second order section
% TimeVector – Time space
function plot_StepResponseAnalogFilter(SA1, SA2, TimeVector)

% Plotting the Magnitude Response for the 4th order digital filter: (b11+b12*z^-1+b13*z^-2)(b21+b22*z^-1+b23*z^-2)/(a11+a12*z^-1+a13*z^-2)(a21+a22*z^-1+a23*z^-2)
% Input parameters:
% ZB1 – numerator 1 of the second order section
% ZA1 – denominator 1 of the second order section
% ZB2 – numerator 2 of the second order section
% ZA2 – denominator 2 of the second order section
% Tsample – sampling time in seconds
function plot_DigitalMagnitudeResponse(ZB1, ZA1, ZB2, ZA2, Tsample)

% Plotting the Phase Response for the 4th order digital filter: (b11+b12*z^-1+b13*z^-2)(b21+b22*z^-1+b23*z^-2)/(a11+a12*z^-1+a13*z^-2)(a21+a22*z^-1+a23*z^-2)
% Input parameters:
% ZB1 – numerator 1 of the second order section
% ZA1 – denominator 1 of the second order section
% ZB2 – numerator 2 of the second order section
% ZA2 – denominator 2 of the second order section
% Tsample – sampling time in seconds
function plot_DigitalPhaseResponse(ZB1, ZA1, ZB2, ZA2, Tsample)

% Plotting the Step Response for the 4th order digital filter: (b11+b12*z^-1+b13*z^-2)(b21+b22*z^-1+b23*z^-2)/(a11+a12*z^-1+a13*z^-2)(a21+a22*z^-1+a23*z^-2)
% Input parameters:
% ZB1 – numerator 1 of the second order section
% ZA1 – denominator 1 of the second order section
% ZB2 – numerator 2 of the second order section
% ZA2 – denominator 2 of the second order section
% SampleNumber – Sample number of the step function
function plot_StepResponseDigitalFilter(ZB1, ZA1, ZB2, ZA2, SampleNumber)

% Plotting the non-causal Magnitude Response for the 4th order digital filter: (b11+b12*z^-1+b13*z^-2)(b21+b22*z^-1+b23*z^-2)/(a11+a12*z^-1+a13*z^-2)(a21+a22*z^-1+a23*z^-2)
% Input parameters:
% ZB1 – numerator 1 of the second order section
% ZA1 – denominator 1 of the second order section
% ZB2 – numerator 2 of the second order section
% ZA2 – denominator 2 of the second order section
% Tsample – sampling time in seconds
function plot_non_causal_DigitalMagnitudeResponse(ZB1, ZA1, ZB2, ZA2, Tsample)

% Plotting the non-causal Magnitude Response for the 4th order digital filter: (b11+b12*z^-1+b13*z^-2)(b21+b22*z^-1+b23*z^-2)/(a11+a12*z^-1+a13*z^-2)(a21+a22*z^-1+a23*z^-2)
% Input parameters:
% ZB – filter numerator
% ZA – filter denominator
% Xinput – input samples
% Nblock – block size
% Nov – overlap size
function Youtput = non_causal_digitalFilter(ZB, ZA, Xinput, Nblock, Nov)

% Integer Test Signal to File
% Support C-Code
% Input parameters:
% SignalAmpl – amplitude of the test signal
% SignalFreqInHz – frequency of the test signal in Hz
% SignalLength – table length
% Tsample – sampling time in s
% FileNameString – output file name
% Output parameters:
% TestSignal – test signal vector
function TestSignal = IntegerTestSignals2file(SignalAmpl, SignalFreqInHz, SignalLength, Tsample, FileNameString)

% Integer Table to File
% Support C-Code
% Input parameters:
% ParameterVector – data array
% FileNameString – output file name
function IntegerParamVector2file(ParameterVector, FileNameString)

% Float Table to File
% Support C-Code
% Input parameters:
% ParameterVector – data array
% FileNameString – output file name
function FloatParamVector2file(ParameterVector, FileNameString)

% IIR Filter Coefficients Quantization
% Input parameters:
% ZA, ZB – IIR coefficients
% ScaleFactor – integer scale
% Output parameters:
% ZA_quantization – integer format
% ZB_quantization – integer format
function [ZA_quantization ZB_quantization] = IirFilterQuantization(ZA, ZB, ScaleFactor)

 

9. C Language non_causal_filter.c/h

There are two variants:
non_causal_filter_integer.c/h – fixed point arithmetic
non_causal_filter_float.c/h – float arithmetic
Base functions:

/*
IIR Non-Causal Filter Initialization
Clear the both histories:
iir_history_first and iir_history_second
Input
void
Return
void
*/
void iir_non_causal_init(void);

/*
IIR Non-Causal Filter
Input
sint16 buffer[N_BLOCK_SIZE+N_OVERLAP_SIZE] – pointer on the input sample buffer
buffer[0…N_OVERLAP_SIZE-1] – from previous processing step: first H(z) was already done
Return
buffer[N_BLOCK_SIZE] – filter output buffer
*/
void iir__non_causal_filter(sint16* buffer)

/* Test */
int main(void); – Check the filter with the step signal and sum of the two sines 500Hz and 5kHz

Note: The code is used the iir functions from the iir_filter.c/h. See the article “IIR Filter and C Implementation Using Octave GNU Tool“

 

10. Download the DigitalChebyshevFilter.m and non_causal_filter.c/h

You can download the files:
DigitalChebyshevFilter.m
non_causal_filter_integer.c/h
non_causal_filter_float.c/h
with the button:
Download
 

 

11. Literature / References

[1] L.Rabiner, B.Gold “Theory and application of digital signal processing“, Prentice-Hall, Ing Englewood Cliffs, New Jercy 1975
[2] Krzysztof Sozanski “A Linear-Phase IIR Filter for Audio Signal Interpolator“, Insitute of Electrical Engineering, University of Zielona Gora 65-246 Zielona Gora, Poland
[3] U. Tietze, Ch. Schenk “Halbleiter-Schaltungstechnik“, Springer-Verlag Berlin, Heidelberg New York, 1980
[4] K.P. Kovacs, I. Racz “Transiente Vorgaenge in Wechselstrommaschinen“, Verlag der ungarischen Akademie der Wissenschaften, Budapest, 1959

 

 

Views: 2,371