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]

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

If I convolute g*g, I get something very beautiful.
If I convolute (g*f)*(g*f), I get something very beautiful.
If I deconvolute g*-1g, I get something very beautiful.
If I deconvolute (g*f)*-1(g*f), I get something very beautiful.
If I convolute g*(g*f), I get something acceptably beautiful.

But why is it that when I deconvolute (g*f)*-1g, 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?

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-FLIM/steadystate/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]

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-FLIM/steadystate/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.

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]

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

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

Not sure what you mean with [tex]H \approx 0[/tex].

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

Deconvolution is the procedure which reverses the image degradation due to convolution (i.e., inverse Fourier transform.) By modeling the microscope optics as a linear and shift-invariant system, the PSF can be used to describe the transformation of any image by the microscope. A typical fluorescence microscope image can be

[measured image] = [PSF] * [desired image]
where '*' symbol represents the mathematical operation of convolution
[and PSF is the pointspread function of the imaging system]. The deconvolution (*-1), the mathematical inverse of convolution can be represented as

[desired image] = [measured image] (*-1) [PSF]
The goal of deconvolution is to solve the equation for the desired image.

when they say

[measured image] = [PSF] * [desired image]

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

when they say

[desired image] = [measured image] (*-1) [PSF]

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].

All the functions invloved surely become 0 when they tend to either infinity.

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

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 :

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]

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

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]

Yeah, I got all of that, but I understand the software I'm using is supposed to take care of the calculation.

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]).

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)

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].

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?

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]

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.

i don't know 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).