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 $c(\pm N, \pm N)$, which are all equal (the Fourier coefficients inherit the symmetries of $F$ and thus 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$.

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$.