


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.


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;