Designing a multi-band FIR filter on MATLAB I am trying to plot the frequency response of a dual recursive running sum (RRS) filter, with the $\mathcal 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.
According to your filter equations, your polynomials for $b1$, $b2$ and $a$ should only contain 1, zeros and -1. Your don't need to evaluate the last coefficient (num_exp, and others) at the complex exponential - That is done inside the freqz()
Take a look at this example here
So you should have:
a=[1, -1]
b1 = [1, 0, 0, ..., 0, -1]
b2 = [1, 0, 0, ..., 0, -1]

So $a$, $b1$, $b2$ contain the coefficients of the $z^{-n}$ terms.