We start from the follwing expression of the Hamiltonian:
$$ H = \sum_{n=0}^{N-1} \left [t \left (c^\dagger_{i+1}c_{i} + c^\dagger_{i}c_{i+1}\right) - \mu f(i)\left ( 2c^\dagger_{i}c_{i}-1\right) + \sum_{l=1}^{N-1}\frac{\Delta}{d^\alpha_l}\left (c^\dagger_{i+l}c^\dagger_{i} + c_{i}c_{i+l}\right)\right], $$where $t$ is the hopping amplitude, $\mu$ is the chemical potential, $\Delta$ is the superconducting pairing amplitude and $c_{i}(c^\dagger_{i})$ are the annihilation (creation) operators at the $i$-th site of the chain. The superconducting pairing is taken to decay as a power law $l^{-\alpha}$, where $l$ denotes the distance between the sites and the scaling exponent $\alpha \in \mathbb{R}$. For a constant chemical potential $f(i)=1$, we recover the long-range Kitaev chain with homogeneous onsite potential. The AAH chemical potential corresponds to $f(i) = \cos\left (2 \pi n \frac{F_{n-1}}{F_{n}} + \phi \right )$, where $F_{n-1},F_n$ are integers from the Fibonacci sequence and $\phi$ is the phase. Therefore, our Hamiltonian has a periodicity of $F_n = q$ sites.
Let's consider the following basis $\chi = \left(c_{0}, c^\dagger_{0}, c_{1}, c^\dagger_{1}, ..., c_{N-1}, c^\dagger_{N-1}\right)^T$. The Hamiltonian can be expressed as follows:
$$ H = \chi^\dagger H_N \chi $$where: $$ H_N = \begin{pmatrix} A_0 & B & C_2 & C_3 & \cdots & C_{N-2} & C_{N-1}\\ B^\dagger & A_1 & B & C_2 & \cdots & C_{N-3} & C_{N-2}\\ C_2^\dagger & B^\dagger & A_2 & B & \cdots & C_{N-2} & C_{N-3}\\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots\\ C_{N-2}^\dagger & C_{N-3}^\dagger & C_{N-4}^\dagger & C_{N-5}^\dagger & \cdots & A_{N-2} & B\\ C_{N-1}^\dagger & C_{N-2}^\dagger & C_{N-3}^\dagger & C_{N-4}^\dagger & \cdots & B^\dagger & A_{N-1}, \end{pmatrix} $$
where
$$ A_{j} = \begin{pmatrix} -\mu f(j) & 0 \\ 0 & +\mu f(j) \end{pmatrix}, $$$$ B = \begin{pmatrix} t/2 & -\Delta \\ +\Delta & -t/2 \end{pmatrix}, $$$$ C_l = \begin{pmatrix} 0 & -\Delta/l^\alpha \\ +\Delta/l^\alpha & 0 \end{pmatrix}. $$We compute the Hamiltonian with OBC both in de BdG and the Majorana basis.
The Hamiltonian with anti-periodic boundary conditions takes the following form: $$ H = \sum_{i=0}^{N-2} \frac{t}{2} \left (c_i^\dagger c_{i+1} - c_{i+1} c_i^\dagger + c_{i+1}^\dagger c_i - c_{i} c_{i+1}^\dagger\right) -\frac{t}{2}\left(c_{N-1}^\dagger c_0 - c_0 c_{N-1}^\dagger + c_0^\dagger c_{N-1} - c_{N-1} c_0^\dagger\right)\\ + \sum_{i=0}^{N-1} \mu f(i)\left (c^\dagger_{i}c_{i}-c_{i}c^\dagger_{i}\right) +\sum_{i=0}^{N-1}\sum_{l=1}^{N-1-i} \frac{\Delta}{d^\alpha_l}\left (c_{i+1}^\dagger c_i^\dagger - c_i^\dagger c_{i+l}^\dagger + c_i c_{i+l} - c_{i+l} c_{i}\right). $$
where $d_l = \text{min}(l, N-l)$.
Since the system has a periodicity of $q$ sites, we can split it into supercells. Let's consider that we have $L = N/q$ of such supercells, and let's consider the spinor $\chi_u = \left(c_{qu}, c^\dagger_{qu}, c_{qu+1}, c^\dagger_{qu+1}, ..., c_{qu+(q-1)}, c^\dagger_{qu+(q-1)}\right)^T$. The Hamiltonian can be expressed as follows:
$$ H = \sum_{u=0}^{L-1} \left[ \chi_u^\dagger H_{\text{local}} \chi_u + ( \chi_u^\dagger H_{\text{hop}} \chi_{u+1} + {\rm h.c.}) + \sum_{l=1}^{L-1}( \chi_u^\dagger H_{\text{l}} \chi_{u+l} + {\rm h.c.}) \right], $$where,
$$ H_{\text{local}} = \begin{pmatrix} A_0 & B & C_2 & \cdots & C_{q-1}\\ B^\dagger & A_1 & B & \cdots & C_{q-2}\\ C_2^\dagger & B^\dagger & A_2 & \cdots & C_{q-3}\\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots\\ C_{q-2}^\dagger & C_{q-3}^\dagger & C_{q-4}^\dagger & \cdots & B\\ C_{q-1}^\dagger & C_{q-2}^\dagger & C_{q-3}^\dagger & \cdots & A_{q-1} \end{pmatrix}, $$$$ H_{\text{hop}} = \begin{pmatrix} 0 & 0 & \cdots & 0\\ 0 & 0 & \cdots & 0\\ \vdots & \vdots & \ddots & \vdots\\ B' & 0 & \cdots & 0 \end{pmatrix}, $$and,
$$ H_{\text{l}} = \begin{pmatrix} C_{l,0,0} & C_{l,0,1} & \cdots & C_{l,0,q-1}\\ C_{l,1,0} & C_{l,1,1} & \cdots & C_{l,1,q-1}\\ C_{l,2,0} & C_{l,2,1} & \cdots & C_{l,2,q-1}\\ \vdots & \vdots & \ddots & \vdots & \vdots\\ C_{l,q-2,0} & C_{l,q-2,1} & \cdots & C_{l,q-2,q-1}\\ C_{l,q-1,0} & C_{l,q-1,1} & \cdots & C_{l,q-1,q-1} \end{pmatrix}, $$where $B' = \frac{t}{2}\sigma_z$ and $C_{l,x,y} = -\frac{\Delta}{2d_{l,x,y}^{\alpha}} i\sigma_y$. Now,
$$ d_{l,x,y} = \text{min}\left(lq-(x-y), N-(lq-(x-y))\right). $$Note that in order to impose anti-periodic boundary conditions, we have used that $c_{uq+x+q} = -c_{uq+x}$ in $H_{\text{local}}$ and $c_{(u+L)q + x} = -c_{uq+x}$ in $H_{\text{hop}}$ and $H_{\text{l}}$.
Now, we want to write the Hamiltonian in momentum space. We transform the spinor $\chi_u$ as follows:
$$ \chi_u = \frac{1}{\sqrt{L}} \sum_{k} e^{iku} \chi_k, \quad \chi_u^\dagger = \frac{1}{\sqrt{L}} \sum_{k} e^{-iku} \chi_k^\dagger, $$where $\chi_k= \left(c_{k,0}, c^\dagger_{-k,0}, c_{k, 1}, c^\dagger_{-k, 1}, ..., c_{k, q-1}, c^\dagger_{-k, q-1})\right)^T$. For finite systems, there is a discrete set of momentums $k$ allowed in the system. Because of the anti-periodic boundary conditions, we want the following to hold: $$ \chi_{u+L} = -\chi_{u}, $$ therefore: $$ \frac{1}{\sqrt{L}} \sum_{k} e^{ik(u+L)} \chi_k = -\frac{1}{\sqrt{L}} \sum_{k} e^{iku} \chi_k, $$ meaning that $kL = (2m + 1)\pi$ where $m$ takes values from $\{0,1,2,...L-1\}$.
Now, if we do the Fourier transform of the Hamiltonian we obtain:
$$ H = \sum_{k} \left[ \chi_k^\dagger H_\text{local} \chi_k + \left(e^{ik} \chi_k^\dagger H_\text{hop} \chi_k + {\rm h.c.}\right) + \sum_{l=1}^{L-1}\left(e^{ikl}\chi_k^\dagger H_{l} \chi_{k} + {\rm h.c.}\right)\right], $$where we have used: $$ \sum_{u=0}^{L-1} e^{i(k-k')u} = L\delta_{kk'}. $$
For the infinite system, $L$ goes to infinity. Then, we need to change $d_l = l$ and $d_{l,x,y} = lq-(x-y))$ because the anti-periodic boundary conditions do not need to be explicitly imposed anymore. The Hamiltonian in momentum space looks as follows:
$$ H = \sum_{k} \left[\chi_k^\dagger H_\text{local} \chi_k + \left(e^{ik}\chi_k^\dagger H_\text{hop} \chi_k + \text{h.c.}\right) + \left(\chi_k^\dagger H_\text{inf} \chi_{k} + \text{h.c.}\right)\right] $$where:
$$ H_{l, inf} = \begin{pmatrix} C_{0,0} & C_{0,1} & C_{0,2} & C_{0,3} & \cdots & C_{0,q-2} & C_{0,q-1}\\ C_{1,0} & C_{1,1} & C_{1,2} & C_{1,3} & \cdots & C_{1,q-2} & C_{1,q-1}\\ C_{2,0} & C_{2,1} & C_{2,2} & C_{2,3} & \cdots & C_{2,q-2} & C_{2,q-1}\\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots\\ C_{q-2,0} & C_{q-2,1} & C_{q-2,2} & C_{q-2,3} & \cdots & C_{q-2,q-2} & C_{q-2,q-1}\\ C_{q-1,0} & C_{q-1,1} & C_{q-1,2} & C_{q-1,3} & \cdots & C_{q-1,q-2} & C_{q-1,q-1} \end{pmatrix}, $$where $C_{x,y} = ig_{xy}\sigma_y$. The function $g_{xy}$ is defined as follows
$$ g_{xy} = -\sum_{l=1}^{\infty} \frac{\Delta e^{ikl}}{[q-(x-y)]^\alpha} = \frac{e^{ik}}{q^{\alpha}} \text{HLP}_{\alpha} \left(k, q, x-y\right), $$where HLP is the Hurwitz Lerth Phi function, also called the Lerch trascendent function.