Deconvolution - signal processing

Gonzolo

Hi. I have a data gaussian g(t) and data that I suspect is a convolution (g*f). I want to find f(t), so I need to deconvolute (using Origin preferably).

If I convolute [tex]g*g[/tex], I get something very beautiful.
If I convolute [tex](g*f)*(g*f)[/tex], I get something very beautiful.
If I deconvolute [tex]g*^-^1g[/tex], I get something very beautiful.
If I deconvolute [tex](g*f)*^-^1(g*f)[/tex], I get something very beautiful.
If I convolute [tex]g*(g*f)[/tex], I get something acceptably beautiful.

But why is it that when I deconvolute [tex](g*f)*^-^1g[/tex], which should recover the f I'm looking for, I get the most horrible noise graph I've ever seen? Roughly a very thick (zigzag) useless straight line. How do I recover f with minimal noise?

Edited with Latex notation after first reply

rbj

may i suggest posting to comp.dsp (USENET) about this? you might have to straighten out your semantics because i know signal processing but i don't know what the hell you mean by "(g*f)*-1g". it's the "*-1" that makes me curious.

since we can really mark this up with real math using: [tex] x(t)=e^{j \omega t} [/tex] in this forum, why not use it so we can know what you're typing about. of course, you can use that form USENET, so you will need to express your math with ASCII somehow if you post to comp.dsp.

r b-j

Gonzolo

I understand that where * can be the symbol for convolution, [tex]*^-^1[/tex] can be the symbol for decovolution (-1 is meant to be an "exponent" to the *, like [tex]2^-^1 = 1/2[/tex]). I just Latex-ed it, thank.

http://www.kcci.virginia.edu/FRET-FL...ate/deconv.htm

My problem seems to simply be that when I deconvolute, I get so much noise, it buries the function I'm looking for.

rbj

even referring to your reference, i think the notation is useless. it is still convolution, but that the Fourier Transform of one of the elements of convolution has been inverted (as in reciprocal).

To say [tex] f = h (*) g [/tex] is the same as saying [tex] F = H G [/tex] where [tex] F [/tex], [tex] G [/tex], and [tex] H [/tex] are the Fourier Transform of [tex] f [/tex], [tex] g [/tex], and [tex] h [/tex]. To recover [tex] G [/tex] from [tex] F [/tex] and [tex] H [/tex], we divide: [tex] G = H^{-1} F [/tex] and then inverse Fourier Transform that to get [tex] g [/tex]. i guess you could say, with some abuse of notation that this is [tex] g = (h^{-1}) (*) f [/tex]. the notation you referred to make no sense.

now if [tex] H \approx 0 [/tex] anywhere, then the inverse (reciprocal) could get very nasty. it is not always possible to invert a "transfer function" (which is the role that i'm giving [tex] H [/tex]. if your "deconvolution" requires that, you got some basic mathematical problems.

r b-j

Gonzolo

The notation is an abuse, but probably the next best thing without going into the Fourier transforms. Anyhow.

Not sure what you mean with [tex]H \approx 0[/tex]. All the functions invloved surely become 0 when they tend to either infinity. My raw data (f*g and g) is really not that noisy, only one hump each, same order of magnitude. Heck, here it is (f*g left, g right), any time scale is good :

0 0
0 0.00339
-0.03315 0
0.04222 0.00339
0.12988 0
0.33524 0.00339
0.63104 0.01017
0.81368 0.02373
0.93687 0.06441
0.99383 0.14915
1 0.34915
0.9216 0.64407
0.70303 0.92542
0.54662 1
0.46232 0.92542
0.40534 0.64407
0.35129 0.34915
0.29124 0.14915
0.25822 0.06441
0.22819 0.02373
0.19514 0.01017
0.18916 0.00339
0.15315 0
0.14712 0.00339
0.12906 0
0.1231 0.00339
0.09908 0
0.09908
0.08407
0.08407
0.08105
0.07804
0.06305
0.0601
0.06322
0.06303
0.06005
0.04513
0.04804
0.05107
0.04199
0.04508
0.048
0.03891
0.04206
0.04502
0.03903
0.03594
0.03622
0.03605
0.02721
0.03305
0.02705
0.03301
0.03301
0.02407
0.03
0.027
0.02709
0.03005
0.02404
0.021
0.02404
0.02109
0.02398
0.02407
0.02712
0
0

rbj

well, i think you make it harder for yourself by not "going into the Fourier transforms."

[tex] H [/tex] is the Fourier Transform of the "PSF" function in the reference you have given:

when they say
i say
[tex] f = h (*) g [/tex]
and
[tex] F = H G [/tex].

when they say
i say
[tex] g = (h^{-1}) (*) f [/tex]
and
[tex] G = H^{-1} F [/tex].

inverting (computing the reciprocal of) [tex] H [/tex] to get [tex] H^{-1} [/tex] is nasty if [tex] H \approx 0 [/tex], even worse if [tex] H=0 [/tex].

well, if the Fourier Transform of the function you're trying to deconvolute is going to zero anywhere, you got a fundamental mathematical problem.

if you have f(*)g and g, and what you want is f, then you have to F.T. both f(*)g and g, divide the F.T. of f(*)g with the F.T. of g, take that result and inverse F.T. got that?

r b-j

Gonzolo

Yeah, I got all of that, but I understand the software I'm using is supposed to take care of the calculation. All I should have to do is click on the data columns and click the "deconvolute" command. When I do that, I don't get anything useful. I roughly get 1, -1, 1.05, -1, 1, -1, 1.... over time, with only a hint of the f I'm looking for (the 1.05)

I understand the transforms. But I didn't know H=0 was problematic. Gaussians and their transforms go to 0, yet they easily convolute and deconvolute, don't they? I can't imagine a useful function that doesn't go to 0. Wouldn't the transform of any pulse, or signal do so?

rbj

you may have heard the adage/acronym: GIGO (garbage-in, garbage-out). your software may be working fine, but when the original convolution (that you're trying to undo) was done, some information was lost (if [tex] H \approx 0 [/tex]). when you try to reverse that loss of information, you are amplifying what was left (after the loss), but you're amplifying it with an obscene gain ([tex] H^{-1} \approx \infty [/tex]).

your gain at high frequencies (around Nyquist [tex] f \approx F_s/2 [/tex]), got very large, probably because [tex] H(F_s/2) \approx 0 [/tex].

some transfer functions for discrete-time linear systems (which is what you're dealing with) settle down to a nonzero constant as [tex] f = F_s/2 [/tex] . those transfer functions invert just fine. the ones that go to zero anywhere, do not invert just fine. you have to fudge them somehow, if you want to invert them.

r b-j

Gonzolo

Ok, I follow. Interesting stuff. So the best way to do this would be to have a higher sampling rate during the measurement, right?

And if that's not possible, I could interpolate between the points I have.

rbj

i dunno how a higher sampling rate or interpolating can help given what you have previously said. you say that both functions are a gaussian pulse, and we both understand that a gaussian pulse is an "eigenfunction" under the Fourier transform operator. that is, the F.T. of [tex] e^{- \pi t^2} [/tex]
is [tex] e^{- \pi f^2} [/tex]. i'm not so sure how that will look with discrete-time functions but it should be very similar if the sampling rate is fast enough and if the sample size is large enough. even though the gaussian pulse never gets to zero, it gets awful damn close in a hurry. so you end up dividing by something close to zero.

perhaps adding a small DC value to the F.T. of your g function will do it. that would be adding a small spike at the g(0) term before transformation.

i dunno.

r b-j

Gonzolo

One function is a gaussian pulse. The one I'm looking should fit with a decreasing exponential [tex]h = ae^-^b^t[/tex]. (more specifically a sum of N decreasing exponentials, such that the greater the N the better the precision, but N=1 is fine for now). Same order of duration as the gaussian. The data I have is for the gaussian, and for the convolution of the two.

The FFT of my discrete gaussian does look like another gaussian (I verified that it would get smoother with a better sampling rate).