How to initialize the filter's step response

diegoferigo commented 3 years ago

Fist off, thanks a lot for the project! It seems nice and clean, and importing it in another CMake projects is as easy as:

FetchContent Usage

```cmake find_package(iir QUIET) option(USE_SYSTEM_IIR "Use system-installed IIR rather than a private copy" ${iir_FOUND}) if(USE_SYSTEM_IIR AND NOT ${iir_FOUND}) message(FATAL_ERROR "Failed to find system iir") endif() if(NOT ${USE_SYSTEM_IIR}) include(FetchContent) FetchContent_Declare(iir GIT_REPOSITORY https://github.com/berndporr/iir1) FetchContent_GetProperties(iir) if(NOT iir_POPULATED) FetchContent_Populate(iir) add_subdirectory(${iir_SOURCE_DIR} ${iir_BINARY_DIR} EXCLUDE_FROM_ALL) endif() endif() add_library(MyLib SHARED MyLib.h MyLib.cpp) target_link_libraries(MyLib PRIVATE iir::iir_static) ```

This being said, I'm facing a problem. A scipy filter like the following:

from scipy import signal

order = 5
fs = 1000.0
nyq = fs / 2.0
cutoff = 100.0

b, a = signal.butter(N=order, Wn=(cutoff / nyq), btype="low")

could be called (simulating real-time data) with:

for x in data:
    out = signal.lfilter(b=b, a=a, x=[x])

Now, if I want to initialize the initial conditions of filter's step response, I can do the following:

z = data[0] * signal.lfilter_zi(b=b, a=a)

for x in data:
    out, z = signal.lfilter(b=b, a=a, x=[x], zi=z)

@berndporr how can I do the same with iir1?

berndporr commented 3 years ago

The filter is a series of 2nd order biqads. The values of the individual states are stored here: https://github.com/berndporr/iir1/blob/master/iir/Cascade.h#L180 but depend which topography you are using and for that reason not exposed to the outside. Bottomline: sorry. No possible.

diegoferigo commented 3 years ago

I see, thanks for the reply! It would be nice having the possibility to access that state somehow, so that from downstream we could implement a custom initialization if needed. Maybe a StateType& getStates(); getter?

If I understood well, that array would match the zi array used in scipy example above, right?

diegoferigo commented 3 years ago

For the moment, another workaround is running the filter to convergence when the first sample is received:

    Iir::Butterworth::LowPass<5> butterworth{};

    try {
        butterworth.setup(/*sampleRate=*/1000, /*cutoffFrequency=*/100.0);
    }
    catch (const std::exception& e) {
        std::cerr << "Wrong parameters" << std::endl;
        return;
    }

    auto initialize = [&butterworth](const double target,
                                     const double threshold = 1e-3,
                                     const size_t convergedSteps = 50) {
        for (size_t i = 0; i < convergedSteps; ++i) {
            if (std::abs(butterworth.filter(target) - target)
                > threshold * std::abs(target)) {
                i = 0;
            }
        }
    };

    const double firstSample = 42.0;
    initialize(firstSample);

berndporr commented 3 years ago

A few points:

there is no 1:1 match from the "zi" to the internal structure as lfilter uses a MATLAB style N-order filter structure. You can switch to SOS though in python which would make it compatible but I've got no idea how the "zi" then map into the 2nd order states. So that would just create confusion and misunderstandings.
Trying to force the filter to a steady state by forcing its internal states to that value creates unexpected behaviour. I'd say you have been just lucky but because it's a recursive system in general there will be oscillations after forcing the delay lines to another const value and the internal values are non-trivially connected to the input. Even more with a chain of 2nd order filters.
So your "workaround" is better and given the frequency of your filter you also know how long it takes to settle. Of course even your pre-seeding with 42 will make the filter oscillate. Imagine the next reading is 999. Then the filter will generate a massive step response. Or imagine the next samples will a sine-wave after 42. Then the filter will slowly settle to filter that sine-wave until it has observed its past.
So not even the "workaround" is a good idea but rather that you take your real signal, feed that into your filter and let the filter settle on that. Discard these samples till it's settled.
Looking at your code this won't work and will throw an exception as your cutoff is way above Nyquist (max cutoff is 0.0005).
You get filter order 6 for free because it's a chain of 2nd order filters.

diegoferigo commented 3 years ago

there is no 1:1 match from the "zi" to the internal structure as lfilter uses a MATLAB style N-order filter structure. You can switch to SOS though in python which would make it compatible but I've got no idea how the "zi" then map into the 2nd order states. So that would just create confusion and misunderstandings.

Trying to force the filter to a steady state by forcing its internal states to that value creates unexpected behaviour. I'd say you have been just lucky but because it's a recursive system in general there will be oscillations after forcing the delay lines to another const value and the internal values are non-trivially connected to the input. Even more with a chain of 2nd order filters.

I see, for the reference, scipy.signal.lfilter_zi is the scipy documentation I originally used in the example of the first post.

So your "workaround" is better and given the frequency of your filter you also know how long it takes to settle. Of course even your pre-seeding with 42 will make the filter oscillate. Imagine the next reading is 999. Then the filter will generate a massive step response. Or imagine the next samples will a sine-wave after 42. Then the filter will slowly settle to filter that sine-wave until it has observed its past.

So not even the "workaround" is a good idea but rather that you take your real signal, feed that into your filter and let the filter settle on that. Discard these samples till it's settled.

Thanks for the clarification. If it was not clear, in the MWE above, I meant 42 being the first data sample of the real signal.

Looking at your code this won't work and will throw an exception as your cutoff is way above Nyquist (max cutoff is 0.0005).

You're right, I copied-and-pasted that snippet and I mistakenly replaced the sample rate with the sample time. I updated the example to be sound.

berndporr commented 3 years ago

As I said before. You cannot force an IIR filter with one "real" sample into a steady state. If the filter looks 1000 samples in the past for example 100Hz cutoff at Fs=1kHz then you need to send in at least 10 samples of your actual signal to settle and not just one sample. What you do is that you make your filter believe it will filter DC forever what I'm sure is not its purpose.

user706 commented 4 months ago

To others coming by here ... I recommend using a different library that can actually have the 1st output set to a desired value, i.e. proper initialisation:

rtfilter, which is also available on debian, ubuntu, etc. See also stackoverflow

To initialize the filter, such that the 1st output is a specific value, you can use the following rtfilter functions:

// filter.h

void rtf_filter_set_lazy_init(hfilter filt, int do_lazy_init);

void rtf_init_filter(hfilter filt, const void* data);

(We can do it in matlab, python etc. so of course we want it in C/C++.)

berndporr commented 4 months ago

It's all done here: https://github.com/mmlabs-mindmaze/rtfilter/blob/master/src/init-filter-func-template.c how about you write a nice pull request for this library? ;)

user706 commented 4 months ago

It's all done here: https://github.com/mmlabs-mindmaze/rtfilter/blob/master/src/init-filter-func-template.c how about you write a nice pull request for this library? ;)

Will that work? Is rtfilter not based on a different method, using a, b instead of sos??

Is iir1 based on sos? And if so, what exactly does it mean when you write

You can switch to SOS though in python which would make it compatible but I've got no idea how the "zi" then map into the 2nd order states. So that would just create confusion and misunderstandings.

?

Does that mean that python can do it with SOS (ref), but you don't understand how they do it??

user706 commented 4 months ago

Here's the python guys working on SOS and initial conditions (a nice read!):

https://github.com/scipy/scipy/issues/2444 (search for initial)

References this: https://github.com/scipy/scipy/pull/3717/files (note sosfilt_zi)

user706 commented 4 months ago

berndporr / iir1

How to initialize the filter's step response #24