The desired inequality should be true iff
$$
c < c_0 := (r - \sqrt{r^2-1})^2
\quad\ \text{where} \quad\
r = \frac{|a|}{2b}
$$
(NB the hypotheses $b>0$ and $a < -2b$ imply $r>1$, so $0 < c_0 < 1$).
Numerical computation suggests that $f(n) \sim A c_0^n \left/ \sqrt{n} \right.$
for some $A>0$. It should be possible to prove this by completing
the following analysis.
It will be convenient to set $N = n+1$, and sum over $j,k$ rather than
$i,j$ because we'll need $i = \sqrt{-1}$.
For $x,y \in ({\bf R}/{\bf Z})^2$ define
$$
F(x,y) = \frac1{2b}\frac
{(1 - \cos(4\pi x)) (1 - \cos(4\pi y))}
{2r - (\cos 2\pi x + \cos 2 \pi y)}.
$$
Then
$$
f(n) = \frac1{N^2}
\sum_{j=1}^{N-1} \sum_{k=1}^{N-1} (-1)^{j+k}
F\Bigl(\frac{j}{2N}, \frac{k}{2N}\Bigr) \, .
$$
Now $F(x,y) = F(-x,y) = F(x,-y)$, and
$$
F(0,y) = F(1/2,y) = F(x,0) = F(x,1/2) = 0
$$
for all $x,y$ thanks to the factors
$1 - \cos(4\pi x)$ and $1 - \cos(4\pi y)$ in the numerator of $F$.
So we can write $f(n)$ as an alternating sum over $\frac1{2N}$-lattice points:
$$
f(n) = \frac1{(2N)^2}
\sum_{j=1}^{2N} \sum_{k=1}^{2N} (-1)^{j+k}
F\Bigl(\frac{j}{2N}, \frac{k}{2N}\Bigr)
$$
(each term in the original sum for $f$ appears four times here,
and the added terms with $j,k \in \{N, 2N\}$ all vanish).
But such a sum can be expressed in terms of the Fourier expansion
$$
F(x,y) = \sum_{r\in\bf Z}\sum_{s\in\bf Z} \phi(r,s) \exp 2\pi i (rx+sy)
$$
of $F$, which converges absolutely because $F$ is smooth.
The alternating sum of $\exp 2 \pi i (rx+sy)$ is $(2N)^2$ if
$(r,s) \equiv (N,N) \bmod 2N$, and zero otherwise. Therefore
$f(n) = \sum\!\sum_{r,s} \phi(r,s)$ with the sum extending over all $(r,s)$
such that $r \equiv s \equiv N \bmod 2N$. The simplest terms
in this sum are the four coefficients $\phi(\pm N, \pm N)$,
which are all equal (since they satisfy $\phi(r,s) = \phi(-r,s) = \phi(r,-s)$ for all $r,s$).
Thus we expect that
$$
f(n) \sim 4 \phi(-N,-N) = 4 \int_{-1/2}^{1/2} \int_{-1/2}^{1/2}
\exp (2\pi i N (x+y)) \, F(x,y) \, dx \, dy
$$
unless $\phi(N,N)$ is unusually small.
We can now explain and quantify the observed exponential cancellation
in the alternating sum that defines $f$. For large $r,s$
the Fourier coefficients $\phi(r,s)$ decay exponentially but no faster,
because $F$ is not just smooth but analytic for
$x,y \in ({\bf C}/{\bf Z})^2$ in some neighborhood of $({\bf R}/{\bf Z})^2$,
but does have singularities for some complex $(x,y)$.
We can shift the integral in the complex direction:
$$
\phi(N,N) = e^{-4\pi N w} \int_{-1/2}^{1/2} \int_{-1/2}^{1/2}
\exp (2\pi i N (x+y)) \, F(x+iw,y+iw) \, dx \, dy
$$
as long as $F$ is analytic for
$\left|\mathop{\rm Im}(x)\right|, \left|\mathop{\rm Im}(y)\right| \leq w$.
This is the case for all
$$
w < w_0 := \frac1{2\pi} \cosh^{-1} r
= \frac1{2\pi} \log \bigl(r + \sqrt{r^2-1}\bigr)
= -\frac1{4\pi} \log c_0.
$$
In fact we can take $w = w_0$: the integrand
$F(x+iw_0,y+iw_0)$ blows up at $(x,y)=(0,0)$, but the integral
still converges absolutely because if $\left| F(x+iw_0,y+iw_0) \right| > M$
then $|x+y| \ll M^{-1}$ and $|x-y| \ll M^{-1/2}$. Thus
$$
\phi(N,N) = c_0^N \int_{-1/2}^{1/2} \int_{-1/2}^{1/2}
\exp (2\pi i N (x+y)) \, F(x+iw_0,y+iw_0) \, dx \, dy.
$$
This shows that $f(n) \ll c_0^N$, and thus also $f(n) \ll c_0^n$;
some more analysis of the double integral near the singularity at $(0,0)$
should show that $f(n) / c^n \to \infty$ for all $c<c_0$.
