Partial code is printed. This I had done two years back in
Octave. In case if you need full listing get in touch on
bharat at arithos dot com.
Regards
Bharat Pathak
//=============================================================
function [w1 w2 w3 yout smin smax] = biquad_3stage(sos, xin)
% This function implements 3 biquads in time domain.
% Basically by means of difference equation.
% This program can be used to 3 reasons.
% 1. perform fractional bits analysis.
% 2. find overflow at critical adder nodes.
% 3. fixing of limit cycle oscillations.
fb_sel1 = 1; % 0 means floor, 1 means round and 2 means fix
fb_sel2 = 1; % 0 means floor, 1 means round and 2 means fix
fb_sel3 = 2; % 0 means floor, 1 means round and 2 means fix
b0c = sos(3,1);
b1c = sos(3,2);
b2c = sos(3,3);
a1c = sos(3,5);
a2c = sos(3,6);
for i = 1 : length(xin)
%------------------ first stage --------------------------
sum1 = b0a*xin(i) + b1a*xp + b2a*xpp - a1a*w1p - a2a*w1pp;
if(sum1 < s1min)
s1min = sum1;
elseif(sum1 > s1max)
s1max = sum1;
else
end
%------------------ second stage --------------------------
sum2 = b0b*w1(i) + b1b*w1p + b2b*w1pp - a1b*w2p - a2b*w2pp;
if(sum2 < s2min)
s2min = sum2;
elseif(sum2 > s2max)
s2max = sum2;
else
end
%------------------ third stage --------------------------
if(fb_sel3 == 0)
w3(i) = floor(sum3 * k)/k;
elseif(fb_sel3 == 1)
w3(i) = round(sum3 * k)/k;
else
w3(i) = fix(sum3 * k)/k;
end
xpp = xp;
xp = xin(i);
w1pp = w1p;
w1p = w1(i);
yout(i) = round(sum3 * k1)/k1;
end
smin = [s1min s2min s3min];
smax = [s1max s2max s3max];
endfunction
Reply by katwalatapan●March 24, 20102010-03-24
Hello,
I'm very new to the implementation of the biquad filter. It would really
help if someone could post links where I can learn how to implement a
biquad filter as well as cascaded biquad filter in octave or matlab. I'd
appreciate if someone could point me as to how do I decide the filter
co-efficient and test them with respect to the desired frequency response.
Thank you.