Impulse response of Butterworth filter

[rafael.guerra]

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")

Thanks in advance.

[mgkuhn]

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

[rafael.guerra]

@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.

[mgkuhn]

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)

[rafael.guerra]

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.