[Introduction]
[Signal
arithmetic] [Signals
and
noise] [Smoothing]
[Differentiation]
[Peak
Sharpening] [Harmonic
analysis] [Fourier
convolution] [Fourier
deconvolution] [Fourier
filter] [Peak
area
measurement] [Linear
Least
Squares] [Multicomponent
Spectroscopy] [Iterative
Curve
Fitting] [Hyperlinear
quantitative
absorption spectrophotometry] [Appendix
and
Case Studies] [Peak
Finding
and Measurement] [iPeak]
[iSignal]
[Peak
Fitters] [iFilter]
[iPower]
[List
of
downloadable software] [Interactive
tools]

[Applications] [Matlab/Octave]

Fourier
deconvolution
is
the converse of Fourier convolution
in
the sense that division is the converse of multiplication. If you
know that * m* times

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

^{}

by the rules for the division
of
complex numbers. Many computer languages will perform this
operation automatically when the two quantities divided are
complex.

**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
*Fourier* deconvolution, because in Fourier deconvolution,
the
underlying peak shape is unknown but the broadening function is
assumed to be known; whereas in iterative least-squares curve
fitting
it's just the reverse: the peak shape must be known but the width
of
the broadening process, which determines the width and shape of
the
peaks in the recorded data, is unknown. Thus the term "spectral
deconvolution" is ambiguous: it might mean the Fourier
deconvolution of a response function from a spectrum, or it might
mean the decomposing of a spectrum into its separate additive peak
components. These are different processes; don't get them
confused.

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. In some cases the
physical
convolution can be measured experimentally by applying a single
spike
impulse ("delta") function to the input of the system, then
that data used as a deconvolution vector. Deconvolution can also
be
used to determine the form of a convolution operation that has
been
previously applied to a signal, by deconvoluting the original and
the
convoluted signals. These two types of application of Fourier
deconvolution are shown in the two figures below.

**Fourier deconvolution*** is used here
to
remove the distorting influence of an exponential tailing
response
function from a recorded signal (Window 1, top left) that is the
result of an unavoidable RC low-pass filter action in the
electronics. The response function (Window 2, top right) must be
known and is usually either calculated on the basis of some
theoretical model or is measured experimentally as the output
signal
produced by applying an impulse (delta) function to the input of
the
system. The response function, with its maximum at x=0, is
deconvoluted from the original signal . The result (bottom,
center)
shows a closer approximation to the real shape of the peaks;
however,
the signal-to-noise ratio is unavoidably degraded compared to
the
recorded signal, because the Fourier deconvolution operation is
simply recovering the original signal before the low-pass
filtering,
noise and all. (Matlab/Octave
script)*

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.

**A different application of Fourier deconvolution***
is to reveal the nature of an unknown data transformation
function
that has been applied to a data set by the measurement
instrument
itself. In this example, the figure in the top left is a
uv-visible
absorption spectrum recorded from a commercial photodiode array
spectrometer (X-axis: nanometers; Y-axis: milliabsorbance). The
figure in the top right is the first
derivative of this spectrum produced by an (unknown)
algorithm in
the software supplied with the spectrometer. The objective here
is to
understand the nature of the differentiation/smoothing
algorithm that the instrument's software uses. The signal
in the
bottom left is the result of deconvoluting the derivative
spectrum
(top right) from the original spectrum (top left). This
therefore
must be the convolution function used by the differentiation
algorithm in the spectrometer's software. Rotating and expanding
it
on the x-axis makes the function easier to see (bottom right).
Expressed in terms of the smallest whole numbers, the
convolution
series is seen to be +2, +1, 0, -1, -2. This simple
example of
"reverse engineering" would make it easier to compare
results from other instruments or to duplicate these result on
other
equipment.*

When
applying Fourier deconvolution to experimental data, for example
to
remove the effect of a known broadening or low-pass filter
operator
caused by the experimental system, there are three serious
problems
that limit the utility of the method: (1) the convolution
occurring
in the physical domain might not be accurately modeled by a
mathematical convolution; (2) the width of the convolution - for
example the time constant of a low-pass filter operator or the
shape
and width of a spectrometer slit function - must be known, or at
least adjusted by the user to get the best results, and (3) a
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.

In practice, the width of the physical convolution is seldom known exactly, so the deconvolution must be adjusted empirically to yield the best results. Similarly, the width of the final smooth operation must also be adjusted for best results. The results will seldom be perfect, but are often an improvement over the original.

As a method for *peak sharpening*, deconvolution can be
compared to the methods
described earlier, in particular the much simpler power
method, in which the raw signal is simply raised to some
positive
power *n.* The Matlab/Octave script PowerMethod.m
compares these two methods for a pair of overlapping peaks,
showing
that power method (with *n*=5) is capable of reducing the
peak
width and increasing the separation-to-width ratio better than
deconvolution, but it does not restore the original peak
positions,
heights, and widths as does the deconvolution method (graphic).

SPECTRUM, the freeware signal-processing application for Mac OS8 and earlier, 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 (*y*) convoluted with a transfer function *c*
before being measured. In line 7, *c* is deconvoluted from *yc*,
in an attempt to recover the original *y*. This requires
that
the transfer function *c* be known. The rectangular signal
pulse
is recovered in the lower right (*ydc*), complete with the
noise
that was present in the original signal. The Fourier
deconvolution reverses not only the signal-distorting effect of
the
convolution by the exponential function, but also its low-pass
noise-filtering effect. As explained above, there is significant
amplification of any noise that is added *after* the
convolution
by the transfer function (line 5). This script demonstrates that
there is a big difference between noise added *before* the
convolution (line 3), which is recovered unmodified by the Fourier
deconvolution along with the signal, and noise added *after*
the
convolution (line 6), which is amplified compared to that in the
original signal. Execution time: 0.03 seconds in Matlab; 0.3
seconds
in Octave. Download this
script.

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 **before
**the
convolution

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 **after
**the
convolution

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

% deconvoluting c from yc

`%
Plot all the steps`

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 *without *using
the
built-in Matlab/Octave "deconv" function by dividing
the Fourier transforms of *yc *and *c* using the
built-in
Matlab/Octave "fft.m" function and inverse transform the
result with the built-in Matlab/Octave "ifft.m" function.
Note that *c* must be zero-filled
to match the size of *yc*. The results are essentially the
same
(except for the numerical floating point precision of the
computer,
which is usually negligible), and it's actually ten times faster
than
using the deconv function:

`ydc=ifft(fft(yc)./fft([c
zeros(1,2000)])).*sum(c);
`

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

DeconvDemo4.m (on the right) shows a Gaussian deconvoluted from a Gaussian function and an attempt to recover the original peak width. This is an example of "self deconvolution", so-called because the shape of the deconvolution function is the same as the shape of the peaks in the signal, in this case both

In all of the above simulations, the deconvolution method works as well as it does because the signal-to-noise ratio of the "observed signal" (upper right quadrant) is quite good; the noise is not even visible on the scale presented here. In the absence of any knowledge of the width of the deconvolution function, finding the right deconvolution width depends upon experimentally minimizing the wiggles that are seen when the deconvolution width is incorrect, and a poor signal-to-noise ratio will make this much more difficult. Of course smoothing can reduce noise, especially high-frequency (blue) noise, but smoothing also slightly increases the width of peaks, which works counter to the point of deconvolution, so it must not be over used. The image on the left shows the widths of the peaks (as full width at half maximum); the width of the deconvoluted peaks (lower right quadrant) are only slightly larger than in the (unobserved) underlying peaks (upper left quadrant) either because of imperfect deconvolution or the broadening effects of the smoothing needed to reduce the high frequency noise. As a rough but practical rule of thumb, if there is any

Here is another example, shown below on the left. (Download this script). 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 expgaussian function with time

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

Noise=.001; % <<<< Change the noise here

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 Gaussian

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

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:

FitResults =

1 25.006 0.99828 10.013 10.599 69.83

GoodnessOfFit =

0.15575 0.99998

tstep=(endt-startt)/NumSegments;

tc=startt:tstep:endt;

**iSignal **version 5.7
has a **Shift-V** keypress that displays the menu of Fourier
convolution and deconvolution operations that allow you to
convolute
a Gaussian or exponential function with the signal, or to
deconvolute
a Gaussian or exponential function from the signal, and asks you
for
the width or the time constant (in X units).

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

Revised June, 2017. This page is part
of "**A
Pragmatic Introduction to Signal Processing**", created
and maintained by Prof.
Tom
O'Haver , Department of Chemistry and Biochemistry, The
University of Maryland at College Park. Comments, suggestions and
questions should be directed to Prof. O'Haver at toh@umd.edu.

Unique visits since May 17, 2008: