Another Approach

I'm assuming here that $A$ is roughly equal to one. If it's a lot bigger or a lot smaller, then the problem does change quite a bit.

A first estimate to getting $x(t)$ would simply to take the mean of the observations

$$x_1(t) = \frac{1}{N}\sum_{n=0}^{N-1}s_n(t)$$

All the $x_0(t)$ terms should add coherently whereas the time and phase shifted terms will be most uncorrelated.

Next we take a look at the transfer functions (I'll use capital letters for frequency domain)

$$S_n(\omega) = X(\omega)\cdot (1 + A\cdot e^{j(\theta-\omega)\tau_n})$$

So transfer function between $S$ and $X$ basically a complex comb filter that has maxima at $(\theta-\omega)\tau_n = 2k\pi$. If $A< 1 that it also has minima at $(\theta-\omega)\tau_n = (2k+1)\pi$.

Now you can use your first estimate $x_1(t)$ to estimate the transfer function,

$$H_{n,1} = \frac{S_n(\omega)}{X_1(\omega)}$$

This should roughly look like a comb filter. You can use the location of the maxima and maybe the maxima to estimate $\tau_{n,1}$ and the height of the maxima to estimate $A_{n,1}$.

Now you use these estimates to refine you estimate for $x(t)$

$$x_2(t) = \frac{1}{N}\sum_{n=0}^{N-1}(s_n(t)-A_{n,1}x_1(t-\tau_{n,1})e^{j\theta \tau_{n,1}})$$

Rinse and repeat until the residual error becomes small or it stops converging.

Another Approach

I'm assuming here that $A$ is roughly equal to one. If it's a lot bigger or a lot smaller, then the problem does change quite a bit.

A first estimate to getting $x(t)$ would simply to take the mean of the observations

$$x_1(t) = \frac{1}{N}\sum_{n=0}^{N-1}s_n(t)$$

All the $x_0(t)$ terms should add coherently whereas the time and phase shifted terms will be most uncorrelated.

Next we take a look at the transfer functions (I'll use capital letters for frequency domain)

$$S_n(\omega) = X(\omega)\cdot (1 + A\cdot e^{j(\theta-\omega)\tau_n})$$

So transfer function between $S$ and $X$ basically a complex comb filter that has maxima at $(\theta-\omega)\tau_n = 2k\pi$. If $A< 1$ that it also has minima at $(\theta-\omega)\tau_n = (2k+1)\pi$.

Now you can use your first estimate $x_1(t)$ to estimate the transfer function,

$$H_{n,1} = \frac{S_n(\omega)}{X_1(\omega)}$$

This should roughly look like a comb filter. You can use the location of the maxima and maybe the maxima to estimate $\tau_{n,1}$ and the height of the maxima to estimate $A_{n,1}$.

Now you use these estimates to refine you estimate for $x(t)$

$$x_2(t) = \frac{1}{N}\sum_{n=0}^{N-1}(s_n(t)-A_{n,1}x_1(t-\tau_{n,1})e^{j\theta \tau_{n,1}})$$

Rinse and repeat until the residual error becomes small or it stops converging.

Another Approach

I'm assuming here that $A$ is roughly equal to one. If it's a lot bigger or a lot smaller, then the problem does change quite a bit.

A first estimate to getting $x(t)$ would simply to take the mean of the observations

$$x_1(t) = \frac{1}{N}\sum_{n=0}^{N-1}s_n(t)$$

All the $x_0(t)$ terms should add coherently whereas the time and phase shifted terms will be most uncorrelated.

Next we take a look at the transfer functions (I'll use capital letters for frequency domain)

$$H_n(\omega) = X(\omega)\cdot (1 + A\cdot e^{j(\theta-\omega)\tau_n})$$

That's basically a complex comb filter that has maxima at $(\theta-\omega)\tau_n = 2k\pi$ and minima at $(\theta-\omega)\tau_n = (2k+1)\pi$

Now you can use your first estimate $x_1(t)$ to normalize the individual spectra,

$$S_{n,1,norm} = \frac{S_n(\omega)}{X_1(\omega)}$$

This should roughly look like a comb filter. You use the location of the minima and maxima to estimate $\tau_{n,1}$ and the height of the maxima (or depth of the minima) to estimate $A_{n,1}$.

Now you use these estimates to refine you estimate for $x(t)$

$$x_2(t) = \frac{1}{N}\sum_{n=0}^{N-1}(s_n(t)-A_{n,1}x_1(t-\tau_{n,1})e^{j\theta \tau_{n,1}})$$

Rinse and repeat until the residual error becomes small or it stops converging.

Another Approach

I'm assuming here that $A$ is roughly equal to one. If it's a lot bigger or a lot smaller, then the problem does change quite a bit.

A first estimate to getting $x(t)$ would simply to take the mean of the observations

$$x_1(t) = \frac{1}{N}\sum_{n=0}^{N-1}s_n(t)$$

All the $x_0(t)$ terms should add coherently whereas the time and phase shifted terms will be most uncorrelated.

Next we take a look at the transfer functions (I'll use capital letters for frequency domain)

$$S_n(\omega) = X(\omega)\cdot (1 + A\cdot e^{j(\theta-\omega)\tau_n})$$

So transfer function between $S$ and $X$ basically a complex comb filter that has maxima at $(\theta-\omega)\tau_n = 2k\pi$. If $A< 1 that it also has minima at $(\theta-\omega)\tau_n = (2k+1)\pi$.

Now you can use your first estimate $x_1(t)$ to estimate the transfer function,

$$H_{n,1} = \frac{S_n(\omega)}{X_1(\omega)}$$

This should roughly look like a comb filter. You can use the location of the maxima and maybe the maxima to estimate $\tau_{n,1}$ and the height of the maxima to estimate $A_{n,1}$.

Now you use these estimates to refine you estimate for $x(t)$

$$x_2(t) = \frac{1}{N}\sum_{n=0}^{N-1}(s_n(t)-A_{n,1}x_1(t-\tau_{n,1})e^{j\theta \tau_{n,1}})$$

Rinse and repeat until the residual error becomes small or it stops converging.