function [h,N,fc] = firdig(fp, fs, Ap, As, fm, window)
%
% function [B,A] = firdig(fp, fs, Ap As, fm, window)
%
% This function calculates an FIR digital filter with the window specified by variable window
%   type can be 'boxcar', 'hamming', 'hanning', 'blackman' or 'bartlett' 
%
% fp is a vector of passband frequencies (Hz)
% fs is a vector of stopband frequencies (Hz)
% Ap is the attenuation in db at the passband frequencies
% As is the attenuation in db at the stopband frequencies
% fm is the sampling frequency (Hz)
%


% Check the parameters fp and fs

close all;

np = size(fp);
ns = size(fs);

if (~isequal(np,ns))
    fprintf(1,'Error: fp and fs are of different sizes\n'); return;
elseif (min(np) == 1)  % fp and fs are vectors of the same size
    if (max(np) == 1)  % fp and fs are scalars (filter must be either a lowpass or highpass)
        if fp < fs
            filter = 'lp';
        elseif fp > fs
            filter = 'hp';
        else
            fprintf(1,'Error: fs and fp are equal\n'); return;
        end
    elseif (max(np)==2) % fp and fs are vectors (filter must be a bandpass or bandstop)
        if (fs(1) < min(fp)) & (fs(2) > max(fp))
            filter = 'bp';
        elseif (fp(1) < min(fs)) & (fp(2) > max(fs))
            filter = 'bs';
        else 
            fprintf(1,'Error: Values of fp and fs does not correspond to a known kind of filter\n'); return;
        end
    else
        fprintf(1,'Error: The length of fs or fp is bigger than 2\n'); return;
    end 
else
    fprintf(1,'Error: fp, fs or both are matrices\n'); return;
end

% Ordering fp and fs from smaller to bigger frequencies

if filter=='bp' | filter=='bs'
    x = min(fp); y = min(fs);
    if (x==2) tmp = fp(1); fp(1) = fp(2); fp(2) = tmp; end
    if (y==2) tmp = fs(1); fs(1) = fs(2); fs(2) = tmp; end
end


% Check the Attenuation parameters

ap = size(Ap);
as = size(As);

if (~isequal(ap,as))
    fprintf(1,'Error: Ap and As are of different sizes\n'); return;
elseif (ap ~= [1 1])
    fprintf(1,'Error: Ap and As are not scalars\n');return;
elseif (Ap <= 0 | As <= 0)
    fprintf(1,'Error: Ap or As or both are less or equal than 0\n'); return;
end

%% Check the sampling frequency
maxfreq = max(max([fp fs]));

if fm < 2*maxfreq
    fprintf(1,'Warning: fm < twice the maximum frequency specified in in fp and fs\n');
    fprintf(1,'         I suggest to increase the sampling frequency\n');
    return;
end

% Convert to normalized frequencies
%
Fp = fp/fm;
Fs = fs/fm;

%% Convert specs to lowpass

switch filter
case {'lp'}
    Fp_lp = Fp;
    Fs_lp = Fs;
    Atp = Ap;
    Ats = As;
    Apfig = Ap;
    Asfig = As;
case {'hp'}
    Fp_lp = 0.5-Fp;
    Fs_lp = 0.5-Fs;
    Atp=Ap;% -20*log10(1-10^(-As/20));
    Ats=As;% -20*log10(1-10^(-Ap/20));
    Apfig = Ap;
    Asfig = As;
case {'bp'}
    [minimum,idx]=min([Fp(1)-Fs(1) Fs(2)-Fp(2)]);
    if idx==1,
        Fs(2)=Fp(2)+minimum;
    else
        Fs(1)=Fp(1)-minimum;
    end
    Fp_lp = (Fp(2)-Fp(1))/2;
    Fs_lp = (Fs(2)-Fs(1))/2;
    Fo = (Fp(1)+Fp(2))/2;
    Atp = Ap;
    Ats = As;
    Apfig = [Ap Ap];
    Asfig = [As As];
case {'bs'}
    [minimum,idx]=min([Fs(1)-Fp(1) Fp(2)-Fs(2)]);
    if idx==1
        Fp(2)=Fs(2)+minimum;
    else
        Fp(1)=Fs(1)-minimum;
    end  
    
    Fp_lp = (Fs(2)-Fs(1))/2;
    Fs_lp = (Fp(2)-Fp(1))/2;
    Fo = (Fp(2)+Fp(1))/2;
    
    Atp=-20*log10(1-10^(-As/20));
    Ats=-20*log10(1-10^(-Ap/20));
    
    Apfig = [Ap Ap];
    Asfig = [As As];
otherwise
    fprintf(1,'Error: unknown window specified');
end

%% Check the filter window

bw = Fs_lp-Fp_lp;
lwindow = lower(window);

switch lwindow
case {'boxcar'}
    N = round(0.81/bw);
case {'hamming'}
    N = round(1.91/bw);
case {'hanning'}
    N = round(1.97/bw);
case {'blackman'}
    N = round(2.82/bw);
case {'bartlett'}
    N = round(1.62/bw);
otherwise
    fprintf(1,'Error: Unknown filter window'); return;
end


%% Design the Lowpass filter
%% Set Fc (cut frequency) to Fp_lp
Fc = Fp_lp;

run = 1;
t_steps=4;
while run==1
    n = 0:N-1;
    eval(['w = ' lwindow '(' num2str(N) ');']);
    hlp = 2*Fc*sinc(2*Fc*(n-(N-1)/2)).* w';
    At = freqz(hlp,1,[Fp_lp Fs_lp]*fm,fm);
    Atdb = -20*log10(abs(At));
    dFc = (Fs_lp-Fp_lp)/t_steps;
    if (Atdb(1) > Atp) | (Atdb(2) < Ats)
        if (Fc > Fs_lp)
            N = N + 1;
	    if rem(N,2)==0 & filter=='bs' %% if filter is bandstop
                                          %the order must be odd
	      N=N+1;
	    end
            Fc = Fp_lp;
        else
	  Fc = Fc + dFc;
	end
	
    elseif (t_steps<256)
      % Go back to previous order. Just in case!!
      if (filter=='bs')
	N=N-2;
      else
	N=N-1;
      end
      Fc = Fc-dFc;  %% Go back to previous freq
      t_steps = t_steps*2;
    else
      run=0;
    end
end

fc = Fc*fm;

[H,f] = freqz(hlp,1,500,1);
Hdb = 20*log10(abs(H));
figure, subplot(2,1,1), plot(f, Hdb, [Fp_lp Fs_lp], -[Atp Ats], '*r');

%% Transform to aproppiate filter
switch filter
case {'hp'}
    delta = zeros(1,N); delta((N+1)/2) = 1;
    h = ((-1).^n).*hlp; 
case {'bp'}
    h = 2*cos(2*pi*Fo*(n-(N-1)/2)).*hlp;
case {'bs'}
    delta = zeros(1,N); delta((N+1)/2) = 1;
    h = delta - 2*cos(2*pi*Fo*(n-(N-1)/2)).*hlp;
case {'lp'}
    h = hlp;
otherwise
end

[H,f] = freqz(h,1,500,fm);
Hdb = 20*log10(abs(H));
subplot(2,1,2), plot(f, Hdb, [fp fs], -[Apfig Asfig], '*r');
axis([-inf inf -As-20 1]);

figure; stem(h);

fprintf(1,'The order of the resulting filter is %d\n', N-1);
fprintf(1,'The cut frequency of the lowpass prototype is %f Hz\n',fc);

return;