absorption spectrophotometry] [Appendix
Case Studies] [Peak
and Measurement] [iPeak]
downloadable software] [Interactive
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
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 representing the sine and cosine 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
by the rules for the division
complex numbers. Many computer languages will perform this
operation automatically when the two quantities divided are
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 four 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;
(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);
(4) if the denominator vector contains even a single zero value, 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 using the "remove zeros" function, rmz.m, or by smoothing the data before convolution and by constraining the Fourier deconvolution to a frequency region where the denominator is sufficiently high.
You can see the 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 many cases, the width of the physical convolution is not 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 result will
seldom be perfect, especially if the original signal is noisy, but
it is often a better approximation to the real underlying signal
than the recorded data without deconvolution.
As a method for peak sharpening, deconvolution can be
compared to the derivative
peak sharpening method described earlier or to the power method, in
which the raw signal is simply raised to some positive power n.
SPECTRUM, the freeware signal-processing application for Mac OS8 and earlier, includes a Fourier deconvolution function.
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
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
% The sum(c2) is included simply to scale the amplitude of the result to match the original y.
% 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')
% Deconvolution demo 2
% Underlying signal with a single peak (Gaussian) of unknown
% height, position, and width.
% 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; % <<<< Change the exponential time constant here
% Guess, or use prior knowledge, or curve fit one peak, to
% determine time constant (tc), then compute transfer function cc
% Use "deconv" to recover original signal uyy by deconvoluting cc
% from yy. It's necessary to zero-pad the observed signal as shown here.
% Plot the signals and results in 4 quadrants
plot(xx,uyy);title('Underlying 4 Gaussian signal, uyy');
plot(xx,cc);title('Exponential transfer function, cc')
plot(xx,yy);title('observed broadened and noisy signal, yy');
plot(xx,yydc);title('After deconvoluting transfer function, yydc')
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 the more complicated exponentially broadened Gaussian model.
If the exponential factor "tc" is not known, it can be determined by iterative curve fitting, manually adjusting the exponential factor ('extra') interactively to get the best fit:
which in this particular case gives a best fit when the exponential factor "tc" is adjusted to about 69.9 (very close the correct value of 70 in this simulation).
Alternatively, you can use peakfit.m with the unconstrained 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,FitError]=peakfit([xx;yy],0,0,1,31,70,10, [20 10 50])
FitResults = Peak# Position Height Width Area tc
1 25.006 0.99828 10.013 10.599 69.83
The value of the exponential factor determined by this method is 69.8, again close to 70. However, if the signal is very noisy, there will be quite a bit of uncertainty in the value of the exponential factor so determined - for example, the value will vary a bit if slightly different regions of the signal are selected for measurement (e.g. by panning or zooming in ipf.m or by changing the center and window arguments in peakfit.m). See Appendix H: Fourier deconvolution vs curve fitting for another example with four overlapping Gaussians.
Segmented deconvolution. If the peak widths or tailing vary substantially across the signal, you can use a segmented deconvolution, which allows the deconvolution vector to adapt to the local conditions in different signal regions. SegExpDeconv(x,y,tc) divides x,y into a number of equal-length segments defined by the length of the vector 'tc', then each segment is deconvoluted with an exponential decay of the form exp(-x./t) where t is corresponding element of the vector tc. Any number and sequence of t values can be used. SegExpDeconvPlot.m is the same except that it plots the original and deconvoluted signals and shows the divisions between the segments by vertical magenta lines to make it easier to adjust the number and values of the segments. SegGaussDeconv.m and SegGaussDeconvPlot.m are the same except that they perform a symmetrical (zero-centered) Gaussian deconvolution. SegDoubleExpDeconv.m and SegDoubleExpDeconvPlot.m perform a symmetrical (zero-centered)exponential deconvolution. If the peak widths increase regularly across the signal, you can calculate a reasonable initial value for the vector 'tc' by giving only the number of segments (NumSegments) , the first value, "startt", and the last value, "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,
courtesy of Natalie Harmann.
Revised December, 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 firstname.lastname@example.org.
Unique visits since May 17, 2008: