Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I've been trying to teach myself about deconvolution through fourier transforms, but I seem to be missing something simple, as my results are garbage.

I start by defining a test function and a window which I will convolve with it

Window[width_, x_] = UnitStep[x + width/2] UnitStep[width/2 - x];
Test[x_] = UnitStep[x] UnitStep[count - x] Sin[0.1 x]^2;

I then sample the convolution across a set of data points

count=100;

MakeData[width_] = 
  Simplify[Table[
    Integrate[Window[width, x - i] Test[x], {x, -Infinity, Infinity}, 
     Assumptions -> {count > width > 0}], {i, 0, count}]];

tempData = MakeData[20];

I make a similar list for my sample window

tempWindow = Table[Window[20, x], {x, -count/2, count/2}];

I use the convolution theorem to deconvolve my window from my data

deconvolved = InverseFourier[Fourier[tempData]/Fourier[tempWindow]];

I was expect that this will return the original function

expected = Table[Test[x], {x, 0, count}];

However, the results are completely different.

Graph of Results

There's obviously something I'm doing wrong in my convolution or deconvolution, but I don't see what it is.

share|improve this question
1  
For starters, you're generating the data by sampling the theoretical result of the convolution but trying to deconvolve some other sampled window... Deconvolutions (noisy) are generally unstable. See this message of mine for an example demonstrating the convolution theorem (a and b are random vectors). –  rm -rf Jan 30 at 0:36
1  
I'm not sure if this question is related to MMA. Anyway, when your window has zeros, as usually happens, the result is not unique and unstable. You can try making your smallest values of Fourier@tempData be really zero, and add some small constant value to the denominator to get zero out of those 0/0. Furthermore, your convolution formula only works because it is even. Otherwise it should say Window[width, x - i]. Furthermore, given that you didn't rotate your window samples, and sampled from -count/2, your expected should be rotated right count/2 –  Rojo Jan 30 at 1:19

1 Answer 1

up vote 7 down vote accepted

While it's tempting to attribute the errors you're observing to floating-point errors due to zeros in the DFT of the window, this is actually not the case here:

Window[width_, x_] := UnitStep[x + width/2] UnitStep[width/2 - x];
Test[x_] := UnitStep[x] UnitStep[count - x] Sin[1/10 x]^2;
tempWindow = Table[Window[20, x], {x, -count/2, count/2}];
{Min[Abs[Fourier[tempWindow]]],Max[Abs[Fourier[tempWindow]]]}

(*Out: {0.00455694, 2.08958}*)

Instead, the primary issue is a misuse of the convolution theorem. The inverse DFT of the product of the DFT's of the samples of two functions is NOT the samples of the convolution of the two functions; rather, it is the cyclic convolution of the samples of the two functions, which is an important distinction.

Your code samples the convolution of the two functions (I modified it to be a bit cleaner):

count = 100;
expr=Integrate[Window[width, x - i] Test[x], {x, -Infinity, Infinity}];
tempData = N@Table[expr/.width->20,{i,0,count}];

It then samples the window function, and then attempts to deconvolve the sampled window from the samples of the convolution of the two functions:

deconvolved = Abs[InverseFourier[Fourier[tempData]/Fourier[tempWindow]]];

The result looks like nonsense because tempData is not actually a cyclic convolution of tempWindow with the samples of Test, so you have no reason to expect that the cyclic deconvolution of tempData with the samples of Window will give you the samples of Test. It's a bit of a subtle error, but it's actually very important.

For reference, to demonstrate the cyclic convolution and deconvolution, try this:

tempTest = Table[Test[x], {x, -count/2, count/2}];
tempConv = Re@InverseFourier[Fourier[tempTest] Fourier[RotateLeft[tempWindow,count/2]]];
tempDeconv = Re@InverseFourier[Fourier[tempConv]/Fourier[RotateLeft[tempWindow, count/2]]];
ListPlot[tempTest]
ListPlot[tempConv]
ListPlot[tempDeconv - tempTest]

enter image description here

enter image description here

enter image description here

share|improve this answer
1  
What a nice answer & explanation. +1 –  rasher Jan 30 at 7:59
    
That was an incredibly clear explanation. Thank you. –  user640078 Jan 30 at 19:50

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.