% m-file file name: dbubsf.m
% Function for designing digital Butterworth bandstop filter.
% Specifications are given in terms of amax, amin, f1bsf, f2bsf, f3bsf, f4bsf.
% Written by Dr. James S. Kang, Professor, 
% Department of Electrical and Computer Engineering,
% California State Polytechnic University, Pomona, 
% amax = maximum loss in dB in the pass band.
% amin = minimum loss in dB in the stop band.
% f1bsf = lower passband cutoff frequency in Hz. f1bsf < f3bsf 
% f2bsf = upper passband cutoff frequency in Hz. f4bsf < f2bsf
% f3bsf = lower stopband cutoff frequency in Hz. f1bsf < f3bsf < f4bsf
% f4bsf = upper stopband cutoff frequency in Hz. f3bsf < f4bsf < f2bsf 
% f1bsf < f3bsf < f4bsf < f2bsf
% fs = sampling rate in Hz.
% dp = location of digital poles.
% dz = location of digital zeros.
% H = Discrete-Time Fourier Transform (DTFT)
% h = impulse response.
% References:
% M. E. Van Valkenburg, Analog Filter Design, Holt, Rinehart and Winston, Inc., 1982.
% Richard W. Daniels, Approximation Methods for Electronic Filter Design, McGraw-Hill, 1974.
% Andreas Antoniou, Digital Filters: Analysis, Design, and Applications, 2nd ed., McGraw-Hill, 1993.
% Leland B. Jackson, Digital Filters and Signal Processing, 3rd ed., Kluwer Academic Publishers, 1996.
% T. W. Parks and C. S. Burrus, Digital Filter Design, John Wiley and Sons, Inc., 1987.
% J. G. Proakis and D. G. Manolakis, Digital Signal Processing, 3rd ed., Prentice-Hall, 1996.
% A. V. Oppenheim and R. W. Shafer, Discrete-Time Signal Processing, Prentice-Hall, 1989.
function[dp,dz,h] = dbubsf(amax,amin,f1bsf,f2bsf,f3bsf,f4bsf,fs)
Ts = 1/fs;     % Ts = sampling interval.
q1bsf = 2*pi*f1bsf*Ts;     % q1bsf = digital lower passband cutoff frequency in radians.
q2bsf = 2*pi*f2bsf*Ts;     % q2bsf = digital upper passband cutoff frequency in radians.
q3bsf = 2*pi*f3bsf*Ts;     % q3bsf = digital lower stopband cutoff frequency in radians.
q4bsf = 2*pi*f4bsf*Ts;     % q4bsf = digital upper stopband cutoff frequency in radians.
w1bsf = 2*fs*tan(q1bsf/2);   % w1bsf = prewarped lower passband cutoff frequency in radians/s.
w2bsf = 2*fs*tan(q2bsf/2);   % w2bsf = prewarped upper passband cutoff frequency in radians/s.
w3bsf = 2*fs*tan(q3bsf/2);   % w3bsf = prewarped lower stopband cutoff frequency in radians/s.
w4bsf = 2*fs*tan(q4bsf/2);   % w4bsf = prewarped upper stopband cutoff frequency in radians/s.
if w1bsf*w2bsf > w3bsf*w4bsf w4bsf = w1bsf*w2bsf/w3bsf; end
if w1bsf*w2bsf < w3bsf*w4bsf w3bsf = w1bsf*w2bsf/w4bsf; end
N = ceil(log((10^(0.1*amin)-1)/(10^(0.1*amax)-1))/(2*log((w2bsf-w1bsf)/(w4bsf-w3bsf)))); % N = order of the filter.
disp('Digital Butterworth Bandstop Filter Design')
disp(' ')
disp('Lower passband cutoff frequency q1bsf =')
disp(q1bsf)
disp('Upper passband cutoff frequency q2bsf =')
disp(q2bsf)
disp('Lower stopband cutoff frequency q3bsf =')
disp(q3bsf)
disp('Upper stopband cutoff frequency q4bsf =')
disp(q4bsf)
disp('Prewarped lower passband cutoff frequency w1bsf =')
disp(w1bsf)
disp('Prewarped upper passband cutoff frequency w2bsf =')
disp(w2bsf)
disp('Prewarped lower stopband cutoff frequency w3bsf =')
disp(w3bsf)
disp('Prewarped upper stopband cutoff frequency w4bsf =')
disp(w4bsf)
disp('The order of the filter is')
disp(N)
%
% Calculating the normalized lowpass pole locations.  
% Vector s is the normalized poles.
%
if rem(N,2) == 1 s(1) = -1; else s(1) = 0; end
if rem(N,2) == 1 min = (N+1)/2; else min = N/2 + 1; end
if N == 1 max = 1; elseif rem(N,2) == 1 max = N-1; else max = N; end
if rem(N,2) == 1 factor = 0; else factor = 1/2; end
if N == 1 maxa = 1; elseif rem(N,2) == 1 maxa = (N-1)/2; else maxa = N/2; end
if rem(N,2) == 1 evena = 0; else evena = 1; end
for m = min:max
   s(m) = cos(pi*(m - factor)/N) + j*sin(pi*(m - factor)/N);
end
for ka = 1:maxa
   s(2*ka-evena) = s((N+1+evena+2*(ka-1))/2);
   s(2*ka+1-evena) = conj(s(2*ka-evena));
end
%
% Plotting the normalized analog lowpass poles.
%
for n = 1:201
   x(n) = (n-101)/100;
   ty(n) = sqrt(1-x(n).*x(n));
   by(n) = - ty(n);
end
figure
plot(x,ty)
hold on
plot(x,by)
plot(real(s),imag(s),'bx'); grid;
title('Normalized Analog Lowpass Pole Locations')
xlabel('Real Part')
ylabel('Imaginary Part')
disp('Normalized Analog Lowpass Pole Locations')
disp(s);
%
% Calculating the normalized analog lowpass zero locations.  
% Vector szb is the normalized analog zeros.
%
for k = 1:N
   sz(k) = 10^50;
end
%
% Calculating the frequency transformed analog poles and zeros.
% Lowpass to bandstop transformation.
%
Wo = 1/(10^(0.1*amax)-1)^(1/(2*N));  % Wo = half-power frequency.
w21=(w2bsf-w1bsf)/2;
w2m1=w2bsf-w1bsf;
w12t4=4*w1bsf*w2bsf;
for k = 1:N
   s(k) = Wo*s(k);
   R1(k) = real(s(k));
   I1(k) = imag(s(k));
   R2(k) = R1(k)^2;
   I2(k) = I1(k)^2;
   mag(k) = R2(k) + I2(k);
   A(k) = (I2(k)-R2(k))*w2m1^2/mag(k)^2 + w12t4;
   B(k) = 2*abs(R1(k))*abs(I1(k))*w2m1^2/mag(k)^2;
   RA(k) = (1/2)*sqrt((sqrt(A(k)^2+B(k)^2)-A(k))/2);
   RB(k) = j*(1/2)*(sqrt((sqrt(A(k)^2+B(k)^2)+A(k))/2)-I1(k)*w2m1/mag(k));
end
for m = 1:2*N
   if rem(m,2) == 1
      if I1((m+1)/2) > 0 sfs(m) = w21*R1((m+1)/2)/mag((m+1)/2) + RA((m+1)/2) + RB((m+1)/2); end
      if I1((m+1)/2) < 0 sfs(m) = w21*R1((m+1)/2)/mag((m+1)/2) - RA((m+1)/2) + RB((m+1)/2); end    
      R1a = (w2m1/R1((m+1)/2))^2;
      if I1((m+1)/2) == 0 & R1a >= w12t4 sfs(m) = w21/R1((m+1)/2) + (1/2)*sqrt(R1a-w12t4); end  
      if I1((m+1)/2) == 0 & R1a < w12t4 sfs(m) = w21/R1((m+1)/2) + j*(1/2)*sqrt(-R1a+w12t4); end     
   else 
      sfs(m) = conj(sfs(m-1)); 
      if I1(m/2) == 0 & (w2m1/R1(m/2))^2 >= w12t4     
         sfs(m) = w21/R1(m/2) - (1/2)*sqrt(w2m1^2/R1(m/2)^2-w12t4); end   
   end
end
if N==1 sfs(1) = w21/R1(1) + (1/2)*sqrt((w2m1/R1(1))^2-w12t4); end
if N==1 sfs(2) = w21/R1(1) - (1/2)*sqrt((w2m1/R1(1))^2-w12t4); end
for m = 1:2*N
   if rem(m,2) == 1 szfs(m) = j*sqrt(w1bsf*w2bsf); else szfs(m) = conj(szfs(m-1)); end
end
%
% Plotting the frequency transformed bandpass poles.
%
wo = sqrt(w1bsf*w2bsf)/Wo;
for na = 1:201
   xa(na) = wo*(na-101)/100;
   tya(na) = sqrt(wo^2-xa(na)^2);
   bya(na) = - tya(na);
end
figure
plot(xa,tya)
hold on
plot(xa,bya)
plot(real(sfs),imag(sfs),'bx'); grid;
title('Frequency Transformed Bandstop Pole Locations')
xlabel('Real Part')
ylabel('Imaginary Part')
disp('Frequency Transformed Bandstop Pole Locations')
disp(sfs)
disp('Frequency Transformed Bandstop Zero Locations')
disp(szfs)
%
% Transforming analog poles and zeros to digital poles and zeros using bilinear transformation.
%
for m = 1:2*N
   dp(m) = (1 + (Ts/2)*sfs(m))/(1 - (Ts/2)*sfs(m));
   dz(m) = (1 + (Ts/2)*szfs(m))/(1 - (Ts/2)*szfs(m));
end
%
% Plotting digital poles and zeros on the complex z-plane
%
figure
plot(x,ty)
hold on
plot(x,by)
plot(real(dp),imag(dp),'bx')
plot(real(dz),imag(dz),'ro'); grid;
title('Digital Poles and Zeros')
xlabel('Real Part')
ylabel('Imaginary Part')
disp('Digital Pole Locations')
disp(dp)
disp('Digital Zero Locations')
disp(dz)
%
for i = 1:N
	d2(i,1)=1; d2(i,2)=-2*real(dp(2*i-1));d2(i,3) = (abs(dp(2*i-1)))^2;
	n2(i,1)=1; n2(i,2)=-2*real(dz(2*i-1));n2(i,3) = (abs(dz(2*i-1)))^2;
end
disp('Numerator Second Order Terms')
disp(n2)
disp('Denominator Second Order Terms')
disp(d2)
%
% Plotting the magnitude response and the phase response with respect to the normalized digital frequency.
%
M = 128;  % M is the number of points the frequency response is calculated.
step = 2*pi/M;
for m = 1:M
   H(m) = (1-dz(1)*exp(j*(-1)*(m-1)*step))*(1-dz(2)*exp(j*(-1)*(m-1)*step))/[(1-dp(1)*exp(j*(-1)*(m-1)*step))*(1-dp(2)*exp(j*(-1)*(m-1)*step))];
end
if N > 1
   for n = 3:2*N
     for m = 1:M
       H(m) = H(m).*[(1-dz(n)*exp(j*(-1)*(m-1)*step))/(1-dp(n)*exp(j*(-1)*(m-1)*step))];
     end
   end
end
max = abs(H(1));
for m = 2:M
   if abs(H(m)) > max max = abs(H(m)); end
end
disp('maximum gain is')
disp(max)
b0=1/max;
disp('Constant Term b0 =')
disp(b0);
%
%Polynomial Form
%
nump(1) = n2(1,1); nump(2) = n2(1,2); nump(3) = n2(1,3); 
denp(1) = d2(1,1); denp(2) = d2(1,2); denp(3) = d2(1,3);
if N > 1
	for i=2:N nump = conv(nump,n2(i,:)); denp = conv(denp,d2(i,:)); end
end
disp('Numerator Polynomial')
disp(b0*nump)
disp('Denominator Polynomial')
disp(denp)
%
for m = 1:M
   f(m) = (m-1)/M; 
   H(m) = H(m)/max;
   mag(m) = abs(H(m));
   if H(m) == 0 dB(m) = -200; else dB(m) = 20*log10(mag(m)); end
   if dB(m) < -200 dB(m) = -200; end
   pha(m) = angle(H(m));
end
figure
plot(f,mag); grid;
title('Magnitude Response')
xlabel('Frequency')
ylabel('Magnitude')
figure
plot(f,dB); grid;
title('Magnitude Response in dB')
xlabel('Frequency')
ylabel('Magnitude in dB')
figure
plot(f,pha); grid;
title('Phase Response')
xlabel('Frequency')
ylabel('Phase in Radians')
%
% Calculating and Plotting the Impulse Response of the Filter.
%
h = ifft(H);
for m = 1:M
   k(m) = m-1; z(m) = 0; Lo(m) = 0; Hi(m) = 0;
   if h(m) >= 0 Hi(m) = h(m); else Lo(m) = h(m); end
end
maxh = h(1); minh = h(1);
for m = 2:M
   if h(m) > maxh maxh = h(m); end
   if h(m) < minh minh = h(m); end
end
figure
%errorbar(k,z,Lo,Hi,'w')
stem(k,h); grid;
title('Impulse Response')
xlabel('Time')
ylabel('Amplitude')