for "A Pragmatic Introduction to Signal Processing"
The is a list of all the
functions, scripts, data files, and spreadsheets used in this
essay. You can right-click and select
"Save link as..." to download them to your
computer. There are approximately 200 Matlab/Octave m-files
(functions and demonstration scripts); place these into the
Matlab or Octave "path" so you can use them just like any other
built-in feature. (Difference
between scripts and functions). To display the built-in help for these functions and
script, type "help <name>" at the command prompt
(where "<name>" is the name of the script or function).
Items marked or on
this page are recently added within the past year or so.
Spreadsheet
files are available in .xls format for Excel, but they can
also be opened in OpenOffice Calc.
If you are unsure whether you have all the latest versions, the
simplest way to update all my functions, scripts, tools,
spreadsheets and documentation files is to download the latest site archive ZIP file (approx. 200
MBytes), then right-click on the zip file and click "Extract
all". Then list the contents of the extracted folder by
date and then
drag and drop any new or newly updated
files into a folder in your
Matlab/Octave path. The ZIP files contains all the
files used by this web site in one directory, so you can
search for them by file name or sort them by date to determine
which ones that have changed since the last time you downloaded
them.
If you try to run one of my scripts or
functions and it gives you a "missing function" error, look
for the missing item here, download it into your path, and try
again. The script testallfunctions.m
is intended to test
for the existence of all/most of the functions in this
collection. If it comes to a function that is not installed on
your system, or if one of them does not run, it will stop with
an error, alerting you of the problem. It takes about 5
minutes to run in Matlab on a contemporary PC (slower in
Octave).
Some of these functions have been requested by users, suggested
by Google search terms, or corrected and expanded based on
extensive user feedback; you could almost consider this an
international "crowd-sourced" software project. I wish to express my thanks
and appreciation for all those who have made useful
suggestions, corrected errors, and especially those who
have sent me data from their work to test my programs
on. These contributions have really helped to correct
bugs and to expand the capabilities of my programs.
Most of these shape functions take three required input arguments: the independent variable ("x") vector, the peak position, "pos", and the peak width, "wid", usually the full width at half maximum. The functions marked 'variable shape' require an additional fourth input argument that determines the exact peak shape. The sigmoidal and exponential shapes (alpha function, exponential pulse, up-sigmoid, down-sigmoid, Gompertz, FourPL, and OneMinusExp) have different input argument variables names.
Gaussian y =
gaussian(x,pos,wid)
exponentially-broadened Gaussian (variable shape)
Triangle-broadened
Gaussian (variable shape)
bifurcated
Gaussian (variable shape)
Flattened Gaussian (variable
shape)
Clipped Gaussian (variable
shape)
Lorentzian (aka 'Cauchy') y
= lorentzian(x,pos,wid)
exponentially-broadened Lorentzian (variable shape)
Clipped
Lorentzian (variable shape)
Gaussian/Lorentzian
blend (variable shape)
Voigt profile (variable
shape)
lognormal
logistic distribution (for logistic function, see up-sigmoid)
Pearson 5 (variable shape)
alpha function
exponential
pulse
sech2, a peak shape intermediate
between Gaussian and Lorentzian. See Sech2ShapeComparison.m
(graphic).
plateau (variable shape, symmetrical
product of sigmoid and down sigmoid, similar to Flattened Gaussian)
Breit-Wigner-Fano resonance (BWF)
(variable shape)
triangle
exponentially-broadened
triangle (variable shape)
Gaussian/Triangle blend (variable shape)
rectanglepeak
tsallis distribution (variable shape,
similar to Pearson 5)
up-sigmoid (logistic function or "S-shaped"). Simple upward going
sigmoid.
down-sigmoid ("Z-shaped") Simple downward going
sigmoid.
Gompertz, 3-parameter logistic, a
variable-shape sigmoidal:
y=Bo*exp(-exp((Kh*exp(1)/Bo)*(L-t)
+1))
FourPL, 4-parameter logistic, y =
maxy+(miny-maxy)./(1+(x./ip).^slope)
OneMinusExp, Asymptotic rise to
flat plateau: g =
1-exp(-wid.*(x-pos))
peakfunction.m, a function that
generates many different peak types specified by number.
modelpeaks, a function that simulates multi-peak time-series signal
data consisting of any number of peaks of the same shape. Syntax
is model= modelpeaks(x, NumPeaks, peakshape, Heights, Positions,
Widths, extra), where 'x' is the independent variable vector,
'NumPeaks' in the number of peaks, 'peakshape' is the peak shape
number, 'Heights' is the vector of peak heights, 'Positions' is
the vector of peak positions,'Widths' is the vector of peak
widths, and 'extra' is the additional shape parameter required
by the exponentially broadened, Pearson, Gaussian/Lorentzian
blend, BiGaussian and BiLorentzian shapes. Type 'help
modelpeaks'. To create noisy peaks, use one of the following
noise functions to create some random noise to add to the
modelpeaks array.
modelpeaks2, a function that simulates multi-peak time-series
signal data consisting of any number of peaks of different
shapes. Syntax is y=modelpeaks2(t, Shape, Height, Position,
Width, extra)where 'shape' is a vector of peak type numbers and
the other input arguments are the same as for modelpeaks.m. Type
'help modelpeaks2'
ShapeDemo demonstrates 16 basic peak shapes graphically, showing
the variable-shape peaks as multiple lines.
SignalGenerator.m is a script to create and plot realistic
computer-generated signal consisting of multiple peaks on a
variable baseline plus variable random noise. You may change the
lines here marked by <<< to modify the character of the
signal peaks, baseline, and noise.
stdev.m Octave and Matlab compatible standard deviation function (because the regular built-in std.m function behaves differently in Matlab and in Octave). rsd.m is the relative standard deviation (the standard deviation divided by the mean). prsd.m is the percent relative standard deviation (100 times the relative standard deviation).
PercentDifference.m A simple
function that calculates the percent difference between two
numbers or vectors, i.e. 100*(b-a)/a, where a and b can be
scalar or vector.
val2ind(x,val)
Returns the index and the value of the element of vector x that
is closest to "val". If more than one element is equally close,
returns vectors of indicies and values, Example: If x=[1 2
4 3 5 9 6 4 5 3 1], then val2ind(x,6)=7 and val2ind(x,5.1)=[5
9]. This is useful for accessing subsets of x,y data sets; for
example, the code sequence x1=7;
x2=8; irange=val2ind(x,x1):val2ind(x,x2); xx=x(irange);
yy=y(irange); plot(xx,yy) will
plot x vs y only over the range of x values from 7 to 8. See
also PeakFindingandMeasurement.htm#UsingP
halfwidth
and tenthwidth: [FWHM,slope1,slope2,hwhm1,hwhm2] = halfwidth(x,y,xo)
uses linear interpolation between points to compute the
approximate FWHM (full width at half maximum) of any smooth peak
whose maximum is at x=xo, has a zero baseline, and falls to
below one-half of the maximum height on both sides. Not accurate
if the peak is noisy or sparsely sampled. If the additional
output arguments are supplied, it also returns the leading and
trailing edge slopes, slope1 and slope2, and the leading and
trailing edge half widths at half maximum, hwhm1 and hwhm2,
respectively. If x0 is omitted, it determines the halfwidth of
the largest peak. Example: xo=500; width=100; x=1:1000;
y=exp(-1.*((x-xo)/(0.60056120439323.*width)).^2);
halfwidth(x,y,xo). The analogous function
[twidth,slope1,slope2,hwhm1,hwhm2] = tenthwidth(x,y,xo) computes
the full width at 1/10 maximum, and just for the heck of it, hundredthwidth, [hwidth,slope1,slope2] = hundredthwidth(x,y,xo),
computes the full width at 1/100 maximum.
MeasuringWidth.m is a script that compares two methods of measuring the
full width at half maximum of a peak: gaussian fitting (using peakfit.m)
and direct interpolation (using halfwidth.m). The two methods
agree exactly for a finely-sampled noiseless Gaussian on a zero
baseline, but give slightly different answers if any of these
conditions are not met. The halfwidth function works well for
any finely-sampled smooth peak shape on a zero baseline, but the
peakfit function is better at resisting random noise and it has
the ability to correct for some types of baseline and it has a
wide selection of peak shapes to use as a model. See the help
file.
IQrange.m, estimates the standard deviation of a set of numbers
by dividing its "interquartile
range" (IQR) by 1.34896,
an alternative to the usual standard deviation calculation that
works better for computing the dispersion (spread) of a data set
that contains outliers. Essentially it's the standard deviation
with outliers removed. Syntax is b = IQrange(a).
rmnan(a) "ReMoveNotANumber".
Removes NaNs ("Not a Number") and Infs ("
Infinite") from vectors, replacing with nearest real numbers and
printing out the number of changes (if any are made). Use this
to prevent subsequent operations from stopping on an error.
rmz(a)
ReMoves Zeros from vectors, replacing with
nearest non-zero numbers
and printing out the number of changes (if any are made).
Use this to remove zeros from vectors that will subsequently be
used as the denominator of a division.
[a,changes]=nht(a,b); "no higher than"
replaces any numbers in vector a that are greater than
the scalar b with b. Optionally "changes"
returns the number of changes.
[a,changes]=nlt(a,b); "no
lower than" replaces any numbers in vector a
that are less than the scalar b with b.
Optionally "changes" returns the number of changes.
makeodd(a):
Makes the elements of vector "a" the next higher odd integers.
This can be useful in computing smooth widths to insure that the
smooth will not shift the maximum of peaks. For example, makeodd([1.1
2 3 4.8 5 6 7.7 8 9]) = [1 3 3 5 5 7 9 9 9]
condense(y,n), function to reduce the length of vector y by
replacing each group of n successive values by their average.
The similar function condensem.m
works for matrices. Use to re-sample an oversampled signal.
Mentioned on Smoothing.html and iSignal.html.
val2ind(x,val),
returns the index and the value of the element of vector x that
is closest to val. Example: if x=[1 2 4 3 5 9 6 4 5 3 1], then
val2ind(x,6)=7 and val2ind(x,5.1)=[5 9]. For some examples of
how this can be used, see PeakFindingandMeasurement.htm#UsingP and InteractivePeakFitter.htm#FitResultsTable. testcondense.m is a script that demonstrates of the effect of boxcar
averaging using the condense.m function to reduce noise without
changing the noise color. Shows that it reduces the measured
noise, removing the high frequency components, resulting in a
faster fitting execution time and a lower fitting error, but no
more accurate measurement of peak parameters.
NumAT(m,threshold): "Numbers Above Threshold":
Counts the number of adjacent elements in the vector 'm' that
are greater than or equal to the scalar value 'threshold'.
Returns a matrix listing each group of adjacent values, their
starting index, the number of elements in that group, and the
sum of that group, and the mean. Type "help NumAT" and try the
example.
isOctave.m Returns 'true' if this code is being executed by Octave. Returns 'false' if this code is being executed by MATLAB, or any other MATLAB variant. Useful in those few cases where there is a small difference between the syntax or operation of Matlab and Octave functions, as for example , , and .m.
Data plotting. The Matlab/Octave scripts plotting.m and plotting2.m show how to plot multiple signals using matrices and
subplots (multiple small plots in a single Figure window). The scripts realtimeplotautoscale.m and
realtimeplotautoscale2.m demonstrating plotting in real time (animated graphic).
plotit, version 2, (previously named 'plotfit'), is an easy-to-use function for plotting x,y data in matrices or in separate vectors. Syntax: [coef,RSquared,StdDevs,BootResults] = plotit(xi,yi,polyorder,datastyle,fitstyle). It can also fit polynomials to the data and compute errors. Click here or type "help plotit" at the Matlab/Octave prompt for examples.
whitenoise, pinknoise, bluenoise propnoise, sqrtnoise, bimodal: different types of random noise that might be encountered in physical measurements. Type "help whitenoise", etc, for help and examples.
noisetest.m is a
self-contained Matlab/Octave function for
demonstrating different noise types. It plots Gaussian
peaks with four different types of added noise with
the same standard deviation: constant white noise;
constant pink (1/f) noise; proportional white noise;
and square-root white noise, then fits a Gaussian
model to each noisy data set and computes the average
and the standard deviation of the peak height,
position, width and area for each noise type. See SignalsAndNoise.html.
Se also NoiseColorTest.m.
SubtractTwoMeasurements.m
is a Matlab/Octave script demonstration of measuring
the noise and signal-to-noise ratio of a stable
waveform by subtracting two measurements of the signal
waveform, m1 and m2 and computing the standard
deviation of the difference The signal must be stable
between measurements (except for the random noise).
The standard deviation of the measured noise is given
by sqrt((std(m1-m2)^{2})/2).
NoiseColorTest.m, a
function that demonstrates the effect of smoothing
white, pink, and blue noise. It displays a graphic of
5 noise color types both before and after smoothing, as
well as their frequency
spectra. All noise samples have a standard
deviation of 1.0 before smoothing. You can change the
smooth width and type in lines 6 and 7.
CurvefitNoiseColorTest.m,
a function that demonstrates the effect of white,
pink, and blue noise on curve fitting a single
Gaussian peak.
RANDtoRANDN.m is
a script that demonstrates how the expression 1.73*(RAND()
- RAND() + RAND() - RAND()) approximates normally-distributed random
numbers with zero mean and a standard deviation of 1.
See SignalsAndNoise.html.
RoundingError.m. A
script that demonstrates digitization (rounding) noise
and shows that adding noise and then ensemble
averaging multiple signals can reduce the overall
noise in the signal. A rare example where adding noise
is actually beneficial. See CaseStudies.html#Digitization.
DigitizedSpeech.m, an
audible/graphic demonstration of rounding error on
digitized speech. It starts with an audio recording of
the spoken phrase "Testing, one, two, three",
previously recorded at 44000 Hz and saved in WAV
format (download link), rounds off the amplitude data
progressively to 8 bits (256 steps), 4 bits (16
steps), and 1 bit (2 steps), and then the same with
random white noise added before the rounding (2 steps
+ noise), plots the waveforms and plays the resulting
sounds, demonstrating both the degrading effect of
rounding and the remarkable improvement caused by
adding noise. See CaseStudies.html#Digitization.
CentralLimitDemo.m,
script that demonstrates that the more independent uniform
random variables are combined, the probability distribution
becomes closer and closer to normal (Gaussian). See SignalsAndNoise.html.
EnsembleAverageDemo.m is a
Matlab/Octave script that demonstrates ensemble averaging to
improved the signal-to-noise ratio of a very noisy signal. Click for graphic. The
script requires the "gaussian.m"
function to be downloaded and placed in the Matlab/Octave
path, or you can use any other peak
shape
function, such as lorentzian.m
or rectanglepulse.m.
EnsembleAverageDemo2.m
is a Matlab/Octave script that demonstrates the effect
of amplitude noise, frequency noise, and phase noise
on the ensemble averaging of a sine waveform.
EnsembleAverageFFT.m
is a Matlab/Octave script that demonstration of the
effect of amplitude noise, frequency noise, and phase noise
on the ensemble averaging of a sine waveform
signal. Shows that: (a) ensemble averaging reduces the
white noise in the signal but not the frequency or
phase noise, (b) ensemble averaging the Fourier
transform has the same effect as ensemble averaging
the signal itself, and (c) the effect of phase noise
is reduced if the power spectra are ensemble averaged.
EnsembleAverageFFTGaussian.m does the same for a
Gaussian peak signal, where variation in peak width is
frequency noise and variation in peak position is
phase noise.
iPeakEnsembleAverageDemo.m ia a self-contained demonstration of the iPeak function. In this example, the signal contains a repeated pattern of two overlapping Gaussian peaks of width 12, with a 2:1 height ratio. These patterns occur a random intervals, and the noise level is about 10% of the average peak height. Using iPeak's ensemble average function (Shift-E), the patterns can be averaged and the signal-to-noise ratio significantly improved. See CaseStudies.html#Ensemble.
PeriodicSignalSNR.m is a Matlab/Octave script demonstrating the estimation of the peak-to-peak and root-mean-square signal amplitude and the signal-to-noise ratio of a periodic waveform, estimating the noise by looking at the time periods where its envelope drops below a threshold. See CaseStudies.html#SNR.
iPeakEnsembleAverageDemo.m is a demonstration of iPeak's ensemble average function. In this example, the signal contains a repeated pattern of two overlapping Gaussian peaks, 12 points apart, both of width 12, with a 2:1 height ratio. These patterns occur a random intervals throughout the recorded signal, and the random noise level is about 10% of the average peak height. Using iPeak's ensemble average function (Shift-E), the patterns can be averaged and the signal-to-noise ratio significantly improved.
LowSNRdemo.m is a script that compares several different methods of peak measurement with very low signal-to-noise ratios. It creates a single peak, with adjustable shape, height, position, and width, adds constant white random noise so the signal-to-noise ratio varies from 0 to 2, then measures the peak height and position by each method and computes the average error. Four methods are compared: (1) the peak-to-peak measure of the smoothed signal and background; (2) a peak finding method based on findpeakG; (3) unconstrained iterative least-squares fitting (INLS) based on the peakfit.m function; and (4) constrained classical least squares fitting (CLS) based on the cls2.m function. See Appendix J: How Low can you Go? Performance with very low signal-to-noise ratios.
RandomWalkBaseline.m
simulates a Gaussian peak with randomly variable position and
width superimposed on a drifting "random walk" baseline.
Compare to WhiteNoiseBaseline.m.
See CaseStudies.html#RandomWalk.
AmplitudeModulation.m is a
Matlab/Octave script simulation of modulation and synchronous
detection, demonstrating the noise reduction capability. See CaseStudies.html#Modulation.
DerivativeNumericalPrecisionDemo.m.
Self-contained function that demonstrates how the numerical
precision limits of the computer effects the first through fourth
derivatives of a smooth ("noiseless") Gaussian band, showing
both the waveforms (in Figure 1) and their frequency spectra
(in Figure 2). The numerical precision limit of the computer
creates random noise at very high frequencies, which is
emphasized by differentiation, and by the fourth derivative
that noise overwhelms the signal frequencies at lower
frequencies. Smoothing with a Gaussian (three passes of a
sliding-average) smooth with a smooth ratio of 0.2 removes
most of the noise. With real experimental data, even the
tiniest amounts of noise in the original data would be much
greater than this. Used in CaseStudies.html#Numerical.
RegressionNumericalPrecisionTest.m
is a Matlab/Octave script that demonstrates how the numerical
precision limits of the computer effects the Classical Least Squares
(multilinear regression) of two very closely-spaced
"noiseless" overlapping Gaussian peaks. This uses three
different mathematical formulation of the least-squares
calculation that give different results when the numerical
precision limits of the computer are reached. But practically,
the difference between these methods is unlikely to be seen;
even the tiniest bit of added random noise (line 15) or signal
instability produces a far greater error. Used in CaseStudies.html#Numerical.
RegressionADCbitsTest.m.
Demonstration of the effect of analog-to-digital converter
resolution (defined by the number of bits in line 9) on
Classical Least Squares (multilinear regression) of two
closely-spaced overlapping Gaussian peaks. Normally, the
random noise (line 10) produces a greater error than the ADC
resolution. Used in CaseStudies.html#Numerical.
fastsmooth, versatile function for fast data smoothing. The syntax is SmoothY=fastsmooth(Y,w,type,ends). See Smoothing.html#Matlab. Note: Greg Pittam has published a modification of the fastsmooth function that tolerates NaNs (Not a Number) in the data file (nanfastsmooth(Y,w,type,tol)) and a version for smoothing angle data (nanfastsmoothAngle(Y,w,type,tol)). Click for animated example.
SegmentedSmooth.m, segmented multiple-width data smoothing function based on the fastsmooth algorithm. The syntax is SmoothY = SegmentedSmooth(Y,smoothwidths,type,ends). This function divides Y into a number of equal-length segments according to the length of the vector 'smoothwidths', then smooths each segment with a smooth of width defined by the sequential elements of vector 'smoothwidths' and smooth type 'type'. Type "help SegmentedSmooth" for examples. DemoSegmentedSmooth.m demonstrates the operation (click for graphic). See CaseStudies.html#Segmented.
medianfilter , median-based filter function for eliminating narrow spike artifacts. The syntax is mY=medianfilter(y,Width), where "Width" is the number of points in the spikes that you wish to eliminate. Type "help medianfilter" at the command prompt.
killspikes.m is a threshold-based filter function for eliminating narrow spike artifacts. The syntax is fy= killspikes(x, y, threshold, width). Each time it finds a positive or negative jump in the data between y(n) and y(n+1) that exceeds "threshold", it replaces the next "width" points of data with a linearly interpolated segment spanning x(n) to x(n+width+1), See killspikesdemo. Type "help killspikes" at the command prompt.
testcondense.m is a script that demonstrates of the effect of boxcar averaging using the condense.m function, which performs a non-overlapping boxcar averaging function, to reduce noise without changing the noise color. Shows that it reduces the measured noise, removing the high frequency components, resulting in a faster fitting execution time and a lower fitting error, but unfortunately no more accurate measurement of peak parameters.
SmoothWidthTest.m is a
Matlab/Octave script that demonstrates the effect of smoothing
on the peak height, random white noise, and signal-to-noise
ratio of a noisy peak signal. Produces an animation showing the
effect of progressively wider smooth widths, then draws a graph
of peak height, noise, and signal-to-noise ratio vs smooth
ratio. Click to see gif animation.
You can change the peak shape and width in line 8 and the smooth type
in line 9: 1=rectangle; 2=triangle;
3=pseudo Gaussian. The script requires the "gaussian.m"
function
to be downloaded and placed in the Matlab/Octave path, or you
can use any other peak shape
function, such as lorentzian.m
or rectanglepulse.m, etc.
SmoothExperiment.m, very simple
script that demonstrates the effect of smoothing on the
position, width, and height of a single Gaussian peak. Requires
that the fastsmooth.m and peakfit.m functions be present in the
path. See Smoothing.html#Matlab
smoothdemo.m compares the performance
and speed of four types of smooth
operations: (1) sliding-average, (2) triangular, (3)
pseudo-Gaussian (equivalent to three passes of a
sliding-average), and (4) Savitzky-Golay. These smooth
operations are applied to a single noisy Gaussian peak. The peak
height of the smoothed peak, the standard deviation of the
smoothed noise, and the signal-to-noise ratio are all measured
as a function of smooth width. See Smoothing.html#Matlab. The modified script, smoothdemoWavelet.m adds wavelet
denoising as method #5 (requires the Wavelet Toolkit).
SmoothOptimization.m, script that shows why you
don't need to smooth data prior to least-squares curve
fitting; it compares the effect of smoothing on the
signal-to-noise ratio of peak height of a noisy Gaussian
peak, using three different measurement methods. Requires that the fitgauss2.m,
gaussfit.m, gaussian.m,
and fminsearch.m functions be present in the path. See CurveFittingC.html#Smoothing.
SmoothVsCurvefit.m, comparison
of peak height measurement by taking the maximum of the smoothed
signal and by curve fitting the original unsmoothed data.
Requires peakfit.m and gaussian.m in path.
DemoSegmentedSmooth.m
demonstrates the operation of SegmentedSmooth.m
with a signal consisting of noisy variable-width peaks that get
progressively wider. Requires SegmentedSmooth.m and gaussian.m in the path.
DeltaTest.m. A simple
Matlab/Octave script that demonstrates the shape of any
smoothing algorithm can be determined by applying that
smooth to a delta function, a signal consisting of all zeros except for one
point.
iSignal performs several different kinds of smoothing,
segmented smoothing, median filtering, and spike removal
(as well as differentiation, peak sharpening,
least-squares measurements of peak position, height,
width, and area, signal and noise amplitudes, frequency
spectra in selected regions of the signal, and
signal-to-noise ratio of peaks). m-file link: isignal.m. Click here to download the ZIP file
"iSignal6.zip". Click
for
animated
example.
The script RealTimeSmoothTest.m
demonstrates real-time smoothing, plotting the raw
unsmoothed data as a black line and the smoothed data in
red. In this case the script pre-calculates simulated data
in line 28 and then accesses the data point-by-point in
the processing loop (lines 30-51). The total number of
data points is controlled by 'maxx' in line 17 (initially
set to 1000) and the smooth width (in points) is
controlled by 'SmoothWidth' in line 20. Animated graphic.
Data Smoothing Tool (download link: DataSmoothing.mlx) is an interactive Live Script that can apply several types of smoothing to experimental data stored on disk. It can perform spike removal, sliding average smooths with up to 5 passes, Savitsky-Golay and Fourier low-pass filtering, and wavelet denoising (which requires the Matlab Wavelet Toolkit).
Differentiation and peak sharpening
deriv, deriv2, deriv3, deriv4, derivxy and secderivxy, simple functions for
computing the derivatives of time-series data. See Differentiation.html#Matlab
SmoothDerivative.m
combines differentiation and smoothing. The syntax is
SmoothedDeriv =
SmoothedDerivative(x,y,DerivativeOrder,w,type,ends) where
'DerivativeOrder' determines the derivative order (0, 1, 2,
3, 4, 5), 'w' is the smooth width, 'type' determines the
smooth mode:
If type=0,
the signal is not smoothed.
If type=1,
rectangular (sliding-average or boxcar)
If type=2,
triangular (2 passes of sliding-average)
If type=3,
pseudo-Gaussian (3 passes of sliding-average)
If type=4,
Savitzky-Golay smooth
and 'ends' controls how the "ends" of the signal (the first
w/2 points and the last w/2 points) are handled
If ends=0,
the ends are zeroed
If ends=1,
the ends are smoothed with progressively smaller smooths the
closer to the end.
Type "help SmoothDerivative" for some examples (graphic).
SlopeAnimation.m is an animated Matlab/Octave demonstration that shows that the first derivative of a signal is the slope of the tangent to the signal at each point.
sharpen,
Resolution enhancement (peak sharpening) by the even-derivative
method. Syntax is SharpenedSignal=
sharpen(signal, factor1, factor2, SmoothWidth). See ResolutionEnhancement.html.
Related demos: SegmentedSharpen.m,
DemoSegmentedSharpen.m (graphic), SharpenedGaussianDemo.m (graphic), SharpenedGaussianDemo4terms.m
(graphic), SharpenedLorentzianDemo.m
(graphic), SharpenedLorentzianDemo4terms.m. The
script SharpenedOverlapDemo.m
(graphic) demonstrates
the effect of sharpening on perpendicular
drop area measurements of two overlapping Gaussians peaks
with adjustable height, separation, and width, calculating the
percent different between the area measured on the overlapping
peak signal compared to the true areas of the isolated peaks. SharpenedOverlapDemo.m is a
script that automatically determines the optimum degree of
even-derivative sharpening that minimizes the errors of
measuring peak areas of two overlapping Gaussians by the
perpendicular drop method using the autopeaks.m
function. (Graphic 1, Graphic
2).
ProcessSignal, a Matlab/Octave command-line function that performs smoothing, See differentiation, peak sharpening, and median filtering on the time-series data set x,y (column or row vectors). Similar to iSignal, without the plotting and interactive keystroke controls. 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).
derivdemo1.m, a function that demonstrates the basic shapes of derivatives. See Differentiation.html#BasicProperties
DerivativeShapeDemo.m is a function that demonstrates the first derivatives of 16 different peak shapes. (graphic)
derivdemo2.m, a function that demonstrates the effect of peak width on the amplitude of derivatives. See Differentiation.html#BasicProperties
derivdemo3.m, a function that demonstrates the effect of smoothing on the first derivative of a noisy signal. See Differentiation.html#Smoothing
derivdemo4.m, a function that demonstrates the effect of smoothing on the second derivative of a noisy signal. See Differentiation.html#Smoothing
DerivativeDemo.m is a
self-contained Matlab/Octave demo function that uses ProcessSignal.m and plotit.m to demonstrate an application
of differentiation to the quantitative analysis of a peak
buried in an unstable background (e.g. as in various forms of
spectroscopy). The object is to derive a measure of peak
amplitude that varies linearly with the actual peak amplitude
and is minimally effected by the background and the noise. To
run it, just type DerivativeDemo at the command prompt. You
can change several of the internal variables (e.g. Noise,
BackgroundAmplitude) to make the problem harder or easier.
Note that, despite the fact that the magnitude of the
derivative is numerically smaller than the original signal
(because it has different units), the signal-to-noise ratio of
the derivative is better, and the derivative signal is
linearly proportional to the actual peak height, despite the
interference of large background variations and random noise.
See Differentiation.html
iSignal is an
interactive function for Matlab that performs differentiation, smoothing,
and peak sharpening for
time-series signals, up to the 5^{th}
derivative, automatically including the
required type of smoothing. Simple keystrokes allow
you to adjust the smoothing and sharpening parameters
while observing the effect on your signal dynamically.
Click
here to download the ZIP file "iSignal7.zip". Click
for
animated
example. Version 7, June 2019, includes the
ability to symmetrize exponentially broadened peaks
peaks by weighted first derivative addition.
demoisignal.m
for Matlab is a self-running script that demonstrates the
features of iSignal (and requires
that the latest version of iSignal, and version 6 of plotit.m, be present in your Matlab
path). Demonstrates panning and zooming, smoothing,
differentiation, frequency spectrum, peak measurement, and
derivative spectroscopy calibration (in conjunction with plotit.m version 6).
iSignalDeltaTest
is a Matlab/Octave script that demonstrates
the frequency response (power spectrum) of the
smoothing and differentiation functions of iSignal by a applying them
to a delta function. Change the
smooth type, smooth width, and derivative order and
see how the power spectrum changes.
The script RealTimeSmoothFirstDerivative.m demonstrates real-time smoothed differentiation, using a simple adjacent-difference algorithm (line 47) and plotting the raw data as a black line and the first derivative data in red. The script RealTimeSmoothSecondDerivative.m computes the smoothed second derivative by using a central difference algorithm (line 47). Both of these scripts pre-calculate the simulated data in line 28 and then accesses the data point-by-point in the processing loop (lines 31-52). In both cases the maximum number of points is set in line 17 and the smooth width is set in line 20.
The
script RealTimePeakSharpening.m
demonstrates real-time peak
sharpening using the second
derivative technique. It uses
pre-calculated simulated data in
line 30 and then accesses the data
point-by-point in the processing
loop (lines 33-55). In both cases
the maximum number of points is set
in line 17 and the smooth width is
set in line 20 and the weighting
factor (K1) is set in line 21. In
this example the smooth width is 101
points, which accounts for the the
delay in the sharpened peak compared
to the original.
FrequencySpectrum.m (syntax fs=FrequencySpectrum(x,y)) returns a matrix containing the real part of the Fourier power spectrum of x,y. To plot it, type "plotit(fs)". Input data can be separate x and y vectors of both in a matrix. Type "help FrequencySpectrum".
PlotFrequencySpectrum.m
plots the frequency spectrum or periodogram of the signal x,y on
linear or log coordinates. The syntax is PowerSpectrum=
PlotFrequencySpectrum(x, y, plotmode, XMODE, LabelPeaks).
Type "help PlotFrequencySpectrum" for details. Try
this example:
x= [0:.01:2*pi]';
y=sin(200*x)+randn(size(x));
subplot(2,1,1); plot(x,y); subplot(2,1,2);
PowerSpectrum=PlotFrequencySpectrum(x,y,1,0,1);
CompareFrequencySpectrum.m.
A script that compares two signals (upper panel) and their
frequency spectra (lower panel) with the original signal shown
in blue and the modified signal in green. plotmode: =1:linear,
=2:semilog X, =3:semilog Y; =4: log-log). XMODE: =0 for
frequency Spectrum (x is frequency); =1 for periodogram (x is
time). Define the signal modification in line 15. You can load a
signal stored in .mat format or create an simulated signal for
testing. Must have PlotFrequencySpectrum.m
in the path.
PlotSegFreqSpect.m plots and segmented Fourier spectrum (syntax PSM=(x,y,NumSegments,MaxHarmonic,LogMode)). It breaks y into 'NumSegments' equal length segments, multiplies each by an apodizing window function, computes the power spectrum of each segment, returns the power spectrum matrix (PSM), and plots the result of the first 'MaxHarmonic' Fourier components as as a contour plot. See HarmonicAnalysis.html#Software
iSignalDeltaTest is a Matlab/Octave script that demonstrates the frequency response (power spectrum) of the smoothing and differentiation functions of iSignal by a applying them to a delta function. Change the smooth type, smooth width, and derivative order and see how the power spectrum changes.SineToDelta.m. A demonstration animation (animated graphic) showing the waveform and the power spectrum of a rectangular pulsed sine wave of variable duration (whose power spectrum is a "sinc" function) changing continuously from a pure sine wave at one extreme (where its power spectrum is a delta function) to a single-point pulse at the other extreme (where its power spectrum is a flat line). GaussianSineToDelta.m is similar, except that it shows a Gaussian pulsed sine wave, whose power spectrum is a Gaussian function, but which is the same at the two extremes of pulse duration (animated graphic).
iSignal is a multi-purpose interactive signal processing tool (for Matlab only) that includes a Frequency Spectrum mode, toggled on and off by the Shift-S key; it computes frequency spectrum of the segment of the signal displayed in the upper window and displays it in the lower window (in red). You can use the pan and zoom keys to adjust the region of the signal to be viewed or press Ctrl-A to select the entire signal. Press Shift-S again to return to the normal mode. See CaseStudies.html#Harmonic for a relevant example. Click for animated example.
iPower, a keyboard-controlled interactive power spectrum
demonstrator, useful for teaching and learning about the power
spectra of different types of signals and the effect of signal
duration and sampling rate. Single keystrokes allow you to
select the type of signal (12 different signals included), the
total duration of the signal, the sampling rate, and the global
variables f1 and f2 which are used in different ways in the
different signals. When the Enter key is pressed, the
signal (y) is sent to the Windows WAVE audio device. Press K to
see a list of all the keyboard commands. (m-file link: ipower.m).
Slideshow of examples.
The
script RealTimeFrequencySpectrumWindow.m
computes and plots the Fourier frequency spectrum of a
signal. It loads
the simulated real-time data from a ".mat file" (in
line 31) and then
accesses that data point-by-point in the processing
'for' loop. A critical variable in this case is
"WindowWidth" (line 37),
the number of data points taken to compute each
frequency spectrum. If the data stream is an
audio signal, it's also possible to play the sound
through the computer's sound system synchronized
with the display of the frequency spectra (set
"PlaySound" to 1).
ExpBroaden, exponential broadening function. Syntax is yb = ExpBroaden(y,t). Convolutes the vector y with an exponential decay of time constant t. Mentioned on SignalsAndNoise.html and InteractivePeakFitter.htm.
GaussConvDemo.m, a script that demonstrates that a Gaussian of unit height, Fourier convoluted with a centered Gaussian of the same width is a shorter, broader Gaussian (with a height of 1/sqrt(2) and a width of sqrt(2) and of equal area to the original Gaussian). The bottom panel shows an attempt to recover the original y from the convoluted result. You can optionally add noise to show how convolution smooths the noise and how Fourier deconvolution restores it. Requires gaussian.m in path.
CombinedDerivativesAndSmooths.txt. Convolution coefficients for computing the first through fourth derivatives, with rectangular, triangular and Gaussian smooths.
Convolution.txt, simple examples of whole-number convolution vectors for smoothing and differentiation.
deconvolutionexample.m, a
simple example script that demonstrates the use of the Matlab
Fourier deconvolution 'deconv' function. See Deconvolution.html.
DeconvDemo.m, a Fourier deconvolution
demo script with a signal containing four Gaussians broadened by
an exponential function (graphic).
DeconvDemo2.m is a similar script
for a single Gaussian (graphic). DeconvDemo3.m demonstrates
deconvolution of a Gaussian convolution function from a rectangular pulse (animated graphic). DeconvDemo4.m (animated graphic) demonstrates "self deconvolution" applied a
signal consisting of a Gaussian peak that is broadened
by the measuring instrument, and an attempt to recover
the original peak width. DeconvDemo5.m (graphic) shows an attempt to resolve two closely-spaced underlying
peaks that are completely unresolved in the observed signal. See CaseStudies.html#deconvolution.
Variation of this include versions with Lorentzian
peaks and one with a triangular
convolution function.
deconvgauss.m. function ydc=deconvgauss(x,y,w) deconvolutes a Gaussian function of width 'w' from vector y, returning the deconvoluted result.
deconvexp.m. function ydc=deconvexp(y,tc) deconvolutes an exponential function of time constant 'tc' from vector y, returning the deconvoluted result.
SegExpDeconv(x,y,tc)
is a segmented version of deconvexp.m;
it 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. Useful when the peak width and/or exponential tailing
of peaks varies across the signal duration. 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. 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.
iSignal
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). Click here to download the ZIP file
"iSignal6.zip"
Data convolution tool. The
interactive Live Script DeconvoluteData.mlx can perform Fourier
self-deconvolution on you own data stored in disk.
FouFilter.m, Fourier filter
function, with variable band-pass, low-pass, high-pass, or
notch (band reject). The syntax is ry=FouFilter(y,
samplingtime, centerfrequency, frequencywidth, shape,
mode. Version 2, March 2019. See FourierFilter.html. The
scripts TestFouFilter.m
and TestFouFilter2.m
demonstrate the use of this function.
SegmentedFouFilter.m
is a segmented version of FouFilter.m, which
applies different center frequencies and widths to
different segments of the signal. The syntax is the same
as FouFilter.m except that the two input arguments
"centerFrequency" and "filterWidth" must be vectors with
the values of centerFrequency of filterWidth for each
segment. The signal is filtered for each value of centerFrequency of
filterWidth, and the
output is constructed by concatenating all the
filtered segments,
divided
equally into a number of segments determined by the
length of centerFrequency and filterWidth, which must be
equal in length. Type help SegmentedFouFilter for
help and examples.
iFilter,
interactive Fourier filter. (m-file link: ifilter.m), which uses the pan
and zoom keys to control the center frequency and the
filter width. Click
here
for
animated example. Select from
low-pass, high-pass, band-pass, band-reject, harmonic
comb-pass, or harmonic comb-reject filters. Click here to watch or download an
mp4 video of iFilter filtering a noisy Morse code
signal, with sound (watch the title of the figure as the
video plays).
MorseCode.m is a script that
uses iFilter to demonstrate the abilities and limitations
of Fourier filtering. It creates a pulsed fixed frequency
sine wave that spells out "SOS" in Morse code
(dit-dit-dit/dah-dah-dah/dit-dit-dit), adds random white
noise so that the SNR is very poor (about 0.1 in this
example), then uses a Fourier bandpass filter tuned to the
signal frequency, in an attempt to isolate the signal from
the noise. As the bandwidth is reduced, the
signal-to-noise ratio begins to improve and the signal
emerges from the noise until it becomes clear, but if the
bandwidth is too narrow, the step response time is too
slow to give distinct "dits" and "dahs". Use the ? and "
keys to adjust the bandwidth. (The step response time is
inversely proportional to the bandwidth). Press 'P' or the
Spacebar to hear the sound. You must install iFilter.m in the Matlab path. Click here to watch or download an
mp4 video of this script in operation, with sound
(watch the title of the figure as the video plays), or watch on YouTube.
TestingOneTwoThree.wav
is a 1.58 sec
duration audio recording of the spoken phrase "Testing,
one, two, three", recorded at a sampling rate of 44000
Hz and saved in WAV format. When loaded into iFilter
(v=wavread('TestingOneTwoThree.wav');) set to bandpass mode and tuned to a narrow
segment that is well above the frequency range of most
of the signal, it might seem as if though this passband
would miss most of the frequency components in the
signal, yet even in this case the speech is
intelligible, demonstrating the remarkable ability of
the ear-brain system to make do with a highly
compromised signal. Press P or space to hear the
filter's output. Different filter settings will change
the timbre
of the sound. See InteractiveFourierFilter.htm#SpeechExample;
Graphic.
IntegrationTest.m is a
digital simulation of the effect of sampling rate (data
density) on the accuracy of peak area measurements
for single isolated sparsely sampled Gaussian peaks. It
shows that the trapezoidal method is surprisingly accurate
and requires a minimum of only 2.5 points in the base width
of the peak as previously predicted by reference 71.
PerpDropAreas.m
[AreaVector]=PerpDropAreas(x,y,startx,endx,MaxVector)
measures he peak areas of the peaks in x, y, starting
an x value of startX and ending at endX, with
specified peak positions in the vector MaxVector,
which can be of any length. Uses the halfway point
method. Returns the areas in the vector PDMeasAreas
and the midpoint indices in the optional second output
argument. See https://terpconnect.umd.edu/~toh/spectrum/Integration.html#Matlab
HeightAndArea.m is a demonstration script that uses measurepeaks.m to measure the peaks in computer-generated signals consisting of a series of Gaussian peaks with gradually increasing widths that are superimposed in a curved baseline plus random white noise. It plots the signal and the individual peaks and compares the actual peak position, heights, and areas of each peak to those measured by measurepeaks.m using the absolute peak height, peak-valley difference, perpendicular drop, and tangent skim methods. Prints out a table of the relative percent difference between the actual and measured values for each peak and the average error for all peaks.
measurepeaks.m automatically
detects peaks in a signal, similar to findpeaksSG. It
returns a table of
peak number, position, absolute peak height, peak-valley
difference, perpendicular drop area, and tangent skim area
of each peak. It can plot
the signal and the individual
peaks if the last (7^{th}) input argument
is 1. Type "help measurepeaks" and try the seven examples
there, or run HeightAndArea.m to run a test of the
accuracy of peak height and area measurement with signals
that have multiple peaks with noise, background, and some
peak overlap. The script testmeasurepeaks.m
will run all of the examples with a 1-second pause between
each (requires measurepeaks.m and gaussian.m in the path).
The
script SharpenedOverlapDemo.m
(graphic) demonstrates
the effect of sharpening on perpendicular
drop area measurements of two overlapping Gaussians peaks
with adjustable height, separation, and width, calculating the
percent different between the area measured on the overlapping
peak signal compared to the true areas of the isolated peaks.
ComparePDAreas.m compares
the
effect of digital processing on the areas of a set of
peaks measured by the perpendicular drop method. Syntax is [P1,P2,coef,R2] =
ComparePDAreas(x,orig,processed,PeakSensitivity), where
x=independent variable (e.g. time); orig = original signal
y values; processed = processed signal y values; P1 = peak
table of original signal; P2 = peak table of processed
signal; PeakSensitivity = approximate number of peaks that
would fit into the entire x-axis range (larger numbers
> more peak detected). Displays a scatter plot of
original areas vs processed areas for each peak and
returns the peak tables, P1 and P2 respectively, and the
slope, intercept, and R2 values, which should ideally be
1,0, and 1, if the processing has no effect at all on peak
area.
SharpenedOverlapCalibrationCurve.m
is a script that simulates the calibration curve of three
overlapping Gaussian peaks. Even-derivative sharpening
(the red line in the signal plots) is used to improve the
resolution of the peaks to allow perpendicular drop area
measurement. A straight line is fit to the calibration
curve and the R2 is calculated, in order to demonstrate
(1) the linearity of the response, and (2) the
independence of the overlapping adjacent peaks. Must have
gaussian.m, derivxy.m, autopeaks.m, val2ind.m,
halfwidth.m, fastsmooth.m, and plotit.m in the path.
iSignal is a downloadable Matlab function that performs various signal processing functions described in this tutorial, including one-at-a-time manual measurement of peak area using Simpson's Rule and the perpendicular drop method. Click to view or right-click > Save link as... here, or you can download the ZIP file with sample data for testing. The animated GIF iSignalAreaAnimation.gif (click to view) shows iSignal applying the perpendicular drop method to a series of four peaks of equal area. (Look at the bottom panel to see how the measurement intervals, marked by the vertical dotted magenta lines, are positioned at the valley minimum on either side of each of the four peaks). It also has a built-in peak fitter, activated by the Shift-F key, based on peakfit.m, that measures the areas of overlapping peak of known shape. There is also an automatic peak finding function based on the autopeaks function, activated by the J or Shift-J keys, which displays a table of peak number, position, absolute peak height, peak-valley difference, perpendicular drop area, and tangent skim area of each peak in the signal.
peakfit, a command-line function for multiple peak fitting by iterative non-linear least-squares.It measures the peak position, height, width, and area of overlapping peaks, and it has several ways to correct for non-zero baselines. For best results, it requires that the peak shape of your peaks be among those listed here.
PeakCalibrationCurve.m is an Matlab/Octave simulation of the calibration of a flow injection or chromatography system that produces signal peaks that are related to an underlying concentration or amplitude ('amp'). The measurepeaks.m function is used to determine the absolute peak height, peak-valley difference, perpendicular drop area, and tangent skim area. The Matlab/Octave script PeakShapeAnalyticalCurve.m shows that, for a single isolated peak whose shape is constant and independent of concentration, if the wrong model shape is used, the peak heights measured by curve fitting will be inaccurate, but that error will be exactly the same for the unknown samples and the known calibration standards, so the error will "cancel out" and the measured concentrations will be accurate, provided you use the same inaccurate model for both the known standards and the unknown samples. See CaseStudies.html#Calibration.
PowerTransformTest.m
is a simple script that demonstrates the power method
of peak sharpening to aid in reducing in peak overlap.
The scripts PowerMethodGaussian.m
and PowerMethodLorentzian.m
compare the power methods to deconvolution, for Gaussian
and Lorentzian peak, respectively. PowerMethodCalibrationCurve
is a variant of PeakCalibrationCurve.m
that evaluates the power method
in the context of a flow injection or chromatography
measurement. The self-contained function PowerMethodDemo.m
demonstrates the power method for measuring the area of
small shouldering peak that is partly overlapped by a
much stronger interfering peak (Graphic). It also
demonstrates the effect of random noise, smoothing, and
any uncorrected background under the peaks.
SumOfAreas.m. Demonstrates
that even drastically non-Gaussian peaks can be fit with
up to five overlapping Gaussian components, and that the
total area of the components approaches the area under
the non-Gaussian peak as the number of components
increases (graphic). In
most cases only a few components are necessary to obtain
a good estimate of the peak area.
TestLinearFit
effect
of
number of points.txt. Effect of sample size on
least-square error estimates by Monte Carlo Simulation,
Algebraic propagation-of-errors, and the bootstrap
method, using the Matlab script TestLinearFit.m.
LeastSquaresCode.txt.
Simple pseudocode for calculating the first-order
least-square fit of y vs x, including the Slope and
Intercept and the predicted standard deviation of the
slope (SDslope) and intercept (SDintercept).
CalibrationQuadraticEquations.txt.
Simple pseudocode for calculating the second-order
least-square fit of y vs x, including the constant, x,
and x2 terms.
plotit, version 2, (previously named 'plotfit'), is a function for plotting x,y data in matrices or in separate vectors. It optionally fits the data with a polynomial of order n if n is included as the third input argument. In version 6 the syntax is [coef, RSquared, StdDevs] = plotit(x,y) or plotit(x,y,n) or optionally plotit(x, y, n, datastyle, fitstyle), where datastyle and fitstyle are optional strings specifying the line and symbol style and color, in standard Matlab convention. For example, plotit(x,y,3,'or','-g') plots the data as red circles and the fit as a green solid line (the default is red dots and a blue line, respectively). Plotit returns the best-fit coefficients 'coeff', in decreasing powers of x, the standard deviations of those coefficients 'StdDevs' in the same order, and the R-squared value. Type "help plotit" at the command prompt for syntax options. See CurveFitting.html#Matlab. This function works in Matlab or Octave and has a built-in bootstrap routine that computes coefficient error estimates (STD and % RSD of each coefficient) by the bootstrap method and returns the results in the matrix "BootResults" (of size 5 x polyorder+1). The calculation is triggered by including a 4^{th} output argument, e.g. [coef, RSquared, StdDevs, BootResults]= plotit(x,y,polyorder). This works for any positive integer polynomial order. The variation plotfita animates the bootstrap process for instructional purposes. The variation logplotfit plots and fits log(x) vs log(y), for data that follows a power law relationship or that covers a very wide numerical range.
RSquared.m Computes the R2 (Rsquared or correlation coefficient) in both Matlab and Octave. Syntax RS=RSquared(polycoeff,x,y).
trypoly(x,y) fits the data in x,y with a series of polynomials of degree 1 through length(x)-1 and returns the coefficients of determination (R^{2}) of each fit as a vector, allowing you to evaluate how polynomials of various orders fit your data. To plot as a bar graph, write bar(trypoly(x,y)); xlabel('Polynomial Order'); ylabel('Coefficient of Determination (R2)'). Here's an example. See related function testnumpeaks.m.
trydatatrans(x,y,polyorder) tries 8 different simple data transformations on the data x,y, fits the transformed data to a polynomial of order 'polyorder', displays results graphically in 3 x 3 array of small plots and returns all the R^{2} values in a vector.
LinearFiMC.m,
a script that compares standard deviation of slope and intercept
for a first-order least-squares fit computed by random-number
simulation of 1000 repeats to predictions made by closed-form
algebraic equations. See CurveFitting.html#Reliability
TestLinearFit.m, a script that
compares standard deviation of slope and intercept for a
first-order least-squares fit computed by random-number
simulation of 1000 repeats to predictions made by closed-form
algebraic equations and to the bootstrap sampling method.
Several different noise models can be selected by
commenting/uncommenting the code in lines 20-26. See CurveFitting.html#Reliability
GaussFitMC.m, a function that
demonstrates Monte Carlo simulation of the measurement of the
peak height, position, and width of a noisy x,y Gaussian peak.
See CurveFitting.html
GaussFitMC2.m, a function that
demonstrates measurement of the peak height, position, and width
of a noisy x,y Gaussian peak, comparing the gaussfit parabolic
fit to the fitgaussian iterative fit. See CurveFitting.html
SandPfrom1950.mat is a MAT file containing the daily value of the S&P 500 stock market index vs time from 1950 through September of 2016. These data are used by FitSandP.m a Matlab/Octave script that performs a least-squares fit of the compound interest equation to the daily value, V, of the S&P 500 stock market index vs time, T, from 1950 through September of 2016, by two methods: (1) the iterative curve fitting method, and (2) by taking the logarithm of the values and fitting those to a straight line. SnPsimulation.m. Matlab/Octave script that simulates the S&P 500 stock market index by adding proportional random noise to data calculated by the compound interest equation with a known annual percent return, then fits the equation to that noisy synthetic data by the two methods above. See CaseStudies.html#StockMarket
gaussfit.m function [Height, Position, Width]=gaussfit(x,y). Takes the natural log of y, fits a parabola (quadratic) to the (x,ln(y)) data, then calculates the position, width, and height of the Gaussian from the three coefficients of the quadratic fit.
lorentzfit.m function [Height, Position, Width]=lorentzfit(x,y). Takes the reciprocal of y, fits a parabola (quadratic) to the (x,1/y) data, then calculates the position, width, and height of the Lorentzian from the three coefficients of the quadratic fit.
OverlappingPeaks.m is a demo script that shows how to use gaussfit.m to measure two overlapping partially gaussian peaks. It requires careful selection of the optimum data regions around the top of each peak (lines 15 and 16). Try changing the relative position and height of the second peak or adding noise (line 3) and see how it effects the accuracy. This function needs the gaussian.m and gaussfit.m functions in the path. Iterative methods work much better in such cases.
allpeaks.m. allpeaks(x,y) A super-simple peak
detector for x,y, data sets that lists every y value that has
lower y values on both
sides; allvalleys.m is the
same for valleys, lists every y value that has higher y values on both
sides.
peaksat.m. (Peaks Above
Threshold) P=peaksat(x,y,threshold) Lists every y value that (a) has lower y
values on both sides and (b) is above the specified
threshold. Returns
a 2 by n matrix P with the x and y values of
each peak, where n is the number of detected peaks.
findpeaksx.m, P=findpeaksx(x,y, SlopeThreshold,
AmpThreshold,SmoothWidth, FitWidth, smoothtype) is a simple
command-line function to locate and count the positive
peaks in a noisy data sets. It's an alternative to the
findpeaks function in the Signal Processing Toolkit. It
detects peaks by looking for downward zero-crossings in
the smoothed first derivative that exceed SlopeThreshold
and peak amplitudes that exceed AmpThreshold, and
returns a list (in matrix P) containing the peak number
and the position and height of each peak. It can find
and count over 10,000 peaks per second in very large
signals. Type "help findpeaksx.m". See PeakFindingandMeasurement.htm.
The
variant findpeaksxw.m
additionally measures the width
of the peaks. See the demonstration script demofindpeaksxw.m.
findpeaksG.m and findvalleys.m, automatically find the peaks or valleys in a signal and and measure their position, height, width and area by curve fitting. The syntax is P= findpeaksG(x, y, SlopeThreshold, AmpThreshold, SmoothWidth, FitWidth, smoothtype). Returns a matrix containing the peak parameters for each detected peak. For peak of Lorentzian shape, use findpeaksL.m instead. See PeakFindingandMeasurement.htm#findpeaks.
findpeaksplot.m is a simple variant of findpeaksG.m that also plots the x,y data and numbers the peaks on the graph (if any are found). Syntax: findpeaksplot(x, y, SlopeThreshold, AmpThreshold, SmoothWidth, FitWidth, smoothtype)
OnePeakOrTwo.m is a demo script that creates a signal that might be interpreted as either one peak at x=3 on a curved baseline or as two peaks at x=.5 and x=3, depending on context. In this demo, the findpeaksG.m function used called twice, with two different values of SlopeThreshold to demonstrate.
iPeak automatically finds and measures multiple peaks in a signal. (m-file link: ipeak.m). Check out the Animated step-by-step instructions. The ZIP file ipeak7.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. testipeak.m is a script that tests for the proper installation and operation of iPeak by running quickly through all eight examples and six demos for iPeak. Assumes that ipeakdata.mat has been loaded into the Matlab workspace. Click for slideshow of examples. The syntax is P=ipeak(DataMatrix, PeakD, AmpT, SlopeT, SmoothW, FitW, xcenter, xrange, MaxError, positions, names)
findpeaksSG.m is a segmented variant
of findpeaksG with the same syntax,
except that the peak detection parameters can be vectors, dividing up the
signal into regions optimized for peaks of different widths. The
syntax is P =
findpeaksSG(x, y, SlopeThreshold, AmpThreshold, smoothwidth,
peakgroup, smoothtype). This works
better than findpeaksG when the peak widths vary greatly over
the duration of the signal. The script TestPrecisionFindpeaksSG.m
demonstrates the application. Graphic. See CaseStudies.html#Segmented.
findpeaksSGw.m
autofindpeaks.m (and autofindpeaksplot.m) are similar to findpeaksSG.m except that you can leave out the peak detection parameters and just write "autofindpeaks(x,y)" or "autofindpeaks(x,y,n)" where n is the peak capacity, roughly the number of peaks that would fit into that signal record (greater n looks for many narrow peaks; smaller n looks for fewer wider peaks). It also prints out the input argument list for use with any of the findpeaks... functions. In version 1.1, you can call autofindpeaks with the output arguments [P,A] and it returns the calculated peak detection parameters as a 4-element row vector A, which you can then pass on to other functions such as measurepeaks, effectively giving that function the ability to calculate the peak detection parameters from a single single number n (for example: x=[0:.1:50];y=5+5.*sin(x)+randn(size(x)); [P,A]=autofindpeaks(x,y,3); P=measurepeaks(x,y,A(1),A(2),A(3),A(4),1);). Type "help autofindpeaks" and run the examples. The script testautofindpeaks.m runs all the examples in the help file, additionally plotting the data and numbering the peaks (like autofindpeaksplot.m). Graphic animation.
[M,A]=autopeaks.m and autopeaksplot.m. Peak detection and height and area measurement for peaks of arbitrary shape in x,y time series data. The syntax is [P,DetectionParameters] = autofindpeaks(x, y, SlopeThreshold, AmpThreshold, smoothwidth, peakgroup, smoothtype), but similar to autofindpeaks.m, the peak detection parameters SlopeThreshold, AmpThreshold, smoothwidth peakgroup, smoothtype can be omitted and the function will calculate estimated initial values. Uses the measurepeaks.m algorithm for measurement, returning a table in the matrix M containing the peak number, position, absolute peak height, peak-valley difference, perpendicular drop area, and tangent skim area of each peak. Optionally returns the peak detection parameters that it calculates in the vector A. Using the simple syntax M=autopeaks(x,y) works well in some cases, but if not try M=autopeaks(x,y,n), using different values of n (roughly the number of peaks that would fit into the signal record) until it detects the peaks that you want to measure. For the most precise control over peak detection, you can specify all the peak detection parameters by typing M=autopeaks(x,y, SlopeThreshold, AmpThreshold, smoothwidth, peakgroup). autopeaksplot.m is the same but it also plots the signal and the individual peaks (in blue) with the maximum (red circles), valley points (magenta), and tangent lines (cyan) marked. The script testautopeaks.m runs all the examples in the autopeaks help file, with a 1-second pause between each one, printing out results in the command window and additionally plotting and numbering the peaks (Figure window 1) and each individual peak (Figure window 2); it requires gaussian.m and fastsmooth.m in the path. iSignal has a peak finding function based on the autopeaks function, activated by the J or Shift-J keys, which displays a table of peak number, position, absolute peak height, peak-valley difference, perpendicular drop area, and tangent skim area of each peak in the signal.
findpeaksG2d.m is a variant of findpeaksSG that can be used to locate the positive peaks and shoulders in a noisy x-y time series data set. Detects peaks in the negative of the second derivative of the signal, by looking for downward slopes in the third derivative that exceed SlopeThreshold. See TestFindpeaksG2d.m. Syntax: P = findpeaksG2d(x, y, SlopeThreshold, AmpThreshold, smoothwidth, peakgroup, smoothtype)
measurepeaks.m automatically detects peaks in a signal, similar to findpeaksSG. M = measurepeaks(x, y, SlopeThreshold, AmpThreshold, SmoothWidth, FitWidth, plots). It returns a table M of peak number, position, absolute peak height, peak-valley difference, perpendicular drop area, and tangent skim area of each peak. It can plot the signal and the individual peaks if the last (7^{th}) input argument is 1. Type "help measurepeaks" and try the seven examples there, or run HeightAndArea.m to run a test of the accuracy of peak height and area measurement with signals that have multiple peaks with noise, background, and some peak overlap. Generally, its values for perpendicular drop area are best for peaks that have no background, even if they are slightly overlapped, whereas its values for tangential skim area are better for isolated peaks on a straight or slightly curved background. Note: this function uses smoothing (specified by the SmoothWidth input argument) only for peak detection; it performs measurements on the raw unsmoothed y data. In some cases it may be beneficial to smooth the y data yourself before calling measurepeaks.m, using any smooth function of your choice. The script testmeasurepeaks.m will run all of the examples in the measurepeaks help file with a 1-second pause between each (requires measurepeaks.m and gaussian.m in the path). Graphic animation.
findpeaksT.m and findpeaksTplot.m are variants of findpeaks that measure the peak parameters by constructing a triangle around each peak with sides tangent to the sides of the peak. Graphic example.
findpeaksb.m is a variant of findpeaksG.m that more accurately measures peak parameters by using iterative least-square curve fitting based on peakfit.m. This yields better peak parameter values that findpeaks alone, because it fits the entire peak, not just the top part, and because it has provision for 33 different peak shapes and for background subtraction (linear or quadratic). Works best with isolated peaks that do not overlap. Syntax is P = findpeaksb(x,y, SlopeThreshold, AmpThreshold, smoothwidth, peakgroup, smoothtype, window, PeakShape, extra, AUTOZERO). The first seven input arguments are exactly the same as for the findpeaksG.m function; if you have been using findpeaks or iPeak to find and measure peaks in your signals, you can use those same input argument values for findpeaksb.m. The demonstration script DemoFindPeaksb.m shows how findpeaksb3 works with multiple overlapping peaks . Type "help findpeaksb" at the command prompt. See PeakFindingandMeasurement.htm. Compare this to the related findpeaksfit.m and findpeaksb3, next. Click for slideshow of examples.
findpeaksSb.m is a segmented variant of findpeaksb.m. It has the same syntax as findpeaksb.m, P = findpeaksb(x, y, SlopeThreshold, AmpThreshold, smoothwidth, peakgroup, smoothtype, window, PeakShape, extra, NumTrials, AUTOZERO), except that the input arguments SlopeThreshold, AmpThreshold, smoothwidth, peakgroup, window, width, PeakShape, extra, NumTrials, autozero, and fixedparameters, can all optionally be scalars or vectors with one entry for each segment, in the same manner as findpeaksSG.m. Returns a matrix P listing the peak number, position, height, width, area, percent fitting error and "R2" of each detected peak. DemoFindPeaksSb.m demonstrates this function by creating a series of Gaussian peaks whose widths increase by a factor of 25-fold and that are superimposed in a curved baseline with random white noise that increases gradually; four segments are used, changing the peak detection and curve fitting values so that all the peaks are measured accurately. Graphic. Printout. See CaseStudies.html#Segmented.
findpeaksb3.m is a variant of findpeaksb.m that fits each detected peak together with the previous and following peaks found by findpeaksG.m. It deals better with overlapping peaks than findpeaksb.m does, and it handles larger numbers of peaks better than findpeaksfit.m, but it fits only those peaks that are found by findpeaks. The syntax is function P=findpeaksb3(x,y, SlopeThreshold, AmpThreshold, smoothwidth, peakgroup, smoothtype, PeakShape, extra, NumTrials, AUTOZERO, ShowPlots). The first seven input arguments are exactly the same as for the findpeaksG.m function; if you have been using findpeaks or iPeak to find and measure peaks in your signals, you can use those same input argument values for findpeaksb3.m. The demonstration script DemoFindPeaksb3.m shows how findpeaksb3 works with multiple overlapping peaks.
findpeaksfit.m is essentially a
serial combination of findpeaksG.m
and peakfit.m. It uses
the number of peaks found by findpeaks and their peak positions
and widths as input for the peakfit.m function, which then fits
the entire signal with the specified peak model. This deals with
non-Gaussian and overlapped peaks better than findpeaks alone.
However, it fits only those peaks that are found by findpeaks.
The syntax is [P,
FitResults, LowestError, BestStart, xi, yi] = findpeaksfit(x,
y, SlopeThreshold, AmpThreshold, smoothwidth, peakgroup,
smoothtype, peakshape, extra, NumTrials, autozero,
fixedparameters, plots). The first
seven input arguments are exactly the same as for the findpeaksG.m
function; if you have been using findpeaks or iPeak to find and
measure peaks in your signals, you can use those same input
argument values for findpeaksfit.m. The remaining six input
arguments of findpeaksfit.m are for the peakfit function; if you
have been using peakfit.m or ipf.m
to fit peaks in your signals, you can use those same input
argument values for findpeaksfit.m. Type "help findpeaksfit" for
more information. See PeakFindingandMeasurement.htm#findpeaksfit.
Click
for
animated
example.
peakstats.m uses the same algorithm as
findpeaksG.m, 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 displaying the x,t data plot with numbered
peaks in figure window 1, the table of peak statistics in the
command window, and the histograms of the peak intervals,
heights, widths, and areas in figure window 2. Type "help
peakstats". See PeakFindingandMeasurement.htm.
Version 2, March 2016, adds median and mode.
tablestats.m (PS=tablestats(P,displayit)) is similar to peakstats.m except that it accepts as
input a peak table P such as generated by findpeaksG.m,
findvalleys.m, findpeaksL.m, findpeaksb.m, findpeaksplot.m,
findpeaksnr.m, findpeaksGSS.m, findpeaksLSS.m, or
findpeaksfit.m, any function that return a table of peaks with
at least 4 columns listing peak number, height, width, and area.
Computes the peak intervals (the x-axis interval between
adjacent detected peaks) and the maximum, minimum, average, and
percent standard deviation of each, and optionally displaying
the histograms of the peak intervals, heights, widths, and areas
in figure window 2. The optional last argument displayit = 1 if
the histograms are to be displayed, otherwise not.
findpeaksnr.m is a variant of findpeaksG.m that additionally computes the signal-to-noise ratio (SNR) of each peak and returns it in the 5^{th} 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 actual data and the least-squares fit over the top part of the peak). See PeakFindingandMeasurement.htm
findpeaksE.m is a variant of findpeaksG.m that additionally estimates the percent relative fitting error of each peak (assuming a Gaussian peak shape) and returns it in the 6^{th} column of the peak table.
findpeaksGSS.m and findpeaksLSS.m, for Gaussian and
Lorentzian peaks respectively, are variants of findpeaksG.m and
findpeaksL.m that additionally compute the 1% start and end
positions return them in the 6^{th} and 7^{th} columns of the peak
table. See PeakFindingandMeasurement.htm
findsquarepulse.m (syntax S=findsquarepulse(t,y,threshold) locates the rectangular pulses in the signal t,y that
exceed a y-value of "threshold" and determines their start time,
average height (relative to the adjacent baseline) and width. DemoFindsquare.m creates a test
signal and calls findsquarepulse.m to demonstrate.
findsteps.m P= findpulses(x, y, SlopeThreshold, AmpThreshold, SmoothWidth, peakgroup) locates positive transient steps in noisy x-y time series data, by computing the first derivative of y that exceed SlopeThreshold, computes the step height as the difference between the maximum and minimum y values over a number of data point equal to "Peakgroup". Returns list (P) with step number, x and y positions of the bottom and top of each step, and the step height of each step detected; "SlopeThreshold" and "AmpThreshold" control step sensitivity; higher values will neglect smaller features. Increasing "SmoothWidth" ignores small sharp false steps caused by random noise or by "glitches" in the data acquisition. See findsteps.png for a real example. findstepsplot.m plots the signal and numbers the peaks.
idpeaks, peak identification function. The syntax is [IdentifiedPeaks,AllPeaks] = idpeaks(DataMatrix, AmpT, SlopeT, sw, fw, maxerror, Positions, Names). Locates and identifies peaks in DataMatrix that match the position of peaks in the array "Positions" with matching names "Names". Type "help idpeaks" for more information. Download and extract idpeaks.zip for a working example, or see Example 8 on PeakFindingandMeasurement.htm
idpeaktable.m [IdentifiedPeaks]=idpeaktable(P, maxerror, Positions, Names). Compares the found peak positions in peak table "P" to a database of known peaks, in the form of an cell 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"), that peak is considered a match and its peak position, name, error, and amplitude are entered into the output cell array "IdentifiedPeaks". The peak table may be one returned by any of my peak finder or peak fitting functions, having one row for each peak and columns for peak number, position, and height as the first three columns.
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.
DemoFindPeak.m, a demonstration
script using the findpeaksG function on noisy synthetic data.
Numbers the peaks and prints out the peak parameters in the
command window. Requires that gaussian.m
and findpeaksG.m be present in the
path. See PeakFindingandMeasurement.htm#findpeaks.
TestFindpeaksG2d.m.
Demonstration script for findpeaks2d.m, which shows that this
function can locate peaks resulting in 'shoulders' that do not
produce a distinct maximum in the original signal. Detects peaks
in the negative of the smoothed second derivative of the signal
(shown as the dotted line in the figure). Requires gaussian.m,
findpeaksG.m, findpeaksG2d.m, fastsmooth.m, and peakfit.m in the
path. Graphic. Also uses the
TestFindpeaksG2d results as the
"start" value for iterative peak fitting using peakfit.m, which takes
longer to compute but gives more accurate results, especially
for width and area:
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 5^{th} column.
triangulationdemo.m is a demo
function (screen graphic)
that compares findpeaksG (which
determines peak parameters by curve-fitting a Gaussian to the
center of each peak) to findpeaksT,
which determines peak parameters by the triangle construction
method (drawing a triangle around each peak with sides that are
tangent to the sides of the peak). Performs the comparison with
4 different peak shapes: plain Gaussian, bifurcated Gaussian,
exponential modified Gaussian, and Breit-Wigner-Fano). In some
cases, the triangle construction method can be more accurate
that the Gaussian method if the peak shape is asymmetric.
findpeaksfitdemo.m, a
demonstration script of findpeaksfit automatically finding and
fitting the peaks in a set of 150 signals, each of which may
have 1 to 3 noisy Lorentzian peaks in variable locations.
Requires the findpeaksfit.m and lorentzian.m functions
installed. This script was used to generate the GIF animation findpeaksfit.gif.
FindpeaksComparison.m. Which
to use: findpeaksG, findpeaksb, findpeaksb3, or findpeaksfit?
This script compares all four functions applied to a
computer-generated signal with multiple peaks with variable types
and amounts of baseline and random noise. (Requires all these functions, plus modelpeaks.m,
findpeaksG, and findpeaksL.m, in the Matlab/Octave path. Type
"help FindpeaksComparison" for details). Results are displayed
graphically in figure windows 1, 2, and 3 and printed out
in a table of parameter accuracy and elapsed time for
each method. You may change
the lines in the script marked by <<< to modify the
number and character and amplitude of the signal peaks,
baseline, and noise. (Make the signal similar to yours to
discover which method works best for your type of signal). The
best method depends mainly on the shape and amplitude of the
baseline and on the extent of peak overlap.
iPeakEnsembleAverageDemo.m is a demonstration script for iPeak's ensemble average function. In this example, the signal contains a repeated pattern of two overlapping Gaussian peaks, 12 points apart, both of width 12, with a 2:1 height ratio. These patterns occur a random intervals throughout the recorded signal, and the random noise level is about 10% of the average peak height. Using iPeak's ensemble average function (Shift-E), the patterns can be averaged and the signal-to-noise ratio significantly improved.
ipeakdata.mat, data set for demonstrating idpeaks.m or the peak identification function of iPeak; includes a high-resolution atomic spectrum and a table of known atomic emission wavelenghs. See PeakFindingandMeasurement.htm#ipeak
Which to use: iPeak or Peakfit? Try these Matlab demo functions 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 sub-functions. Just place them in your path and type their name at the command prompt. You can download all these demos together in idemos.zip. They require no input or output arguments.
SpikeDemo1.m and SpikeDemo2.m are Matlab/Octave scripts that demonstrate how two measure spikes (very narrow peaks) in the presence of serious interfering signals. See CaseStudies.html#spikes.
PowerTransformTest.m is a simple script that demonstrates the power method of peak sharpening to aid in reducing in peak overlap. PowerMethodCalibrationCurve is a variant of PeakCalibrationCurve.m that evaluates the power method in the context of a flow injection or chromatography measurement. powertest2 is a self-contained function that demonstrates the power method for measuring the area of small shouldering peak (Graphic).
The script realtimepeak.m demonstrates simple real-time peak detection based on derivative zero-crossing, using mouse clicks to simulate data. Each time your mouse clicks form a peak (that is, go up and then down again), the program will register and label the peak on the graph (as illustrated on the right) and print out its x and y values. In this case, a peak is defined as any data point that has lower amplitude points adjacent to it on both sides, which is determined by the nested 'for' loops in lines 31-36. The more sophisticated script RealTimeSmoothedPeakDetectionGauss.m uses the technique described on PeakFindingandMeasurement.htm#findpeaks that locates the positive peaks in a noisy data set that rise above a set amplitude threshold, performs a least-squares curve-fit of a Gaussian function to the top part of the raw data peak, computes the position, height, and width (FWHM) of each peak from that least-squares fit and prints out each peak found in the command window. (Animated graphic).cls.m is a classical least-squares function that you can use to fit a computer-generated model, consisting of any number of peaks of known shape, width, and position, but of unknown height, to a noisy x,y signal. The syntax is heights= cls(x,y, NumPeaks, PeakShape, Positions, Widths, extra) where x and y are the vectors of measured signal (e.g. x might be wavelength and y might be the absorbance at each wavelength), 'NumPeaks' is the number of peaks, 'PeakShape' is the peak shape number (1=Gaussian, 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=BiGaussian, 15=BiLorentzian), 'Positions' is the vector of peak positions on the x axis (one entry per peak), 'Widths' is the vector of peak widths in x units (one entry per peak), and 'extra' is the additional shape parameter required by the exponentially broadened, Pearson, Gaussian/Lorentzian blend, BiGaussian and BiLorentzian shapes. cls.m returns a vector of measured peak heights for each peak. See clsdemo.m. (Note: this method is now included in the non-linear iterative peak fitter peakfit.m as peak shape 50. See the demonstration script peakfit9demo.m)
The cls2.m function is similar to cls.m, except that it also measures the baseline (assumed to be flat) and returns a vector containing the background B and measured peak heights H for each peak , e.g. [B H1 H2 H3...].
RegressionDemo.m,
script that demonstrates the classical least squares
procedure for a simulated absorption spectrum of a
5-component mixture at 100 wavelengths. Requires that gaussian.m be present in the path.
See CurveFittingB.html
clsdemo.m is a demonstration
script that creates a noisy signal, fits it using the
Classical Least Squares method with cls.m,
computes the accuracy of the measured heights, then
repeats the calculation using iterative least-squares
(using peakfit.m) for comparison. (This script requires cls.m, modelpeaks.m,
and peakfit.m in the
Matlab/Octave path).
CLSvsINLS.m is a script that compares
the classical least-squares (CLS) method with three
different variations of the iterative method (INLS)
method for measuring the peak heights of three Gaussian
peaks in a noisy test signal, demonstrating that the
fewer the number of unknown parameters, the faster and
more accurate is the peak height calculation.
Non-linear iterative curve fitting and peak fitting
gaussfit, function that performs least-squares fit of a single Gaussian function to an x,y data set, returning the height, position, and width of the best-fit Gaussian. Syntax is [Height,Position,Width] = gaussfit(x,y). The similar function lorentzfit.m performs the calculation for a Lorentzian peak shape. See CurveFitting.html#Transforming. The similar function plotgaussfit does the same thing as gaussfit.m but also plots the data and the fit. The data set can not contain any zero or negative values.
bootgaussfit is an expanded version of gaussfit that provides optional plotting and error estimation. The syntax is [Height, Position, Width, BootResults] = bootgaussfit(x, y, plots). If plots=1, plots the raw data as red dots and the best-fit Gaussian as a line. If the 4^{th} output argument (BootResults) is supplied, computes peak parameter error estimates by the bootstrap method.
fitshape2.m, syntax [Positions, Heights, Widths, FittingError] = fitshape2(x, y, start), is a simplified general-purpose Matlab/Octave function for fitting multiple overlapping model shapes to the data contained in the vector variables x and y. The model is linear combination of any number of basic peak shapes which are defined mathematically as a function of x, with two variables that the program will independently determine for each peak, positions and widths, in addition to the peak heights (i.e. the weights of the weighted sum). You must provide the first-guess starting vector 'start', in the form [position1 width1 position2 width2 ...etc], which specifies the first-guess position and width of each component (one pair of position and width for each peak in the model). The function returns the parameters of the best-fit model in the vectors Positions,Heights,Widths, and computes the percent error between the data and the model in FittingError. It also plots the data as dots and the fitted model as a line. The interesting thing about this function is that the only part that defines the shape of the model is the last line. In fitshape2.m, that line contains the expression for a Gaussian peak of unit height, but you could change that to any other expression or algorithm that computes g as a function of x with two unknown parameters 'pos' and 'wid' (position and width, respectively, for peak-type shapes, but they could represent anything for other function types, such as the exponential pulse, sigmoidal, etc.); everything else in the fitshape.m function can remain the same. This makes fitshape a good platform for experimenting with different mathematical expression as proposed models to fit data. There are also two other variations of this function for models with one iterated variable plus peak height (fitshape1.m) and three iterated variables plus peak height (fitshape3.m). Each has illustrative examples.
peakfit, a versatile
command-line function for multiple peak fitting by
iterative non-linear least-squares. A Matlab File Exchange "Pick
of
the
Week". The full syntax is [FitResults, GOF, baseline, coeff, BestStart,
xi, yi, BootResults] = peakfit(signal, center, window,
NumPeaks, peakshape, extra, NumTrials, start, AUTOZERO,
fixedwidth, plots, bipolar, minwidth). Type "help peakfit". See InteractivePeakFitter.htm#command.
Compared
to the fitshape.m function described previously, peakfit.m has a large
number of built-in peak shapes selected by number, it does not require
(although it can be given) the first-guess position and width of
each component, and it has features for background correction
and other useful features to improve the quality and estimate
the reliability of fits. Test the installation on your computer by
running the autotestpeakfit.m script,
which runs through the whole gauntlet of fitting tests
without pause, printing out what it's doing and the
results, checking to see if the fitting error is greater
than expected and printing out a WARNING if it is. See the Version
history, a brief description of the new features of each
version of peakfit.m from 3.7 to the present.
testnumpeaks(x,y,peakshape,extra,NumTrials,MaxPeaks).
Simple
test
to estimate the number of model peaks required to fit an x,y
data set. Fits data x,y, with shape "peakshape", with optional
extra shape factor "extra", with NumTrials repeat per fit, up to
a maximum of "MaxPeaks" model peaks, displays each fit and
graphs fitting error vs number of model peaks. If two or more
numbers of peaks give about the same error, its best to take the
smaller number.
SmoothVsFit.m is a demonstration script that compares iterative least-square fitting to two simpler methods
for the measurement of the peak height of a single peak of
uncertain width and position and with a very poor
signal-to-noise ratio of 1. The accuracy and precision of the
methods are compared. SmoothVsFitArea.m
does the same thing for the measurement of peak area. See CurveFittingC.html#Fitting_peaks.
ipf.m, is
an interactive multiple peak fitter (m-file link: ipf.m). It uses iterative nonlinear
least-squares to fit any number of overlapping peaks of
the same or different peak shapes to x-y data sets. Demoipf.m is a demonstration
script for ipf.m, with a built-in
simulated signal generator. The true values of the
simulated peak positions, heights, and widths are
displayed in the Matlab command window, for comparison
to the FitResults obtained by peak fitting. Animated step-by-step
instructions . You can also download a ZIP file containing ipf.m plus some
examples and demos. Click for animated
example. Version
13.2, June 2019, added "Working..." display on graph when
fitting is in progress. Previous version 13.1:
December, 2018, has modified the action of the "d" key to save
the model to disc as SavedModel.mat.
SmallPeak.m
is a demonstration of several curve-fitting techniques
applied to the challenging problem of measuring the
height of a small peak that is closely overlapped with
and completely obscured by a much larger peak. It
compares iterative fits by unconstrained, equal-width, and fixed-position
models (using peakfit.m) to a
classical least squares fit in which only the peak heights are unknown (using cls.m). Spread out the four figure
windows so you can observe the dramatic difference in
stability of the different methods. A final table of
relative percent peak height errors shows that the more
the constraints, the better the results (but only if the constraints are
justified). See CaseStudies.html#BuriedPeak.
BlackbodyDataFit.m,
a script that
demonstrates iterative
least-squares
fitting
of the blackbody
equation to a measured spectrum of an incandescent body, for the
purpose of estimating its color temperature. See CurveFittingC.html#Matlab
Demofitgauss.m a script that demonstrates iterative fitting a single Gaussian function to a set of data, using the
fminsearch function. Requires that gaussian.m
and fmsearch.m (in the "Optim 1.2.1" package) be
installed. Demofitgaussb.m
and fitgauss2b.m illustrate
a modification of this technique to accommodate shifting baseline (Demofitlorentzianb.m
and fitlorentzianb.m for
Lorentzian peaks). This modification is now incorporated
to peakfit.m (version 4.2 and later), ipf.m (version 9.7
amd later), findpeaksb.m (version 3 and later), and
findpeaksfit, (version 3 and later). See CurveFittingC.html#Matlab
Demofitgauss2.m a script that demonstrates iterative fitting of two overlapping Gaussian functions to a set of data, using the fminsearch function. Requires that gaussian.m and fmsearch.m (in
the "Optim 1.2.1" package) be installed. Demofitgauss2b.m is the baseline-corrected extension.
See CurveFittingC.html#Matlab
VoigtFixedAlpha.m and VoigtVariableAlpha.m
demonstrate two different ways to fit peaks with variable shapes, such as
the Voigt profile, Pearson, Gauss-Lorentz blend, and the
bifurcated and exponentially-broadened shapes, which are defined
not only by a peak position, height, and width, but also by an
additional "shape" parameter that fine-tunes the shape of the
peak. If that parameter is equal for all peaks in a group, it can be passed as an
additional input argument to the shape function, as shown in VoigtFixedAlpha.m
for the Voigt shape (which is calculated as a convolution of Gaussian
and Lorentzian components with different widths). If
the shape parameter is allowed to be different
for each peak in the group and is to
be determined by iteration (just as is position and width), then
the routine must be modified to accommodate three, rather than two, iterated variables,
as shown in VoigtVariableAlpha.m.
Although the fitting error is lower with variable alphas, the execution time is longer and
the alphas values so determined are
not very stable, with respect to
noise in the data and the starting guess values, especially for
multiple peaks. The script
VoigtShapeFittingDemonstration.m
uses the general-purpose
Matlab/Octave function peakfit.m
to fit both fixed and variable shape Voigt peaks and to
calculate the Gaussian width component, Lorentzian width
component, alpha, and the peak width. It computes the
theoretical Voigt profile and adds random noise for realism.
(Requires voigt.m, halfwidth.m,
and peakfit.m in the path). See CurveFittingC.html#Fitting_peaks
Demofitmultiple.m. Demonstrates
an iterative fit to sets of computer-generated noisy peaks of
different types, knowing only the shape types and variable shape
parameters of each peak. Iterated parameters are shape, height,
position, and width of all peaks. Requires the fitmultiple.m and peakfunction.m functions. View screen shot. See CurveFittingC.html#Fitting_peaks
BootstrapIterativeFit.m, a function that demonstrates bootstrap
estimation of the variability of an iterative
least-squares fit to a single noisy Gaussian peak. Form
is: BootstrapIterativeFit(TrueHeight,
TruePosition,
TrueWidth,
NumPoints, Noise, NumTrials). See
CurveFitting.html#Bootstrap
BootstrapIterativeFit2.m, a function that demonstrates bootstrap
estimation of the variability of an iterative
least-squares fit to two noisy Gaussian peaks. Form is: BootstrapIterativeFit2(TrueHeight1,
TruePosition1,
TrueWidth1,
TrueHeight2, TruePosition2, TrueWidth2, NumPoints, Noise,
NumTrials). See CurveFitting.html#Bootstrap
DemoPeakfitBootstrap.m.
Self-contained demonstration function for peakfit.m, with
built-in signal generator. Demonstrates bootstrap error
estimation. See CurveFitting.html#Bootstrap
DemoPeakfit.m, Demonstration script (for peakfit.m) that generates
an overlapping peak signal, adds noise, fits it with peakfit.m,
then computes the accuracy and precision of peak parameter
measurements. Requires that peakfit.m be present in the path. See InteractivePeakFitter.htm#peakfitdemos
peakfit9demo.
Demonstrates multilinear regression (shape 50) available in peakfit.m version 9
(Requires modelpeaks.m and peakfit.m in the Matlab path).
Creates a noisy model signal of three peaks of known shapes,
positions, and widths, but unknown heights. Compares multilinear
regression in Figure 1 with unconstrained iterative non-linear
least squares in Figure 2. For shape 50, the 10th input argument
fixedparameters must be a matrix listing the peak shape
(column 1), position (column 2), and width (column 3) of each
peak, one row per peak. peakfit9demoL
is similar but uses Lorentzian peaks (specified in the fixedparameters matrix and in
the PeakShape vector).
DemoPeakFitTime.m is a simple script that demonstrates how to use peakfit.m to apply multiple curve fits to a signal that is changing with time. The signal contains two noisy Gaussian peaks in which the peak position of the second peak increases with time and the other parameters remain constant (except for the noise). The script creates a set of 100 noisy signals (on line 5) containing two Gaussian peaks where the position of the second peak changes with time (from x=6 to 8) and the first peak remains the same. Then it fits a 2-Gaussian model to each of those signals (on line 8), displays the signals and the fits graphically with time as a kind of animation (click to play animation), then plots the measured peak position of the two peaks vs time on line 12.
isignal, can be used as a command-line function in Octave, but its interactive features currently work only in Matlab. The syntax is isignal(DataMatrix, xcenter, xrange, SmoothMode, SmoothWidth, ends, DerivativeMode, Sharpen, Sharp1, Sharp2, SlewRate, MedianWidth). See iSignal.
testpeakfit.m, a test script that demonstrates 36 different examples on InteractivePeakFitter.htm#Examples. Use for testing that peakfit and related functions are present in the path. Updated for peakfit 6. autotestpeakfit.m does the same same without pausing between functions and waiting for a keypress.
Multiple peak fits with different
profiles. ShapeTestS.m and ShapeTestA.m tests the data in its
input arguments x,y, assumed to be a single isolated peak,
fits it with different candidate model peak
shapes using peakfit.m, plots each fit
in a separate figure window, and prints out a table of fitting
errors in the command window. ShapeTestS.m
tries seven different candidate symmetrical model peaks, and ShapeTestA.m tries six different
candidate asymmetrical model peaks. The one with the lowest
fitting error (and R2 closest to 1.000) is presumably the best
candidate. Try the examples in their help
files. But beware: if there is too much
noise in your data, the results can be misleading. For
example, a multiple Gaussians model is likely to fit best
because it has more degrees of freedom and can "fit the
noise", even if the actual peak shape is
something other than a Gaussian. (The function peakfit.m has
many more built-in shapes to choose from, but still it is a finite list and there is
always the possibility that the actual underlying peak shape
is not available in the software you are using or that it is
simply not describable by a single mathematical function).
WidthTest.m is a script that
demonstrates that constraining some of the peak parameters of
a fitting model to fixed values, if those values are
accurately known, improves that accuracy of measurement of the
other parameters, even though it increases the fitting error.
Requires installation of the GL.m and peakfit.m functions (version 7.6 or
later) in the Matlab/Octave path.
The script NumPeaksDemo.m
demonstrates one way to attempt to estimate the minimum number
of model peaks needed to fit a set of data, plotting the
fitting error vs the number of model peaks, and looking for
the point at which the fitting error reaches a minimum. This
script creates a noisy computer-generated signal containing a
user-selected 3, 4, 5 or 6 underlying Lorentzian peaks and
uses peakfit.m to fit the data to a series of models
containing 1 to 10 model peaks. The correct number of
underlying peaks is either the fit with the lowest fitting
error, or, if two or more fits have about the same fitting
error, the fit with the least number of peaks, as in this example, which actually has
4 underlying peaks. If the data are very noisy, however, the
determination becomes unreliable. (To make this demo closer to
your type of data, you could change Lorentzian to Gaussian or
any other model shape, or change the peak width, number of
data points, or the noise level). This script requires that peakfit.m and the appropriate shape functions
(gaussian.m, lorentzian.m, etc) be present in the path. The
function testnumpeaks.m does
this for your own x,y data.
Peakfit Time Tests. These are a series of
scripts that demonstrate how the execution time of the
peakfit.m function varies with
the peak shape (PeakfitTimeTest2.m
and PeakfitTimeTest2a.m, with number
of peaks in the model (PeakfitTimeTest.m), and with
the number of data points in the fitted region (PeakfitTimeTest3.m). This issue
is discussed on InteractivePeakFitter.htm#time.
TwoPeaks.m is a simple 8-line
script that compares findpeaksG.m and peakfit.m with a
signal consisting to two noisy peaks. findpeaksG.m and
peakfit.m must be in the Matlab/Octave path.
peakfitVSfindpeaks.m
performs a direct comparison of the accuracy of findpeaksG vs
peakfit. This script generates four very noisy peaks of
different heights and widths, then applies findpeaksG.m and
peakfit.m to measure the peaks and compares the results. The
peaks detected by findpeaks are labeled "Peak 1", "Peak 2",
etc. If you run this script several times, you'll find that
both methods work well most of the time, with peakfit giving
smaller errors in most cases, but occasionally findpeaks will
miss the first (lowest) peak and rarely it will detect an
extra peak that is not there if the signal is very noisy.
CaseStudyC.m is a
self-contained Matlab/Octave demo function that demonstrates
the application of several techniques described on this site
to the quantitative measurement of a peak buried in an
unstable background, a situation that can occur in the
quantitative analysis applications of various forms of
spectroscopy and remote sensing. See Case Studies C.
GaussVsExpGauss.m Comparison of
alternative models for the unconstrained exponentially
broadened Gaussians, shapes 31 and 39. Shape 31 (expgaussian.m) creates the
shape by performing a Fourier convolution of a
specified Gaussian by an exponential decay of
specified time constant, whereas shape 39 (expgaussian2.m) uses a
mathematical expression for the final shape so
produced. Both result in the same shape
but are parameterized differently. Shape 31
reports the peak height and position as that of the
original Gaussian before broadening, whereas shape 39
reports the peak height of the broadened result. Shape
31 reports the width as the FWHM of the original
Gaussian and shape 39 reports the standard deviation
(sigma) of that Gaussian. Shape 31 reports the
exponential factor an the number of data points
and shape 39 reports the reciprocal of time constant
in time units. See
Figure windows 2
and 3. You must
have peakfit.m (version 8.4)
gaussian.m, expgaussian.m, expgaussian2.m, findpeaksG.m, and halfwidth.m in the
Matlab/Octave path. DemoExpgaussian.m
is a script that gives another, more detailed,
exploration of the effect of exponential broadening on
a Gaussian peak (requires gaussian.m,
expgaussian.m, halfwidth.m, val2ind.m, and peakfit.m to be in the
Matlab/Octave path).
AsymmetricalOverlappingPeaks.m is a multi-step script that demonstrates the use of a combination of first-derivative symmetrization before curve fitting to analyze a complex mystery peak.
Keystroke-operated interactive signal processing tools (for Matlab only)
The three interactive functions described above, iPeak, iSignal, and ipf, all have several keystroke commands in common: all share the same set of pan and zoom keys (cursor arrow keys, < and > keys, [ and ] keys, etc) to adjust the portion of the signal that is displayed in the upper panel. All use the K key to display the list of keystroke commands. All use the T key to cycle through the four baseline connection modes. All use the Shift-Ctrl-S, Shift-Ctrl-F, and Shift-Ctrl-P keys to transfer the current signal between iSignal, ipf, and iPeak, respectively. To make it easier to transfer settings from one of these functions to other related functions, all use the W key to print out the syntax of other related functions, with the pan and zoom settings and other numerical input arguments specified, ready for you to Copy, Paste and edit into your own scripts or back into the command window. For example, you can convert a curve fitting operation performed in ipf.m into the command-line peakfit.m function; or you can convert a peak finding operation performed in ipeak.m into a command-line findpeaksG.m or findpeaksb.m functions.
Hyperlinear Quantitative Absorption Spectrophotometry
tfit.m, a self-contained
command-line Matlab/Octave function that demonstrates a computational method for quantitative
analysis by multiwavelength absorption spectroscopy which uses
convolution and iterative curve fitting to correct for
spectroscopic non-linearity. The syntax is tfit(TrueAbsorbance). TFitStats.m is a script that
demonstrates the reproducibility of the method. TFitCalCurve.m compares the
calibration curves for single-wavelength, simple regression,
weighted regression, and TFit methods. TFit3.m
is a demo function for a mixture of 3 absorbing components;
the syntax is TFit3(TrueAbsorbanceVector), e.g. TFit3([3 .2 5]). Download all these as
a ZIP file. Click
for
animated
example.
TFitDemo.m is a keypress-operated
interactive explorer for
the Tfit method (for Matlab only), applied to the
measurement of a single component with a Lorentzian
(or Gaussian) absorption peak, with controls that
allow you to adjust the true absorbance (Peak A),
spectral width of the absorption peak (AbsWidth),
spectral width of the instrument function (InstWidth),
stray light, and the noise level (Noise) continuously
while observing the effects graphically and
numerically. See TFit.html. Click
for
animated
example.
MAT files (for Matlab and Octave) and Text files (.txt)
DataMatrix2
is a computer-generated test signal consisting of 16
symmetrical Gaussian peaks with random white noise added. Can
be used to test the peakfit.m function. See CurveFittingC.html#Exponential_broadening
DataMatrix3 is a
computer-generated test signal consisting of 16 Gaussian peaks
with random white noise that have been exponentially broadened
with a time constant of 33 x-axis units. See CurveFittingC.html#Exponential_broadening
udx.txt: a text file containing the 2 x
1091 matrix that consists of two Gaussian peaks with different
sampling intervals. It's used as an example in Smoothing and in Curve Fitting.
TimeTrial.txt
is a text file comparing the speed of several different signal
processing tasks running on 5 different systems: (1) Windows
10, 64 bit, 3.6 GHz core i7, with 16 GBytes RAM, using Matlab
9.9 (R2020b) Update 3 academic, (2) Matlab 2017b Home, and (3)
Matlab Online R2018b in Chrome, both running on older desktop
PCs, (4) Matlab Mobile on an iPad, and (5) Octave 3.6.4. The
Matlab/Octave code that generated this is TimeTrial.m, which runs all the tasks
one after the other and prints out the elapsed times for your
particular machine, in addition to the times previously
recorded for each task on each of the five software systems.
Generally speaking, the more recent version of Matlab or of
Octave are faster than their older versions. Matlab versions
are almost always faster than Octave, which is free.
Readability.txt. English
language readability analysis of IntroToSignalProcessing.pdf
performed by http://www.online-utility.org/english/readability_test_and_improve.jsp
Note 1: If you see a yellow bar at the top of the spreadsheet window, click the "Enable Editing" button.
Note 2: If your browser changes the file extension of these spreadsheets to .zip when they are downloaded, rename the files to their original file extensions (.ods, .xls, or .xlsx) before running them.
Note: These spreadsheets have no protected cells, so there is nothing stopping you from changing the formulas accidentally. If you mess up, just use the Undo function (Ctrl-Z) or download another copy. But this also means you can modify any aspect of these spreadsheets for your own purposes, which you are invited to do.
Random numbers and noise. The spreadsheets RandomNumbers.xls (for Excel)
and RandomNumbers.ods (for
OpenOffice) demonstrate how to create a column of
normally-distributed random numbers (like white noise) in a
spreadsheet that has only a uniformly-distributed random
number function. Also shows how to compute the interquartile
range and the peak-to-peak value and how they compare to the
standard deviation. See SignalsAndNoise.html.
The same technique is used in the spreadsheet SimulatedSignal6Gaussian.xlsx,
which computes and plots a simulated signal consisting of up
to 6 overlapping Gaussian bands plus random white noise.
Smoothing. The spreadsheets smoothing.ods (for Open office
Calc) and smoothing.xls (for
Microsoft Excel) demonstrate a 7-point rectangular (sliding
average) in column C and a 7-point
triangular smooth in column E, applied to the data in column
A. You can type in (or Copy and Paste) any data you like into
column A. You can extend the
spreadsheet to longer columns of data by dragging the last row
of columns A, C, and E down as needed. You
can change the smooth width by changing the equations in
columns C or E. The spreadsheet MultipleSmoothing.xls for
Excel or Calc demonstrates a more flexible method that allows
you to define various types of smooths by typing a few integer
numbers. The spreadsheets UnitGainSmooths.xls
and UnitGainSmooths.ods
contain a collection of unit-gain convolution coefficients for
rectangular, triangular, and Gaussian smooths of width 3 to 29
in both vertical (column) and horizontal (row) format. You can
Copy and Paste these into your own spreadsheets. Convolution.txt lists some simple
whole-number coefficient sets for performing single and
multi-pass smoothing.
VariableSmooth.xlsx demonstrates an even more powerful
and flexible technique, especially for very large and
variable smooth widths, that uses the spreadsheet AVERAGE
and INDIRECT functions. It allows you to change the smooth
width simply by changing the value of a single cell. See
Smoothing.html for details.
SegmentedSmoothTemplate.xlsx
is a segmented multiple-width data smoothing
spreadsheet template, which can apply individually specified different smooth
widths to different regions of the signal, especially useful
if the widths of the peaks or the noise level varies
substantially across the signal. In this version there are
20 segments. SegmentedSmoothExample.xlsx
is an example with data (graphic).
A related sheet GradientSmoothTemplate.xlsx
(graphic) performs a
linearly increasing (or decreasing) smooth width across the
entire signal, given only the start and end values,
automatically generating as many segments are are necessary.
Differentiation. DerivativeSmoothingOO.ods
(for OpenOffice Calc) and DerivativeSmoothing.xls
(for Excel) demonstrate the application of differentiation for
measuring the amplitude of a peak that is buried in a broad
curved background. Differentiation and smoothing are both
performed together. Higher order derivatives are computed by
taking the derivatives of previously computed derivatives. DerivativeSmoothingWithNoise.xlsx
is a related spreadsheet that demonstrates the
dramatic effect of smoothing on the signal-to-noise
ratio of derivatives on a noisy signal. It uses the
same signal as DerivativeSmoothing.xls,
but adds simulated white noise to the Y data. You can
control the amount of added noise. SecondDerivativeXY2.xlsx,
demonstrates locating and measuring changes in the second
derivative (a measure of curvature or acceleration) of a
time-changing signal, showing the apparent increase in noise
caused by differentiation and the extent to which the noise
can be reduced by smoothing (in this case by two passes of a
5-point triangular smooth). The smoothed second derivative
shows a large peak at the point where the acceleration changes
and plateaus on either side showing the magnitude of the
acceleration before and after the change (2 and 4,
respectively). Convolution.txt
lists simple whole-number coefficient sets for performing
differentiation and smoothing. CombinedDerivativesAndSmooths.txt
lists the sets of unit-gain coefficients that perform 1^{st} through 4^{th} derivatives with various degrees
of smoothing. See Differentiation.html
Convolution. Spreadsheets can be used to perform "shift-and-multiply" convolution for small data sets (for example, MultipleConvolution.xls or MultipleConvolution.xlsx for Excel and MultipleConvolutionOO.ods for Calc), which is essentially the same technique as the above spreadsheets for smoothing and differentiation. Use this spreadsheet to investigate convolution, smoothing, differentiation, and the effect of those operations on noise and signal-to-noise ratio. (For larger data sets the performance is slower that Fourier convolution, which is much easier done in Matlab or Octave than in spreadsheets). Convolution.txt lists simple whole-number coefficient sets for performing differentiation and smoothing.
Peak Area Measurement. EffectOfDx.xlsx demonstrates that the simple equation sum(y)*dx accurately measures the peak area of an isolated Gaussian peak if there are at least 4 or 5 points visibly above the baseline. EffectOfNoiseAndBaseline.xlsx demonstrates the effect of random noise and non-zero baseline, showing that the area is more sensitive to non-zero baseline that the same amount of random noise. PeakSharpeningAreaMeasurementDemo.xlsm (screen image) demonstrates the effect of derivative peak sharpening on perpendicular drop area measurements of two overlapping Gaussian peaks. Sharpening the peaks reduces the degree of overlap and can greatly reduce the peak area measurement error errors made by the perpendicular drop method.
Curve Fitting. LeastSquares.xls and LeastSquares.odt perform
polynomial least-squares fits to a straight-line model and QuadraticLeastSquares.xls
and QuadraticLeastSquares.ods
does the same for a quadratic (parabolic model. There are
specific versions of these spreadsheets that also calculate
the concentrations of the unknowns (download complete set as
CalibrationSpreadsheets.zip).
Multi-component spectroscopy. RegressionTemplate.xls and RegressionTemplate.ods (graphic
with example data) perform multicomponent analysis using
the matrix method for a fixed
5-component, 100 wavelength data set. RegressionTemplate2.xls
uses a more advanced spreadsheet
technique that allows the template to automatically adjust
to different numbers of components and wavelengths. Two examples
show the same template with data entered for a mixture
of 5 components measured at 100 wavelengths (RegressionTemplate2Example.xls)
and for 2 components at 59 wavelengths (RegressionTemplate3Example.xls).
Peak fitting. A set of spreadsheets using the Solver function to perform iterative non-linear peak fitting for multiple overlapping peak models is described here. There are versions for Gaussian and for Lorentzian peak shapes, with and without baseline, for 2-6 peak models and 100 wavelengths (with instructions for modification). All of these have file names beginning with "CurveFitter...".
Peak detection and measurement. The spreadsheet PeakDetection.xls and PeakDetection.xlsx implement the simple derivative zero-crossing peak detection method. The input x,y data are contained in Sheet1, column A and B, rows 9 to 1200. (You can paste your own data there). The amplitude threshold and slope threshold are set in cells B4 and E4, respectively. Smoothing and differentiation are performed by the convolution technique used by the spreadsheet DerivativeSmoothing.xls described previously. The Smooth Width and the Fit Width are both controlled by the number of non-zero convolution coefficients in row 6, 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. The peak index numbers, X-axis positions, and peak heights are listed in columns AC to AF. Peak heights are computed two ways: "Height" is based on slightly smoothed Y values (more accurate if the data are noisy) and "Max" is the highest individual Y value near the peak (more accurate if the data are smooth or if the peaks are very narrow). See PeakDetectionExample.xlsx/.xls) for an example with data already pasted in. PeakDetectionDemo2.xls/xlsx is a demonstration with a user-controlled computer-generated series of peaks. PeakDetectionSineExample.xls is a demo that generates a sinusoidal signal with an adjustable number of peaks.
An extension
of that method is made in PeakDetectionAndMeasurement.xlsx
(screen image),
which
makes the assumption that the peaks are Gaussian and
measures their height, position, and width more
precisely using a least-squares technique, just like "findpeaksG.m". 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. See PeakFindingandMeasurement.htm#Spreadsheet.
If the peaks in the data are too much overlapped, they may not
make sufficiently distinct maxima to be detected reliably. If
the noise level is low, the peaks can be artificially sharped
by the derivative
sharpening technique described previously. This is
implemented by PeakDetectionAndMeasurementPS.xlsx
and its demo version PeakDetectionAndMeasurementDemoPS.xlsx.
Spreadsheets for the TFit Method: Hyperlinear Quantitative Absorption Spectrophotometry. TransmissionFittingTemplate.xls (screen image) is an empty template for a single isolated peak; TransmissionFittingTemplateExample.xls (screen image) is the same template with example data entered. TransmissionFittingDemoGaussian.xls (screen image) is a demonstration with a simulated Gaussian absorption peak with variable peak position, width, and height, plus added stray light, photon noise, and detector noise, as viewed by a spectrometer with a triangular slit function. You can vary all the parameters and compare the best-fit absorbance to the true peak height and to the conventional log(1/T) absorbance.
TransmissionFittingCalibrationCurve.xls
(screen
image) includes an Excel macro
(see CaseStudies.html#Using_Macros)
that automatically constructs a calibration curve comparing
the TFit and conventional log(1/T) methods, for a series of 9
standard concentrations that you can specify. To create a
calibration curve, enter the standard concentrations in AF10 -
AF18 (or just use the ones already there, which cover a
10,000-fold concentration range), enable macros, then press Ctrl-f (or click Developer tab, click Macros, select macro2, and click Run). This macro
constructs and plots the calibration curve for both the TFit
(blue dots) and conventional (red dots) methods and computes
the R2 value for the TFit calibration curve, in the upper
right corner of graph. (Note: you can also use this
spreadsheet to compare the precision and reproducibility of
the two methods by entering the same concentration 9 times
in AF10 - AF18. The result should be a straight flat line with
zero slope).
Advanced spreadsheet techniques: "SpecialFunctions.xlsx" (Graphic) demonstrates the applications of the MATCH, INDIRECT, COUNT, IF, and AND functions when dealing with data arrays of variable size. "IndirectLINEST.xls" (Graphic link) demonstrates the particular benefit of using the INDIRECT function in array functions such as INV and LINEST.
Copyright (c) 2021, 2024 Thomas C. O'Haver
Permission is hereby granted, free of charge, to any person
obtaining a copy of this software and associated documentation
files (the "Software"), to deal in the Software without
restriction, including without limitation the rights to use,
copy, modify, merge, publish, distribute, sublicense, and/or
sell copies of the Software, and to permit persons to whom the
Software is furnished to do so, subject to the following
conditions:
The above copyright notice and this permission notice shall be
included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
OTHER DEALINGS IN THE SOFTWARE.
This page is part of "A Pragmatic Introduction to Signal Processing", a retirement
project and international
community service, 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.
Updated May, 2023