[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 x equals n, where m and n are known but x is unknown, then x equals n divided by m. Conversely if you know that m convoluted with x equals n, where m and n are known but x is unknown, then x equals m deconvoluted from n.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. 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.

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, 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 (

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

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

% 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 April 2017. This page is part of "

Unique visits since May 17, 2008: