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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%