DSP Filter application

[roble]

Hi, I am trying to use the filt(f, x[, si]) function from the DSP.jl module.

My raw data looks like this:

On this data I apply a lowpass filter:

df.data_raw  # my raw data
w0 = 1/120 # cutoff frequency
df.data_lp= filt(digitalfilter(Lowpass(w0, fs=1/5), Butterworth(4)), df.data_raw )

The filtered data now looks like this:

Now my problem is that the filter assumes a step between the beginning and the first value of my data_raw vector. This results in a slowly rising line right in the beginning of the filtered data. I have tried using the si argument to define the initial filter value. Unfortunately it resulted in an error:

w0 = 1/120
df.data_lp = filt(digitalfilter(Lowpass(w0, fs=1/5), Butterworth(4)), df.data_raw, [df.data_raw[1]])

MethodError: no method matching filt(::DSP.Filters.ZeroPoleGain{Complex{Float64},Complex{Float64},Float64}, ::Array{Float64,1}, ::Array{Float64,1})

Closest candidates are:

filt(!Matched::Union{Number, AbstractArray{T,1} where T}, ::Union{Number, AbstractArray{T,1} where T}, ::AbstractArray{T,N} where N) where T at C:\Users\roble\.julia\packages\DSP\q9iEF\src\dspbase.jl:16

filt(!Matched::Union{Number, AbstractArray{T,1} where T}, ::Union{Number, AbstractArray{T,1} where T}, ::AbstractArray{T,N} where N, !Matched::AbstractArray{S,N} where N) where {T, S} at C:\Users\roble\.julia\packages\DSP\q9iEF\src\dspbase.jl:16

filt(!Matched::DSP.Filters.PolynomialRatio, ::Any, ::Any) at C:\Users\roble\.julia\packages\DSP\q9iEF\src\Filters\filt.jl:35

...

Does anyone know how to correctly use the filt function to achieve this?

Many thanks in advance!

[rafael.guerra]

Using filtfilt helps for this problem:

using DSP, Plots
t = 1:1000; fs = 1/5; f0 = 1/750   # sampling and cutoff frequencies
rawdata = 500*(exp.(-t/500) .+ 0.1*rand(length(t)))  # noisy input data
data_filt = filtfilt(digitalfilter(Lowpass(f0, fs=fs), Butterworth(4)), rawdata)
plot(rawdata,label="input")
plot!(data_filt, lc=:blue, lw=2,label="filtered using filtfilt")

[roble]

Wow, that worked well. Many thanks for your help!

[sairus7]

Not necessarily use filtfilt - you can setup initial state equal to filter output for repeated first point, like this:

Assuming Transposed-Direct-Form-II filter:
for lowpass or bandstop - for constant input signal filter outputs the same constant (x = 1, y = 1):

s0[order] = b[order + 1] - a[order + 1]
for k in order:-1:2
    s0[k - 1] = b[k] - a[k] + s0[k]
end

for highpass or bandpass - for constant input signal filter outputs zero (x = 1, y = 0):

s0[order] = b[order + 1]
for k in order:-1:2
    s0[k - 1] = b[k] + s0[k]
end

Then, before applying filter you should multiply s0 by first point amplitude (x = x[1], y = x[1] or y = 0):

for k in 1:order
    s0[k] = s0[k] * x[1]
end

[roble]

So this is how i can find the initial filter state. Thank you for your explanations

As of now, i have an understanding of analog transfer functions. However, i am new to digital filtering, hence i do not know what a transposed-direct-form-II filter is…

I’d like to learn more about this. Can you recommend a book which covers the theory?

[rafael.guerra]

@sairus7, following your initial state recipe the following result was obtained:

using DSP, Plots
N = 1000; t = 1:N; fs = 1/5; f0 = 1/750   # sampling and cutoff frequencies
rawdata = 500*(exp.(-t/500) .+ 0.1*rand(length(t)))  # noisy input data
lpfilt = digitalfilter(Lowpass(f0, fs=fs), Butterworth(4))
data_filt1 = filtfilt(lpfilt, rawdata)

b, a = coefb(lpfilt), coefa(lpfilt)
order = 4
s0 = zeros(order)
s0[order] = b[order + 1] - a[order + 1]
for k in order:-1:2
    s0[k - 1] = b[k] - a[k] + s0[k]
end
s0 *= rawdata[1]
data_filt0 = filt(b, a, rawdata, s0) 

plot(rawdata,label="input")
plot!(data_filt1, lc=:blue, lw=2,label="filtfilt")
plot!(data_filt0, lc=:red, lw=2,label="filt with initial state")

Could you please comment on this result. Thanks in advance.

[sairus7]

That’s because it is assumed that before the first point there was a constant level equal to first point. This assumption is better than zero, but is not ideal.

So in your case it abruptly begins to slope down from initial constant level, and a lowpass filter slowly trying to catch up (phase delay).

In this case, zero-phase filter is better, if you can use it. That means you have the whole chunk of signal beforehand, so you can filter it both in forward and backward directions.

[RoyiAvital]

One thing to remember when using filtfilt() is that the applied filter isn’t exactly the designed filter (It is its Squared Magnitude). So beside having all samples at once it also changes the ripple property of the filter.
In most cases it doesn’t make any difference yet in some it does (Or at least must be taken into account).

By the way, instead of using the general form of the filtering (Having both b and a) one could use convolution with boundary conditions. For this task, instead of replicate boundary one might use mirror and get similar effect to filtfilt().

[rafael.guerra]

@RoyiAvital, I came out with the following alternative to filtfilt by making the input signal symmetric and by zero’ing the phase of the filter in the Fourier domain. I believe this produces a single-pass 0-phase filter and so, its spectral amplitude is not squared.

PS: edited the FFT 0-phase filter code. Depending on the actual signal & noise, the result (see updated figure) may be better than with filtfilt.

using DSP, FFTW, Plots
import DSP.freqz

N = 1000; t = 1:N; fs = 1/5; f0 = 1/750   # sampling and cutoff frequencies
y = 500*(exp.(-t/500) .+ 0.1*rand(length(t)))  # noisy input y
lpfilt = digitalfilter(Lowpass(f0, fs=fs), Butterworth(4))

# (a) filt with initial state by @sairus7
b, a = coefb(lpfilt), coefa(lpfilt)
order = 4
s0 = zeros(order)
s0[order] = b[order + 1] - a[order + 1]
for k in order:-1:2
    s0[k - 1] = b[k] - a[k] + s0[k]
end
s0 *= y[1]
y_filt0 = filt(b, a, y, s0)

# (b) filtfilt
y_filt1 = filtfilt(lpfilt,y)

# (c) rg's FFT 0-phase filter
y2 = [y[end:-1:1]; y]   # make real input symmetric around t=0
f0 = fftshift(LinRange(-fs/2,fs/2,2N+1)[1:2N])
Hf = abs.(freqz(lpfilt,f0,fs)) .* fft(y2)
y_filt2 = real.((ifft(Hf)))[N+1:end]

plot(y2[N+1:end], lw=0.5, lc=:grey, xl="Time", label="Input signal with random noise")
plot!(y_filt0, lc=:red, lw=2,label="(a) filt with initial state (Butterworth LP 4th order)")
plot!(y_filt1, lc=:blue, lw=3, yl="Amplitude", label="(b) filtfilt (Butterworth LP 8th order)")
plot!(y_filt2, lc=:lime, lw=2,label="(c) FFT 0-phase filt (Butterworth LP 4th order)")

The new result is given by the green curve in plot below:

[RoyiAvital]

It looks decent.
Another approach, for this case at least, can be a parametric approach.
It looks like you will be able to easily create a model for this kind of a signal and just estimate the parameters.

[rafael.guerra]

The input signal examples in this post seem extreme, in the sense that they do not look like the typical ones found in DSP problems. There is often a data acquisition anti-alias filter and we deal with continuous, band-limited, oscillatory signals. More appropriate smoothing or fitting techniques exist for the type of data posted. The goal (for me) was to test the DSP limits/envelope.

[roble]

The last solution is really neat! I did not know that you could omit the phase and just use the amplitude in the frequency domain. The advantage of this one is that it applies the filter only once, wheras (I believe) filtfilt results in a squared filter amplitude.

In my case, the differences in performance is not that important. Hence, I do not need to do benchmarks.

To continue, I chose the filtfilt solution since it has the null-phase property and is simple to implement. Since most solutions are nicely presented and compared, I will set this post as the final solution.

[rafael.guerra]

Please note that the filter was constructed so that its amplitude spectrum is symmetric and has only real numbers. Its phase should be zero.

[roble]

They shure do seem a bit extreme. The data shows the ambient light measurement during an afternoon. By nature, the signal is noisy and can have fluctuations (in case of a shadow or cloud) faster than the nyquist frequency of my measurement. As proposed by @RoyiAvital I could have tried to fit a model. However with a model, It is difficult to regard the effects of dawn and mist.

The main purpose of filtering the data is to remove fluctuations faster than a certain frequency. Thanks to your help, I could do it. For me, it was a pleasure to learn more about digital filters!