Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

transfer function coefficients and v.s. Matlab #148

Open
BinDong314 opened this issue Feb 17, 2022 · 3 comments
Open

transfer function coefficients and v.s. Matlab #148

BinDong314 opened this issue Feb 17, 2022 · 3 comments

Comments

@BinDong314
Copy link

BinDong314 commented Feb 17, 2022

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?

  1. 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
Copy link

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
Copy link
Contributor

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
Copy link
Contributor

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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
3 participants