# DSP Filter application

**URL:** <https://discourse.julialang.org/t/dsp-filter-application/55156>\
**Category:** Signal and Image Processing\
**Tags:** dsp\
**Created:** [February 12, 2021, 3:19pm UTC](https://discourse.julialang.org/t/dsp-filter-application/55156 "2021-02-12T15:19:10Z")\
**Posts on this page:** 14\
**Page:** 1

<div class="post-metadata">

**Author:** ![roble](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/roble/32/12303_2.png) [@roble](https://discourse.julialang.org/u/roble)\
**Post date:** [February 12, 2021, 3:19pm UTC](https://discourse.julialang.org/t/dsp-filter-application/55156/1 "2021-02-12T15:19:10Z")

</div>

Hi, I am trying to use the `filt(f, x[, si])` function from the [DSP.jl](https://docs.juliadsp.org/stable/filters/#DSP.filt) module.

My raw data looks like this:

 ![unfiltered](https://global.discourse-cdn.com/julialang/original/3X/e/b/eb90abd693819ed73c8b47715e0b57ed0bfe21f6.png)

On this data I apply a lowpass filter:

```julia
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:

 ![filtered](https://global.discourse-cdn.com/julialang/original/3X/3/0/3086730b2ec80ed2fbd52abd9267e4904d8981e7.png)

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:

```julia
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!

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [February 12, 2021, 3:50pm UTC](https://discourse.julialang.org/t/dsp-filter-application/55156/3 "2021-02-12T15:50:09Z")

</div>

Using `filtfilt` helps for this problem:

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

```

![filtfilt_filtered_input](https://global.discourse-cdn.com/julialang/original/3X/b/a/bafc6c223d251b05ec25e17426d452871d025cc2.png)

---

<div class="post-metadata">

**Author:** ![roble](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/roble/32/12303_2.png) [@roble](https://discourse.julialang.org/u/roble)\
**Post date:** [February 12, 2021, 4:03pm UTC](https://discourse.julialang.org/t/dsp-filter-application/55156/4 "2021-02-12T16:03:53Z")

</div>

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

 ![filtered_correct](https://global.discourse-cdn.com/julialang/original/3X/5/b/5bb97f2f77a5bb5ca2f7ad5f3092d0a1c1bc8f23.png)

---

<div class="post-metadata">

**Author:** ![sairus7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sairus7/32/10816_2.png) [@sairus7](https://discourse.julialang.org/u/sairus7)\
**Post date:** [February 12, 2021, 7:52pm UTC](https://discourse.julialang.org/t/dsp-filter-application/55156/5 "2021-02-12T19:52:16Z")

</div>

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

```julia
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):

```julia
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):

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

```

---

<div class="post-metadata">

**Author:** ![roble](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/roble/32/12303_2.png) [@roble](https://discourse.julialang.org/u/roble)\
**Post date:** [February 13, 2021, 10:32pm UTC](https://discourse.julialang.org/t/dsp-filter-application/55156/6 "2021-02-13T22:32:15Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [February 14, 2021, 3:00pm UTC](https://discourse.julialang.org/t/dsp-filter-application/55156/7 "2021-02-14T15:00:44Z")

</div>

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

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

```

![filtfilt_and_init_state_filt_results](https://global.discourse-cdn.com/julialang/original/3X/8/4/840b1f77837dd37539e1ed7ca3a1c81ea87b31b0.png)

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

---

<div class="post-metadata">

**Author:** ![sairus7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sairus7/32/10816_2.png) [@sairus7](https://discourse.julialang.org/u/sairus7)\
**Post date:** [February 15, 2021, 5:25am UTC](https://discourse.julialang.org/t/dsp-filter-application/55156/8 "2021-02-15T05:25:45Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![RoyiAvital](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/royiavital/32/571_2.png) [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)\
**Post date:** [February 15, 2021, 6:08am UTC](https://discourse.julialang.org/t/dsp-filter-application/55156/9 "2021-02-15T06:08:33Z")

</div>

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()`.

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [February 15, 2021, 10:15pm UTC](https://discourse.julialang.org/t/dsp-filter-application/55156/10 "2021-02-15T22:15:53Z")

</div>

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

```julia
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:  
 ![FFT_0phase_filt_filtfilt_and_init_state_filt](https://global.discourse-cdn.com/julialang/original/3X/e/2/e25af5108127fae0b44fd3708a63afcfdcc80c1c.png)

---

<div class="post-metadata">

**Author:** ![RoyiAvital](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/royiavital/32/571_2.png) [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)\
**Post date:** [February 16, 2021, 5:38am UTC](https://discourse.julialang.org/t/dsp-filter-application/55156/11 "2021-02-16T05:38:09Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [February 16, 2021, 6:48am UTC](https://discourse.julialang.org/t/dsp-filter-application/55156/12 "2021-02-16T06:48:28Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![roble](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/roble/32/12303_2.png) [@roble](https://discourse.julialang.org/u/roble)\
**Post date:** [February 17, 2021, 4:07pm UTC](https://discourse.julialang.org/t/dsp-filter-application/55156/13 "2021-02-17T16:07:49Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [February 17, 2021, 4:26pm UTC](https://discourse.julialang.org/t/dsp-filter-application/55156/14 "2021-02-17T16:26:02Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![roble](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/roble/32/12303_2.png) [@roble](https://discourse.julialang.org/u/roble)\
**Post date:** [February 17, 2021, 4:28pm UTC](https://discourse.julialang.org/t/dsp-filter-application/55156/15 "2021-02-17T16:28:45Z")

</div>

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!