index previous next

Peak Finding and Measurement

[findpeaks]   [findvalleys]   [Signal-to-noise ratio]   [iPeak]  [iPeak demo scripts]  [Spreadsheets]  [Video Demonstration]  [Have a question? Email me]

A common requirement in signal processing is to detect peaks in a signal and to measure their positions, heights, and widths. A common way to do this is to make use of the fact that the first derivative of a peak has a downward-going zero-crossing at the peak maximum.  But the presence of random noise in real experimental signal will cause many false zero-crossing simply due to the noise. To avoid this problem, the technique described here first smooths the first derivative of the signal, before looking for downward-going zero-crossings, and then it takes only those zero crossings whose slope exceeds a certain predetermined minimum (called the "slope threshold") at a point where the original signal exceeds a certain minimum (called the "amplitude threshold").  By carefully adjusting the smoothwidth, slope threshold, and amplitude threshold, it's possible to detect only the desired peaks and ignore peaks that are too small, too wide, or too narrow.  Moreover, because smoothing can distort peak signals, reducing peak heights and increasing peak widths, this technique estimates the position, height, and width of each peak by least-squares curve-fitting of a segment of the original unsmoothed signal in the vicinity of the zero-crossing. Thus, even if heavy smoothing of the first derivative is necessary to provide reliable discrimination against noise peaks, the peak parameters extracted by curve fitting are not distorted. (This is useful primarily for signals that have several data points in each peak, not for spikes that have only one or two points).  This technique is capable of measuring peak positions and heights quite accurately, but the measurements of peak widths and areas is accurate only if the peaks are approximately Gaussian in shape.  (For the most accurate measurement of highly overlapped peaks, iPeak version 3 and above is capable of utilizing non-linear iterative curve fitting).

The routine is now available in four different versions:
(1) the basic command-line findpeaks function (findpeaks.m and findvalleys.m) for Matlab or Octave (not to be confused with the findpeaks function in the Signal Processing Toolbox).
(2) an interactive keypress-operated function, called iPeak (ipeak.m). for adjusting the peak detection criteria in real-time to optimize for any particular peak type (Matlab only).
(3) an older script using mouse-controlled sliders for interactive control, which works only in Matlab 6.5 but not in more recent versions.
(4) A spreadsheet, available in Excel and in OpenOffice open document formats.

Click here to download the ZIP file "PeakFinder.zip", which includes findpeaks.m, ipeak.m, and a sample data file and demo scripts for testing. You can also download iPeak and other programs of mine from the Matlab File Exchange.

The basic command-line functions (findpeaks.m and findvalleys.m) for Matlab or Octave.

function P=findpeaks(x,y,SlopeThreshold,AmpThreshold,SmoothWidth,FitWidth,smoothtype)
This is my basic command-line function to locate and measure the positive peaks in a noisy data sets. It detects peaks by looking for downward zero-crossings in the smoothed first derivative that exceed SlopeThreshold and peak amplitudes that exceed AmpThreshold, and determines the position, height, and approximate width of each peak by least-squares curve-fitting the top part of the peak. Returns a list (in matrix P) containing the peak number and the estimated position, height, and width of each peak.
It can find and measure over 1000 peaks per second in very large signals.

Version
5.0 (December, 2012) has improved error catching for very noisy data.

The data are passed to the findpeaks function in the vectors x and y (x = independent variable, y = dependent variable). The other parameters are user-adjustable:

There is also a similar function for finding valleys (minima), called findvalleys.m, which works the same way as findpeaks.m, except that it locates minima instead of maxima.  Right-click the link and select "Save link as...."

Examples:

>> x=[0:.01:50];y=cos(x);P=findpeaks(x,y,0,-1,5,5)
P =
            1       6.2832            1       2.3569
            2       12.566            1       2.3569
            3        18.85            1       2.3568
            4       25.133            1       2.3569
            5       31.416            1       2.3569
            6       37.699            1       2.3568
            7       43.982            1       2.3569

>> x=[0:.01:50];y=cos(x);P=findvalleys(x,y,0,-1,5,5)
P =
            1       3.1416           -1       2.3571
            2       9.4248           -1       2.3571
            3       15.708           -1       2.3571
            4       21.991           -1       2.3571
            5       28.274           -1       2.3571
            6       34.558           -1       2.3571
            7       40.841           -1      
2.3571


Note: The accuracy of the measurements of  peak position, height, and width by the findpeaks function depends on the shape of the peaks, the extent of peak overlap, and signal-to-noise ratio.  The width measurements particularly are strongly influenced by peak overlap, noise, and the choice of FitWidth. Isolated peaks of Gaussian shape are measured most accurately.  See "ipeakdemo" (below) for an accuracy trial. For the most accurate measurements of highly overlapping peaks of various shapes, use the Interactive Peak Fitter (ipf.m)..

The function findpeaksplot.m is a simple variant of findpeaks.m that also plots the x,y data and numbers the peaks on the graph (if any are found).  An example is shown on the right.

The function peakstats.m uses the same algorithm as findpeaks, but it computes and returns a table of summary statistics of the peak intervals (the x-axis interval between adjacent detected peaks), heights, widths, and areas, listing the maximum, minimum, average, and percent standard deviation of each, and optionally plotting the x,y data with numbered peaks in figure window 1, printing the table of peak statistics in the command window, and plotting the histograms of the peak intervals, heights, widths, and areas in figure window 2. Type "help peakstats". The syntax is the same as findpeaks, with the addition of a 8th input argument to control the display and plotting. Example:

>> x=[0:.1:1000];y=5+5.*cos(x)+randn(size(x));PS=peakstats(x,y,0,-1,15,23,3,1);


Peak Summary Statistics
158 peaks detected
          Interval      Height      Width          Area
Maximum    6.6195       10.8991       4.7749       47.3653
Minimum    5.9086       9.1022       2.4123       26.3111
Mean       6.2833       10.0183       3.2346       34.4027
% STD      1.7362       3.2154       13.028       10.9314


With the last input argument omitted or equal to zero, the plotting and printing in the command window are omitted; the numerical values of the peak statistics table are returned as a 4x4 array, in the same order as the example above.

Findpeaksnr.m is a variant of findpeaks.m that additionally computes the signal-to-noise ratio (SNR) of each peak and returns it in the 5th column of the peak table. The SNR is computed as the ratio of the peak height to the root-mean-square residual (difference between the original unsmoothed data and the least-squares fit over the top part of the peak).  Example:

>> x=[0:.01:5];P=findpeaksnr(x,x.*sin(x.^2).^2+.1*whitenoise(x),.001,1,15,10)
P =
            1       1.3224       1.3021      0.40005       32.829
            2       2.1867       2.1641      0.46717       49.023
            3       2.7996       2.7787      0.36404       49.856
            4       3.3198       3.3233      0.23066       61.675
            5        3.763       3.7765      0.19827       96.878
            6       4.1584       4.1666       0.1956       94.435
            7        4.519       4.5264      0.18169       116.62
            8       4.8547       4.8518      0.16985       198.42


Using the peak table. All these functions return the peak table as a matrix (with one row and for each peak detected and columns listing the peak number, position, height, and width). You can assign this matrix to a variable (e.g. P, in the examples above) and then use Matlab/Octave functions to extract specific information from that matrix. For example:

mean(P(:,3)) returns the average peak height (because peak height is in the 3rd column);
max(P(:,3)) returns the maximum peak height;
hist(P(:,3)) displays the histogram of peak heights;
std(P(:,4))/mean(P(:,4)) returns the relative standard deviation of the peak widths (column 4)
for n=1:length(P)-1;d(n)=max(P(n+1,2)-P(n,2));end  creates "d" as the vector of x-axis distance between peaks.

Using my val2ind.m function:
val2ind(P(:,3),7.5) returns the peak number of the peak whose height is closest to 7.5;
P(val2ind(P(:,2),7.5),3) returns the peak height (3rd column) of the peak whose position (2nd column) is closest to 7.5;
and so on.

DemoFindPeak (DemoFindPeak.m)

Click to view enlarged
                figure
Click to view enlarged figure
A simple demonstration script using the findpeaks function on noisy synthetic data. Numbers the peaks and prints out the P matrix in the Matlab command window: 
    Peak #   Position   Height    Width
P =
1.0000 10.0235 0.1814 1.7506
2.0000 29.9964 0.6412 1.9931
3.0000 50.0205 0.4315 2.1433
4.0000 70.0120 1.3935 2.0202
5.0000 79.9567 0.2580 1.8267
6.0000 100.0086 1.0862 2.0306
......
DemoFindPeakSNR is a variant of DemoFindPeak.m that uses findpeaksnr.m to compute the signal-to-noise ratio (SNR) of each peak and returns it in the 5th column.  

Peak Identification  

The command line function idpeaks.m is used for identifying peaks according to their x-axis maximum positions. The syntax is
[IdentifiedPeaks,AllPeaks]=idpeaks(DataMatrix,AmpT,SlopeT,sw,fw,maxerror,Positions,Names)
It finds peaks in the signal "DataMatrix" (x-values in column 1 and y-values in column 2), according to the peak detection parameters AmpT, SlopeT, sw, fw (see the "findpeaks" function above), then compares the found peak positions (x-values) to a database of known peaks, in the form of an array of known peak maximum positions ('Positions') and matching cell array of names ('Names'). If the position of a peak found in the signal is closer to one of the known peaks by less than the specified maximun error ('maxerror'), that peak is considered a match and its peak position, name, error, and peak amplitude (height) are entered into the output cell array "IdentifiedPeaks". The full list of detected peaks, identified or not, is returned in "AllPeaks". Use "cell2mat" to access numeric elements of IdentifiedPeaks,e.g. cell2mat(IdentifiedPeaks(2,1)) returns the position of the first identified peak, cell2mat(IdentifiedPeaks(2,2)) returns its name, etc. Obviously for your own applications, it's up to you to provide your own array of known peak maximum positions ('Positions') and matching cell array of names ('Names') for your particular types of signals.

Example:  Download idpeaks.zip, extract it, and place the extracted files in the Matlab or Octave path. This contains a high-resolution atomic emission spectrum of copper ('spectrum') and a data table of known Cu I and II atomic lines ('DataTable') containing their positions and names.

>> load DataTable
>> load spectrum
>> idpeaks(Cu,0.01,.001,5,5,.01,Positions,Names)

ans=
'Position' 'Name' 'Error' 'Amplitude'
[ 221.02] 'Cu II 221.027' [ -0.0025773] [ 0.019536]
[ 221.46] 'Cu I 221.458' [ -0.0014301] [ 0.4615]
[ 221.56] 'Cu I 221.565' [-0.00093125] [ 0.13191]
.................... etc...

The Matlab iPeak Function (ipeak.m), Version 5.3  Download as a ZIP file with sample data.


Example 1


Peaks in annual sunspot numbers from 1700 to 2009 (download the datafile)*.


Press U key to switch between peak and valley mode.




                     Example 2




ECG data (normal and inverted)

ipeak

ipeak4
Example 3


Example 5. Pressing 'L' toggles ON and OFF the peak labels in the upper window.


Example 8. The peak identification function applied to a high-resolution atomic spectrum.


Example 8. Three peaks near 296 nm isolated and identified. Press the I key to display the peak ID names.



Normal Peak Fit (N key)



iPeak 3.2 and later, showing log scale (Y key) and Autozero mode (T key)



Click to view larger figures.
iPeak is a keyboard-operated Interactive Peak Finder for time series data, based on the "findpeaks.m" function, for Matlab only. (What's new in Version 5.3?) It accepts data in a single vector, a pair of vectors, or a matrix with the independent variable in the first column and the dependent variable in the second column:
 
Example 1:   One input argument; data in single vector
              
>> y=cos(.1:.1:100);ipeak(y) 

Example 2:   One input argument; data in two columns of a matrix
              
>> x=[0:.01:5]';y=x.*sin(x.^2).^2;M=[x y];ipeak(M)

Example 3:   Two input arguments; data in separate x and y vectors
              
>> x=[0:.1:100];y=(x.*sin(x)).^2;ipeak(x,y);

Example 4:   Additional input argument (after the data) to control peak sensitivity.

       >> x=[0:.1:100];y=5+5.*cos(x)+randn(size(x));ipeak(x,y,10);
   or  >> ipeak([x;y],10);
   or  >> ipeak(humps(0:.01:2),3)
   or  >> x=[0:.1:10];y=exp(-(x-5).^2);ipeak([x' y'],1)

The additional numeric argument 
is an estimate of the ratio of the typical peak width to the length of the entire data record (PeakD). Small values detect fewer peaks; larger values detect more peaks. It effects only the starting values for the peak detection parameters. (This is just a quick way to set reasonable initial values of the peak detection parameters, rather than specifying each one individually as in the following example).

iPeak displays the entire signal in the lower half of the Figure window and an adjustable zoomed-in section in the upper window. Pan and zoom the portion in the upper window using the cursor arrow keys. (Press Shift-C to change the plot color.)  Adjust the peak detection parameters AmpThreshold (A/Z keys), SlopeThreshold (S/X), SmoothWidth (D/C), FitWidth (F/V) so that it detects the desired peaks and ignores those that are too small, too broad, or too narrow to be of interest. Detected peaks are numbered from left to right. Press P to display the peak table of all the detected peaks (Peak #, Position, Height, Width, Area). Press Shift-P to save peak table as disc file. Press U to switch between peak and valley mode (version 3.7 and above). Note: in version 5.1 and later, to speed up the operation for signals over 100,000 points in length, the lower window is refreshed only when the number of detected peaks changes or if the Enter key is pressed. Press K to see all the keystroke commands.

Peak Summary Statistics. Version 5.3 (April, 2013) adds the E key command to print a table of summary statistics of
the peak intervals (the x-axis interval between adjacent detected peaks), heights, widths, and areas, listing the maximum, minimum, average, and percent standard deviation, and displaying the histograms of the peak intervals, heights, widths, and areas in figure windows 2 through 5.
 
Peak Summary Statistics
15 peaks detected
Autozero OFF
          Interval      Height      Width          Area
Maximum    6.3795       10.5308      3.2354       34.943
Minimum    6.1649        9.7355      2.6671       29.9008
Mean       6.291        10.1559      3.0149       32.5771
% STD      0.91178       1.904       5.2584        4.3022


Example 5: Six input arguments. As above, but input arguments 3 to 6 directly specifies initial values of AmpThreshold (AmpT), SlopeThreshold (SlopeT), SmoothWidth (SmoothW), FitWidth (FitW) . PeakD is ignored in this case, so just type a '0' as the second argument after the data matrix).

                 >> ipeak(datamatrix,0,.5,.0001,20,20);

Keystrokes allow you to pan and zoom the upper window, to inspect each peak in detail if desired.  You can set the initial values of pan and zoom in optional input arguments 7 ('xcenter') and 8 ('xrange').  See example 6 below.

The peak cloest to the center of the upper window is labeled in the upper left of the top window, and it peak position, height, and width are listed. The Spacebar/Tab keys
jump to the next/previous detected peak and displays it in the upper window at the current zoom setting (use the up and down cursor arrow keys to adjust the zoom range) The Y key toggles between linear and log y-axis scale in the lower window (a log axis is good for inspecting signals with high dynamic range; it effects only the lower window display and has no effect on the peak detection or measurements). 
 
Example 6: Eight input arguments. As above, but input arguments 7 and 8 specifiy the initial pan and zoom settings, 'xcenter' and 'xrange', respectively.  In this example, the x-axis data are wavelengths in nanometers (nm), and the upper window zooms in on a very small 0.4 nm region centered on 249.7 nm. (These data, provided in the ZIP file, are from a high-resolution atomic spectrum).

        >> load ipeakdata.mat
       
>> ipeak(Sample1,0,100,0.05,3,4,249.7,0.4);

Autozero mode. The T key toggles the autozero mode off and on.  When autozero is OFF, peak heights are measured relative to zero.  (If the peaks are superimposed on a background, use the baseline subtract keys - B and G - first to subtract the background).  When autozero is ON, peak heights are automatically measured relative to the local baseline on either side of the peak; use the zoom controls to isolate the peaks so that the signal returns to the local baseline between the peaks as displayed in the upper windowWhen autozero is ON, the peak heights, widths, and areas in the peak table (R or P keys) will be automatically corrected for the baseline.  (Autozero OFF will give better results when the baseline is zero, or has been subtracted using the B key, even if the peaks are partly overlapped. Autozero ON will work best if the peaks are well separated so that the signal returns to the local baseline between the peaks.  If the peaks are highly overlapped, or if they are not Gaussian in shape, the best results will be obtained by using the curve fitting function - the N or M keys). 
Example 7:  Nine input arguments. As example 6, but the 9th input argument turns ON the autozero mode (equivalent to pressing the T key).  If not specified, autozero is initially OFF.

       
>> ipeak(Sample1,0,100,0.00,3,4,249.7,0.4,1);

Normal and Multiple Peak fitting: The N key applies iterative curve fitting to the detected peaks that are displayed in the upper window (referred to here as "Normal" curve fitting).  The use of the iterative least-squares function can result in more accurate peak parameter measurements that the normal peak table (R or P keys), especially if the peaks are non-Gaussian in shape or are highly overlapped. If the peaks are superimposed on a background, turn on the Autozero mode (T key). Then use the pan and zoom keys to select a peak or a group of overlapping peaks in the upper window, with the signal returning all the way to the local baseline at the ends of the upper window.  Make sure that AmpThreshold, SlopeThreshold, SmoothWidth are adjusted so that each peak is numbered once.  Then press the N key, type a number for the desired peak shape at the prompt in the Command window and press Enter (1=Gaussian (default), 2=Lorentzian, 3=logistic, 4=Pearson, 5=exponentially broadened Gaussian; 6=equal-width Gaussians; 7=Equal-width Lorentzians; 8=exponentially broadened equal-width Gaussian, 9=exponential pulse, 10=sigmoid, 11=Fixed-width Gaussian, 12=Fixed-width Lorentzian; 13=Gaussian/Lorentzian blend; 14=bifurcated Gaussian, 15=bifurcated Lorentzian), then type in a number of repeat trial fits and press Enter (the default is 1; start with that and then increase if necessary).  If you have selected a variable-shape peak (numbers 4, 5, 8 ,13, 14, or 15), the program will ask you to type in a number that controls the shape ("extra" in the peakfit input arguments). The program will then perform the fit, display the results graphically in Figure window 2, and print out a table of results in the command window, e.g.:

Peak shape (1-8): 2
Number of trials: 1

Least-squares fit to Lorentzian peak model
Fitting Error 1.1581e-006%
          Peak#     Position     Height      Width         Area 
            1          100            1           50       71.652
            2          350            1          100       146.13
            3          700            1          200       267.77

There is also a "Multiple" peak fit function (M key) that will attempt to apply
iterative curve fitting to all the detected peaks in the signal simultaneously. Before using this function, it's best to turn off the Autozero (T key) and use the multi-segment baseline correction function (B key) to remove the background (because the autozero function will probably not be able to subtract the baseline from the entire signal). Then press M and proceed as for the normal curve fit. A multiple curve fit may take a minute or so to complete if the number of peaks is large, possibly longer than the Normal curve fitting function on each group of peaks separately.  Make sure all the peaks are detected and numbered before activating this function, because it depends on the peak table for the number of peaks and the starting values, and it will fit only those peaks that it has already located and numbered.

The N and M key fitting functions perform non-linear iterative curve fitting using the peakfit.m function. The number of peaks and the starting values of peak positions and widths for the curve fit function are automatically supplied by the the findpeaks function, so it is essential that the peak detection variables in iPeak be adjust so that all the peaks in the selected region are detected and numbered once.  (For more flexible curve fitting, use ipf.m, which allows manual optimization of peak groupings and start positions). 

Note 1: If the peaks are too overlapped to be detected and numbered separately, try pressing the H key to activate the sharpen function before pressing M (version 4.0 and above only).

Note 2
If you plan to use a variable-shape peak (numbers 4, 5, 8 ,13, 14, or 15) for the Multiple peak fit, it's a good idea to obtain a reasonable value for the requested "extra" shape parameter by performing a Normal peak fit on an isolated single peak (or small group of partly-overlapping peaks) of the same shape, then use that value for the Multiple curve fit of the entire signal.

Note 3: Although the peak widths can vary from peak to peak, the curve fit routines assume that the  peak shape is the same for all peaks in one fitting operation, so if the peak shape varies accross the signal, use the Normal peak fit to fit each section with a different shape rather than the Multiple peak fit.

Peak identification. There is an optional "peak identification" function if optional input arguments 9 ('MaxError'), 10 ('Positions'), and 11 ('Names') are included. The "I" key toggles this function ON and OFF. This function compares the found peak positions (maximum x-values) to a database of known peaks, in the form of an array of known peak maximum positions ('Positions') and matching cell array of names ('Names'). If the position of a found peak in the signal is closer to one of the known peaks by less than the specified maximun error ('MaxError'), then that peak is considered a match and its name is displayed next to the peak in the upper window. When when the 'O' key is pressed (the letter 'O'), the peak positions, names, errors, and amplitudes are printed out in a table in the command window.
 
Example 8: Eleven input arguments. As above, but also specifies 'MaxError', 'Positions', and 'Names' in optional input arguments 9, 10, and 11, for peak identification function. Pressing the 'I' key toggles off and on the peak identification labels in the upper window.These data (provided in the ZIP file) are from a high-resolution atomic spectrum (x-axis in nanometers).

       
>> load ipeakdata.mat
       
>> ipeak(Sample1,0,100,0.05,3,6,296,5,0.1,Positions,Names);

Pressing "
O" prints the peak positions, names, errors, and amplitudes in a table in the command window.
    Name            Position    Error        Amplitude
    'Mg I 295.2'     [295.2]    [0.058545]    [129.27]
    'Cu 296.1 nm'    [296.1]    [0.045368]    [124.6]
    'Hg 297.6 nm'    [297.6]    [0.023142]    [143.95]

Note: The ZIP file contains the latest version of the iPeak function as well as some sample data to demonstrate peak identification (Example 8).  Obviously for your own applications, it's up to you to provide your own array of known peak maximum positions ('Positions') and matching cell array of names ('Names') for your particular types of signals.

Keyboard Controls: 
 
 Pan signal left and right...Coarse pan: < or >                            
                             Fine pan: left or right cursor arrow keys

                             Nudge
one point left or right: [ and ]

 Zoom in and out.............Coarse zoom: / or '  
                             Fine zoom: up or down cursor arrow keys
 Resets pan and zoom.........ESC
 Refresh entire plot.........Enter (Updates cursor position in lower plot)
 Change plot color...........Shift-C (cycles through standard colors)
 Adjust AmpThreshold.........A,Z  (Larger values ignore short peaks)
 Adjust SlopeThreshold.......S,X  (Larger values ignore broad peaks)
 Adjust SmoothWidth..........D,C  (Larger values ignore sharp peaks}
 Adjust FitWidth.............F,V  (Adjust to cover just top part of peaks)
 Toggle sharpen mode ........H  Helps detect overlapped peaks.
 Baseline correction.........B, then click baseline at multiple points
 Restore original signal.....G  to cancel previous background subtraction
 Invert signal...............-  Invert (negate) the signal (flip + and -)
 Set minimum to zero.........0 (Zero) Sets minimum signal to zero
 Toggle log y mode OFF/ON....Y  Plot log Y axis on lower graph 
 Toggle autozero OFF/ON......T  Auto background subtraction on upper graph
 Toggle valley mode OFF/ON...U  Switch to valley mode
 Print peak table............P  Prints Peak #, Position, Height, Width
 Save peak table.............Shift-P  Saves peak table as disc file
 Step through peaks..........Space/Tab  Jumps to next/previous peak
 Normal peak fit.............N  Fit peaks in upper window with peakfit.m
 Multiple peak fit...........M  Fit all peaks in signal with peakfit.m
 Print keyboard commands.....K  Prints this list
 Print findpeaks arguments...Q  Prints findpeaks function with arguments.
 Print ipeak arguments.......W  Prints ipeak function with all arguments.
 Print report................R  Prints Peak table and parameters
 Print peak statistics.......E  prints mean, std of peak intervals, heights, etc.
 Peak labels ON/OFF......... L  Label all peaks detected in upper window.
 Peak ID ON/OFF..............I  Identifies closest peaks in 'Names' database.
 Print peak IDs..............O  Prints table of peaks IDs


*Sunspot data downloaded from NOAA: ftp://ftp.ngdc.noaa.gov/STP/SOLAR_DATA/SUNSPOT_NUMBERS/INTERNATIONAL/yearly/YEARLY.PLT

iPeak Demo scripts

demoipeak.m is a simple demo script that generates a noisy signal with peaks, calls iPeak, and then prints out a table of the actual peak parameters and a list of the peaks detected and measured by iPeak for comparison. Before running this demo, ipeak.m must be downloaded and placed in the Matlab path. The ZIP file at http://terpconnect.umd.edu/~toh/spectrum/ipeak5.zip contains several demo scripts (ipeakdemo.m, ipeakdemo1.m, etc) that illustrate various aspects of the iPeak function and how it can be used effectively. Download the zip file, right-click and select "Extract all", then put the resulting files in the Matlab path and run them by typing their names at the Matlab command window prompt or by opening them in the Matlab editor and clicking the green "Run" button on the editor toolbar. To test for the proper installation and operation of iPeak, you can run this test script.


ipeakdemo: effect of the peak detection parameters

Four Gaussian peaks with the same heights but different widths (10, 30, 50 and 70 units) This demonstrates the effect of SlopeThreshold and SmoothWidth on peak detection. Increasing SlopeThreshold (S key) will discriminate against the broader peaks. Increasing SmoothWidth (D key) will discriminate against the narrower peaks and noise. FitWidth (F/V keys) controls the number of points around the "top part" of the (unsmoothed) peak that are taken to estimate the peak heights, positions, and widths. A reasonable value is ordinarily about equal to 1/2 of the number of data points in the half-width of the peaks. In this case, where the peak widths are quite different, set it to about 1/2 of the number of data points in the narrowest peak.
ipeakdemo1: the autozero mode

Demonstration of autozero background correction, for separated, narrow peaks on a large baseline. Each time you run this demo, you will get a different set of peaks and noise. A table of the actual peak positions, heights, widths, and areas is printed out in the command window. Jump to the next/previous peaks using the Spacebar/Tab keys.

Hint: Turn on the Autozero mode (T key), adjust the zoom setting so that the peaks are shown one at a time in the upper window, then press the P key to display the peak table. 

ipeakdemo2: peak overlap and the curve fitting functions.

Demonstration of error caused by overlapping peaks on a large offset baseline. Each time you run this demo, you will get a different set of peaks and noise. A table of the actual peak positions, heights, widths, and areas is printed out in the command window. Jump to the next/previous peaks using the Spacebar/Tab keys.

Hint: Use the B key and click on the baseline points, then press the P key to display the peak table. Or turn on the Autozero mode (T key) and use the Normal curve fit (N key) or Multiple curve fit (M key) with peak shape 1 (Gaussian). 

ipeakdemo3 Non-Gaussian peak shapes

Demonstration of overlapping Lorentzian peaks, without an added background. Overlap of peaks causes significant errors in peak height, width, and area. Jump to the next/previous peaks using the Spacebar/Tab keys. Each time you run this demo, you will get a different set of noise. A table of the actual peak positions, heights, widths, and areas is printed out in the command window.

Hint: turn OFF the Autozero mode (T key) and use the Normal curve fit (N key) or Multiple curve fit (M key) with peak shape 2 (Lorentzian). 

ipeakdemo4: dealing with very noisy signals

Detection and measurement of four peaks in a very noisy signal. The signal-to-noise ratio of first peak is 2. Each time you run this demo, you will get a different set of noise. A table of the actual peak positions, heights, widths, and areas is printed out in the command window. Jump to the next/previous peaks using the Spacebar/Tab keys. The peak at x=100 is usually detected, but the accuracy of peak parameter measurement is poor because of the low signal-to-noise ratio.

Hint: With very noisy signals it is usually best to increase SmoothWidth and FitWidth to help reduce the effect of the noise.

ipeakdemo5: dealing with highly overlapped peaks

In this demo the peaks are so highly overlapped that only one or two of the highest peaks yield distinct maxima that are detected by iPeak. The height, width, and area estimates are highly inaccurate because of this overlap. The normal peak fit function ('N' key) would be useful in this case, but it depends on iPeak for the number of peaks and for the initial guesses, and so it would fit only the peaks that were found and numbered.

To help in this case, pressing the 'H' key will activate the peak sharpen function that decreases peak width and increases peak height of all the peaks, making it easier to detect and number all the peaks for use by the peakfit function (N key). Note: peakfit fits the original unmodified peaks; the sharpening is used only to help locate the peaks to provide peakfit with suitable starting values..

Which to use: iPeak or Peakfit?  Download these Matlab demos that compare iPeak.m with Peakfit.m for signals with a few peaks and signals with many peaks and that shows how to adjust iPeak to detect broad or narrow peaks. These are self-contained demos that include all required Matlab functions. Just place them in your path and click Run or type their name at the command prompt.  Or you can download all these demos together in idemos.zip.


FindPeakSliders, for Matlab 6.5 

This is an older version that does not work reliably in more recent versions of Matlab. Development was halted in 2008 and it therefore lacks all the newer features of iPeak.  

Click to view enlarged figure
Chromatogram, before baseline correction

Click to view enlarged
                  figure

After baseline correction

Interactive findpeaks script using mouse-controlled sliders for interactive control. Requires Matlab 6.5. This can be used to determine what values of the parameters give the most reliable peak detection. Load your data set into the vectors x and y (x = independent variable, y = dependent variable), then run this m-file and adjust the sliders to change the four arguments of the peakfind function. The four sliders correspond to the arguments of the findpeaks function described above: SlopeThreshold (SlopeT), AmpThreshold (AmpT), SmoothWidth (Smooth), and FitWidth (Fit). The range of these sliders is easily changed in lines 65 - 70 to suit a wide range of data types. The BG button is used for baseline (background) subtraction: click once on the BG button, then click on the baseline at five points starting to the left of the lowest x-value and ending to the right of the highest x-value. The background will be subtracted. You can repeat as needed. (You can also change the number of baseline points by changing BaselinePoints in the function BG.m). Peak number and the estimated position, height, and width of each peak is returned in the matrix P:
>> P
P =
1.0000 45.0095 0.1328 4.0102
2.0000 90.5301 0.3655 4.7127
3.0000 94.9673 1.4075 4.0967
......

DemoFindPeakSliders

Click to view
                enlarged figure
Click to view enlarged figure
A self-contained interactive demonstration of FindPeakSliders applied to noisy synthetic data set consisting of a random number of narrow peaks superimposed on a gently curved background. Requires Matlab 6.5. Use the sliders to explore the effect of the variables SlopeThreshold (SlopeT), AmpThreshold (AmpT), SmoothWidth (Smooth), and FitWidth (Fit), and the baseline correct (BG) button. Peak number and the estimated position, height, and width of each peak is returned in the matrix P.

DemoFindPeakSliders2


Click to view enlarged figure
A self-contained interactive demonstration of FindPeakSliders applied to a data set containing four simple peaks with increasing peak height and peak width. Use this to understand the difference between the variables SlopeThreshold (SlopeT), which discriminates on the basis of peak width, and AmpThreshold (AmpT), which discriminates on the basis of peak amplitude. Peak number and the estimated position, height, and width of each peak is returned in the matrix P.
>> P
P =
1.0000 49.9096 1.9895 20.1960
2.0000 126.1034 3.0585 42.3508
3.0000 201.3947 4.0752 63.3354
4.0000 299.9242 4.9925 82.2000


FindPeakSlidersG

Click to view
                enlarged figure
Click to view enlarged figure
Click to view
                enlarged figure
Click to view enlarged figure

Interactive findpeaks script with peak-zoom feature, for predefined data in x,y. Requires Matlab 6.5. Load a typical data set into the vectors x,y, then run this m-file and adjust the sliders to determine what values of the parameters give the most reliable peak detection. Peak number and position, height, and width of each peak is returned in the matrix P. This version has an additional slider, Peak, which zooms in on any one of the detected peaks and shows the fitted top of the peak (FitWidth) as a blue curve. To see all the peaks, set the Peaks slider to zero (all the way down).
Note: If the slider ranges are not appropriate for your data, change them in lines 63-69.

DemoFindPeakSlidersG

A self-contained interactive demonstration of FindPeakSlidersG. Requires Matlab 6.5.

Video Demonstration of FindPeakSliders

This is a 26-second, 1.4 MByte video (PeakFind2.wmv) of FindPeakSliders (version 1.1), showing the effect of varying the four parameters SlopeThreshold (SlopeT), AmpThreshold (AmpT), SmoothWidth, and FitWidth (labeled "PeakWidth" in this video).

Peak Finding and Measurement Spreadsheet

The spreadsheet pictured above implements the "findpeaks" peak finding and measurement method in a spreadsheet.  The input x,y data are contained in Sheet1, column A and B, rows 8 to 263. You can Copy and Paste your own data there. The amplitude threshold and slope threshold are set in cells B3 and E3, respectively. 

Smoothing and differentiation are performed by the convolution technique used by the spreadsheets DerivativeSmoothing.xls described previously. The Smooth Width and the Fit Width are both controlled by the number of non-zero convolution coefficients in row 5, columns J through Z. (In order to compute a symmetrical first derivative, the coefficients in columns J to Q must be the negatives of  the positive coefficients in columns S to Z). The original data and the smoothed derivative are shown in the two charts on Sheet1.

To detect peaks in the data, a series of three conditions are tested for each data point in columns F, G, and H, corresponding to the three nested loops in findpeaks.m:  

  1. Is the signal greater than Amplitude Threshold? (line 45 of findpeaks.m; column F in the spreadsheet)
  2. Is there a downward directed zero crossing in the smoothed first derivative? (line 43 of findpeaks.m; column G in the spreadsheet)
  3. Is the slope of the derivative at that point greater than the Slope Threshold? (line 44 of findpeaks.m; column H in the spreadsheet)

If the answer to all three questions is yes (highlighted by blue cell coloring), a peak is registered at that point (column I), counted in column J, and assigned an index number in column K.  

For the first 10 peaks found, the x,y original unsmoothed data are copied to Sheets 2 through 11, respectively, where that segment of data is subjected to a Gaussian least-squares fit, using the same technique as GaussianLeastSquares.xls.  The best-fit Gaussian parameter results are copied back to Sheet1, in the table in columns AH-AK.  (In its present form. the spreadsheet is limited to measuring 10 peaks, although it can detect any number of peaks.  Also it is limited in Smooth Width and Fit Width by the 17-point convolution coefficients).  

The spreadsheet is available in OpenOffice (ods) and in Excel (xls) and (xlsx) formats. They are functionally equivalent and differ only in minor cosmetic aspects. An example spreadsheet, with data, is available.  A demo version, with a calculated noisy waveform that you can modify, is also available. 

A comparison of this spreadsheet to its Matlab/Octave equivalent findpeaksplot.m is instructive. On the positive side, the spreadsheet literally "spreads out" the data and the calculations spatially over a large number of cells and sheets, breaking down the discrete steps in a very graphic way.  In particular, the use of conditional formatting in columns F through K make the peak detection decision process more evident for each peak, and the least-squares sheets 2 through 11 lay out every detail of those calculations. Spreadsheet programs have many flexible and easy-to-use formatting options to make displays more attractive. On the down side, a spreadsheet as complicated as this one is far more difficult to construct than its Matlab/Octave equivalent. Much more serious, the spreadsheet is less flexible and harder to expand to larger signals and larger number of peaks. In contrast, the Matlab/Octave equivalent, while requiring some understanding of programming, is faster, much more flexible, and can easily handle signals and smooth/fit widths of any size. Moreover, because it is written as a Matlab/Octave function, it can be readily employed as an element in your own custom Matlab/Octave programs to perform even larger tasks. 

Tom O'Haver
Professor Emeritus
Department of Chemistry and Biochemistry
The University of Maryland at College Park
toh@umd.edu
http://terpconnect.umd.edu/~toh/

Last updated, May, 2013

Website created with Document made with
          KompoZer


index previous next

Number of unique visits since May 17, 2008:

counter create hit