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

transfer function coefficients and v.s. Matlab #148

Open BinDong314 opened 2 years ago

BinDong314 commented 2 years ago

Hi, KFR Community, Thanks for the help in advance for two issues.

1) how to get the transfer function coefficients of IIR, such as "the 8th-order Chebyshev type I filter, lowpass" below

   plot_save("chebyshev1_lowpass8", output,
              options + ", title='8th-order Chebyshev type I filter, lowpass'");

    {
        zpk<fbase> filt                       = iir_lowpass(chebyshev2<fbase>(8, 80), 0.09);
        std::vector<biquad_params<fbase>> bqs = to_sos(filt);
        output                                = biquad<maxorder>(bqs, unitimpulse());
    }

Code is from : https://github.com/kfrlib/kfr/blob/master/examples/iir.cpp

I understand that "filt" contains zeor-pole-gain forms, any help function inside to get transfer function coefficients?

2) I print out some information of "filt" for the " the 12th-order Butterworth filter, bandpass' below. They are different from Matlab code. Any hint to get the same results?

   zpk<fbase> filt = iir_bandstop(butterworth<fbase>(2), 0.0200, 0.6400);
    std::cout << "filt.k = " << filt.k << "\n";
    univector<complex<fbase>> z = filt.z;
    univector<complex<fbase>> p = filt.p;

    std::cout << "filt.z =  \n";

    for (int i = 0; i < z.size(); i++)
    {
        std::cout << z[i].real() << " + i " << z[i].imag() << "\n";
    }

    std::cout << "filt.p =  \n";

    for (int i = 0; i < z.size(); i++)
    {
        std::cout << p[i].real() << " + i " << p[i].imag() << "\n";
    }

    std::vector<biquad_params<fbase>> bqs = to_sos(filt);

    for (int i = 0; i < bqs.size(); i++)
    {
        biquad_params<fbase> norm = bqs[i].normalized_all();
        std::cout << "a =" << norm.a0 << ", " << norm.a1 << ", " << norm.a2 << " b=" << norm.b0 << ", "
                  << norm.b1 << ", " << norm.b2 << "\n";
    }
`

 Output from KFR: 
```cpp
filt.k = 0.190617
filt.z =
0.905633 + i 0.424062
0.905633 + i 0.424062
0.905633 + i -0.424062
0.905633 + i -0.424062
filt.p =
0.955593 + i -0.0442401
0.955593 + i 0.0442401
-0.251103 + i 0.403472
-0.251103 + i -0.403472

Output from Matlab:

>> fbands=[0.0200    0.6400]
>> [bfz1,bfp1,bfk1]=butter(2,fbands,'bandpass')

bfz1 =

     1
     1
    -1
    -1

bfp1 =

  -0.2511 + 0.4035i
  -0.2511 - 0.4035i
   0.9556 + 0.0442i
   0.9556 - 0.0442i

bfk1 =

    0.4127

Bests, Bin

mhagh3448 commented 1 year ago

Hello, I downloaded the file and put the headers in the include path, but when using the filter, I encountered a return 1 error.. Please help me... I have been facing this error for several days and I am new

Jalmenara commented 3 months ago

Building upon this issue, now that the input parameters of the iir_lowpass filter are clearer.

I am getting a mismatch between matlab and kfr-cpp at the zeros, but not at the poles and the gain. Here it is:

Here is my kfr-cpp code:

// Filter design:
using namespace kfr;
int order = 8;
double Fcutoff_Hz = 400;
double Fsampling_Hz = 2500;
filt = iir_lowpass(butterworth<double>(order),Fcutoff_Hz, Fsampling_Hz);
// Printing parameters
univector<std::complex<double>> z = filt.z;
univector<std::complex<double>> p = filt.p;

std::cout << "filt.z =  \n";
for (int i = 0; i < z.size(); i++){
    std::cout << z[i].real() << " + i " << z[i].imag() << "\n";
}
std::cout << "filt.p =  \n";
for (int i = 0; i < p.size(); i++){
    std::cout << p[i].real() << " + i " << p[i].imag() << "\n";
}
std::cout << "filt.k = " << filt.k << "\n";

Output:

filt.z =  
-1 + i 0
-1 + i 0
-1 + i 0
-1 + i 0
-1 + i 0
-1 + i 0
-1 + i 0
-1 + i 0
filt.p =  
0.460048 + i -0.71099
0.364735 + i -0.477871
0.314816 + i -0.275602
0.293105 + i -0.0901044
0.293105 + i 0.0901044
0.314816 + i 0.275602
0.364735 + i 0.477871
0.460048 + i 0.71099
filt.k = 0.000544958 

And this is the matlab counterpart:

order = 8;
Fcutoff_Hz = 400;
Fsampling_Hz = 2500;
Fnyquist_Hz = Fsampling_Hz/2;
Wn = Fcutoff_Hz/Fnyquist_Hz;
[b, a] = butter(order, Wn, 'low')
[z,p,k] = tf2zpk(b,a)

Output:

z =

  -1.0217 + 0.0000i
  -1.0152 + 0.0154i
  -1.0152 - 0.0154i
  -0.9998 + 0.0215i
  -0.9998 - 0.0215i
  -0.9848 + 0.0151i
  -0.9848 - 0.0151i
  -0.9786 + 0.0000i

p =

   0.4600 + 0.7110i
   0.4600 - 0.7110i
   0.3647 + 0.4779i
   0.3647 - 0.4779i
   0.3148 + 0.2756i
   0.3148 - 0.2756i
   0.2931 + 0.0901i
   0.2931 - 0.0901i

k =

   5.4496e-04

Note that the order is different but the values are the same. Is this a bug, or am I missing something...?

Thanks @BinDong314 for opening this issue. It helped me getting ob track for my own purposes with kfr.

Jalmenara commented 3 months ago

I add some info. Interestingly enough, the matlab result is not the same when requesting directly the zpk representation. With the same input parameters:

[z,p,k] = butter(order, Wn, 'low')

Output (which is, indeed, the same as in kfr-cpp):

z =

    -1
    -1
    -1
    -1
    -1
    -1
    -1
    -1

p =

   0.4600 + 0.7110i
   0.4600 - 0.7110i
   0.3647 + 0.4779i
   0.3647 - 0.4779i
   0.2931 + 0.0901i
   0.2931 - 0.0901i
   0.3148 + 0.2756i
   0.3148 - 0.2756i

k =

   5.4496e-04