%**************************************************************************
% Filer parameters
%**************************************************************************
f_s    = 44100;                      % Sample rate
f_low  =    80;
f_high =   200;
sb_att =    80;
pb_att =     0.001;

%**************************************************************************
% Related parameters
%**************************************************************************
f_h    = [f_low f_high];                % Cutoff frequencies
a_h    = [1 0];                         % Desired amplitudes
dev    = [(10^(pb_att/20)-1)/(10^(pb_att/20)+1) 10^(-sb_att/20)];

%**************************************************************************
% Filter order estimation
%**************************************************************************
[N_lp,fo,ao,w] = firpmord(f_h,a_h,dev,f_s);
if (mod(N_lp,2) == 1)
    N_lp = N_lp + 1;
end

%**************************************************************************
% Lowpass filter design
%**************************************************************************
h_lp        = firpm(N_lp,fo,ao,w);

%**************************************************************************
% Corresponding highpass filter
%**************************************************************************
h_hp           = zeros(1,N_lp+1);
h_hp(N_lp/2+1) = 1;
h_hp           = h_hp - h_lp;

%**************************************************************************
% Analysis
%**************************************************************************
fig = figure(1);
set(fig,'Units','Pixel');
set(fig,'Position',[200 100 1000 600]);

N_FFT = 4096*16;
lw    =       1;
f_max =     500;

H_lp    = freqz(h_lp, 1, N_FFT);
H_hp    = freqz(h_hp, 1, N_FFT);
H_ent   = freqz(h_lp+h_hp, 1, N_FFT);
f_help  = linspace(0, f_s/2, length(H_lp));

subplot('Position',[0.07 0.55 0.9 0.38]);
plot(f_help, 20*log10(abs(H_lp)+eps),'b', ...
     f_help, 20*log10(abs(H_hp)+eps),'r', ...
     f_help, 20*log10(abs(H_ent)+eps),'k', ...
     'LineWidth',lw);
legend('Lowpass', 'Highpass', 'Sum / allpass')
ylabel('Magnitude in dB')
set(gca,'XTickLabel','');
grid on
ylim([-100, 5])
xlim([0, f_max])
title('Magnitude responses');

g_lp  = grpdelay(h_lp, 1, N_FFT);
g_hp  = grpdelay(h_hp, 1, N_FFT);
f_help               = linspace(0, f_s/2, length(H_lp));

subplot('Position',[0.07 0.10 0.9 0.38]);
plot(f_help, g_lp,'b', ...
     f_help, g_hp+10,'r', ...
     'LineWidth',lw);
legend('Lowpass', 'Highpass (+10 samples)')
ylabel('Group delay in samples')
xlabel('Frequency in Hz')
grid on
ylim([length(h_lp)/2-100, length(h_lp)/2+100])
xlim([0, f_max])
title('Group delays');