% Frequency response for Butterworth and Chebyshev 1 filters
%
clf
filt_type = input('Enter filter type; 1 = Butterworth; 2 = Chebyshev type 1 ');
n_max = input('Enter maximum order of filter ');
fc = input('Enter cutoff frequency (3-dB for Butterworth) in Hz ');
if filt_type == 2
   R = input('Enter Chebyshev filter ripple in dB ');
end
W = logspace(0, 3, 1000);						% Set up frequency axis; hertz assumed
for n = 1:n_max
   if filt_type == 1
      [num,den]=butter(n, 2*pi*fc, 's');	    % Generate numerator and denominator polynomials for Butterworth
   elseif filt_type == 2
      [num,den]=cheby1(n, R, 2*pi*fc, 's');      % Generate numerator and denominator polynomials for Chebyshev
   end
   H = freqs(num, den, W);						% Generate complex frequency response of chosen filter
   [phase, mag] = cart2pol(real(H),imag(H));	% Convert to polar coordinates
   phase2 = unwrap(phase);
   % Do plot of the amplitude response; frequency in Hz
   subplot(2,1,1),semilogx(W/(2*pi),20*log10(mag)),axis([min(W/(2*pi)) max(W/(2*pi)) -40 0]),...
   if n == 1									% Put on plot labels and title
   	  grid on
   	  ylabel('|H| in dB')
      hold on									% Hold plot to do all orders wanted;
      if filt_type == 1
         title(['Butterworth filter responses: order 1 - ',num2str(n_max),'; cutoff freq = ',num2str(fc),' Hz'])
      elseif filt_type == 2
         title(['Chebyshev filter responses: order 1 - ',num2str(n_max),'; ripple = ',num2str(R),' dB; cutoff freq = ',num2str(fc),' Hz'])
      end
   end
   % Do plot of phase response
   subplot(2,1,2),semilogx(W/(2*pi), 180*phase2/pi),axis([min(W/(2*pi)) max(W/(2*pi)) 180*min(phase2)/pi 0]),...
   if n == 1
      grid on
      hold on
      xlabel('f, Hz'),ylabel('Phase in degrees')
   end
end