I have some function $I(x)$ which I want to reconstruct as an infinite sum of 'simpler' basis functions $d(x)$. Each of these basis functions will have an (unknown) weighting $W(x)$ which I seek to find.
If my working is correct, I believe this reconstruction results in a convolution integral. As a result $W(x)$ can be found using the Fourier and inverse Fourier transform.
However when I test this numerically, it's wrong. Where am I making my mistake? Thanks.
Trying to use the same function shifted by various amounts as a basis is generally not a good idea. When you convolve a Gaussian with some $W(x)$, you get "$W(x)$ with Gaussian blur applied", so you cannot represent functions that have sharper features than the Gaussian width (at least if you restrict yourself to "normal'' functions, not derivatives of the Dirac delta).
Besides, the Fourier transform of $I(x) = x^2$ is the second derivative of the delta function—definitely not good for numerical evaluation.
Although in your particular case, when $I(x)$ does not have sharp features, you still can get a sensible answer analytically. If you insist on numerics, you have to make $I(x)$ smoothly go to zero at large $|x|$ (e. g., multiply by a wide Gaussian) to avoid integration up to infinity.