“NonCausal“ 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:
The (1) means that all the filter poles must be inverted:
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 noncausal 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 21 shows a noncausal filter with a zero phase (linear phase).
Figure 21: NonCausal IIR Filter with Null Phase/Delay Response
We describe the signal spectra at the output of each block of the noncausal filter using the signal property of the y(n) and y ( n) which have conjugate spectrum

Spectrum of the g(n):

Spectrum of the g(n):

Spectrum of the y(n):

Spectrum of the output signal y(n):
(7) matches the filter shown in Figure 22.
Figure 22: Equivalent filter for Figure 21
Notes:

The filter (7) has a zero phase (linear phase). It is not possible to implement a zerolatency 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):
Noncausal filter:
Ripples of the noncausal filter is verdoppelt, and also reduces the bandwidth compared to H(z). Usually the filter bandwidth is defined at 3dB. The noncausal filter has a passband of the original filter 6dB (verdoppelt).
3. Time Reverse (TR) Implementation from [2]
See the Figure 21. 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 31.
Figure 31: 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 32.
Figure 32: Reverse Time for g(n) Blocks
At the output second H(z) the y ( n) is got . See the Figure 33.
Figure 33: Reverse Time for y(n) Blocks
The subsequent time reversion for y( n) blocks gets the output signal y(n). See the Figure 34.
Figure 34: 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 21 can be started. Thus the filter has the delay:
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:
[2] is defined another time delay of the filter:

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:
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: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 21 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 21. 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:

I think that the method on Figure 21 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:
Figure 41: RC Analog Low Pass Filter
For example: the RC circuit (low pass filter) is described with the first order differential equation:
or in the Laplace form:
Now let’s feed the Heaviside function 1(t) to the RC circuit input
Laplace transform of the 1(t):
Then using (16) and (18), we get the filter response at the output:
This corresponds to the signal:
The transition process takes approximately 5τ:
The transition process is almost complete. The exponent values (21) reach about zero.
[4] is defined a generalized time constant as follows:
The area under the graph g(t) – g (+∞) is given according to the area of the rectangle [g(t) – g(+∞)] τ.
(22) is correct if
Let’s check the definition (22) for the RC filter:
We will calculate the generalized time constant (22) for the circuit:
Let’s use two properties of Laplace Transform:
Now the step function Heaviside 1(t) will be fed to the input (24). Then the signal at the filter output is:
Using (25) and (26) can be written:
Note: (28) is true for the case n>m. If n=m then H (+∞)=bn/an
Now we’ll calculate:
Suppose t → +∞, which corresponds to s → 0 in (30), then (22) can be written using (26):
Calculating (31) for (24):
Using the 5τ rule for (32):
Let’s check the formula (33) for (16):
Then:
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:
5. Chebyshev Analog Filter Prototype
Now we will develop as an example:
Analog Low Pass Chebyshev Filter with the Fband=1kHz, fourthorder, 0.5dB Ripple in frequency band.
Using the [3] the normalized Laplace transform of the Chebyshev filter can be written:
Using the (38), (39)
Then
Analog Chebyshev Filter (40) has a magnitude response in Figure 5.1, a phase response in Figure 52 and a step response in Figure 53.
Figure 51: 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
Figure 52: Phase Response of the Analog Low Pass Chebyshev Filter, Fband=1kHz, 0.5dB Ripple, 4th order
The phase response is nonlinear.
Figure 53: 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):
Figure 53 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:
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:
Choose the sampling frequency:
To calculate the bilinear transformation, we use the readymade functions of the Octave GNU Tool:
1. Second Order Section (SOS):
[ZB1, ZA1] = bilinear ([1], [8.6987e8, 4.1829e4, 1], 83.3e6)
ZB1 = [0.016343 0.032686 0.016343]
ZA1 = [1.00000 1.60636 0.67173]
2. Second Order Section (SOS):
[ZB2, ZA2] = bilinear ([1], [2.9153e8, 5.8060e5, 1], 83.3e6)
ZB2 = [0.052085 0.104169 0.052085]
ZA2 = [1.00000 1.64645 0.85479]
Then digital Chebyshev Filter has follow zTransform function with second order sections:
Digital Chebyshev Filter (45) has a magnitude response in Figure 6.1, a phase response in Figure 62 and a step response in Figure 63.
Figure 61: Magnitude Response of the Digital Low Pass Chebyshev Filter, Fband=1kHz, 0.5dB Ripple, 4th order, fs=12kHz (Ts=83.3μs)
Figure 62: Phase Response of the Digital Low Pass Chebyshev Filter, Fband=1kHz, 0.5dB Ripple, 4th order, fs=12kHz (Ts=83.3μs)
Figure 63: Step Response of the Digital Low Pass Chebyshev Filter, Fband=1kHz, 0.5dB Ripple, 4th order
7. “NonCausal“ Digital Chebyshev Filter with the Linear Phase
Now let’s analyze and implement the filter shown in Figure 21. Using (9), we get the Magnitude Response of the filter. See the Figure 71.
Figure 71: Magnitude Response of the noncausal Digital Low Pass Chebyshev Filter, Fband=1kHz, 0.5dB Ripple, 4th order, fs=12kHz (Ts=83.3μs). See the Figure 21
Comparing Figure 71 and Figure 61, 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):
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 noncausal filter from Figure 21 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 72 shows the step response of the noncausal filter. The filter response is very similar to the response of an ideal lowpass filter with a linear phase (syncfilter). The filter response appears before the edge of the step signal. However this digital filter is physically realizable.
Figure 72: Step Response of the NonCausal 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 73. Only the 500Hz signal remains at the filter output. See the Figure 74
Figure 73: Sum of the two sines 500Hz and 5kHz. Filter Input Signal.
Figure 74: Response on the Sum of the two sines 500Hz and 5kHz. NonCausal Low Pass Chebyshev Filter 4th Order. Filter Output Signal.
Comparing Figure 73 and Figure 74 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:
8. Octave GNU Script DigitalChebyshevFilter.m
The script supports all plots of the article, generates the help files for the C code, calculates the noncausal 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 noncausal 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 noncausal 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 CCode
% 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 CCode
% Input parameters:
% ParameterVector – data array
% FileNameString – output file name
function IntegerParamVector2file(ParameterVector, FileNameString)
% Float Table to File
% Support CCode
% 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 NonCausal Filter Initialization
Clear the both histories:
iir_history_first and iir_history_second
Input
void
Return
void
*/
void iir_non_causal_init(void);
/*
IIR NonCausal Filter
Input
sint16 buffer[N_BLOCK_SIZE+N_OVERLAP_SIZE] – pointer on the input sample buffer
buffer[0…N_OVERLAP_SIZE1] – 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:
11. Literature / References
[1] L.Rabiner, B.Gold “Theory and application of digital signal processing“, PrenticeHall, Ing Englewood Cliffs, New Jercy 1975
[2] Krzysztof Sozanski “A LinearPhase IIR Filter for Audio Signal Interpolator“, Insitute of Electrical Engineering, University of Zielona Gora 65246 Zielona Gora, Poland
[3] U. Tietze, Ch. Schenk “HalbleiterSchaltungstechnik“, SpringerVerlag Berlin, Heidelberg New York, 1980
[4] K.P. Kovacs, I. Racz “Transiente Vorgaenge in Wechselstrommaschinen“, Verlag der ungarischen Akademie der Wissenschaften, Budapest, 1959