[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: Convolve and Deconvolve not working right??
From: |
Robert A. Macy |
Subject: |
Re: Convolve and Deconvolve not working right?? |
Date: |
Sun, 25 Mar 2007 21:58:41 -0700 |
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