Re: Convolve and Deconvolve not working right??

[Robert A. Macy]

Chris,

Thank you for response and the paper

quoting from Chris:
> the convolution theorem for DFT/FFT is 'circular'
> or 'periodic' convolution; so zero-padding makes the
> two types of convolution equivalent.  Perhaps the
> single most important thing to keep in mind about
> DFT is that it is a periodized approximation to FT.

THIS IS EXACTLY what I thought was going wrong here.  I'm
now trying to figure out how to remove this effect. 

I think I got closer when I 'centered' the function and did
not use an offset.  where gt(1)=1;and dt(1)=1, elsewhere =
0

The fft of dt then gave me 1 over all ranges as expected,
thought I was on the right track.  

>> gf=fft(gt);  %where gt is centered at the origin
then 
>> rf=1./gf;rt_half=ifft(rf);
to find rt whole I shifted it out
>> rt(51:101)=rf_half;rt(51:-1:1)=rf_half;
or should that have been conjugate?, maybe not

but that produced a function that the sum(rt)=.89, which is
a lot closer to 1, however, it inverted the signal ??

 Robert 

On Mon, 26 Mar 2007 08:15:11 -0400
 "Chris Zarowski" <address@hidden> wrote:
> Hi again.  Comments ...
> 
> The Fourier transform of delta(t) (Dirac) is unity
> (depending on the definition of FT - there are variants).
> So r(t) * g(t) = delta(t) ('*' is linear convolution)
> yields
>    R(f) = 1/G(f)
> in the Fourier domain. (g(t) = Gaussian, r(t) = 'magic.')
> 
> The Fourier transform (FT) of a Gaussian pulse is
> a Gaussian pulse (the Gaussian pulse is an eigenfunction
> of the Fourier transform operator). So, if G(f) is
> small (the tails making trouble) then R(f) will blow up. 
> I have attached notes I wrote over 10 years ago
> on wavelets and wavelet packets. Lemma 5.7 on p. 77
> derives the FT of a Gaussian pulse. This result
> may be found elsewhere (e.g., math tables).
> 
> An ad hoc fix for the small tail problem
> (crude attempt at regularization) is to use
>   R(f) = 1/(G(f) + epsilon)  , epsilon is 'small'
> 
> When approximating the FT with discrete Fourier
> transform (DFT) (as implemented with FFT) keep
> in mind that if you do not sample the underlying
> analog signals rapidly enough there will be aliasing
> errors.  
> Also, wrap-around error can occur in FFT-based
> convolution if there is not adequate zero-padding.
> Recall that LTI filtering is 'linear' convolution,
> whereas
> the convolution theorem for DFT/FFT is 'circular'
> or 'periodic' convolution; so zero-padding makes the
> two types of convolution equivalent.  Perhaps the
> single most important thing to keep in mind about
> DFT is that it is a periodized approximation to FT.
> 
> If you are trying to approximate the Dirac generalized
> function keep in mind it has no amplitude, but only
> a strength (its area).  The zero-mean Gaussian pulse
> will approximate Dirac in the limit as variance goes to
> zero. However, the pulse height increases without
> bound as variance decreases.
>   
> If you use a square pulse then cutting the
> width of it in half means you must double the height
> to keep the area (strength) a constant. Maybe this
> explains sum(rt) not being what you expect.
>  
> Christopher J. Zarowski, Ph. D., P. Eng. (Ontario)
>      IEEE Senior Member
>      NanoDotTek
>      http://www.nanodottek.net
> 
> ----- Original Message ----- From: "Robert A. Macy"
> <address@hidden>
> To: "Miroslaw Kwasniak" <address@hidden>
> Cc: <address@hidden>
> Sent: Monday, March 26, 2007 12:58 AM
> Subject: Re: Convolve and Deconvolve not working right??
> 
> 
> > 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
> >