Generalization of Wigner overlap formula

I want to generalize the Wigner overlap formula,
$Tr( F G ) = 2 pi int_{-infty}^{infty} dq int_{-infty}^{infty} dq W_F(q,p) W_G(q,p)$,
where $W_F(q,p)$ and $W_G(q,p)$ are the Wigner functions of the operators $F$ and $G$, respectively.

This formula is stated in literature for two operators $F,G$ (see e.g. Measuring the quantum states of light by Ulf Leonhardt) and some sources state that it is easy to generalize it to more than two operators.

My approach was to proof the statement for two operators and try to adapt the proof for three operators, hoping to find some pattern.

Thus, first my approach for two operators $F$ and $G$:
begin{align}
Tr[FG] &= int_{-infty}^{infty} dq_1 langle q_1 | F G |q_1 rangle = int_{-infty}^{infty}dq_1 int_{-infty}^{infty}dq_2 langle q_1 | F | q_2 rangle langle q_2 |G |q_1 rangle
&= int_{-infty}^{infty}dq int_{-infty}^{infty}dx_1 langle q- frac{x_1}{2} | F | q + frac{x_1}{2} rangle langle q + frac{x_1}{2} |G |q -frac{x_1}{2} rangle
&= int_{-infty}^{infty}dq int_{-infty}^{infty}dx_1 int_{-infty}^{infty}dx_2 langle q- frac{x_1}{2} | F | q + frac{x_1}{2} rangle langle q – frac{x_2}{2} |G |q +frac{x_2}{2} rangle delta(x_1+x_2)
&= int_{-infty}^{infty}dq int_{-infty}^{infty}dx_1 int_{-infty}^{infty}dx_2 int_{-infty}^{infty}dp frac{1}{2pi} e^{ip(x_1+x_2)} langle q- frac{x_1}{2} | F | q + frac{x_1}{2} rangle langle q – frac{x_2}{2} |G |q +frac{x_2}{2} rangle
&=2 piint_{-infty}^{infty}dq int_{-infty}^{infty}dp int_{-infty}^{infty}dx_1 frac{1}{2pi} e^{ipx_1} langle q- frac{x_1}{2} | F | q + frac{x_1}{2} rangle frac{1}{2pi} int_{-infty}^{infty}dx_2 e^{ip x_2}langle q – frac{x_2}{2} |G |q +frac{x_2}{2} rangle
&= 2 pi int_{-infty}^{infty} dq int_{-infty}^{infty} W_F(q,p) W_G(q,p)
end{align}

Now I try to generalize this as similar as possible for three operators $F, G, H$:
begin{align}
Tr[FGH] &= int_{-infty}^{infty} dq_1 langle q_1|FGH |q_1rangle = int_{-infty}^{infty} dq_1 int_{-infty}^{infty} dq_2 int_{-infty}^{infty} dq_3 langle q_1|F|q_2 rangle langle q_2|G|q_3 rangle langle q_3| H |q_1rangle
&= int_{-infty}^{infty} dq int_{-infty}^{infty} dx_1 int_{-infty}^{infty} dx_3 langle q – frac{x_1}{2}|F|q + frac{x_1}{2} rangle langle q + frac{x_1}{2} |G|q – frac{x_3}{2} rangle langle q – frac{x_3}{2} | H |q – frac{x_1}{2}rangle
&= …
end{align}

Well, one observes that the trick from the proof for two operators doesn’t work here, because if I choose $x_2$ to be $-x_1$ for the middle part, we require $x_3 = – x_2$, hence $x_3 = x_1$. So, the last part would have the form $|q-frac{x_1}{2} | H |q – frac{x_1}{2}rangle$ and not those we require to proceed.

Is there anything (stupid?) that I oversee?
Does anyone have an idea how one can prove the formula for three operators?

Thank you in advance for your help!