kfrlib / kfr

Fast, modern C++ DSP framework, FFT, Sample Rate Conversion, FIR/IIR/Biquad Filters (SSE, AVX, AVX-512, ARM NEON)
https://www.kfrlib.com
GNU General Public License v2.0
1.62k stars 248 forks source link

biquad_filter packet processing #177

Closed netlcod closed 1 year ago

netlcod commented 1 year ago

My task is to filter the incoming values generated in a batch. For example, the packet size is 100 samples. Figure "generated1" - visualization of the signal under test (1 Hz sine wave, sampling rate 100 Hz). Figure "generated2" - visualization of 100 packets.

generated1: generated1 generated2: generated2

I use an IIR filter of second order (in this case a highpass filter with fc = 20 Hz). Figure "filtered1" is the filter's response to one packet. Figure "filtered2" - processing of 100 packets with the filter. The "total" picture is a cropping for several packets (the filter's response is amplified for clarity). filtered1: filtered1 filtered2: filtered2 total: total

The problem is that the filter does not keep its internal state between calls when processing incoming packets. In the picture the "expected" is what I expect from my filter. So the question is how do I keep the filter state between calls? expected:
expected

My test code is attached below.

#include <QDebug>

#include <kfr/base.hpp>
#include <kfr/dsp.hpp>
#include <kfr/base/generators.hpp>

using namespace kfr;

int main()
{
    double fs = 100;
    const int samples = 100;
    double dt = 1.0 / fs;

    double frequency = 1;
    double amplitude = 1;
    double phase = 0;
    double offset = 0;

    std::unique_ptr<biquad_filter<fbase, 32>> bqfilter;
    double fc = 20;
    std::vector<biquad_params<fbase>> bqs = to_sos(iir_highpass(butterworth<fbase>(12), fc, fs));
    bqfilter.reset(new biquad_filter<fbase, 32>(bqs));

    univector<fbase> generated;
    univector<fbase> filtered;

    int iterations = 100;
    double last = 0;
    for (int i = 0; i < iterations; i++) {
        kfr::univector<kfr::fbase> time = kfr::linspace(last, last+samples*dt, samples, false, kfr::ctrue);
        kfr::univector<kfr::fbase> batch = amplitude * kfr::sin(2*kfr::c_pi<double>*frequency*time + phase) + offset;
        foreach(fbase val, batch) {
            generated.emplace_back(val);
        }

        bqfilter->apply(batch);
        last += samples * dt;

        foreach(fbase val, batch) {
            filtered.emplace_back(val);
        }
    }
    qDebug() << "end";
}
netlcod commented 1 year ago

This example also doesn't work as I expect https://github.com/kfrlib/kfr/issues/59#issuecomment-584579088

netlcod commented 1 year ago

Hello! I solved the problem. The problem was not in the filter code, but in the signal generator code. When creating time samples: kfr::univector<kfr::fbase> time = kfr::linspace(m_last, m_last+m_batchSize*m_dt, m_batchSize, true, kfr::ctrue); 4 the argument of linspace (endpoint)should have been set to true