%**************************************************************************
% Parameter
%**************************************************************************
f_s          = 44100;
f_crossover  =   150;
N_butter_pre =     2;

%**************************************************************************
% Design (preliminary) Butterworth filter 
%**************************************************************************
[b_hp, a_hp] = butter(N_butter_pre, f_crossover/f_s*2, 'high');            
[b_lp, a_lp] = butter(N_butter_pre, f_crossover/f_s*2, 'low');

%**************************************************************************
% "Double" filters
%**************************************************************************
b_lp_ent = conv(b_lp,b_lp);
a_lp_ent = conv(a_lp,a_lp);
b_hp_ent = conv(b_hp,b_hp);
a_hp_ent = conv(a_hp,a_hp);

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

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

H_lp_ent = freqz(b_lp_ent, a_lp_ent, N_FFT);
H_hp_ent = freqz(b_hp_ent, a_hp_ent, N_FFT);
f_help               = linspace(0, f_s/2, length(H_lp_ent));

subplot('Position',[0.07 0.55 0.9 0.38]);
plot(f_help, 20*log10(abs(H_lp_ent)+eps),'b', ...
     f_help, 20*log10(abs(H_hp_ent)+eps),'r', ...
     f_help, 20*log10(abs(H_hp_ent + H_lp_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_ent = grpdelay(b_lp_ent, a_lp_ent, N_FFT);
g_hp_ent = grpdelay(b_hp_ent, a_hp_ent, N_FFT);
f_help               = linspace(0, f_s/2, length(H_lp_ent));

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