Below summarizes efficient phase estimators for this application updated to include both a phase range of +/- 30 degrees and +/- 60 degrees. This is given in two parts, estimators for a real IF (intermediate frequency) signal, and estimators for a baseband complex signal. At the end are additional considerations related to acquisition.
For an additional estimator provided by Richard Lyons, please see his answer at this other post.
Efficient Phase Estimators for Real IF Signals
Product detector: For real signals, a common phase estimator (detector) is a multiplier followed by a low pass. For this application where sensitivity is desired over $\pm 30$ degrees, the signals are nominally 90 degrees in phase resulting is an estimate that is proportional to the sine of the phase between two signals:
$$y(t,\phi) =A_1\cos(\omega_ct)A_2\sin(\omega_ct+\phi) = \frac{A_1 A_2\sin(\phi) + A_1A_2\sin(\omega_ct+\phi)}{2} $$
Where when followed with a low pass filter removes the time varying component resulting in:
$$y(\phi) = \text{LPF}[y(t,\phi)] =\frac{A_1A_2}{2}\sin(\phi) $$
Where $\text{LPF}[\cdot]$ is the time average provided by a low pass filter.
As a demodulator, this can be implemented either in a coherent receiver where $A_1\cos(\omega_ct)$ is estimated during acquisition and provided as an NCO in a digital implementation (or VCO in an analog), or in a non-coherent receiver (where the interest is in the phase difference between two successive symbols) where the demodulation is done by multiplying the received signal with a time delayed copy of itself, delayed by one symbol duration plus the duration of a quarter cycle of the carrier to convert $\cos$ to $\sin$ (when the IF carrier is sufficiently larger than the symbol rate):
$$y(\phi_2-\phi_1) = \text{LPF}\bigg[A\cos(\omega_ct+\phi_1)A\cos(\omega_c(t-T_s-T_c)+\phi_2)\bigg]$$
Where $T_s$ is the symbol duration in seconds and $T_c = 1/(4f_c)$ is a quarter cycle of the IF carrier in seconds with the IF carrier frequency as $f_c$ in Hz. Resulting in:
$$y(\phi_2-\phi_1) = \text{LPF}\bigg[A\cos(\omega_ct+\phi_1)A\sin(\omega_c(t-T_s)+\phi_2)\bigg]$$
$$ y(\phi_2-\phi_1) = \frac{A^2}{2}\sin(\phi_2-\phi_1)$$
$$ y(\Delta\phi) = \frac{A^2}{2}\sin(\Delta\phi)$$
A very efficient way to implementing either of the above approaches digitally is by hard limiting the input signal which reduces the above to a simple XOR of the most significant bit of the waveform. For accuracy, this requires ensuring the inputs to the XOR operation are a 50% duty cycle, but the result is linearly proportional to phase! It is usable over a $\pm 90°$ range with a linear phase result. Further, hard limiting a phase modulated waveform provides a 3 dB SNR improvement in positive SNR conditions (since all AM noise is removed), but can be more susceptible to jamming an interference (3 dB loss in negative SNR conditions). This is an approach to be considered due to the simplicity and high phase linearity.
As above, the X-or phase detector could be used in a coherent receiver where the NCO is also simplified to a 1 bit output (basically the MSB of a counter and you increment the count rate to adjust the frequency as part of a carrier tracking loop), or non-coherently where the MSB of the received signal is XOR'd with a delayed copy. As with the multiplier the two input signals would be in quadrature to center the detector over its unambiguous range.
Efficient Phase Estimators for Complex Baseband Signals
Given a generic signal as $Ae^{j\phi}= I + jQ$, the actual phase is given by $\phi =\tan^{-1}(Q/I)$ or $\phi =\sin^{-1}(Q/A)$.
The following is a summary of efficient phase demodulation approximations for a variation over $\pm30$ degrees and $\pm60$ degrees, assuming carrier recovery and timing is established during the $0 / 180$ acquisition period. Initial thoughts on approaches to efficient acquisition are also included at the bottom of this post.
Summary of Results
Below is a table summarizing the peak and rms phase error for various estimators. Estimators that were included in earlier versions of this post that offer no advantage to those listed below have been removed. As Ben suggests in the comments, the Q/A estimators are attractive for FPGA implementation since A is assumed to be constant over the duration of the packet.
Plots showing the relative performance are included below:
Detailed Descriptions
The estimators that are scaled by the envelope magnitude $A$ (Q/A, Q/A Juha and Q/Est(A)) are preferred since $A$ can be readily determined during acquisition of the 0/180 signal, and only needs to be determined once for relatively short packets, or is a parameter from the AGC otherwise. In a constant envelope phase modulated signal such as this, the received signal can be simply hard limited if there isn't a concern with potential 3dB loss with stronger out of band interference (or the complete loss from hard-limiting in the presence of a coherent jammer). Further, there is no need to actually divide by $A$, assuming $A$ is maintained to be constant over the packet duration the result will be linearly proportional to the phase and the decision thresholds can be set accordingly.
Q/A
$$\phi =\sin^{-1}\bigg(\frac{Q}{A}\bigg)$$
$$\frac{Q}{A} = sin(\phi)$$
for small $\phi$, $sin(\phi) \approx \phi$ for $\phi$ in radians:
$$\phi \approx \frac{Q}{A}$$
Q/A Juha
Similar to @JuhaP's suggestion in the comments of removing the linear slope error for the $Q/I$ estimator, here applied to the Q/A Estimator. The coefficient is found from the linear portion of the remaining terms in the Taylor Series expansion that weren't used, minimizing the error:
For ±30° Operation:
$$\phi \approx 1.0475\frac{Q}{A}$$
For ±60° Operation:
$$\phi \approx 1.150\frac{Q}{A}$$
Q/Est(A)
A fast and very efficient approach for estimating magnitude is the $\alpha$ max plus $\beta$ min algorithm where the maximum between $|I|$ and $|Q|$ scaled by coefficient $\alpha$ is added to the minimum scaled by coefficient $\beta$. At 30° range, $Q$ would always be the minimum and $I$ always positive so this would simplify to $\alpha I + \beta|Q|$. A common choice for FPGA implementation is $\alpha = 1$ and $\beta =1/2$ since this minimized the error over all phases with bit shift divisions, but in this case $\alpha = 1$ and $\beta =1/4$ is a better choice given the narrowed phase range of $±30°$. If multipliers were acceptable, the optimized coefficients are $\alpha = 0.961$ and $\beta =0.239$. The plot below summarizes the two choices:
$$\phi \approx \frac{Q}{\alpha I + \beta |Q|}$$
option 1: $\alpha =1$, $\beta = 0.25$
option 2: $\alpha =0.961$, $\beta = 0.239$
Also not plotted below but shown above is the option optimized for use over ±60°:
$\alpha =0.85$, $\beta = 0.45$
Note these are not optimized for the estimate of $A$, but to minimize the phase estimation error.
Q/I Phase Approximation
$$\phi =\tan^{-1}\bigg(\frac{Q}{I}\bigg)$$
$$\frac{Q}{I} = tan(\phi)$$
for small $\phi$, $tan(\phi) \approx \phi$ for $\phi$ in radians:
$$\phi \approx \frac{Q}{I}$$
As @JuhaP mentioned in the comment, the linear slope component of the error could be removed by multiplying by 0.9289 resulting in (This one is labeld Q/A JuhaP in the plot). The coefficient below is slightly different than his suggestion but minimizes the error as it was found from the linear portion of the remaining terms in the Taylor Series expansion that weren't used rather than his approach of a first order terms of a polynomial fit to arctan:
$$\phi \approx 0.9289\frac{Q}{I}$$
Taylor Series Phase Approximations
The first term is the Q/A and Q/I approximations covered above for $\sin^{-1}$ and $\tan^{-1}$ respectively. Going beyond that is NOT recommended if efficiency is paramount but included for accuracy comparison.
arcsin
$$sin^{-1}(n) = \sum_{n=0}^\infty \frac{2n!}{2^{2n}(n!)^2}\frac{x^{2n+1}}{2n+1} \text{ for } |n|\le1$$
$$sin^{-1}\bigg(\frac{Q}{A}\bigg)= \frac{Q}{A} +\frac{1}{6} \bigg(\frac{Q}{A}\bigg)^3 +\frac{3}{40}\bigg(\frac{Q}{A}\bigg)^5 ... \text{ for } |Q/A|\le1$$
arctan
$$tan^{-1}(n) = \sum_{n=0}^\infty (-1)^n\frac{x^{2n+1}}{2n+1} \text{ for } |n|\le1$$
$$tan^{-1}\bigg(\frac{Q}{I}\bigg) = \frac{Q}{I} -\frac{1}{3} \bigg(\frac{Q}{I}\bigg)^3 +\frac{1}{5}\bigg(\frac{Q}{I}\bigg)^5 ... \text{ for } |Q/I|\le1$$
Using the first two terms for each results in:
$$\phi \approx \frac{Q}{A} +\frac{1}{6} \bigg(\frac{Q}{A}\bigg)^3$$.
$$\phi \approx \frac{Q}{I} -\frac{1}{3} \bigg(\frac{Q}{I}\bigg)^3$$.
A linear slope can also be removed from either of these with a gain constant multiplication as was done with the Q/A and Q/I estimators to further minimize the error.
Other Estimators
Juha Squared
@JuhaP offered this interesting estimator in the comments. Not very efficient but highly accurate with square terms:
$$\phi \approx \frac{3QI}{Q^2 + 3I^2}$$.
Acquisition
Efficient Acquisition for 0/180 preamble:
One idea that comes to mind for acquisition during the $0/180$ transitions is to use $\text{sign}(I_2)Q_1-\text{sign}(I_1)Q_2$ to get the change in phase between two symbols, which can the be corrected in a fast coverging and simple loop by derotating the incoming signal. This appraoch would work well if the frequency offset is such that the phase does not rotate more than $\pm \pi/2$ between successive signals, otherwise a course FLL can be used first to get the offset within this acquisition range.
For a coherent receiver approach a PLL would be used to lock/track an NCO or PLL to the carrier and my squaring the received signal a reference tone at twice the carrier can be tracked for all modulations presented here (both the bi-phase acquisition interval and the 30 degree phase modulation when doubled will produce a distinct tone at 2x the carrier). Similarly a Costas Loop would track both signals while providing the reference signal that is nominally 90 degrees in phase with the carrier, thus providing both carrier recovery and phase demodulation.
Sources:
Taylor Series Expansions for arcsin and arctan:
https://proofwiki.org/wiki/Book:Murray_R._Spiegel/Mathematical_Handbook_of_Formulas_and_Tables
