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 mx equals n, where m and n are known but x is unknown, then x equals n/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 nIn 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

by the usual 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 underlying peak shape must be known but the broadening process, which determines the width and shape of the peaks in the recorded data, is unknown.  

The practical significance of Fourier deconvolution in signal processing is that it can be used as an artificial (i.e. 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.


Figure 12. 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 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.

Figure 13. A different application of the deconvolution function 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. 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.
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.

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
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 (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. Even with ten times the previous noise level (Noise=.01), the values of peak parameters determined by curve fitting are still quite good, and even with 100x more noise (Noise=.1) the peak parameters are more accurate than you might expect for that amount of noise (because that noise is blue).  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 interactive curve fitting, using ipf.m, 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
xx=5:.1:65;
% Underlying signal with a single peak (Gaussian) of unknown height, position, and width.
uyy=gaussian(xx,25,10);
% Observed signal yy, with noise added AFTER the broadening convolution (ExpG)
Noise=.001;
tc=70; % Time Constant
yy=expgaussian(xx,25,10,-tc)'+Noise.*randn(size(xx));
% Guess, or use prior knowledge, or curve fit one peak to ExpG to determine
% time constant (20), then compute transfer function cc
cc=exp(-(1:length(yy))./tc);  
% Attempt to recover original signal uyy by deconvoluting cc from yy
yydc=deconv([yy zeros(1,length(yy)-1)],cc).*sum(cc); % It's necessary to zero-pad the observed signal as shown here.
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 recovered signal overlaid with underlying signal
%  plot(xx,uyy,xx,yy) % plot observed signal overlaid with  with underlying signal
% Curve fit recovered signal to a Gaussian to determine peak parameters
% [FitResults,FitError]=peakfit([xx;yydc],26,42,1,1,0,10)
% Alternatively curve fit observed signal with an exponentially broadened
% Gaussian
% [FitResults,FitError]=peakfit([xx;yy],26,50,1,5,tc,10)


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

index previous next
Revised October, 2014. This page is 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: