Re: Convolve and Deconvolve not working right??

[Robert A. Macy]

Thank you for your suggestions and the paper, too.

I tried putting in all kinds of nonzero values, even as
large as .001 versus max of 1 and the thing still blows up


I tried a series of things that don't work.

My first idea was to convolve the guassiancurve [called
gt]with some 'magical' function should then recreate the
dirac function, or impluse function [called dt].

Now, how to find this magical function:

the idea is that rt [magical function] convolved with gt
equals dt.

since convolution is identical to multiplying in the
frequency domain..

Find the Fourier Transform for both gt and dt:
then rf times gf equals df.  

dividing through finds rf = df/gf

Well I first tried this with the offset where the peaks
were at 51:

for gaussian curve...
>> x=(-50:50);
>> gt=exp( -(x.*x/200) );
which makes a nice shape
>> gf=fft(gt);

and for dirac
>> dt=zeros(1,101);dt(51)=1;
>> df=fft(dt);

then the spectrum of my magical form should be
>> rf=df./gf;rt=ifft(rf);
note rt comes back only real terms as expected
inspection looks about right.

but I think sum(rt) should equal 1, but is always .03...
so that's not right.  

next 
remember the step function that was convolved with gt?
well,
>> out=conv(rt,stept);
does not recreate the step function at all!

ok, well may be it didn't work because time equations
should be centered in time

so I tried half the gt and half the dirac:

interestingly, the df was 1 for all 101 as it should be
and the spectrum of fft(gt) looked about right.

yet this approach [with much fudging] yielded similar
results to the first approach, just different.

I think what Im missing is EXACTLY what the effect of the
fft over a limited range is doing to things.

Need some help here.

            Robert

On Sun, 25 Mar 2007 23:28:15 +0200
 Miroslaw Kwasniak <address@hidden> wrote:
> 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