How to extract IIR coefficients w/Matlab?

Hello.

1. Is there an easy way in Matlab to take a z-domain transfer function
and extract coefficients for an IIR filter type implementation?

I first designed my filter in the s-domain and tweaked it to get my
desired response. I then used Matlab's c2dm with the Tustin method to
get a z-domain transfer function.

I am aware of "fdatool" but my filter started out life in s-domain and
it is more than a simple low-pass or high-pass filter.

2. Is there an associated facility in Matlab that can take the
coefficients for the IIR filter and provide a sample implementation in
C?

Currently, I am manually examining the transfer function, putting it
into an acceptable form and then copying down the coefficients and
writing the recursive implementation in C. However, whenever the order
of the filter increases or decreases, it's very annoying to have to
repeat the entire process.

Regards.
Mark.

Mark wrote:

> Hello.
> 
> 1. Is there an easy way in Matlab to take a z-domain transfer function
> and extract coefficients for an IIR filter type implementation?
> 
> I first designed my filter in the s-domain and tweaked it to get my
> desired response. I then used Matlab's c2dm with the Tustin method to
> get a z-domain transfer function.
> 
> I am aware of "fdatool" but my filter started out life in s-domain and
> it is more than a simple low-pass or high-pass filter.
> 
> 2. Is there an associated facility in Matlab that can take the
> coefficients for the IIR filter and provide a sample implementation in
> C?
> 
> Currently, I am manually examining the transfer function, putting it
> into an acceptable form and then copying down the coefficients and
> writing the recursive implementation in C. However, whenever the order
> of the filter increases or decreases, it's very annoying to have to
> repeat the entire process.
> 
> Regards.
> Mark.
> 

Seems like this would be a fairly easy script to write for a competent C 
programmer -- get the roots and poles of the transfer function, factor 
it into a bunch of 1st- and 2nd- order ones, then generate your code for 
each.

If you can survive a bit of overhead you could just make a generic 
2nd-order filter function with an associated definition structure (or 
use C++) and fill in the blanks on your filter definitions.

-- 

Tim Wescott
Wescott Design Services
http://www.wescottdesign.com

Mark wrote:

> Hello.
> 
> 1. Is there an easy way in Matlab to take a z-domain transfer function
> and extract coefficients for an IIR filter type implementation?
> 
> I first designed my filter in the s-domain and tweaked it to get my
> desired response. I then used Matlab's c2dm with the Tustin method to
> get a z-domain transfer function.
> 
> I am aware of "fdatool" but my filter started out life in s-domain and
> it is more than a simple low-pass or high-pass filter.
> 
> 2. Is there an associated facility in Matlab that can take the
> coefficients for the IIR filter and provide a sample implementation in
> C?
> 
> Currently, I am manually examining the transfer function, putting it
> into an acceptable form and then copying down the coefficients and
> writing the recursive implementation in C. However, whenever the order
> of the filter increases or decreases, it's very annoying to have to
> repeat the entire process.
> 
> Regards.
> Mark.
> 
Hi Mark,

If you say that you have designed your filter in the s-domain and need
IIR filter coeffients, why not using the bilinear function ? Or am I
missing the point. (Maybe bilinear is only with the sptoobox)

Mark <mark@abuse.hotmail> wrote in message news:<hcpnd01bdi00siuslctmc0qf161t16ndtd@4ax.com>...
> Hello.
> 
> 1. Is there an easy way in Matlab to take a z-domain transfer function
> and extract coefficients for an IIR filter type implementation?
> 
> I first designed my filter in the s-domain and tweaked it to get my
> desired response. I then used Matlab's c2dm with the Tustin method to
> get a z-domain transfer function.
> 

Maybe I'm missing something, but aren't the coefficients of your
transfer function the IIR coefficients by definition?

Tom

Tom Seim wrote:

   ...

> Maybe I'm missing something, but aren't the coefficients of your
> transfer function the IIR coefficients by definition?

That would be dandy! Can you show me how?

Jerry
-- 
Engineering is the art of making what you want from things you can get.
�����������������������������������������������������������������������

Jerry Avins wrote:
> Tom Seim wrote:
> 
>   ...
> 
>> Maybe I'm missing something, but aren't the coefficients of your
>> transfer function the IIR coefficients by definition?
> 
> 
> That would be dandy! Can you show me how?
> 
> Jerry

Simple:

            b_n*z^n + ... + b_0
T(z) = ---------------------------
         z^n + a_{n-1} + ... + a_0

x_k = -(a_{k-1}*x_{k-1} + ... + a_0*x{k-n}) + b_n*u_k + ... + b_0*u_{k-n}.

Of course the requirement for precision in both your accumulator and 
your coefficients goes up exponentially with n (i.e. your bit count goes 
up linearly), so if you go much beyond a second order your precision 
requirements are absurd.

It's much better to break the transfer function down into 1st and 
2nd-order ones and cascade, and even then you can have problems if you 
have any highly resonant pole pairs.

-- 

Tim Wescott
Wescott Design Services
http://www.wescottdesign.com

Suodatin Pussi wrote:
> Mark wrote:
> 
>> Hello.
>> 
>> 1. Is there an easy way in Matlab to take a z-domain transfer function
>> and extract coefficients for an IIR filter type implementation?
>> 
>> I first designed my filter in the s-domain and tweaked it to get my
>> desired response. I then used Matlab's c2dm with the Tustin method to
>> get a z-domain transfer function.
>> 
>> I am aware of "fdatool" but my filter started out life in s-domain and
>> it is more than a simple low-pass or high-pass filter.
>> 
>> 2. Is there an associated facility in Matlab that can take the
>> coefficients for the IIR filter and provide a sample implementation in
>> C?
>> 
>> Currently, I am manually examining the transfer function, putting it
>> into an acceptable form and then copying down the coefficients and
>> writing the recursive implementation in C. However, whenever the order
>> of the filter increases or decreases, it's very annoying to have to
>> repeat the entire process.
>> 
>> Regards.
>> Mark.
>> 
> Hi Mark,
> 
> If you say that you have designed your filter in the s-domain and need
> IIR filter coeffients, why not using the bilinear function ? Or am I
> missing the point. (Maybe bilinear is only with the sptoobox)
> 
Mark,
If you need to split up an nth order z polynome in into 2nd (+1nd) order
sections, here's some code I wrote some years ago in Matlab. The idea is
to find the roots of your polynome and combine complex conjungate pairs
to get real filter coeffs. It's old code and there are probably better
ways to do it. Anyway second order section shouldn't be to hard to program.

Suo

function        coef=fil2scnd(bx,ax,fs,startfreq,endfreq);
%function        coef=fil2scnd(bx,ax,fs,startfreq,endfreq);
%
%       program: FIL2SCND - split filter coefficients into
%               second order sections, like
%
%               [a2 b2 a1 b1 b0]
%               [a2 b2 a1 b1 b0]
%               [a2 b2 a1 b1 b0]
%               [a2 b2 a1 b1 b0]
%
%               uses b,a and fs as an input and
%               plots original filter and 2nd order cascade
%               startfreq and endfreq can be specified for plot
%
%               example: get 4th order octave filter at 8kHz
%
%               fs=44100;
%               fc=8000;
%               ripple=0.2;
%               order=4;
%               [b,a]=cheby1(order,ripple,[fc/sqrt(2)/(fs/2)
fc*sqrt(2)/(fs/2)]);
%               fil2scnd(b,a,fs);
%
%       	    300798 v1.0
%

if ~exist('startfreq')
        startfreq=10;
end
if ~exist('endfreq')
        endfreq=fs/2;
end

f=logspace(log10(startfreq),log10(endfreq),1000);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%
%        find roots of polynome to split in n 2nd order sections
%
ar0=roots(ax);

%
%       split into real and imaginary part
%       sort complex conjungates
%
ar0_real=ar0(imag(ar0)==0);
ar0_imag=sort(ar0(~imag(ar0)==0));

%
%       add zero to ar0_real for oddness
%
ar0_real=[ar0_real; 0];

%
% make second order sections of real values
%

for i=1:length(ar0_real)/2
        a=[a;poly([ar0_real(i*2-1) ar0_real(i*2)])];
end;


%
%       make second order sections of complex conjungates
%
for i=1:length(ar0_imag)/2
        a=[a;real(poly([ ar0_imag(i*2-1) ar0_imag(i*2) ]))];
end



%
%       make sure that the first coeff of numerator of each second order
%       section is equal to 1
%
GainA=1;
for i=1:size(a,1)
        GainA=GainA*a(i,1);
        a(i,:)=a(i,:)/a(i,1);
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%
%        find roots of polynome to split in n 2nd order sections
%
br0=roots(bx);
GainB=bx(1);

%
%       split into real and imaginary part
%       and sort the complex conjungates
%
br0_real=br0(imag(br0)==0);
br0_imag=sort(br0(~imag(br0)==0));


%
%       add zero to br0_real for oddness
%
br0_real=[br0_real; 0];


%
% make second order sections of real values
%

for i=1:length(br0_real)/2
        b=[b;poly([br0_real(i*2-1) br0_real(i*2)])];
end;


%
%       make second order sections of complex conjungates
%
for i=1:length(br0_imag)/2
        b=[b;real(poly([ br0_imag(i*2-1) br0_imag(i*2) ]))];
end


%
%       divide gain over all sections
%

Gain=exp(log(GainA*GainB)/size(b,1));

b=b*Gain;


%
%       check filter by convolving second order sections
%


bc=b(1,:);
for i=2:size(b,1)
        bc = conv(bc,b(i,:));
end


ac=a(1,:);
for i=2:size(a,1)
        ac = conv(ac,a(i,:));
end

%
%       plot both the original and the calculated filters
%


[h1,f1]=freqz(bx,ax,f,fs);
[h2,f2]=freqz(bc,ac,f,fs);

semilogx(f,20*log10(abs(h1)),'o',f,20*log10(abs(h2)))
axis([startfreq endfreq -90 10])
grid



%
%       write a and b coeffs for second order sections to screen
%

%b
%a

%
%       order=a2 b2 a1 b1 a0
%
for i=1:size(b,1)
        coef=[coef;a(i,3) b(i,3) a(i,2) b(i,2) b(i,1)];
end




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Tim Wescott wrote:

> Jerry Avins wrote:
> 
>> Tom Seim wrote:
>>
>>   ...
>>
>>> Maybe I'm missing something, but aren't the coefficients of your
>>> transfer function the IIR coefficients by definition?
>>
>>
>>
>> That would be dandy! Can you show me how?
>>
>> Jerry
> 
> 
> Simple:
> 
>            b_n*z^n + ... + b_0
> T(z) = ---------------------------
>         z^n + a_{n-1} + ... + a_0
> 
> x_k = -(a_{k-1}*x_{k-1} + ... + a_0*x{k-n}) + b_n*u_k + ... + b_0*u_{k-n}.
> 
> Of course the requirement for precision in both your accumulator and 
> your coefficients goes up exponentially with n (i.e. your bit count goes 
> up linearly), so if you go much beyond a second order your precision 
> requirements are absurd.
> 
> It's much better to break the transfer function down into 1st and 
> 2nd-order ones and cascade, and even then you can have problems if you 
> have any highly resonant pole pairs.

Sorry to bother you. Tim. I was thinking s-domain transfer function
coefficients to filter to biquad coefficients in one step. I know the
long way around.

Jerry
-- 
Engineering is the art of making what you want from things you can get.
�����������������������������������������������������������������������

Jerry Avins wrote:
> Tim Wescott wrote:
> 
>> Jerry Avins wrote:
>>
>>> Tom Seim wrote:
>>>
>>>   ...
>>>
>>>> Maybe I'm missing something, but aren't the coefficients of your
>>>> transfer function the IIR coefficients by definition?
>>>
>>>
>>>
>>>
>>> That would be dandy! Can you show me how?
>>>
>>> Jerry
>>
>>
>>
>> Simple:
>>
>>            b_n*z^n + ... + b_0
>> T(z) = ---------------------------
>>         z^n + a_{n-1} + ... + a_0
>>
>> x_k = -(a_{k-1}*x_{k-1} + ... + a_0*x{k-n}) + b_n*u_k + ... + 
>> b_0*u_{k-n}.
>>
>> Of course the requirement for precision in both your accumulator and 
>> your coefficients goes up exponentially with n (i.e. your bit count 
>> goes up linearly), so if you go much beyond a second order your 
>> precision requirements are absurd.
>>
>> It's much better to break the transfer function down into 1st and 
>> 2nd-order ones and cascade, and even then you can have problems if you 
>> have any highly resonant pole pairs.
> 
> 
> Sorry to bother you. Tim. I was thinking s-domain transfer function
> coefficients to filter to biquad coefficients in one step. I know the
> long way around.
> 
> Jerry

The OP is converting to the z domain, he just wants a C code generator 
to go from there to C.

I was wondering why you needed to ask!

-- 

Tim Wescott
Wescott Design Services
http://www.wescottdesign.com