I have an output signal $y$ which is an input signal $x$ convolved $\star$ with an impulse response function $h$ with some added noise $n$ :

$$y(t) = h(t) \star x(t) + n(t)$$

I know the input signal $x$ and output signal $y$ and would like to calculate $h$ the impulse response function. I found that deconvolution is not as straightforward as convolution because the input signal contains zeros and then division in the frequency domain would not be defined. Looking around the internet for ways to "deconvolve" if found two methods: Wiener deconvolution and regularized deconvolution. The Wiener deconvolution seemed easier to understand so I wanted to try and implement it in Matlab (the Matlab function deconv gives me errors about the input signal having a zero at the first entry and if I read the help file it only seems to work correctly for polynomials?).

So per the [wikipedia][1] explanation you want to find $g$ so that:

$$\hat{x}(t) = g(t) \star y(t)$$

But then in the definition of how to calculate $G$ they use all variables in the original equation. Also they only show how to find $x$ but probably $x$ and $h$ can be exchanged because convolution is commutative, but I'm unsure about the correct length of both vectors. Currently they are the same length.

$$G(f) = \frac{H^*(f) S(f)}{|H(f)|^2 S(f) + N(f)}$$

where:

 - $H = {\tt fft}(h)$
 - $G = {\tt fft}(g)$
 - $S =$ power spectral density of $x$ ? is this ${\tt fft}(x)$?
 - $N =  $mean power spectral density of $n$, don't really understand what this is

My question is how do I get the impulse response function $h$ without already knowing it (as it is both in the definition of $G$ and in the original equation)? Since I know both input and output it should not be very different from finding the input with a known impulse response function. I want to do this in Matlab. 

  [1]: http://en.wikipedia.org/wiki/Wiener_deconvolution
  [2]: https://i.sstatic.net/ScdWy.png