Home > chronux_2_00 > spikesort > utility > datatools > buttersimple.m

buttersimple

PURPOSE ^

BUTTERSIMPLE Fixed order Butterworth digital filter design.

SYNOPSIS ^

function [B, A] = buttersimple(Wn, option)

DESCRIPTION ^

BUTTERSIMPLE      Fixed order Butterworth digital filter design.
   [B, A] = BUTTERSIMPLE(Wn), where 0.0 < Wn < 1.0, computes the
   numerator/denominator coefficients for a 6th order lowpass digital
   Butterworth filter.

   [B, A] = BUTTERSIMPLE(Wn, 'high'), where 0.0 < Wn < 1.0, computes the
   coefficients for a 6th order high-pass digital Butterworth filter.

   [B, A] = BUTTERSIMPLE([W1 W2]), with 0.0 < W1 <1.0 and 0.0 < W2 < 1.0,
   computes the coefficients for a 6th order bandpass digital filter with
   passband W1 < W < W2.

   Based on TMW Signal Processing Toolbox function BUTTER, but without
   any of the generality and does not require any toolboxes.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [B, A] = buttersimple(Wn, option)
0002 %BUTTERSIMPLE      Fixed order Butterworth digital filter design.
0003 %   [B, A] = BUTTERSIMPLE(Wn), where 0.0 < Wn < 1.0, computes the
0004 %   numerator/denominator coefficients for a 6th order lowpass digital
0005 %   Butterworth filter.
0006 %
0007 %   [B, A] = BUTTERSIMPLE(Wn, 'high'), where 0.0 < Wn < 1.0, computes the
0008 %   coefficients for a 6th order high-pass digital Butterworth filter.
0009 %
0010 %   [B, A] = BUTTERSIMPLE([W1 W2]), with 0.0 < W1 <1.0 and 0.0 < W2 < 1.0,
0011 %   computes the coefficients for a 6th order bandpass digital filter with
0012 %   passband W1 < W < W2.
0013 %
0014 %   Based on TMW Signal Processing Toolbox function BUTTER, but without
0015 %   any of the generality and does not require any toolboxes.
0016 
0017 % Step 1: Compute analog frequencies
0018 if (~all(Wn>0 & Wn<1.0)),  
0019     error('Cutoff frequencies must be in the interval (0,1).');  
0020 end
0021 Wn = 4*tan(pi*Wn/2);
0022 if (length(Wn) == 2)
0023     mode = 'B';
0024     Bw = diff(Wn);          % bandwidth
0025     Wn = sqrt(prod(Wn));    % center frequency
0026 elseif (length(Wn) == 1)
0027     if (nargin > 1)
0028         if (~strcmpi(option, 'high')),  error('Unknown option.');   end;
0029         mode = 'H';
0030     else  mode = 'L';
0031     end
0032 else
0033     error('Wn can not have more than two elements.');
0034 end
0035 
0036 % Step 2: Make Butterworth lowpass prototype
0037 if (mode == 'B')  % if bandpass, start w/ 3rd order prototype
0038     a = [-1 0 0; 1 -1 -1; 0 1 0];
0039     b = [1 0 0]';    
0040     c = [0 0 1];
0041     d = 0;
0042 else  % otherwise 6th order lowpass prototype
0043     a = diag([1 1 1 1 1], -1) + diag([-1 0 -1 0 -1],1);
0044     for j = [1,3,5]
0045         a(j,j) = 2*cos((1-(j/12))*pi);
0046     end
0047     b = [1 0 0 0 0 0]';
0048     c = [0 0 0 0 0 1];
0049     d = 0;
0050 end
0051         
0052 % Step 3: Convert to lowpass/bandpass/highpass & shift frequency
0053 switch(mode)
0054     case 'L',
0055         a = Wn*a;
0056         b = Wn*b;
0057 %       c = c;
0058 %       d = d;
0059     case 'B',
0060         a = Wn*[a/(Wn/Bw) eye(3); -eye(3) zeros(3)];
0061         b = [b*Bw; 0; 0; 0];
0062         c = [c 0 0 0];
0063         d = 0;
0064     case 'H',
0065         d = d - c/a*b;        
0066         c = c/a;
0067         b = -Wn*(a\b);       
0068         a =  Wn*inv(a);
0069 end
0070 
0071 % Step 4: Find discrete equivalent
0072 t1 = eye(6) + a/4;    t2 = eye(6) - a/4;
0073 d = c/t2*b/4 + d;
0074 c = c/t2/sqrt(2);
0075 b = t2\b/sqrt(2);
0076 a = t2\t1;
0077 
0078 % Step 5: Convert to polynomial form
0079 A = poly(a);
0080 B = poly(a-b*c)+(d-1)*A;

Generated on Fri 15-Aug-2008 11:35:42 by m2html © 2003