I Wrote a bessel filter function

[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