I Wrote a bessel filter function

--Boundary-00=_M4MrGX3Zvma0O70
Content-Type: text/plain;
  charset="us-ascii"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline

Dear all,

I used octave to write some code to compute the coefficients of the 
transferfunction of a bessel filter, because I couldn't find these 
coefficients in the web. The reverse bessel polynomial is splitted in 
quadratic terms using the roots() function, I also normalize the function so 
that it is useful for filter design (3dB cutoff at frequency=1*i). Quadratic 
terms are good for designing signal filters with operational amplifiers in 
electrical engineering.

As I'm new to octave I don't know if I did it in the optimal way. On 
http://www.hw-loesungen.de/besselfilter.pdf you can find the mathematical 
background I used to write that code. Here is the code, you can include it in 
Octave.

With this code you can generate the coefficients for bessel filters with as 
many poles as you like, only restricted by numerical accuracy.

Further I would like to thank all people who contributed to octave, it is 
great!

Best,
Wolfgang
-- 
Wolfgang A. Heidrich

--Boundary-00=_M4MrGX3Zvma0O70
Content-Type: text/x-objcsrc;
  charset="us-ascii";
  name="bessel.m"
Content-Transfer-Encoding: 7bit
Content-Disposition: attachment;
	filename="bessel.m"

## Copyright (C) 2007 Wolfgang A. Heidrich
##
## This is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 2, or (at your option)
## any later version.
##
## This is distributed in the hope that it will be useful, but
## WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
## General Public License for more details.
##

## -*- texinfo -*-
## @deftypefn {Function File} {} bessel (@var{n})
##
## For a scalar @var{n}, return
## the coefficients of the quadratic terms of the
## factorized reverse bessel polynomial B(s), normalized so that B(i)=sqrt(2).
## @iftex
## @tex
## $$
## B(s)= (1+a_{1}x+b_{1}x^{2})(1+a_{2}x+b_{2}x^{2})(1+a_{3}x+b_{3}x^{2})\ldots
## $$
## @end tex
## @end iftex
## @ifinfo
##
## @example
## B(s)= (1 + a(1) * x + b(1)*x^2) * (1 + a(2) * x + b(2)*x^2) * ...
## @end example
## @end ifinfo
## @end deftypefn

## Author: Wolfgang Arne Heidrich
## Created: 29 July 2007
## Read http://www.hw-loesungen.de/besselfilter.pdf
## to understand how this code works.


function [a,b] = bessel (arg)

	if ( (nargout!=2) || (nargin != 1) || (size(arg)(1)!=1)
		       	|| (size(arg)(2)!=1)
	  )
		usage ("[a,b]=bessel(n) where n is an integer number");
	endif
	n=arg(1);

	printf("Bessel polynomial number of poles: %d\n", n );
%	printf(" %2d */ ", n );
	% First we generate the bessel polynomial using the recursive formula
	% v(n-ind+1) = 2*(n-ind+1)*cf/ind/(2*n-ind+1);
	cz=1;cn=1; % to be calculated coefficients as integernumbers
	cf=1; % to be calculated coefficients as floating numbers
	puts("c0=1\nc1=1\n");
	v(n)=1;
	v(n+1)=1;
	for ind=2:n
		cz = 2*(n-ind+1)*cz;
		cn = ind*(2*n-ind+1)*cn;
		cf = 2*(n-ind+1)*cf/ind/(2*n-ind+1);
		% just for fun we also calculate the fraction with integer
		% numerator cz and demoninator cn

		if( cz<cn ) min=cz; else min=cn; end;
		w=sqrt(min)+1;
		for k=2:w
			while( (rem(cz,k)==0) && (rem(cn,k)==0) )
				cz=cz/k;
				cn=cn/k;
			end
		end
		% Here I print out the integer numerator cz and demoninator cn,
		% that aren't used further.
		printf("c%i=%i/%i=%g\n", ind, cz, cn,cf );
		%v(n-ind+1)=cz/cn;
		v(n-ind+1)=cf;
	end
	puts("\n");
	r=roots(v); % these are the roots of the polynomial
	% Now we want to find the 3dB-Frequency omega of this filter.
	% i must be here the imaginary unit sqrt(-1) !
	omega=0;
	while( abs(polyval( v,i*omega ) ) < 2 )
		omega++;
	end
	w2 = sqrt(2);
	xa=0;xb=omega;
	diff = 1;
	count=0;
	while( abs(diff) > 1e-10 )
		omega = (xa+xb)/2;
		diff = w2 - abs(polyval( v,i*omega ) );
		if( diff > 0 )
			xa = (xa+xb)/2;
		else
			xb = (xa+xb)/2;
		end
		count++;
	end
%	printf("3dB-cutoff frequency:%g count:%i erg:%g\n",
	%	omega, count, abs(polyval( v,i*omega )) );
	%
	% OK, now we got the n roots of the bessel polynomial. Most of the
	% roots are complex numbers, but for even orders we always have one
	% root with imaginary part = zero. Now we want to determine the quadratic
	% terms: We know: (x-r)(x-~r) = x^2 - 2*x*real(r) + |r|^2
	% And we want this to divide
	% by |r|^2 : (x-r)(x-~r)/(|r|^2) = x^2/(|r|^2) - 2*x*real(r)/(|r|^2) + 1
	% We only take the poles with imgaginary part > 0
	% I assume for conjugated roots in r=roots(v) :
	% 	first comes the root with imaginary part > 0, then 
	%	comes the root with imaginary part < 0.
	% For other roots we know that the imaginary part is always zero.
	% I want to check this, to be sure that we extract only the wanted roots,
	% because I don't know if the accuracy of octave or the computer is good
	% enough, so that for example for a noncomplex root the calculated
	% imaginary part really is exact zero.
	cnt=1;
	for ind=1:n
		if( imag(r(ind)) > 0 )
			if( imag(r(ind+1))!=-imag(r(ind)) )
				error("Expected exact conjugated pole");
			end
			% it is fantastic that octave works so accurately,
			% that imag(r(ind+1)) == -imag(r(ind)) exactly !!!
			a(cnt)=-omega*2.0*real(r(ind))/(abs(r(ind))^2);
			b(cnt)=omega^2*1.0/(abs(r(ind))^2);
			%printf("a=%.10f b=%.10f\n", a(cnt), b(cnt) );
			%printf("%.10f,", a(cnt) );
			%printf("%.10f,", b(cnt) );
			cnt++;
		else
			if( imag(r(ind)) == 0 )
				% (x-r)/(-r) = -x/r + 1
				a(cnt) = -omega*1.0/r(ind);
				b(cnt)=0;
				%printf("a=%.10f b=%.10f\n", a(cnt), 0 );
				%printf("%.10f,", a(cnt) );
				%printf("%.10f,", 0 );
				cnt++;
			end
		end
	end
	puts("\n");
end

--Boundary-00=_M4MrGX3Zvma0O70
Content-Type: text/plain; charset="us-ascii"
MIME-Version: 1.0
Content-Transfer-Encoding: 7bit
Content-Disposition: inline

_______________________________________________
Octave-sources mailing list
[email protected]
https://www.cae.wisc.edu/mailman/listinfo/octave-sources

--Boundary-00=_M4MrGX3Zvma0O70--