function drawPeakResponse(f, separation, order, ftype, fplot)
%drawPeakResponse(f,seperation, order, ftype, fplot)
%takes in vector f and scalar seperation as in distort
%f is a vector of frequencies- the frequency bounds
%there will be length(f) - 1 filters and thus frequency bins
%seperation is a number less than 1.
%Each filter will be centered at the average of the 2 frequency bounds and
%its 3dB width will be the distance between the 2 bounds * seperation
%order is the order of the filters.  For passband must be a multiple of 2.
%ftype is the type of filter to use and can be:
%  'b' for butterworth
%  'c(R)' (no parentheses) for chebychev with ripple R
%fplot is a vector of frequencies for which we will compute the frequency
%responses and plot
%we draw the corresponding frequency response of all filters AS A PEAK
%(rather than power.  PEAK = sqrt(power))


%the following is copied from distort to create the filters:

fs = 44100;

numbins = length(f)-1;

f1 = f(1:(length(f)-1));
f2 = f(2:length(f));
freq = f1+(f2 - f1)/2; %center frequencies needed for synthesis later

%normalize these values so that we can make our filters.  Butter expects
%values between 0 and 1 where 1 corresponds to fs/2.  So when f = fs/2 we
%should get an fnorm of 1.
fnorm = 2* f / fs;

%now make the filter bank
f1norm = fnorm(1:(length(f)-1)); %left norms
f2norm = fnorm(2:length(f)); %right norms
fcnorm = 2*freq/fs; %normalized center frequencies

%these are the new bin 3dB cutoffs
fleft = fcnorm - (f2norm-f1norm) * separation / 2;
fright = fcnorm + (f2norm-f1norm) * separation / 2;

%need to make this generic if we want filters to be a different order
A = zeros(numbins, order + 1);
B = zeros(numbins, order + 1);

%figure out what type of filter to make, not bothering to check the string:
type = ftype(1);

if type ~= 'b'
    R = str2double(ftype(2:end));
end

%make array of H to compute freqz
H = zeros(numbins, length(fplot));

for i = 1:numbins
    if type == 'b'
        [B(i,:), A(i,:)] = butter(order/2, [fleft(i), fright(i)]);
    else
        [B(i,:), A(i,:)] = cheby1(order/2,R, [fleft(i), fright(i)]);
    end
    
    H(i,:) = abs(freqz(B(i,:), A(i,:),fplot,fs));
end

%H is in power ratios, we want the value of peak that will be created.
%CORRECTION: H is a peak response, leave it

for i = 1:numbins
    plot(fplot, H(i,:))
    hold on
end

axis([fplot(1) fplot(end) -.1 1.1])