title: Comparing Naive Inverse Filter to Wiener Filter for Deconvolution in Matlab
body: <p>I am currently trying to compare a simple inverse filter to the wiener filter for deconvolution using matlab.  My starting signal is <code>exp(-t^2)</code> and this is to be convolved with a rect that is nonzero for times -.5 to .5.  I am introducing noise with amplitude in the range -.5 to .5.</p>
  
  <p>Defining my time domain to frequency domain mapping:</p>
  
  <pre><code>f = exp(-t^2) =&gt; F
  
  s = rect =&gt; R
  
  c = f*s =&gt; C
  
  r = noise (see above) =&gt; R
  
  with noise c becomes: c = f*s + n =&gt; C = FxS + N
  </code></pre>
  
  <p>For the first approach I am simply taking the FT of <code>c</code> and dividing it by the FT of <code>f</code> and then doing an inverse FT. This is amounts to <code>s = (approx.) ifft((FxS + N)/F)</code>.</p>
  
  <p>For the second approach I am taking the wiener filter, <code>W</code>,  and multiplying it against <code>C/R</code> and then doing an inverse FT.  This amounts to <code>S = (approx.) ifft(CxW/R)</code>.  </p>
  
  <p>The wiener filter is <code>W = mag_squared(FxS)/(mag_squared(FxS) + mag_squared(N))</code>.</p>
  
  <p>I have used '*' to mean convolution and 'x' to mean multiplication.</p>
  
  <p>I am trying to compare the two deconvolutions of the rect over the time interval -3 to 3.
  Right now, my resulting graphs of the deconvolved rect look nothing like the original.<br>
  Could someone point me in the right direction as to what I'm doing to wrong?  I have tried using ifftshift and different scalings in many different orderings but nothing seems to work.  </p>
  
  <p>Thanks</p>
  
  <p>My matlab code is below:</p>
  
  <pre><code>%%using simple inverse filter
  dt = 1/1000;
  t = linspace(-3,3,1/dt); %time
  s = zeros(1,length(t)); 
  s(t&gt;=-0.5 &amp; t&lt;=0.5) = 1; %rect
  f = exp(-(t.^2)); %function
  r = -.5 + rand(1,length(t)); %noise
  
  S = fft(s);
  F = fft(f);
  R = fft(r);
  C = F.*S + R;
  S_temp = C./F;
  s_recovered_1 = real(ifft(S_temp));  %correct?...works for signal without R (noise)
  
  figure();
  plot(t,s + r);
  title('rect plus noise');
  
  figure();
  hold on;
  plot(t,s,'r');
  plot(t,f,'b');
  legend('rect input','function');
  title('inpute rect and exponential functions');
  hold off;
  
  figure();
  plot(t,s_recovered_1,'black');
  legend('recovered rect');
  title('recovered rect using naive filter');
  
  
  %% using wiener filter
  N = length(s);
  I_mag = abs(I).^2;
  R_mag = abs(R).^2;
  W = I_mag./(I_mag + R_mag);
  S_temp = (C.*W)./F;
  s_recovered_2 = abs(ifft(S_temp));  
  
  figure();
  freq = -fs/2:fs/N:fs/2 - fs/N;
  hold on;
  plot(freq,10*log10(I_mag),'r');
  plot(freq,10*log10(R_mag),'b');
  grid on
  legend('I_mag','R_mag');
  title('Periodogram Using FFT')
  xlabel('Frequency (Hz)')
  ylabel('Power/Frequency (dB/Hz)')
  
  figure();
  plot(t,s_recovered_2);
  legend('recovered rect');
  title('recovered rect using wiener filter');
  </code></pre>