index previous next

Fourier Deconvolution

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 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 of complex numbers.

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 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. Two examples of the 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 deconvolution operation is simply recovering the original signal before the low-pass filtering, noise and all.

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 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 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. Even the floating-point precision limitations of the computer can be enough to result in 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 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 deconvolution function.
Click to enlargeMatlab 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 (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 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 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.

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 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 Gaussian convolution and deconvolution of the same rectangular pulse, utilizing either the fft/ifft formulation immediately above or the custom function deconvgauss.m. (click for graphic). 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 left. (Download this script). In this case, the underlying signal (uyy) is a Gaussian, but in the observed signal (yy) the peak is broadened exponentially resulting in a shifted, shorter, and wider peak, and then a little constant white noise is added after the broadening convolution (cc). Assuming that the exponential broadening time constant ('tc') is known, or can be guessed or measured, the deconvolution of cc from yy successfully removes the broadening (yydc), but at the expense of considerable noise increase. However, the noise remaining in the deconvoluted signal is "blue" (high-frequency weighted) and so is easily reduced by smoothing and has less effect on least-square fits than does white noise. (For a greater challenge, try more noise in line 6 or a bad guess of time constant ('tc') in line 7). To plot the recovered signal overlaid with underlying signal: plot(xx,uyy,xx,yydc). To plot the observed signal overlaid with with underlying signal: plot(xx,uyy,xx,yy). To curve fit the recovered signal to a Gaussian to determine peak parameters: [FitResults,FitError]=peakfit([xx;yydc],26,42,1,1,0,10), which yields excellent values for the original peak positions, heights, and widths. You can demonstrate to yourself that with ten times the previous noise level (Noise=.01 in line 6), the values of peak parameters determined by curve fitting are still quite good, and even with 100x more noise (Noise=.1 in line 6) the peak parameters are more accurate than you might expect for that amount of noise (because that noise is blue). Remember, there is no need to smooth the results of the deconvolution before curve fitting, as seen previously.

  An alternative to the above deconvolution approach is to use iterative curve fitting to fit the observed signal directly with an exponentially broadened Gaussian (shape number 5)
: [FitResults,FitError]=peakfit([xx;yy], 26, 50, 1, 5, 40, 10). Both methods give good values of the peak parameters, but the deconvolution method is faster, because fitting the deconvoluted signal with a Gaussian model is faster than 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.  (See Appendix H: Fourier deconvolution vs curve fitting for an example with four overlapping Gaussians).

% 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)
  % <<<< Change the noise 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.
yydc=deconv([yy zeros(1,length(yy)-1)],cc).*sum(cc);

% 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')

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

% plot recovered signal overlaid with underlying signal:
% plot(xx,uyy,xx,yydc)

% plot observed signal overlaid with underlying signal:
% 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
% [FitResults,FitError]=peakfit([xx;yy],26,50,1,5,tc,10)

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, courtesy of Natalie Harmann.

index previous next
Revised May 2016. 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
Unique visits since May 17, 2008: