In practice, the deconvolution of one signal from another is usually performed by point-by-point division of the two signals in the Fourier domain, that is, dividing the Fourier transforms of the two signals point-by-point and then inverse-transforming the result. Fourier transforms are usually expressed in terms of complex numbers, with real and imaginary parts. If the Fourier transform of the first signal is a + ib, and the Fourier transform of the second signal is c + id, then the ratio of the two Fourier transforms is

Note: The word "deconvolution" appears in the Oxford dictionary, where its meaning is "A process of resolving something into its constituent elements or removing complication in order to clarify it". 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 putative 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 Figures 12 and 13.

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

When
applying 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
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 Figure 12 above.
However, this effect is not observed in the example in Figure
13 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 deconvolution function.

Matlab and Octave have a built-in function for 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')

Here is another example, shown on the left. (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 broadened

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

Revised January, 2015. This page is part of "

Unique visits since May 17, 2008: