I want to estimate the parameters $\mu_i$ and $\sigma^2_i $ of a countable mixture of Gaussians with assumed equal weights, variance and identically spaced means. I intially thought that the Fourier transforms of an iid sample from the mixture mentioned would give above another Gaussian in phase space, but when computed in matlab I get a sharp peak spectral distribution. I used a Gaussian windows function, so I don't understand why the power spectral density tends to be unbounded at zero phase.

If the Gaussians are constrained to be of equal shapes and distances, then you have a scaled theta function or the convolution of a Gaussian with a sha (Ш) function or scaled Dirac comb. Ш is a sum of equally spaced delta functions, and the Fourier transform of Ш is another Ш function, usually with different spacings and amplitudes depending on your conventions. Since the Fourier transform of a Gaussian is a Gaussian, the Fourier transform of your function is a sum of equally spaced delta functions whose amplitudes sample a Gaussian density, something like $$\beta \sum_{n=\infty}^\infty e^{\pi^2 z^2/\sigma  2 \pi i x_0 z}\delta(zn/\alpha)$$ where the parameters $\alpha, \beta, \sigma,$ and $x_0$ are determined by the parameters of your distribution and your Fourier conventions. In particular, there is a delta function at $0$. You should expect this to happen when you have a periodic function whose average over a period isn't $0$. You can read off the spacing of the Gaussians from the reciprocal of the spacing between the delta functions of the Fourier transform, and the way the amplitude drops off determines the shape of the Gaussian. That said, finding the parameters of a finite Gaussian mixture is an extremely common problem when the components are not contrained to have the same shapes and spacings. This is because Gaussian mixtures are one of the standard models used in machine learning. Gaussians are localized, so there isn't a big difference between the sum of the Guassians with nearby centers and the whole function. So, in practice, if you have (noisy) values from your function then instead of taking the Fourier transform, you may want to use some of these tools, and then project the best fit onto the space of scaled theta functions. The EM (expectation maximization) algorithm is one common method, and you should be able to find many Matlab implementations. 

