Forums

Deconvolution by Fourier method

Started by aries44 June 22, 2005
In order to convolve two functions 'a' and 'b', we can take their Fourier
Transform(FT) and multiply them in Fourier domain i.e.
C= FT(a) * FT(b)
c = IFT(C)
and then Inverse Ft(IFT) of 'C' gives us the convolution of 'a' and 'b'.
Now if we want to deconvolve 'a' from 'c' to get 'b' we can do
B = FT(c)/FT(a)
b = IFT(B)

Now if 'a' is a square wave i face the problem of 'DIVIDE by ZERO', and
similar problem with many other shapes. Any solution for that? or any
better idea to perform convolution or deconvolution?
		
This message was sent using the Comp.DSP web interface on
www.DSPRelated.com
"aries44" <omar_shahid2@hotmail.com> wrote in message 
news:ZI2dnWySq-EHFiTfRVn-tA@giganews.com...
> > In order to convolve two functions 'a' and 'b', we can take their Fourier > Transform(FT) and multiply them in Fourier domain i.e. > C= FT(a) * FT(b) > c = IFT(C) > and then Inverse Ft(IFT) of 'C' gives us the convolution of 'a' and 'b'. > Now if we want to deconvolve 'a' from 'c' to get 'b' we can do > B = FT(c)/FT(a) > b = IFT(B) > > Now if 'a' is a square wave i face the problem of 'DIVIDE by ZERO', and > similar problem with many other shapes. Any solution for that? or any > better idea to perform convolution or deconvolution? > > This message was sent using the Comp.DSP web interface on > www.DSPRelated.com
If "a" is a square wave, do you know that it is? If you know that "a" has zeros in its FT then C must have those same zeros and you have a situation where you are dividing zero by zero - which seems like it might be handled. Otherwise, look at "system identification" and "deconvolution" I found this one with Google: http://www.ece.ucsb.edu/~hespanha/ece147c/web/identif.pdf Fred
aries44 wrote:
...
> Now if we want to deconvolve 'a' from 'c' to get 'b' we can do > B = FT(c)/FT(a) > b = IFT(B)
So you have two time domain signals, "a" and "c", and you want to find the LTI system "b" such that c = a * b (* denotes convolution).
> Now if 'a' is a square wave i face the problem of 'DIVIDE by ZERO', and > similar problem with many other shapes. Any solution for that? or any > better idea to perform convolution or deconvolution?
If you don't have control over the signal "a" (if you do, then you should Google for "system identification"), all you can do is specify the transfer function "B" at the postion of the harmonics of "a" - these are simply the amplitude and phase of the harmonics of the transformed square wave in "c" relative to "a". In other words, you can specify "B" only at a set of discrete frequencies f_0, f_1 = 2 f_0, ... (the number of harmonics will depend on your sampling rate). You end up with a bunch of pairs (f_k, B(f_k))_k, k = 0, 1, ..., and now you can fit the FIR or IIR filter according to your requirements. If you can control "a", then perhaps a Gaussian white noise signal might be appropriate. You can generate several "a"'s and measure several "c"'s and average the results together to get a more reliable estimate for the transfer function B. Regards, Andor
A way to get a working solution is to add a small amount of "noise" to
your FT(a) so that you no longer divide by zero

Also approximately, you can band-pass filter your data to stay away of
the zeroes in FT(a)


aries44 wrote:
> In order to convolve two functions 'a' and 'b', we can take their Fourier > Transform(FT) and multiply them in Fourier domain i.e. > C= FT(a) * FT(b) > c = IFT(C) > and then Inverse Ft(IFT) of 'C' gives us the convolution of 'a' and 'b'. > Now if we want to deconvolve 'a' from 'c' to get 'b' we can do > B = FT(c)/FT(a) > b = IFT(B) > > Now if 'a' is a square wave i face the problem of 'DIVIDE by ZERO', and > similar problem with many other shapes. Any solution for that? or any > better idea to perform convolution or deconvolution?
Congratulations! You have discovered the very reason why spectrum division is not used for signal deconvolution. Zeros on (or very close to) the init circle makes the whole thing break down. This is why there have been so much research done on inverse filters and system identification. Deconvolution is a can of worms that might keep you occupied for a long time to come. One really has to look very closely at eny particular problem, and use and exploit whatever idiosyncracities ther might be. Rune
Fred Marshall wrote:
> "aries44" <omar_shahid2@hotmail.com> wrote in message > news:ZI2dnWySq-EHFiTfRVn-tA@giganews.com... > >>In order to convolve two functions 'a' and 'b', we can take their Fourier >>Transform(FT) and multiply them in Fourier domain i.e. >>C= FT(a) * FT(b) >>c = IFT(C) >>and then Inverse Ft(IFT) of 'C' gives us the convolution of 'a' and 'b'. >>Now if we want to deconvolve 'a' from 'c' to get 'b' we can do >>B = FT(c)/FT(a) >>b = IFT(B) >> >>Now if 'a' is a square wave i face the problem of 'DIVIDE by ZERO', and >>similar problem with many other shapes. Any solution for that? or any >>better idea to perform convolution or deconvolution? >> >>This message was sent using the Comp.DSP web interface on >>www.DSPRelated.com > > > If "a" is a square wave, do you know that it is? > > If you know that "a" has zeros in its FT then C must have those same zeros > and you have a situation where you are dividing zero by zero - which seems > like it might be handled. > > Otherwise, look at "system identification" and "deconvolution" > I found this one with Google: > http://www.ece.ucsb.edu/~hespanha/ece147c/web/identif.pdf
There is also a good discussion starting on page 300 of S.W. Smith's "The Scientist and Engineer's Guide to Digital Signal Processing" available on line at http://www.dspguide.com in Chapter 17. Basically, deconvolution is an art. There is no simple formula that works whenever you turn the crank. If you ask for too much, you get nothing useful. You need to decide how much -- or how little -- you can settle for, then tweak the math to avoid the potholes without excessively blurring the result. All engineering is compromise. Deconvolution drives that truth home. Jerry -- Engineering is the art of making what you want from things you can get. &#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;