Without loss of generality assume that $O=(0,0)$ and the ABCD circle has radius 1.
Then
$$\begin{cases} A=(\cos \alpha,\sin\alpha),\\
B = (\cos\beta,\sin\beta),\\
C = (\cos\gamma,\sin\gamma),\\
D = (\cos\delta,\sin\delta),\\
AB = 2\sin\frac{\alpha-\beta}2,\\
BC = 2\sin\frac{\beta-\gamma}2,\\
CD = 2\sin\frac{\gamma-\delta}2,\\
AD = 2\sin\frac{\alpha-\delta}2,\\
\end{cases}
$$
where we assume $2\pi\geq \alpha\geq \beta\geq \gamma\geq\delta\geq 0$, and by Pitot theorem:
$$(\star)\qquad\sin\frac{\alpha-\beta}2 + \sin\frac{\gamma-\delta}2 = \sin\frac{\beta-\gamma}2 + \sin\frac{\alpha-\delta}2.$$
For the incircles coordinates we have
$$\begin{cases}
E = \tan(\frac{\pi}4-\frac{\alpha-\beta}4)\left(\cos\frac{\alpha+\beta}2,\sin\frac{\alpha+\beta}2 \right),\\
F = \tan(\frac{\pi}4-\frac{\beta-\gamma}4)\left(\cos\frac{\beta+\gamma}2,\sin\frac{\beta+\gamma}2 \right),\\
G = \tan(\frac{\pi}4-\frac{\gamma-\delta}4)\left(\cos\frac{\gamma+\delta}2,\sin\frac{\gamma+\delta}2 \right),\\
H = \tan(\frac{\pi}4-\frac{\alpha-\delta}4)\left(\cos\frac{\delta+\alpha}2,\sin\frac{\delta+\alpha}2 \right)
\end{cases}
$$
The points $E,F,G,H$ are concyclic iff
$$\det \begin{bmatrix}
OE^2 & E_x & E_y & 1\\
OF^2 & F_x & F_y & 1\\
OG^2 & G_x & G_y & 1\\
OH^2 & H_x & H_y & 1
\end{bmatrix}
= 0.
$$
In our case it takes the form:
$$\det \begin{bmatrix}
\tan(\frac{\pi}4-\frac{\alpha-\beta}4) & \cos\frac{\alpha+\beta}2 & \sin\frac{\alpha+\beta}2 & \cot(\frac{\pi}4-\frac{\alpha-\beta}4)\\
\tan(\frac{\pi}4-\frac{\beta-\gamma}4) & \cos\frac{\beta+\gamma}2 & \sin\frac{\beta+\gamma}2 & \cot(\frac{\pi}4-\frac{\beta-\alpha}4)\\
\tan(\frac{\pi}4-\frac{\gamma-\delta}4) & \cos\frac{\gamma+\delta}2 & \sin\frac{\gamma+\delta}2 & \cot(\frac{\pi}4-\frac{\gamma-\delta}4)\\
\tan(\frac{\pi}4-\frac{\alpha-\delta}4) & \cos\frac{\delta+\alpha}2 & \sin\frac{\delta+\alpha}2 & \cot(\frac{\pi}4-\frac{\alpha-\delta}4)
\end{bmatrix}
= 0,
$$
which can be verified routinely under the condition $(\star)$.
To verify the identity we can express everything in terms of $X:=e^{I\frac{\alpha-\beta}2}$, $Y:=e^{I\frac{\beta-\gamma}2}$, $Z:=e^{I\frac{\gamma-\delta}2}$, and $T:=e^{I\frac{\delta}2}$. In particular, we have $\sin\frac{\alpha-\beta}2 = \frac{X-X^{-1}}{2I}$, $\tan(\frac{\pi}4-\frac{\alpha-\beta}4) = \frac{(X-I)I}{X+I}$, $\cos\frac{\alpha+\beta}2 = \frac{X(YZT)^2 + X^{-1}(YZT)^{-2}}2$, and so on.
The following SageMath code verifies that the determinant as a rational function over variables $X,Y,Z,T$ reduces to $0$ w.r.t. the polynomial ideal defined by $(\star)$:
def row(t,z):
return [ (t-I)/(t+I)*I, (z + 1/z)/2, (z - 1/z)/2/I, (t+I)/(t-I)/I ]
R.<X,Y,Z,T> = PolynomialRing(QQ[I])
J = ideal( numerator( (X - 1/X) + (Z - 1/Z) - (Y - 1/Y) - (X*Y*Z - 1/(X*Y*Z)) ) )
M = matrix(Frac(R), 4, 4, [row(X,X*(Y*Z*T)^2), row(Y,Y*(Z*T)^2), row(Z,Z*T^2), row(X*Y*Z,X*Y*Z*T^2)] )
print( J.reduce( numerator(det(M)) ) )
The code prints "0" (run it online), thus establishing the determinant identity.
