implement biquad filter in octave or matlab

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.

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