refinepeaks
Syntax
Description
[
refines peaks based on a signal with amplitudes yRPeaks
,xRPeaks
] = refinepeaks(y
,xPeaksIdx
)y
, and indices of the
initial peak estimates xPeaksIdx
, returning the refined peak values
yRPeaks
and location estimates xRPeaks
. This
syntax uses quadratic least squares processing to refine the peaks based on
xRPeaks
and two surrounding points for every selected peak.
[___] = refinepeaks(___,
specifies additional options, such as the peak refinement method and method specifications,
using name-value arguments. Use this syntax with any of the input or output arguments in
preceding syntaxes.Name=Value
)
refinepeaks(___)
with no output arguments draws a
figure that compares the refined peak amplitudes and locations to the initial peak estimates
and their surrounding samples. The figure also plots the fitting curves used to refine the
peaks. The refinepeaks
function plots the three-point sets with
unfilled circles and plots the refined peaks with filled circles.
Examples
Refine Known Peak from Signal
Refine a known peak from a signal specified as a vector or as a MATLAB® timetable.
Create a five-element vector. Set the fourth sample so that it is a local maximum. Plot the function and observe that the fourth sample is a signal peak with an amplitude of 102.4.
Intensity = [100;98.7;97.2;102.4;99.1]; plot(Intensity)
Refine the known peak and compute its refined value and location. The refined peak is located slightly after the original peak on the fourth sample of the signal.
[YRpeak,XRpeak] = refinepeaks(Intensity,4)
YRpeak = 102.4531
XRpeak = 4.1118
Create a timetable from a signal, using a sample rate of 1 Hz.
TT = timetable(Intensity,SampleRate=1)
TT=5×1 timetable
Time Intensity
_____ _________
0 sec 100
1 sec 98.7
2 sec 97.2
3 sec 102.4
4 sec 99.1
Refine the peak and plot the result. Observe that the refined peak occurs 0.11 seconds after the initial estimate, and the refined peak value of 102.45 shows as a filled circle. The original peak is at 0 seconds for reference, while the surrounding points are at -1 second and 1 second from the reference.
refinepeaks(TT,4)
Quadratic Least Squares Peak Refinement
Generate a Rayleigh distribution $\mathit{P}=\mathit{f}\left(\mathit{x}|\mathit{b}\right)=\frac{\mathit{x}}{{\mathit{b}}^{2}}{\mathit{e}}^{-\frac{{\mathit{x}}^{2}}{2{\mathit{b}}^{2}}}$ with a scale parameter b
set to $\pi $.
P = @(x,b) (x/b^2).*exp(-0.5*(x/b).^2); x = 0:2.5:20; b = pi; p = P(x,b); plot(x,p)
Calculate the theoretical and estimated peak amplitude and location. The peak is located about half of a sample from the expected peak scaled location, $\pi $.
yPeakTheo = P(b,b); xPeakTheo = b; yPeakEst = max(p); xPeakIdx = find(p==max(p)); xPeakEst = x(xPeakIdx); disp(table(["Estimated";"Theoretical"], ... [yPeakEst;yPeakTheo],[xPeakEst;xPeakTheo], ... VariableNames=["Peak" "Amplitude" "Location"]))
Peak Amplitude Location _____________ _________ ________ "Estimated" 0.18456 2.5 "Theoretical" 0.19306 3.1416
Refine the peak estimation using the quadratic least squares (QLS) method with default specifications.
[yRPeakDef,xRPeakDef] = refinepeaks(p,xPeakIdx,x);
Refine the peak estimation using the QLS method by specifying an exponential kernel with an exponent value of three.
[yPeakRef,xPeakRef] = refinepeaks(p,xPeakIdx,x, ... Method="QLS",QLSKernel="exp",KernelPower=3);
Compare the peak refinement results. The refined peak amplitude and location with both QLS-method variants yield a better approximation of the theoretical values than do the initial estimates. The exponential kernel yields the best approximation of the peak theoretical values for the specified Rayleigh distribution.
disp(table(["Initial Estimation";"Refined (QLS default)"; ... "Refined (QLS customized)";"Theoretical"], ... [yPeakEst;yRPeakDef;yPeakRef;yPeakTheo], ... [xPeakEst;xRPeakDef;xPeakRef;xPeakTheo], ... VariableNames=["Peak" "Amplitude" "Location"]))
Peak Amplitude Location __________________________ _________ ________ "Initial Estimation" 0.18456 2.5 "Refined (QLS default)" 0.19581 3.2884 "Refined (QLS customized)" 0.19315 3.2398 "Theoretical" 0.19306 3.1416
Nonlinear Least Squares Peak Refinement
Obtain a refined peak location estimate for the main two peaks in a signal using the nonlinear least squares method with a sinc function kernel.
Generate Signal
Radar pulse compression of a linear FM waveform produces a sinc-shaped spectrum, where the peaks frequency locations are proportional to the distance between the radar and the detected object. You can first estimate the peak locations and amplitudes with findpeaks
and then enhance your estimates with refinepeaks
. This example recreates a noiseless pulse-compression signal, finds and refines the signal peak amplitudes and locations.
Generate a signal composed of two sinc-shaped waveforms with peaks of 1 and 1.5 at 4.76 kHz and 35.8 kHz, respectively. The frequency spacing is 2.5 Hz.
aTarget = [1 1.5]; fTarget = 1e3*[4.76 35.8]; freqkHzFull = (0:0.0025:50)'; waveFull = abs(sinc([1 0.5].*(freqkHzFull-fTarget/1e3)))*aTarget';
Downsample the signal by a factor of 200 so the frequency spacing between samples is 0.5 kHz. This example refines the amplitude and location estimates of the downsampled signal peaks and compares the improved estimates to the values in the original signals.
freq = downsample(freqkHzFull,200); wave = downsample(waveFull,200); plot(freqkHzFull,waveFull,freq,wave,"*") legend(["Full signal" "Selected samples"],Location="northwest") xlabel("Frequency (kHz)") ylabel("Magnitude")
Refine Peaks Using Nonlinear Least Squares
Use findpeaks
to make initial estimates of the amplitudes, locations, and half-height widths of the two highest peaks of the signal.
[PV,PL,PW] = findpeaks(wave,NPeaks=2, ... SortStr="descend",WidthReference="halfheight");
Use refinepeaks
to enhance the peak estimation using the nonlinear least squares (NLS) method. Specify the frequency points of the signal and the peak widths. The peak values are significantly closer to the expected values of 1.5 and 1 while the frequency locations approximate well to 35.8 kHz and 4.76 kHz, respectively.
LW = max(PW,2);
[Ypk,Xpk] = refinepeaks(wave,PL,freq,Method="NLS",LobeWidth=LW)
Ypk = 2×1
1.5063
1.0163
Xpk = 2×1
35.8001
4.7628
Plot the amplitudes of refined peaks on the y-axis and the updated peak locations compared with the initial peak estimates on the x-axis. The two initially estimated peaks and their two surrounding samples are separated 0.5 kHz between each other. The refined peaks, shown with filled circles, indicate a repositioning of the actual peak locations compared with the initially estimated peak locations, while approximating the amplitudes to the expected values.
refinepeaks(wave,PL,freq,Method="NLS",LobeWidth=LW) yline(aTarget) % Theoretical peak amplitudes errorBounds = aTarget.*(1+0.03*[-1;1]); yline(errorBounds(:),":") % ±3% error bounds legend("Peak "+[1 2])
Find and Refine Peaks from Three-Dimensional Data
Identify and refine peaks from hourly temperature data across rows, columns, and pages.
Plot Data
Load a MAT file containing a set of temperature readings in degrees Celsius taken every hour at Logan Airport in Boston for the entire month of January, 2011. Reshape the data from the first 28 days into a three-dimensional array with 24 rows (hours of the day), 7 columns (week days), and 4 pages (week numbers). Plot the temperatures across weeks.
load bostemp tempsJan = (tempC(1:24*28))'; timesJan = 1 + (0:numel(tempsJan)-1)/24; % Reshape matrix to 3D arrays tempsJan3D = reshape(tempsJan,[24 7 4]); timesJan3D = reshape(timesJan,[24 7 4]); % Plot temperatures tempsJan2D = reshape(tempsJan3D,[24*7 4]); timesJan2D = reshape(timesJan3D,[24*7 4]); stem3((1:4).*ones(size(timesJan2D)), ... 1+mod(timesJan2D-1,7),tempsJan2D,".") view([-20 25]) set(gca,Ydir="reverse") grid on axis tight xticks(1:4) xlabel("Week Number") ylabel("Day of Week") zlabel("Temperature (°C)")
Peaks Across Day Hours (Rows)
Find the temperature peaks along all day hours, on the third day of the fourth week. By default, refinepeaks
assumes that the first dimension, which has more than three samples, is the dimension in which refine peaks. Specify the indices of the initially estimated peaks as a three-column matrix. Plot the estimated and refined peaks.
j = 3; % third day of week, third column k = 4; % Fourth week, fourth page [~,xPeaks] = findpeaks(tempsJan3D(:,j,k)); tiledlayout("horizontal") nexttile % Find and plot peaks across rows findpeaks(tempsJan3D(:,j,k),timesJan3D(:,j,k)) xlabel("Day number") ylim([-5 5]) nexttile % Refine peaks across rows, use xPeaks2D numeric matrix xPeaks2D = [0 j k] + [1 0 0].*xPeaks; refinepeaks(tempsJan3D,xPeaks2D) ylim([-5 5])
Define a logical array specifying a value of true
at the indices of the initially estimated peaks. Specify the logical array as peak location index input to refine the peaks. Plot the refined peaks. Both ways of specifying the location indices of the initially estimated peaks, either as numeric matrix or logical array, yield the same result.
nexttile % Refine peaks across rows use xPeaks3D logical matrix
xPeaks3D = false(size(tempsJan3D));
xPeaks3D(xPeaks,j,k) = true;
refinepeaks(tempsJan3D,xPeaks3D)
ylim([-5 5])
Peaks Across Week Days (Columns)
Find the temperature peaks at midnight and at 8 AM along all week days of the third week. Specify a logical array of initially estimated peaks to refine the peaks. Plot the estimated and refined peaks.
d = 2; % Second dimension (columns), across week days i1 = 1; % Midnight, first row i2 = 9; % 8 AM, ninth row k = 3; % Third week, third page [~,xPeaks1] = findpeaks(tempsJan3D(i1,:,k)); [~,xPeaks2] = findpeaks(tempsJan3D(i2,:,k)); tiledlayout("horizontal") nexttile % Find and plot peaks across columns findpeaks(tempsJan3D(i1,:,k),timesJan3D(i1,:,k)) hold on findpeaks(tempsJan3D(i2,:,k),timesJan3D(i2,:,k)) hold off ylim([-1.5 3.5]) xlabel("Day number") nexttile % Refine peaks across columns use xPeaks3D logical matrix xPeaks3D = false(size(tempsJan3D)); xPeaks3D(i1,xPeaks1,k) = true; xPeaks3D(i2,xPeaks2,k) = true; refinepeaks(tempsJan3D,xPeaks3D,Dimension=d) ylim([-1.5 3.5])
Peaks Across Weeks (Pages)
Find the temperature peaks at 7 AM during the second day of all weeks. Specify the indices of the initially estimated peaks as a three-column matrix. Plot the estimated and refined peaks.
d = 3; % Third dimension (pages), across weeks i = 8; % 7 AM, seventh row j = 2; % second day of week, second column tempsJanF8A = tempsJan3D(i,j,:); timesJanF8A = timesJan3D(i,j,:); [~,xPeaks] = findpeaks(tempsJanF8A(:)); tiledlayout("horizontal") nexttile % Find and plot peaks across pages findpeaks(tempsJanF8A(:),timesJanF8A(:)); ylim([-8 1]) xlabel("Day number") nexttile % Refine peaks across pages, use xPeaks2D numeric matrix xPeaks2D = [i j 0] + [0 0 1].*xPeaks; refinepeaks(tempsJan3D,xPeaks2D,Dimension=d) ylim([-8 1])
Input Arguments
y
— Signal amplitude data
vector | matrix | N-D array | timetable
Signal amplitude data, specified as a real-valued vector, matrix,
N-dimensional array, or MATLAB^{®} timetable with one- or two-dimensional data representing the input signal
amplitudes. You must specify y
so its representing data contains at
least a three-element vector. The refinepeaks
function uses three
points to refine each peak: the selected peak and two surrounding points.
Example: y = [100 77; 98.7 82.6; 95.2 86.5; 101.4 84.7; 99.1 84.3]
specifies a matrix representing two five-element signal amplitude
vectors.
Example: y = timetable([100;98.7;95.2;101.4;99.1],SampleRate=100)
specifies a timetable with one-dimensional data represented by a five-element column
vector.
Data Types: single
| double
| timetable
xPeaksIdx
— Location indices of initially estimated peaks
scalar, vector, or matrix of positive integers | N-D array of logical values
Location indices of the initially estimated peaks, specified as one of the following:
Scalar, vector or matrix of positive integers — The function assumes that
xPeaksIdx
represents the location indices of the initial peak estimates ofy
. The dimensionality ofxPeaksIdx
depends on how you specifyy
.If you specify
y
as a vector or timetable with one-column data,xPeaksIdx
must be either a scalar representing the index of one peak or a vector representing the indices of multiple peaks.If you specify
y
as a matrix or as an N-dimensional array,xPeaksIdx
must be either a row vector representing the index of one peak or a two-dimensional matrix representing the indices of multiple peaks.
xPeaksIdx
must be of size of M-by-N, wherexPeaksIdx
points at the indices of the M peaks of interest in the N-dimensional inputy
.N-dimensional array of logical values with the same size as
y
— The function assumes that any location index inxPeaksIdx
with a logical value of1
(true
) points to a peak iny
.
If any element of xPeaksIdx
that you specify selects the first
or last element of y
for the Dimension
of
interest, or if it selects a value that is not a peak, then
refinepeaks
does not perform peak processing. Instead,
refinepeaks
returns the same original (unrefined) peak values and
locations in yRPeaks
and xRPeaks
,
respectively.
Example: xPeaksIdx = 4
specifies that the fourth sample of a
vector y
is an initially estimated peak.
Example: xPeaksIdx = [4 1; 7 6; 5 8]
specifies three peaks to
refine: the fourth sample of the first column of y
, the seventh
sample of the sixth column of y
, and the fifth sample of the eighth
column of y
.
Data Types: single
| double
| logical
x
— Scaled locations of signal amplitude data
real-valued vector | datetime
vector | duration
vector
Scaled locations of the signal amplitude data, specified as a real-valued,
datetime
, or duration
vector.
The vector
x
must increase or decrease monotonically and uniformly. It must also have the same length as the number of elements iny
along theDimension
of interest.If you do not specify
x
,refinepeaks
generatesx
from the indices ofy
along theDimension
of interest. In other words,refinepeaks
setsx
to the vector1:Q
, whereQ
is the number of samples iny
along theDimension
of interest.You cannot specify
x
if you specifyy
as a timetable. When you specifyy
as a timetable,refinepeaks
assigns the time values fromy
tox
.
The refinepeaks
function uses x
to adjust xRPeaks
to match the desired signal scale.
xRPeaks
can represent time, frequency, angle, or any other
metric.
Example: x = seconds(0:0.01:1)
specifies the
duration
locations of the signal amplitude data points in a 101-element
vector y
. If you do not specify x
,
refinepeaks
assumes x = 1:101
, given the same
vector y
.
Example: x = 501:699
specifies the locations of the signal
amplitude data points in an N-D array y
with 200
samples across the dimension of interest. If you do not specify x
,
refinepeaks
assumes x = 1:200
, given the same
N-D array y
.
Data Types: single
| double
| datetime
| duration
Name-Value Arguments
Specify optional pairs of arguments as
Name1=Value1,...,NameN=ValueN
, where Name
is
the argument name and Value
is the corresponding value.
Name-value arguments must appear after other arguments, but the order of the
pairs does not matter.
Example: [yRPeaks,xRPeaks] =
refinepeaks(y,xPeaksIdx,x,Dimension=2,Method="NLS",MaxIterations=7)
refines
the peaks with locations in xPeaksIdx
across the second dimension of
the N-D array y
using scaled locations
x
, the nonlinear squares (NLS) method, and seven iterations at
maximum.
Method
— Peak refinement method
"QLS"
(default) | "NLS"
Peak refinement method, specified as one of these values:
"QLS"
—refinepeaks
uses the quadratic least squares (QLS) method to refine the peaks, representing each peak and the corresponding two surrounding points with a curve-fitting quadratic function. See QLS Method for additional input arguments."NLS"
—refinepeaks
uses the nonlinear least squares (NLS) fractional peak estimation method, representing each peak and the corresponding two surrounding points with a curve-fitting sinc function. Therefinepeaks
function uses the Gauss-Newton implementation and a sinc-function kernel to perform the curve-fitting operation. See NLS Method for additional input arguments.
For more information about these peak refinement methods, see Algorithms.
Data Types: char
| string
Dimension
— Dimension along which to refine peaks
positive integer scalar
Dimension along which to refine the peaks in y
, specified as
a positive integer scalar.
If you specify
y
as an N-dimensional matrix,Dimension
must be an integer scalar between 1 and N.If you do not specify
Dimension
,refinepeaks
sets the value to the first dimension with at least three elements. For example, if you specify a signal with amplitudesy
of size 2-by-1-by-2-by-6-by-7,refinepeaks
setsDimension
to4
.
Assume an N-dimensional array y
of size [s_{1}
s_{2} ⋯
s_{n – 1}
s_{n}
s_{n + 1} ⋯
s_{N – 1}
s_{N}] and you intend to refine peaks across the nth
dimension. When you specify Dimension=n
,
refinepeaks
refines each peak of y
assuming the following location indices.
Sample | Location index |
---|---|
Peak |
[p_{1} p_{2} ⋯ p_{n – 1} p_{n} p_{n + 1} ⋯ p_{N – 1} p_{N}] |
Left surrounding |
[p_{1} p_{2} ⋯ p_{n – 1} p_{n}–1 p_{n + 1} ⋯ p_{N – 1} p_{N}] |
Right surrounding |
[p_{1} p_{2} ⋯ p_{n – 1} p_{n}+1 p_{n + 1} ⋯ p_{N – 1} p_{N}] |
Data Types: single
| double
QLSKernel
— Amplitude preprocessing kernel
"none"
(default) | "log"
| "exp"
Amplitude preprocessing kernel for the QLS method, specified as one of these values:
"none"
—refinepeaks
usesy
to refine peaks."log"
—refinepeaks
useslog10(y)
to refine peaks."exp"
—refinepeaks
usesy.^k
to refine peaks, wherek
is the exponent specified inKernelPower
.
Dependencies
You must set Method
to "QLS"
to
specify QLSKernel
.
Data Types: char
| string
KernelPower
— Exponent of exponential QLS kernel
2
(default) | positive scalar up to 100
Exponent of the exponential QLS kernel, specified as a positive scalar.
Dependencies
You must set QLSKernel
to "exp"
to
specify KernelPower
.
Data Types: single
| double
LobeWidth
— Width of interpolating sinc curves on peaks
2
(default) | scalar | column vector
Width of the interpolating sinc curves where the peaks reside, specified as a
scalar or column vector with as many rows as peaks you specified to refine in
xPeaksIdx
. The minimum value of LobeWidth
must be equal or greater than 2.
LobeWidth
represents the estimate of the width of the lobes where the selected peaks are located in terms of number of samples.If you specify
LobeWidth
as a scalar,refinepeaks
applies the same lobe width to all the peak location entries inxPeaksIdx
.If you specify
LobeWidth
as a column vector,refinepeaks
applies the lobe width from each row ofLobeWidth
to its corresponding peak location entry inxPeaksIdx
.
The
refinepeaks
function usesLobeWidth
as an initial estimate for the function kernel in the iterative computation.You can obtain an initial estimate of
LobeWidth
from the initially estimated peaks usingfindpeaks
function. See Tips for more information.Although you do not need to specify exact values for
LobeWidth
, a good initial estimate can help the NLS method converge faster. Conversely, if the values inLobeWidth
are far from the actual values, the NLS method might not converge to the correct values.
Dependencies
You must set Method
to "NLS"
to
specify LobeWidth
.
Data Types: single
| double
MaxIterations
— Maximum number of iterations
8
(default) | positive integer scalar greater than or equal to 3
Maximum number of iterations, specified as a positive integer scalar greater than or equal to 3. For more information about the iteration process, see Peak Iterative Refinement.
Generally, between 5 and 10 iterations are needed to reach the final refined peak estimates.
Dependencies
You must set Method
to "NLS"
to
specify MaxIterations
.
Data Types: single
| double
SquaredErrorThreshold
— Maximum squared error threshold
eps
(default) | positive scalar
Maximum squared error threshold, specified as a positive scalar.
SquaredErrorThreshold
represents the threshold tolerance for
the sum of squared differences in amplitude between each three-point set (peak and two
surroundings) from y
and the corresponding amplitude evaluation
of the refined curve of y
that refinepeaks
fits using the NLS method.
For more information about the iteration process, see Peak Iterative Refinement.
Dependencies
You must set Method
to "NLS"
to
specify MaxIterations
.
Data Types: single
| double
Output Arguments
Tips
You can initially estimate signal peaks with findpeaks
, and then enhance their amplitudes and locations with
refinepeaks
.
Assume you have a signal with amplitudes y
and locations
x
. The following code snippet shows how you can estimate and refine
peaks from y
and x
.
[yPeaks,xPeaksIdx] = findpeaks(y); [yRPeaks,xRPeaks] = refinepeaks(y,xPeaksIdx,x)
Add name-value arguments for further customization. For example, you can specify
LobeWidth
from the width of the initially estimated
peaks.
[yPeaks,xPeaksIdx,xWidths] = findpeaks(y); [yRPeaks,xRPeaks] = refinepeaks(y,xPeaksIdx,x,Method="NLS", ... LobeWidth=max(xWidths,2))
Algorithms
For each selected peak with location index x_{PI}
from the vector of indices xPeaksIdx
, refinepeaks
identifies the indices and amplitudes of the peak y_{1}
and its two surrounding points y_{0} and
y_{2} from y
along
Dimension
. These three samples form the vector Y such that
$$Y=\left[\begin{array}{c}{y}_{0}\\ {y}_{1}\\ {y}_{2}\end{array}\right]=\left[\begin{array}{c}y[{x}_{\text{PI}}-1]\\ y[{x}_{\text{PI}}]\\ y[{x}_{\text{PI}}+1]\end{array}\right].$$
The refinepeaks
function uses the quadratic least square (QLS) or
nonlinear least square (NLS) fractional peak estimation algorithm to refine initially
estimated peaks from signal amplitude data.
QLS Algorithm
The QLS algorithm [2] preprocesses the signal amplitudes of Y, refines the peak location index by fitting and differentiating a quadratic polynomial, and then calculates the refined peak amplitude and scaled location.
Unless you specify QLSKernel
as "none"
, in
which case refinepeaks
skips the data preprocessing step, the QLS
algorithm starts by preprocessing the signal amplitude data Y using the QLSKernel
that you specify. While
logarithmic preprocessing ("log"
QLS kernel) is commonly used for
phased and radar applications, exponential preprocessing ("exp"
QLS
kernel) is used in other signal applications.
The function transforms the signal amplitudes of the peak and surrounding points for peak refinement, such that y_{0}, y_{1}, and y_{2} have these values:
Logarithmic Processing | Exponential Processing |
---|---|
y_{0} = log_{10}|y[x_{PI} – 1]| |
y_{0} = |y[x_{PI} – 1]|^{k} |
y_{1} = log_{10}|y[x_{PI}]| |
y_{1} = |y[x_{PI}]|^{k} |
y_{2} = log_{10}|y[x_{PI} + 1]| |
y_{2} = |y[x_{PI} + 1]|^{k} |
For logarithmic processing, specify the exponentk using
KernelPower
.
The function refines the peak location index by performing quadratic least square estimation using y_{0}, y_{1}, and y_{2} and differentiates the least-square quadratic polynomial to obtain a subsample estimate of the peak location, x_{RPI}, where
$${x}_{\text{RPI}}={x}_{\text{PI}}+\Delta {x}_{\text{PI}}={x}_{\text{PI}}-\frac{{y}_{2}-{y}_{0}}{2{y}_{0}-4{y}_{1}+2{y}_{2}}\text{.}$$
The function reverts the values of y_{0}, y_{1}, and y_{2} to their original values before preprocessing.
The function calculates the scaled location of the refined peak
x_{RPeak} from
x_{RPI} and from the scaled locations of signal
amplitude data x
, if you specify one.
If you specify the scaled locations of signal amplitude data
x
, then the function scales x_{RPI} for each peak and sets each element ofxRPeaks
from x_{RP}, such that$${x}_{\text{RPeak}}={x}_{1}+\Delta {x}_{\text{PI}}\left(\frac{{x}_{2}-{x}_{0}}{2}\right)={x}_{1}-\frac{({y}_{2}-{y}_{0})({x}_{2}-{x}_{0})}{4{y}_{0}-8{y}_{1}+4{y}_{2}},$$
where x_{0} = x[x_{RPI} − 1], x_{1} = x[x_{RPI}], and x_{2} = x[x_{RPI} + 1].
If you do not specify
x
, then the function sets each element ofxRPeaks
such that x_{RPeak} = x_{RPI}.
The function calculates the refined peak amplitude y_{RPeak} from the least-square quadratic polynomial such that
$${y}_{\text{RPeak}}=a{({x}_{\text{RPeak}})}^{2}+b({x}_{\text{RPeak}})+c,$$
where a = (y_{0} – 2y_{1} + y_{2})/2, b = (y_{2} – y_{0})/2, and c = y_{1}.
NLS Algorithm
The NLS algorithm [4] uses a sinc-function kernel to fit the curve formed by the initially estimated peak and its two surrounding points, followed by an minimum-error optimization, and then calculates the refined peak amplitude.
Assume the function f(X;λ) depends on the input column vector X and on the curve-fitting parameter vector λ such that f produces the column vector Y.
$$Y\cong f(X;\lambda )={\lambda}_{\text{\hspace{0.05em}}0}\text{sinc((}X-{\lambda}_{\text{\hspace{0.05em}}1}\text{)}{\lambda}_{\text{\hspace{0.05em}}2}\text{)s}$$
where X = [–1 0 1]^{T} and λ = [λ_{0} λ_{1} λ_{2}]^{T}.
The refinepeaks
function estimates the parameter vector λ to refine the peak location and amplitude such that λ minimizes the residual error S,
$$S={\displaystyle \sum _{i=0}^{2}{({y}_{i}-f({x}_{i},\lambda ))}^{2}}.$$
The refinepeaks
function uses Gauss–Newton optimization, for
which it finds the gradients with respect to λ_{0},
λ_{1}, and
λ_{2} such that
$$\begin{array}{l}\frac{\partial f(X;\lambda )}{\partial {\lambda}_{\text{\hspace{0.05em}}0}}=\text{sinc(}{\lambda}_{\text{\hspace{0.05em}}2}\text{(}X-{\lambda}_{\text{\hspace{0.05em}}1}\text{))}\\ \frac{\partial f(X;\lambda )}{\partial {\lambda}_{\text{\hspace{0.05em}}1}}=\frac{{\lambda}_{\text{\hspace{0.05em}}0}[\text{sinc(}{\lambda}_{\text{\hspace{0.05em}}2}\text{(}X-{\lambda}_{\text{\hspace{0.05em}}1}\text{))}-\text{cos(}\pi \text{\hspace{0.05em}}{\lambda}_{\text{\hspace{0.05em}}2}\text{(}X-{\lambda}_{\text{\hspace{0.05em}}1}\text{))}]}{X-{\lambda}_{\text{\hspace{0.05em}}1}}\\ \frac{\partial f(X;\lambda )}{\partial {\lambda}_{\text{\hspace{0.05em}}2}}=\frac{{\lambda}_{\text{\hspace{0.05em}}0}[\text{cos(}\pi \text{\hspace{0.05em}}{\lambda}_{\text{\hspace{0.05em}}2}\text{(}X-{\lambda}_{\text{\hspace{0.05em}}1}\text{))}-\text{sinc(}{\lambda}_{\text{\hspace{0.05em}}2}\text{(}X-{\lambda}_{\text{\hspace{0.05em}}1}\text{))}]}{{\lambda}_{\text{\hspace{0.05em}}2}}.\end{array}$$
Then, refinepeaks
performs the optimization as an
iterative process, with initial parameter values λ_{0} =
|y[x_{PI}]|, λ_{1} =
eps
, and λ_{2} =
1/L_{W}, where L_{W} is the peak width specified in LobeWidth
.
At the mth iteration, where m varies from 1 to
MaxIterations
, the function calculates the Jacobian matrix J and the residual error vector ΔY.
$$\begin{array}{l}J=\left[\begin{array}{ccc}\frac{\partial f(X;{\lambda}_{m})}{\partial {\lambda}_{\text{\hspace{0.05em}}0}}& \frac{\partial f(X;{\lambda}_{m})}{\partial {\lambda}_{\text{\hspace{0.05em}}1}}& \frac{\partial f(X;{\lambda}_{m})}{\partial {\lambda}_{\text{\hspace{0.05em}}2}}\end{array}\right],\\ \Delta Y=Y-f(X;{\lambda}_{m}).\end{array}$$
Since JΔλ = ΔY, solving the equation system derives in the calculation of the variation of the parameter vector λ such that
$$\begin{array}{l}\Delta \lambda ={(}^{{J}^{\text{T}}}{J}^{\text{T}}\Delta Y,\\ {\lambda}_{m+1}={\lambda}_{m}+\Delta \lambda .\end{array}$$
The function calculates the residual error S for the mth iteration from the squared sum of the elements in ΔY.
The iteration process stops as soon as one of these occurs:
The
MaxIterations
-th iteration is complete.The residual error S is less or equal to
SquaredErrorThreshold
.
The refinepeaks
function first calculates the subsample estimate
of the peak location x_{RPI} using the final value
of λ_{1} from the Peak Iterative Refinement process, where
$${x}_{\text{RPI}}={x}_{\text{PI}}+\Delta {x}_{\text{PI}}={x}_{\text{PI}}+{\lambda}_{1}\text{.}$$
Then, refinepeaks
calculates the scaled location of the refined
peak x_{RPeak} from
x_{RPI} and from the scaled locations of signal
amplitude data x
, if you specify one.
If you specify the scaled locations of signal amplitude data
x
, then the function scales x_{RPI} for each peak and sets each element ofxRPeaks
from x_{RP} such that$${x}_{\text{RPeak}}={x}_{1}+\Delta {x}_{\text{PI}}\left(\frac{{x}_{2}-{x}_{0}}{2}\right)={x}_{1}+{\lambda}_{1}\left(\frac{{x}_{2}-{x}_{0}}{2}\right),$$
where x_{0} = x[x_{RPI} − 1], x_{1} = x[x_{RPI}], and x_{2} = x[x_{RPI} + 1].
If you do not specify
x
, then the function sets each element ofxRPeaks
such that x_{RPeak} = x_{RPI}.
The refinepeaks
function calculates the refined peak amplitude
y_{RPeak} from the fitted sinc function such that
$${y}_{\text{RPeak}}={\lambda}_{\text{\hspace{0.05em}}0}\text{sinc((}{x}_{\text{RPeak}}-{\lambda}_{\text{\hspace{0.05em}}1}\text{)}{\lambda}_{\text{\hspace{0.05em}}2}\text{)},$$
where the coefficients λ_{0}, λ_{1}, and λ_{2} are the final values from the Peak Iterative Refinement process.
References
[1] Richards, Mark A., James A. Scheer, and William A. Holm, eds. Principles of Modern Radar: Basic Principles. Institution of Engineering and Technology, 2010.
[2] Moddemeijer, R. “On the Determination of the Position of Extrema of Sampled Correlators.” IEEE^{®} Transactions on Signal Processing 39, no. 1 (January 1991): 216–19.
[3] Sharp, I., K. Yu, and Y. J. Guo. “Peak and Leading Edge Detection for Time-of-Arrival Estimation in Band-Limited Positioning Systems.” IET Communications 3, no. 10 (October 1, 2009): 1616–27.
[4] Prager, Samuel, Mark S. Haynes, and Mahta Moghaddam. “Wireless Subnanosecond RF Synchronization for Distributed Ultrawideband Software-Defined Radar Networks.” IEEE Transactions on Microwave Theory and Techniques 68, no. 11 (November 2020): 4787–4804.
Extended Capabilities
C/C++ Code Generation
Generate C and C++ code using MATLAB® Coder™.
Version History
Introduced in R2024b
See Also
findpeaks
| islocalmax
| max
MATLAB Command
You clicked a link that corresponds to this MATLAB command:
Run the command by entering it in the MATLAB Command Window. Web browsers do not support MATLAB commands.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list:
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)