[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]

index previous next

 Peak Sharpening (Resolution enhancement)

[The derivative method]   [The power law method]   [Spreadsheets]   [Matlab]   [Interactive tools]


The figure below shows a spectrum on the left that consists of several poorly-resolved (that is, partly overlapping) bands. The extensive overlap of the bands makes the accurate measurement of their intensities and positions impossible, even though the signal-to-noise ratio is very good. Things would be easier if the bands were more completely resolved, that is, if the bands were narrower.

A peak sharpening algorithm has been applied to the signal on the left to artificially improve the apparent resolution of the peaks. In the resulting signal, right, the component bands are narrowed so that the intensities and positions can be measured.

Here use can be made of peak sharpening (or resolution enhancement) algorithms to artificially improve the apparent resolution of the peaks. One of the simplest such algorithms is based on the weighted sum of the original signal and the negative of its second derivative:
  Rj = Yj - k1Y''
where Rj is the resolution-enhanced signal, Y is the original signal, Y'' is the second derivative of Y, and k is a user-selected weighting factor. It is left to the user to select the weighting factor k which gives the best trade-off between resolution enhancement, signal-to-noise degradation, and baseline flatness. The optimum choice depends upon the width, shape, and digitization interval of the signal. The result of the application of this algorithm is shown on the right. The component bands have been artificially narrowed so that the intensities and positions can be measured. The signal-to-noise ratio is degraded, but this can be moderated by smoothing. Nevertheless, this technique will be useful only if the overlap of peaks rather than the signal-to-noise ratio is the limiting factor.

Here's how it works. The figure below shows, in Window 1, a computer-generated peak (with a Lorentzian shape) in red, superimposed on the negative of its second derivative in green). (Click on the figure to see a full-size figure).

Click to view
          enlarged figure
Click to enlarge

The second derivative is amplified (by multiplying it by an adjustable constant) so that the negative sides of the inverted second derivative (from approximately X = 0 to 100 and from X = 150 to 250) are a mirror image of the sides of the original peak over those regions. In this way, when the original peak is added to the inverted second derivative, the two signals will approximately cancel out in the two side regions but will reinforce each other in the central region (from X = 100 to 150). The result, shown in Window 2, is a substantial (about 50%) reduction in the width, and a corresponding increase in height, of the peak. This effect is most dramatic with Lorentzian-shaped peaks; with Gaussian-shaped peaks, the resolution enhancement is less dramatic (only about 20 - 30%). Note that the baseline of either side of the resolution-enhanced peak is not quite flat, especially for a Lorentzian peak, because the cancellation of the original peak and the inverted second derivative is only approximate; the adjustable weighting factor k is selected to minimize this effect. Peak sharpening will have little or no effect on the baseline, because if the baseline is linear, its derivative will be zero, and if it is gradually curved, its second derivative will be very small compared to that of the peak.

Mathematically, this technique is a simplified version of a converging Taylor series expansion, in which only the even order derivative terms in the expansion are taken and for which their coefficients alternate in sign. The above example is the simplest possible version that includes only the first two terms - the original peak and its negative second derivative. Somewhat better results can be obtained by adding a fourth derivative term, with two adjustable factors k1 and k2:

     Rj = Yj - k1Y'' + k2Y''''

The result is a 26% reduction in width for a Gaussian peak, as shown in the figure on the left (Matlab/Octave script), and a 60% reduction for a Lorentzian peak (script). (It's possible to add a sixth derivative term, but the series converges quickly and the results are only slightly improved, at the cost of increased complexity of three adjustable factors). There is no universal optimum value for the derivative weighting factors; it depends on what you consider the best trade-off between peak sharpening and baseline flatness. However, reasonable values for a Gaussian peak are k1 = ((W/dx)2)/25 and k2 = ((W/dx)4)/800, and for Lorentzian peaks,  k1=((W/dx)2)/8 and k2 = ((W/dx)4)/700, where W is the FWHM of the peak and dx is the interval between x-axis data points. Note that the K factors for the 2nd and 4th derivatives vary with the width raised to the 2nd and 4th power respectively, so they can vary over a very wide numerical range for peaks of different width.

An interesting property of this procedure is that it does not change the total peak area (that is, the area under the peaks) because the total area under the curve of any derivative of a peak-shaped signal is zero (the area under the negatives lobes cancels the area under the positive lobes). Therefore this technique can be helpful in measuring the areas under overlapped peaks. However, this technique is not perfect because the baseline on either side of the sharpened peak is not completely flat, leaving some interference from nearby peaks. For a Lorentzian peak, about 80% of the area of the peak is contained in the central maximum. For a Gaussian peak, over 99% of the area of the peak is contained in the central maximum.  (The Matlab/Octave scripts just mentioned compute both the percent reduction in width and area and the percent baseline error for each set of weighting factors chosen).

Because differentiation and smoothing are both linear techniques, the amplitude of a sharpened signal is therefore proportional to the amplitude of the original signal, which allows quantitative analysis applications employing any of the standard calibration techniques). As long as you apply the same signal-processing techniques to the standards as well as to the samples, everything works.

Peak sharpening can be especially useful in automated peak detection and measurement to increase the ability to detect weak overlapping peaks that appear only as shoulders in the original signal. Click for animated example.


The Power Law Method. An even simpler method for peak sharpening involves simply raising each data point to a power n greater than 1 (reference 61). The effect of this is to change the peak shapes, essentially stretching out the highest center region of the peak to greater amplitudes and placing more weight on the points near the peak, resulting in more nearly Gaussian peak shapes (because most peak shapes are locally Gaussian near the peak maximum) and smaller peak widths (full width at half maximum by a factor of typically 2 or 3. Since the positions of the peaks are not moved by this process, the peak resolution (defined as the ratio of peak separation to peak width) is increased. The technique is demonstrated by the Matlab/Octave scripts PowerMethodGaussian.m and PowerMethodLorentzian.m, and by the self-contained function PowerLawDemo.m. In the figure on the left, the blue line shows two slightly overlapping peaks. The other lines are the result of raising the data to the power of n = 2, 3, and 4, and normalizing each to a height of 1.00. The peak widths, measured with the halfwidth.m function, are 19.2, 12.4, 9.9, and 8.4 units for powers 1 through 4, respectively. This method is independent of, and can be used in conjunction with, the derivative method discussed above. However, there are some limitations of the power law method: (a) it only works if the peaks of interest make a distinct maximum (it's not effective for side peaks that are so small that they only form shoulders; there must be a valley between the peaks); (b) the baseline must be zero for best results; (c) for noisy signals there is a slight decrease in signal-to-noise ratio because the smaller width means fewer data points are contributing to the measurement.

  Of course, this method introduces severe non-linearity into the signal, changing the ratios between peak heights (as is evident in the figure) and complicating further processing, especially quantitative measurement calibration. To correct this, after the raw data have been raised to the power n and peaks heights and/or areas have been measured, the resulting peak measures can be simply raised to the power 1/n, restoring the original linearity (but, notably, not the slope) of the calibration curves used in quantitative analytical measurements. (This works because the peak area is proportional to the height times width, and peak height of the power transformed peaks is proportional to the nth power of the original height, but the width of the peak is not a function of peak height at constant n, thus the area of the transformed peaks remains proportional to nth power of the original height). This is demonstrated quantitatively by the simple Matlab/Octave script PowerLawTest.m (graphic) and by the self-contained Matlab/Octave function PowerLawCalibrationCurve.m which takes the power, n - any positive number greater than 1 - as its single input argument (e.g. PowerLawCalibrationCurve(3) will demonstrate the 3rd power ).

Note: Another signal processing technique that can increase the resolution of overlapping peaks is deconvolution, which is more powerful but also more difficult to apply. Deconvolution of the broadening function from the broadened peaks is in principle capable of extracting the underlying peaks shapes, whereas the resolution enhancement technique described here can not be expected to do that.


SPECTRUM, the freeware signal-processing application  for Mac OS8 and earlier, includes this resolution-enhancement algorithm, with adjustable weighting factor and derivative smoothing width.

Resolution Enhancement for Excel and Calc. 
The derivative sharpening method with two derivative terms (2nd and 4th) is available for Excel and Calc (screen image) in the form of an empty template (PeakSharpeningDeriv.xlsx and PeakSharpeningDeriv.ods) or with example data entered (PeakSharpeningDerivWithData.xlsx and PeakSharpeningDerivWithData.ods). You can either type in the values of the derivative weighting factors K1 and K2 directly into cells J3 and J4, or you can enter the estimated peak width (FWHM in number of data points) in cell H4 and the spreadsheet will calculate K1 and K2. There is also a demonstration version with adjustable simulated peaks which you can experiment with (PeakSharpeningDemo.xlsx and PeakSharpeningDemo.ods). If the signal is noisy, adjust the smoothing just as with the smoothing spreadsheets using the 17 coefficients in row 5 columns K through AA. There is also a 20-segment version where the sharpening constants can be specified for each of 20 signal segments (SegmentedPeakSharpeningDeriv.xlsx). For applications where the peak widths gradually increase (or decrease) with time, there is also a gradient peak sharpening template (GradientPeakSharpeningDeriv.xlsx) and an example with data (GradientPeakSharpeningDerivExample.xlsx); you need only set the starting and ending peak widths and the spreadsheet will apply the required sharpening factors K1 and K2.

Resolution Enhancement for Matlab and Octave
The custom Matlab/Octave function enhance.m has the form Enhancedsignal = enhance(signal,k1,k2,SmoothWidth), where "signal" is the original signal vector, the arguments k2 and k4 are 2nd and 4th derivative weighting factors, and SmoothWidth is the width of the built-in smooth. The resolution-enhanced signal is returned in the vector "Enhancedsignal". Click on this link to inspect the code, or right-click to download for use within Matlab or Octave. The values of k1 and k2 determine the trade-off between peak sharpness and baseline flatness; the values vary with the peak shape and width and should be adjusted for your own needs. For peaks of Gaussian shape, a reasonable value for k2 is PeakWidth2/25 and for k4 is PeakWidth4/800 (or PeakWidth2/6 and PeakWidth4/700 for Lorentzian peaks), where PeakWidth is the full-width at half maximum of the peaks expressed in number of data points. 

 Here's a simple example that creates a signal consisting of four Gaussian peaks of equal width, applies the "enhance" function, and compares a plot (shown on the left) of the original signal (in blue) to the resolution-enhanced version (in red): 

   x=[0:.01:18];
  y=exp(-(x-4).^2)+exp(-(x-9).^2)+exp(-(x-12).^2)+exp(-(x-13.7).^2);
  Enhancedsignal=enhance(y,1000,1000000,3);
  plot(x,y,x,Enhancedsignal,'r')


     Four Gaussian peaks of equal width. Original: blue. After sharpening: red.

Peak sharpening can be a useful pre-process before measuring the areas under overlapping peaks, because it's easier and more accurate to measure the areas of peaks that are more completely separated, and because, as mentioned above, this method of sharpening does not change the total area under each peak.

Segmented peak sharpening. If the peak widths vary substantially across the signal, you can use the segmented version SegmentedSharpen.m, for which the input arguments factor1, factor2, and SmoothWidth are vectors. The script DemoSegmentedSharpen.m, shown on the right, uses this function to sharpen four Gaussian peaks with gradually increasing peak widths from left to right with increasing degrees of sharpening, showing that the peak width is reduced by 20% to 22% compared to the original. DemoSegmentedSharpen2.m shows four peaks of the same width sharpened to increasing degrees.

ProcessSignal, a Matlab/Octave command-line function that performs smoothing, differentiation, and peak sharpening on the time-series data set x,y (column or row vectors). Type "help ProcessSignal". Returns the processed signal as a vector that has the same shape as x, regardless of the shape of y. The syntax is

Processed=ProcessSignal(x,y, DerivativeMode, w, type, ends,Sharpen, factor1, factor2, SlewRate, MedianWidth)

iSignal is a Matlab function that performs resolution enhancement for time-series signals, using the above enhance function, with keystrokes that allow you to adjust the 2nd and 4th derivative weighting factors and the smoothing continuously while observing the effect on your signal dynamically. The E key turns the resolution enhancement function on and off. View the code here or download the ZIP file with sample data for testing. iSignal calculates the sharpening and smoothing settings for Gaussian and for Lorentzian peak shapes using the Y and U keys, respectively, using the expression given above. Just isolate a single typical peak in the upper window using the pan and zoom keys, press P to your on the peak measurement mode, then press Y for Gaussian or U for Lorentzian peaks. You can fine-tune the sharpening with the F/V and G/B keys and the smoothing with the A/Z keys. (The optimum settings depends on the width of the peak, so if your signal has peaks of widely different widths, one setting will not be optimum for all the peaks. In such cases, you can use the segmented sharpen function, SegmentedSharpen.m).

      

 Before Peak Sharpening in iSignal     After peak sharpening in iSignal


iSignal 5.95 can also use the power transform method (press the ^ key, enter the power, n (any positive number greater than 1.00) and press Enter.  To reverse this, simply raise to the 1/n power.

iPeak, a Matlab interactive peak detection and measurement program, has a built-in peak sharpening mode that is based on this technique. See ipeakdemo5.


Video Demonstration. This 15-second, 1.7 MByte video (ResEnhance3.wmv ) demonstrates the Matlab interactive resolution enhancement function InteractiveResEnhance.m. The signal consists of four overlapping, poorly-resolved Lorentzian bands. First, the 2nd derivative factor (Factor 1) is adjusted, then the 4th derivative factor (Factor 2) is adjusted, then the smooth width (Smooth) is adjusted, and finally the Factor 2 is tweaked again.

 
index previous next
Last updated July, 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 toh@umd.edu.
Unique visits since May 17, 2008: