In practice, the deconvolution of one signal from another is usually performed by point-by-point

Note: The word "deconvolution" can have two meanings, which can lead to confusion. The Oxford dictionary defines it as "A process of resolving something into its constituent elements or removing complication in order to clarify it", which in one sense applies to Fourier deconvolution. But the same word is also sometimes used for the process of resolving or decomposing a set of overlapping peaks into their separate additive components by the technique of iterative least-squares curve fitting of a proposed peak model to the data set. However, that process is actually conceptually distinct from

The practical significance of Fourier deconvolution in signal processing is that it can be used as a computational way to reverse the result of a convolution occurring in the physical domain, for example, to reverse the signal distortion effect of an electrical filter or of the finite resolution of a spectrometer. Two examples of the application of Fourier deconvolution are shown in the two figures below.

Note
that this process has an effect that is *visually *similar
to resolution enhancement,
although the later is done without specific knowledge of the
broadening function that caused the peaks to overlap.

When
applying Fourier deconvolution to experimental data, to remove the
effect of a known broadening or low-pass filter operator caused by
the experimental system, a very serious signal-to-noise
degradation commonly occurs. Any noise added to the signal by the
system *after* the convolution by the broadening or low-pass
filter operator will be greatly amplified when the Fourier
transform of the signal is divided by the Fourier transform of the
broadening operator, because the high frequency components of the
broadening operator (the denominator in the division of the
Fourier transforms) are typically very small, resulting in a great
amplification of high frequency noise in the resulting
deconvoluted signal. See the Matlab/Octave code example at the
bottom of this page. (Of course, if the denominator contains any *zeros*,
the result will be a "divide-by-zero" error and the whole
operation fails). The problem of low values or zeros in the
denominator can be reduced by smoothing before convolution and by
constraining the Fourier deconvolution to a frequency region where
the denominator is sufficiently high. You can see this
amplification of high frequency noise happening in the example in
the first example above. However, this effect is *not *observed in
the second example, because in that case the noise was present in
the original signal, *before *the convolution performed by
the spectrometer's derivative algorithm. The high frequency
components of the denominator in the division of the Fourier
transforms are typically much *larger *than in the previous
example, avoiding the noise amplification and divide-by-zero
errors, and the only post-convolution noise comes from numerical
round-off errors in the math computations performed by the
derivative and smoothing operation, which is always much smaller
than the noise in the original experimental signal.

SPECTRUM, the freeware signal-processing application that accompanies this tutorial, includes a Fourier deconvolution function.

Matlab and Octave have a built-in function for Fourier deconvolution: deconv. An example of its application is shown below: the vector yc (line 6) represents a noisy rectangular pulse (

x=0:.01:20;y=zeros(size(x));

y(900:1100)=1; % Create a rectangular function y,

% 200 points wide

y=y+.01.*randn(size(y)); % Noise added

c=exp(-(1:length(y))./30); % exponential trailing convolution

% function, c

yc=conv(y,c,'full')./sum(c); % Create exponential trailing rectangular

% function, yc

% yc=yc+.01.*randn(size(yc)); % Noise added

ydc=deconv(yc,c).*sum(c); % Attempt to recover y by

% deconvoluting c from yc

subplot(2,2,1); plot(x,y); title('original y');subplot(2,2,2); plot(x,c);title('c'); subplot(2,2,3); plot(x,yc(1:2001)); title('yc'); subplot(2,2,4); plot(x,ydc);title('recovered y')

Alternatively, you could perform the Fourier deconvolution yourself

The script DeconvDemo3.m is similar to the above, except that it demonstrates

GaussConvDemo.m shows a Gaussian peak convoluted with a Gaussian function and an attempt to recover the original peak, in this case by using the deconvgauss.m function.

Here is another example, shown on the right. (Download this script). In this case, the underlying signal (

xx=5:.1:65;

% Underlying signal with a single peak (Gaussian) of unknown

% height, position, and width.

uyy=gaussian(xx,25,10);

% Compute observed signal yy, using the

% constant tc, adding noise added AFTER the broadening convolution (ExpG)

Noise=.001;

tc=70;

yy=expgaussian(xx,25,10,-tc)'+Noise.*randn(size(xx));

% Guess, or use prior knowledge, or curve fit one peak, to

% determine time constant (tc), then compute transfer function cc

cc=exp(-(1:length(yy))./tc);

% Use "deconv" to recover original signal uyy by deconvoluting cc

yydc=deconv([yy zeros(1,length(yy)-1)],cc).*sum(cc);

% Plot the signals and results in 4 quadrants

subplot(2,2,1);

plot(xx,uyy);title('Underlying 4 Gaussian signal, uyy');

subplot(2,2,2);

plot(xx,cc);title('Exponential transfer function, cc')

subplot(2,2,3);

plot(xx,yy);title('observed broadened and noisy signal, yy');

subplot(2,2,4);

plot(xx,yydc);title('After deconvoluting transfer function, yydc')

% Try more Noise (line 6) or a bad guess of time constant (line 11)

% plot(xx,uyy,xx,yydc)

% plot(xx,uyy,xx,yy)

% Curve fit recovered signal to a Gaussian to determine peak parameters:

% [FitResults,FitError]=peakfit([xx;yydc],26,42,1,1,0,10)

% Or curve fit observed signal with an exponentially broadene

An alternative to the above deconvolution approach is to use iterative curve fitting (which is covered in a further section) to fit the observed signal directly with an exponentially broadened Gaussian (shape number 5)

>> [FitResults,FitError]=peakfit([xx;yy], 26, 50, 1, 5, 70, 10)

Both methods give good values of the peak parameters, but the Fourier deconvolution method is faster, because fitting the deconvoluted signal with a simple Gaussian model is faster than iteratively curve fitting the observed signal with a exponentially broadened Gaussian model. Also, if the exponential factor "tc" is not known, it can be determined by iterative curve fitting, adjusting the exponential factor ('extra') interactively to get the best fit.

Alternatively, you can use peakfit.m with the constrained variable exponentially broadened Gaussian (shape 31), which will automatically find the best value of "tc", but in that case the best results will be obtained if you give it a rough starting guess ("start") at least within a factor of two or so of the correct values:

iSignal version 5.7 has a

This page is also available in French, at http://www.besteonderdelen.nl/blog/?p=41, courtesy of Natalie Harmann.

Revised February 2017. This page is part of "

Unique visits since May 17, 2008: