I am trying to plot the frequency response of a dual recursive running sum (RRS) filter, with the z-transform of each given by:

$\frac{1 - z^{-L_n}}{1 - z^{-1}}$

This is the code I've written:

% (1 - z^-L1) / (1 - z^-1) * (1- z^-L2) / (1 - z^-1) 
L1 = 75;
L2 = 53;
fc = 500;
fs = 44100;

den_exp = -exp(-2j*pi*(fc/fs));
num_exp = -exp(-2j*pi*(fc/fs)*L1);
num2_exp = -exp(-2j*pi*(fc/fs)*L2);

%z transform
b1 = [1, zeros(1, L1-2), num_exp];
b2 = [1, zeros(1, L2-2), num2_exp];
a = [1, den_exp];

[hA, w] = freqz(b1, a);
[hB, w] = freqz(b2, a);
hAB = hA .* hB;

plot((w*fs)/(2*pi), 20*log10(abs(hAB)))
xlabel('Frequency')
ylabel('Magnitude (dB)')
title('Recursive filter with cut-off = 500Hz') 

The frequency response appears to be off, and I can't seem to figure out why.

This is part of a project where I am trying to build a tone control system on MATLAB. I have read a couple of papers on the RRS, and I understand the idea of combining the transfer functions of high, low, and band pass with selected gains and delays, but I am struggling to implement this on MATLAB. For example, how does one add a delay block to a transfer function, as shown here? And how are these transfer functions added together?

I would be grateful if someone could provide some details on how to implement this.