Thread Subject: How to get [b,a] from fdesign

Hi,
  
   Consider,

Hf = fdesign.lowpass('Fp,Fst,Ap,Ast', Fp, Fst, Ap, Ast, Fs);
He = design(Hf,'ellip');

   How can I get the [b,a] from He?


Thanks,
cfy30

Hi cfy30,

 > How can I get the [b,a] from He?
[b,a] = tf(He);

Notice that it is good practice to always keep an IIR design in
second-order sections and avoid computing the transfer function explicitly.

The reason is that forming the transfer function involves polynomial
multiplication, or equivalently convolution of the polynomial
coefficients. Convolution can introduce significant round-off error even
with double precision floating-point arithmetic. This is mainly due to
additions of large numbers with small numbers (that become negligible
when using finite precision arithmetic).

Of course the effect gets worse as the filter order increases since we
are performing more and more convolutions. For low-order designs, it is
not absolutely necessary to use second-order sections, but there still
may be good reasons to do so from an implementation perspective.

As a rather extreme example consider the following design where forming
the transfer function completely distorts the response of the filter and
even makes the filter unstable:

Hf = fdesign.lowpass('Fp,Fst,Ap,Ast',.47,.48,.05,120);
Hc = design(Hf,'cheby1');
[b,a] = tf(Hc);
hfvt = fvtool(Hc,b,a);
legend(hfvt,'Second-order sections','Transfer function')

hth,

-Vincent

---
CFY30 CFY30 wrote:
> Hi,
>
> Consider,
>
> Hf = fdesign.lowpass('Fp,Fst,Ap,Ast', Fp, Fst, Ap, Ast, Fs);
> He = design(Hf,'ellip');
>
> How can I get the [b,a] from He?
>
>
> Thanks,
> cfy30

Thanks Vincent! You must be reading my mind, the problem I
am now running into is, consider the code attached. When I
plot filter1(b1, a1) and filter(b2, a2) separately, they
are fine. But it messed up when I plot the convolution
product. Filter1(b1, a1) is generated by iirlpnorm, I am
trying to find the best fit of a filter measured on bench.
Filter2(b2, a2) is a 5th butterworth. Could you advise how
to make filter1 and filter2 in terms of second order
sections and avoid dealing with their transfer function
explicity?

Thanks,
cfy30

a1 = [ 1
         -6.41419043322137
          17.6015355703759
         -26.7828346053453
          24.4006169931508
         -13.3073107931446
          4.02157413370539
        -0.519390864551612];
        
b1 = [0.000398558959342733
       -0.00141046292089994
        0.00185609288512552
       -0.00107470934235147
       0.00023052105125844];
   
[h1, f1] = freqz(b1, a1, 1024);
figure;
plot(f1, 20*log10(h1), 'r');
hold on;

a2 = [ 1
          -4.6872228997567
          8.79724893069547
         -8.26385275520034
          3.88511080048537
        -0.731276830775834];
               
b2 = [ 2.26420248883308e-007
       1.13210124441654e-006
       2.26420248883308e-006
       2.26420248883308e-006
       1.13210124441654e-006
       2.26420248883308e-007];
             
[h2, f2] = freqz(b2, a2, 1024);
plot(f2, 20*log10(h2), 'b');
grid on;
set(gca ,'xscale', 'log');

a3 = conv(a1, a2);
b3 = conv(b1, b2);

figure;
[h3, f3] = freqz(b3, a3, 1024);
plot(f3, 20*log10(h3), 'b');
grid on;
set(gca ,'xscale', 'log');


Vincent Pellissier <vincentp@mathworks.com> wrote in
message <g91si6$323$1@fred.mathworks.com>...
> Hi cfy30,
>
> > How can I get the [b,a] from He?
> [b,a] = tf(He);
>
> Notice that it is good practice to always keep an IIR
design in
> second-order sections and avoid computing the transfer
function explicitly.
>
> The reason is that forming the transfer function involves
polynomial
> multiplication, or equivalently convolution of the
polynomial
> coefficients. Convolution can introduce significant round-
off error even
> with double precision floating-point arithmetic. This is
mainly due to
> additions of large numbers with small numbers (that
become negligible
> when using finite precision arithmetic).
>
> Of course the effect gets worse as the filter order
increases since we
> are performing more and more convolutions. For low-order
designs, it is
> not absolutely necessary to use second-order sections,
but there still
> may be good reasons to do so from an implementation
perspective.
>
> As a rather extreme example consider the following design
where forming
> the transfer function completely distorts the response of
the filter and
> even makes the filter unstable:
>
> Hf = fdesign.lowpass('Fp,Fst,Ap,Ast',.47,.48,.05,120);
> Hc = design(Hf,'cheby1');
> [b,a] = tf(Hc);
> hfvt = fvtool(Hc,b,a);
> legend(hfvt,'Second-order sections','Transfer function')
>
> hth,
>
> -Vincent
>
> ---
> CFY30 CFY30 wrote:
> > Hi,
> >
> > Consider,
> >
> > Hf = fdesign.lowpass('Fp,Fst,Ap,Ast', Fp, Fst, Ap, Ast,
Fs);
> > He = design(Hf,'ellip');
> >
> > How can I get the [b,a] from He?
> >
> >
> > Thanks,
> > cfy30

Hi Vincent,

   I figured it out already. My solution is to use tf2sos
and df2sos to create the objects. Thanks for your help.


cfy30

"CFY30 CFY30" <cfy30@yahoo.com> wrote in message <g924v2$i70
$1@fred.mathworks.com>...
> Thanks Vincent! You must be reading my mind, the problem
I
> am now running into is, consider the code attached. When
I
> plot filter1(b1, a1) and filter(b2, a2) separately, they
> are fine. But it messed up when I plot the convolution
> product. Filter1(b1, a1) is generated by iirlpnorm, I am
> trying to find the best fit of a filter measured on
bench.
> Filter2(b2, a2) is a 5th butterworth. Could you advise
how
> to make filter1 and filter2 in terms of second order
> sections and avoid dealing with their transfer function
> explicity?
>
> Thanks,
> cfy30
>
> a1 = [ 1
> -6.41419043322137
> 17.6015355703759
> -26.7828346053453
> 24.4006169931508
> -13.3073107931446
> 4.02157413370539
> -0.519390864551612];
>
> b1 = [0.000398558959342733
> -0.00141046292089994
> 0.00185609288512552
> -0.00107470934235147
> 0.00023052105125844];
>
> [h1, f1] = freqz(b1, a1, 1024);
> figure;
> plot(f1, 20*log10(h1), 'r');
> hold on;
>
> a2 = [ 1
> -4.6872228997567
> 8.79724893069547
> -8.26385275520034
> 3.88511080048537
> -0.731276830775834];
>
> b2 = [ 2.26420248883308e-007
> 1.13210124441654e-006
> 2.26420248883308e-006
> 2.26420248883308e-006
> 1.13210124441654e-006
> 2.26420248883308e-007];
>
> [h2, f2] = freqz(b2, a2, 1024);
> plot(f2, 20*log10(h2), 'b');
> grid on;
> set(gca ,'xscale', 'log');
>
> a3 = conv(a1, a2);
> b3 = conv(b1, b2);
>
> figure;
> [h3, f3] = freqz(b3, a3, 1024);
> plot(f3, 20*log10(h3), 'b');
> grid on;
> set(gca ,'xscale', 'log');
>
>
> Vincent Pellissier <vincentp@mathworks.com> wrote in
> message <g91si6$323$1@fred.mathworks.com>...
> > Hi cfy30,
> >
> > > How can I get the [b,a] from He?
> > [b,a] = tf(He);
> >
> > Notice that it is good practice to always keep an IIR
> design in
> > second-order sections and avoid computing the transfer
> function explicitly.
> >
> > The reason is that forming the transfer function
involves
> polynomial
> > multiplication, or equivalently convolution of the
> polynomial
> > coefficients. Convolution can introduce significant
round-
> off error even
> > with double precision floating-point arithmetic. This
is
> mainly due to
> > additions of large numbers with small numbers (that
> become negligible
> > when using finite precision arithmetic).
> >
> > Of course the effect gets worse as the filter order
> increases since we
> > are performing more and more convolutions. For low-
order
> designs, it is
> > not absolutely necessary to use second-order sections,
> but there still
> > may be good reasons to do so from an implementation
> perspective.
> >
> > As a rather extreme example consider the following
design
> where forming
> > the transfer function completely distorts the response
of
> the filter and
> > even makes the filter unstable:
> >
> > Hf = fdesign.lowpass('Fp,Fst,Ap,Ast',.47,.48,.05,120);
> > Hc = design(Hf,'cheby1');
> > [b,a] = tf(Hc);
> > hfvt = fvtool(Hc,b,a);
> > legend(hfvt,'Second-order sections','Transfer function')
> >
> > hth,
> >
> > -Vincent
> >
> > ---
> > CFY30 CFY30 wrote:
> > > Hi,
> > >
> > > Consider,
> > >
> > > Hf = fdesign.lowpass('Fp,Fst,Ap,Ast', Fp, Fst, Ap,
Ast,
> Fs);
> > > He = design(Hf,'ellip');
> > >
> > > How can I get the [b,a] from He?
> > >
> > >
> > > Thanks,
> > > cfy30
>