Re: Convolve and Deconvolve not working right??

[Miroslaw Kwasniak]

On Sat, Mar 24, 2007 at 06:51:53PM -0700, Robert A. Macy wrote:
> Anybody here help with this?
> 
> I tried using conv and deconv to try something, but it does
> not render the expected solution.
> 
> If I use:
> 
> >> stept=zeros(1,501);stept(1:251)=1;
> then smooth it with a guassian shaped curve
> >> smoothedstept=conv(sept,guassiancurve);
> 
> now recover the step by trying to deconvolve doesn't work
> >> out=deconv(smoothedstept,guassiancurve);
> produces an error that the first term is supposed to be
> non-zero.

Octave uses deconv(x,w) == filter(x,w,1). This method is algebraically OK
but not with real floating point computing :(

Errors are in both stages: conv and deconv, but second (inverse filtering)
suffers much more when filter is unstable -> roots of the window are outside
an unit circle.

As I see in case of an gaussian window w=gausswin(n,a) the value
rmax=max(abs(roots(w))) is always above 1. When (rmax-1) is below some small
value (like 1e-14) the deconv is usable IMHO. For narrow windows you have
problems:

n=50
a=2 -> (rmax-1) = 1.3323e-15
a=3 -> (rmax-1) = 0.13907

For some cases I did with succes the deconvolution with this recipe
(deconvolution via FFT with setting values of an window no less then E and
padding to the length of signal):

    E=.1e-4; A=3; % tuning parameters - as I see min(w1) should be above 1e-22

    lw=length(w);
    lx=length(xs);
    w1=linspace(0,max(min(w),E),lx-lw+2);
    w1=[ max(E,w'), w1(end-1:-1:2).^A ];
    xd=real(ifft(fft(xs)./fft(w1))));

You can deconvolute in LS sense with:

   xd=pinv(convmtx(w,length(xs)))*xs;

this gives good results but CPU time and memory are eaten.

About other methods you can read in
E. Pantin, J.-L. Starck, and F. Murtagh,
"Deconvolution and Blind Deconvolution in Astronomy"
http://jstarck.free.fr/Blind07.pdf