Re: Convolve and Deconvolve not working right??

[Fredrik Lingvall]

The problem is that even when you only have the smallest error
in  your  "out" signal a direct inverse (or an LS estimate) may result in
an unstable solution. This is often handled by introducing an uncertainty
parameter (eg., a noise term) and then minimize the mean squared error
which results in the Wiener filter. Essentially, the convolution operation
destroys information that never can be recovered perfectly, but the Wiener
filter does the best that you do (when using a linear method).

/Fredrik


Robert A. Macy wrote:

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
    

_______________________________________________
Help-octave mailing list
address@hidden
https://www.cae.wisc.edu/mailman/listinfo/help-octave
  

-- 
***********************************************************************
Fredrik Lingvall                      
Postdoctoral Research Fellow          address@hidden
                                      Phone: +47 22 84 00 68
                                      Fax:   +47 22 85 24 01
Department of Informatics
Group for Digital Signal Processing
and Image Analysis
Gaustadalléen 23
P.O.Box  1080 Blindern
NO-0316 Oslo
Norway
***********************************************************************