# Impulse response of Butterworth filter

In the following code example, a 6th-order Butterworth filter is defined over 4-50 Hz and the frequency response plotted. How to use DSP.jl to compute the corresponding zero-phase filter time response?
Using `filt(fil,ht1)` in code example below, it seems to provide the minimum-phase response, while the workaround using `ifft` and `fftshift` outputs what I expect but it is rather convoluted.

``````using DSP, FFTW, Plots

fs = 1000; Δt = 1.0/fs; # sampling frequency [Hz] amd period [s]
fc1, fc2 = 4, 50;  # low-cut / high-cut frequencies [Hz]
order = 6; # Butterworth filter order
fil = digitalfilter(Bandpass(fc1,fc2,fs=fs), Butterworth(order));
Hf = freqz(fil,0:1:200,fs);  # 0-200 Hz response
plot(20*log10.(abs.(Hf)),title="Frequency Response", label="Butterworth 6th order (4-50 Hz)",
xlabel="Frequency [Hz]", ylabel="Amplitude [dB]")

n = 500;
t = Δt * vec(-n÷2:n÷2-1);

# minimum phase impulse response?
ht1 = [1.0; zeros(n-1,1)];  # impulse at time 0
ht1 = filt(fil,ht1)
ht1 = [zeros(n÷2,1); ht1[1:n÷2]]
plot(t, ht1, title="Impulse Response", label="Min-phase (4-50 Hz)")

# 0-phase filter time response (using ifft)
ht0 = fftshift(ifft(fftshift(abs.(freqz(fil,-n÷2:n÷2-1,fs)))));  # over -250 to +250 Hz
plot(t, real.(ht0), title="0-phase filter Time Response", label="Zero-phase (4-50 Hz)",
xlabel="Time [s]", ylabel="Amplitude")
``````

2 Likes

Try constructing an order-3 filter and then apply it with filtfilt. This doubles the order while making the filter zero phase.

@mgkuhn, thank you for the feedback. I did not succeed in making DSP `filtfilt` to work in this case. However, the same principle of filtering in two-passes and flipping in time the first filter output, works if done manually:

``````t = Δt * vec(-(n-1):n-1);
ht = [zeros(n-1,1); 1.0; zeros(n-1,1)];  # impulse at time 0
fil3 = digitalfilter(Bandpass(fc1,fc2,fs=fs), Butterworth(3));
ht3 = filt(fil3,ht)
ht6 = filt(fil3, reverse(ht3, dims=1))
``````

This approach seems to be more restrictive (even filter order) than the workaround using `ifft` and `fftshift`. But I am not a `DSP` expert, just an interested user. The documentation for DSP.jl would benefit from more examples.

Works fine for me (Julia v1.5.3, DSP v0.6.9):

``````using Plots
using DSP

fs = 1000; Δt = 1.0/fs; # sampling frequency [Hz] amd period [s]
fc1, fc2 = 4, 50;  # low-cut / high-cut frequencies [Hz]
fil3 = digitalfilter(Bandpass(fc1,fc2,fs=fs), Butterworth(3));
n = 500;
ht = [zeros(n-1,1); 1.0; zeros(n-1,1)];  # impulse at time 0
y = filtfilt(fil3, ht);
plot(y)
``````

1 Like

Looks good. I had only tested `filtfilt` with an impulse on first sample: `ht = [1.0; zeros(n-1,1)]`, which seems to chop the signal.