There are two methods to tackling this problem:
- As pointed out in Yen-Ta Huang's answer and also in this Everett You's answer (EY16) to this related question we can split the creation and annihilation operators into a real and an imaginary part.
- As hinted at in (Capri, 2002; pg448) we can generalize the Bogoliubov transform to work with complex Hamiltonians.
Here I will do a simple example with the following fermionic Hamiltonian:
$$H=\varepsilon c_1^\dagger c_1+\varepsilon c_2^\dagger c_2+\lambda i(c_1^\dagger c_2^\dagger-c_2c_1)\tag{1}$$
Method 1
We let:
$$c_j=a_j+i b_j\quad \text{for}\quad j=1,2 \tag{2}$$
where $a_j^\dagger=a_j$ and $b_j^\dagger=b_j$. As shown in EY16 for $a_j$ and $b_j$ we have the following commutation relations
$$\{a_j,a_j\}=\{b_j,b_j\}=1$$
$$\{a_1,a_2\}=\{b_1,b_2\}=\{a_i,b_j\}=0$$
Thus subbing (2) into (1) we get that (after some algebra):
$$H=2i (\varepsilon a_1b_1+\varepsilon a_2 b_2+\lambda a_1 a_2-\lambda b_1 b_2)$$
$$=2i\begin{pmatrix} a_1 &b_2 \end{pmatrix}\begin{pmatrix} \varepsilon & \lambda \\ \lambda &\varepsilon \end{pmatrix}\begin{pmatrix} b_1 \\ a_2\end{pmatrix}$$
As explained in EY16 a Bogoliubov transformation of $a_j$ and $b_j$ is an orthogonal transformation in the case of fermions. Thus if we let:
$$\begin{pmatrix} b_1 \\ a_2\end{pmatrix}=\begin{pmatrix} \cos(\theta) & \sin(\theta)\\ -\sin(\theta) & \cos(\theta)\end{pmatrix} \begin{pmatrix} e_1 \\ d_2\end{pmatrix}$$
$$\begin{pmatrix} a_1 \\ b_2\end{pmatrix}=\begin{pmatrix} \cos(\theta) & \sin(\theta)\\ -\sin(\theta) & \cos(\theta)\end{pmatrix} \begin{pmatrix} d_1 \\ e_2\end{pmatrix}$$
with the new fermionic creation and annihilation operators being given by $f_j=d_j+ie_j$ with an appropriate choice of $\theta$ this will diagonalize the Hamiltonian
Method 2
In method 2 we simply generalize the Bogoliubov transformation. Consider the transformation:
$$f_j=u_jc_j+v_j c_j^\dagger$$
we are needing to enforce the conditions that:
$$\{f_i,f_j\}=0, \quad \{ f_i,f_j^\dagger\}=\delta_{ij}$$
If we do this we get that we need:
$$u_1v_2+u_2v_1=0\tag{3}$$
and
$$|u_j|^2+|v_j|^2=1\tag{4}$$
(4) implies that we have:
$$u_j=\cos(\theta_j) e^{i\phi_j^u}\quad v_j=\sin(\theta_j) e^{i\phi_j^v}$$
whilst with these (3) implies that:
$$\cos(\theta_1)\sin(\theta_2)=-\cos(\theta_2) \sin(\theta_1),\quad \phi_1^u+\phi_2^v=\phi_2^u+\phi_1^v$$
Putting these together the general Bololiubov transformation of fermionic operators is:
$$e^{i\tilde \phi_1} \begin{pmatrix}e^{i\tilde \phi_2} \cos(\theta_p) & e^{i\tilde \phi_3}\sin(\theta_p)\\ -e^{-i\tilde \phi_3}\sin(\theta_p) & e^{-i\tilde \phi_2}\cos(\theta_p) \end{pmatrix}$$
The standard method of the Bololiubov transformation can then be followed with this.
For reference the general Bololiubov transformation for bosons is (according to my calcuations:
$$e^{i\tilde \phi_1} \begin{pmatrix}e^{i\tilde \phi_2} \cosh(\theta_p) & e^{i\tilde \phi_3}\sinh(\theta_p)\\ e^{-i\tilde \phi_3}\sinh(\theta_p) & e^{-i\tilde \phi_2}\cosh(\theta_p) \end{pmatrix}$$