function [y,h,G,fG] = bandpass(freqs,freqe,pbG,sbG,x,N,Fs,fhp,in)
%[y,h,G] = dsp_eq(sG,x,N,Fs,fhp,in)
%inputs:
%freqs -> the frequency in which to start the bandpass
%freqe -> the frequency in which to stop the bandpass
%pbG   -> the desired passband gain in dB
%sbG   -> the desired stopband gain in dB
%x    -> input vector to filter using the h coeficients
%N    -> desired taps for linear phase fir filter
%Fs   -> sampling frequency Hz
%fhp  -> figure handel for ploting
%in   -> index from 0-3 for a 2x2 subplot region to plot frequency resp.
%outputs:
%y   -> filtered data
%h   -> filter coefecients
%G   -> magnituded response points
%bandpass.m
%author: Nicholas Kottenstette
%date:   12/5/2004
  
%Technique design linear phase filter using frequency sampling
%choose coef. to be symetric with alpha equal to 0
%see p.633 Table 8.3 for details in Proakis
%Keep spacing the same if lower sampled file is equalized
  f = 0;
  df = Fs/N;
  freq = zeros(2*length(freqs)+1,1);
  freq(1:2:length(freq)-1) = freqs;
  freq(2:2:length(freq)-1) = freqe;
  freq(length(freq)) = Fs/2;
  sG = zeros(length(freq),1);
  sG(1:2:length(sG)-1) = sbG;
  sG(2:2:length(sG)-1) = pbG;
  sG(length(sG)) = sbG;
  fl = length(freq);
  sGl = length(sG);
  if fl ~= sGl
    error('bandpass:invalid_input',...
          'freq and sG vectors unequal lengths fl = %d, sGl = %d\n"',...
          fl,sGl);
  end
  if floor(N/2) == (N/2)
    U = N/2 - 1;
  else
    U = (N-1)/2;
  end
  %l = k+1;
  l = 1;
  flag = true;
  G=zeros(1,U);
  fG=zeros(1,U);
  for j = 1:length(freq)
    while f < freq(j)
      if flag == true
        sg = sG(j);
      else
        %sg = (sG(j) + sG(j-1))/2;
        sg = sG(j);
        flag = true;
      end
      G(l) = ((-1)^(l-1))*10^(sg/20);
      fG(l) = f;
      f = f+df;
      l = l+1;
    end
    if j > 2
      flag = false;
    end
  end
  h = zeros(1,N);
  for n = 0:N-1
    sum = 0;
    for k = 1:U
      sum = sum + G(k+1)*cos(2*pi*k*(n+0.5)/N);
    end
    h(n+1) = (1/N)*(G(1)+2.0*sum);
    %Apply a hamming window to min-ripples and keep reasonable aprox with
    %lower freq adjust.
    h(n+1) = (0.54 - 0.46*cos(2*pi*n/(N-1)))*h(n+1);
    %h(n+1) = (0.5 - 0.5*cos(2*pi*n/(N-1)))*h(n+1);
  end

  if fhp ~= -1
    %Plot resulting filter fft-response based on lever settings
    H = dsp_fft(h,N*2^4,[],[]);
    H = H(1:length(H)/2)';
    dfH = Fs/((N*2^4)-1);
    frH = [0:((N*2^4)-1)/2]*dfH;
    frG = [0:df:(Fs/2)-1];
    figure(fhp);subplot(2,2,in);plot(frH,20*log10(abs(H)),frG,20*log10(abs(G)),'.');grid;
    xlabel('freq (Hz)');ylabel('DFT (dB)');title('Bandpass Filter Response');
    axis([0,Fs/2,-40,10]);
  end
  %Implement equalizer filter:
  y = dsp_filter(h,x,N);
  if fhp ~= -1
    %Plot resulting spectrogram
    spec_plot(y,Fs,fhp,'Bandpass Spectrogram',in+1);
    %Plot original spectrogram
    spec_plot(x,Fs,fhp,'Original Spectrogram',in+3);
  end
  
function y = dsp_filter(h,x,N)
  lx = length(x);
  [m,n]=size(x);
  if (m > n)
    x = x';
  end
  x = [zeros(1,N-1) x zeros(1,N)];
  x1i = [0:-1:-(N/2)+1] + N;
  x2i = [-(N-1):1:-N/2] + N;
  h = h(1+[0:(N/2)-1])';
  y = zeros(1,lx+N-1);
  for n=1:(lx+N-1)
    x1 = x(n + x1i);
    x2 = x(n + x2i);
    y(n) = (x1+x2)*h;
  end
  y = y*(2/(max(y)-min(y)));
