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.66k stars 253 forks source link

biquad_filter processing #186

Closed netlcod closed 1 year ago

netlcod commented 1 year ago

I am trying to apply a high-pass filter as shown in the code. The problem occurs when the signal has DC offset. The last filtered sample is an outlier.

#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 = 20;
    double amplitude = 1;
    double phase = 0;
    double offset = -10;

    std::unique_ptr<biquad_filter<fbase, 32>> bqfilter;
    double fc = 1;
    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 = 10;
    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);
        }
    }
}

image image

netlcod commented 1 year ago

In the line bqfilter->apply(batch); filtering gives wrong result. But if you replace it with bqfilter->apply(batch.data(), samples+1), the result becomes adequate. Result with samples+1: image

What i do wring? Or it is library bug?

netlcod commented 1 year ago

It seems to me that the library is the problem. When get_element gets a invalid index, it returns 0. So for the last sample it returns 0 instead of a valid value.

netlcod commented 1 year ago

In the file kfr/base/expression.cpp at line 726 begin_pass(in, inshape.adapt(start), inshape.adapt(stop)); the array size is reduced. If you replace it to the line begin_pass(in, start, stop); the filters work normally.

So I have a question, I need to replace begin_pass(in, inshape.adapt(start), inshape.adapt(stop)); to the line begin_pass(in, start, stop); or i need replace in file shape.hpp in line 283 return other.template trim<dims>()->min(**this - 1); to the return other.template trim<dims>()->min(**this);

netlcod commented 1 year ago

@dancazarin ?

dancazarin commented 1 year ago

@netlcod Thank you for reporting this issue. The adapt function is required here for multidimensional arrays (tensors), so it isn't an option to remove it. But the adapt function assumes a valid index passed to it which isn't a case for the stop index (past-the-end index). This also affected all expressions that relied on start and stop arguments in their begin_pass and end_pass functions. I've pushed a fix, could you try to pull dev branch and rerun your tests?

netlcod commented 1 year ago

@dancazarin yes, it works, thanks!