close all; clear all; clc;

% Set filter sizes
M0 = 63;
M1 = 31;
M2 = 15;

% Calculate the theoretically expected group delay
expected = (M0-1)/2 + 4*(M1-1)/2 + 8*(M2-1)/2 + (M0-1)/2 + 4*(M1-1)/2;

% Design some example filters
function h = design_filter(N, Fp, Ap, Ast)
    Rp = (10^(Ap/20) - 1)/(10^(Ap/20) + 1); 
    Rst = 10^(-Ast/20);
    h = firceqrip(N-1, Fp, [Rp Rst], 'passedge');
end

h0 = design_filter(M0, 1/4, 0.1, 60);
h1 = design_filter(M1, 1/2, 0.1, 60);
h2 = design_filter(M2, 1/3, 0.1, 60);

% Sanity check: M is filter length
assert(length(h0) == M0);
assert(length(h1) == M1);
assert(length(h2) == M2);

% Generate impulse input
Nx = 2*expected;
x = zeros(Nx,1);
x(1) = 1;

% Apply the multi-rate filter
yd0 = downsample(filter(h0, 1, x), 4);
yd1 = downsample(filter(h1, 1, yd0), 2);
yf = filter(h2, 1, yd1);
yi0 = filter(h1, 1, upsample(yf, 2));
yi1 = filter(h0, 1, upsample(yi0, 4));

% Do the same, but rearranged using the Noble identities
h01 = conv(h0, upsample(h1,4));
zd01 = downsample(filter(h01, 1, x), 8);
zf = filter(h2, 1, zd01);
zi01 = filter(h01, 1, upsample(zf, 8));

% Plot
figure(); hold on; grid on;
plot(0:length(yi1)-1, yi1, 'b');
plot(0:length(zi01)-1, zi01, '--r');
xline(expected, '--k');
axis tight;
xlabel('Delay');
ylabel('Amplitude');
legend('Impulse response (direct)', 'Impulse response (Noble)', 'Theoretical group delay');
