refinepeaks

Refine peak value and location estimates

Since R2024b

collapse all in page

Syntax

[yRPeaks,xRPeaks] = refinepeaks(y,xPeaksIdx)

[yRPeaks,xRPeaks] = refinepeaks(y,xPeaksIdx,x)

[___] = refinepeaks(___,Name=Value)

refinepeaks(___)

Description

[yRPeaks,xRPeaks] = refinepeaks(y,xPeaksIdx) refines peaks based on a signal with amplitudes 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.

example

[yRPeaks,xRPeaks] = refinepeaks(y,xPeaksIdx,x) specifies the scaled locations x of the signal amplitude data. The refinepeaks function refines the peaks using quadratic least squares processing and adjusts xRpeaks to match the desired signal scale from x.

example

[___] = refinepeaks(___,Name=Value) 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.

example

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.

example

Examples

collapse all

Refine Known Peak from Signal

Open Live Script

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)

Figure contains an axes object. The axes object contains an object of type line.

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)

Figure contains an axes object. The axes object with title Refined Peaks, xlabel Update To The Peak Location, ylabel Amplitude contains 4 objects of type scatter, line. This object represents Peak 1.

Quadratic Least Squares Peak Refinement

Open Live Script

Generate a Rayleigh distribution $P = f (x | b) = \frac{x}{b^{2}} e^{- \frac{x^{2}}{2 b^{2}}}$ with a scale parameter b set to $π$ .

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)

Figure contains an axes object. The axes object contains an object of type line.

Calculate the theoretical and estimated peak amplitude and location. The peak is located about half of a sample from the expected peak scaled location, $π$ .

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

Open Live Script

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")

Figure contains an axes object. The axes object with xlabel Frequency (kHz), ylabel Magnitude contains 2 objects of type line. One or more of the lines displays its values using only markers These objects represent Full signal, Selected samples.

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

Figure contains an axes object. The axes object with title Refined Peaks, xlabel Update To The Peak Location, ylabel Amplitude contains 13 objects of type scatter, line, constantline. These objects represent Peak 1, Peak 2.

Find and Refine Peaks from Three-Dimensional Data

Open Live Script

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)")

Figure contains an axes object. The axes object with xlabel Week Number, ylabel Day of Week contains an object of type stem.

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

Figure contains 2 axes objects. Axes object 1 with xlabel Day number contains 4 objects of type line. One or more of the lines displays its values using only markers Axes object 2 with title Refined Peaks, xlabel Update To The Peak Location, ylabel Amplitude contains 11 objects of type scatter, line. These objects represent Peak 1, Peak 2, Peak 3, Peak 4.

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

Figure contains 2 axes objects. Axes object 1 with xlabel Day number contains 2 objects of type line. One or more of the lines displays its values using only markers Axes object 2 with title Refined Peaks, xlabel Update To The Peak Location, ylabel Amplitude contains 4 objects of type scatter, line. This object represents Peak 1.

Input Arguments

collapse all

`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 of y. The dimensionality of xPeaksIdx depends on how you specify y.
- 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, where xPeaksIdx points at the indices of the M peaks of interest in the N-dimensional input y.
N-dimensional array of logical values with the same size as y — The function assumes that any location index in xPeaksIdx with a logical value of 1 (true) points to a peak in y.

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 in y along the Dimension of interest.
If you do not specify x, refinepeaks generates x from the indices of y along the Dimension of interest. In other words, refinepeaks sets x to the vector 1:Q, where Q is the number of samples in y along the Dimension of interest.
You cannot specify x if you specify y as a timetable. When you specify y as a timetable, refinepeaks assigns the time values from y to x.

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.

Peak Refinement

collapse all

`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. The refinepeaks 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 amplitudes y of size 2-by-1-by-2-by-6-by-7, refinepeaks sets Dimension to 4.

Assume an N-dimensional array y of size [s₁ s₂ ⋯ 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₁ p₂ ⋯ p_{n – 1} p_n p_{n + 1} ⋯ p_{N – 1} p_N]
Left surrounding	[p₁ p₂ ⋯ p_{n – 1} p_n–1 p_{n + 1} ⋯ p_{N – 1} p_N]
Right surrounding	[p₁ p₂ ⋯ p_{n – 1} p_n+1 p_{n + 1} ⋯ p_{N – 1} p_N]

Data Types: single | double

QLS Method

collapse all

`QLSKernel` — Amplitude preprocessing kernel
`"none"` (default) | `"log"` | `"exp"`

Amplitude preprocessing kernel for the QLS method, specified as one of these values:

"none" — refinepeaks uses y to refine peaks.
"log" — refinepeaks uses log10(y) to refine peaks.
"exp" — refinepeaks uses y.^k to refine peaks, where k is the exponent specified in KernelPower.

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

NLS Method

collapse all

`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 in xPeaksIdx.
- If you specify LobeWidth as a column vector, refinepeaks applies the lobe width from each row of LobeWidth to its corresponding peak location entry in xPeaksIdx.
The refinepeaks function uses LobeWidth 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 using findpeaks 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 in LobeWidth 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

collapse all

`yRPeaks` — Amplitudes of refined peaks
column vector

Amplitudes of the refined peaks, returned as a column vector.

If you specify y or xPeaksIdx in single precision, then refinepeaks performs the peak refinement using single-precision arithmetic and returns yRPeaks in single precision.

`xRPeaks` — Scaled locations of refined peaks
column vector

Scaled locations of the refined peaks, returned as a column vector.

If you specify y or xPeaksIdx in single precision, then refinepeaks performs the peak refinement using single-precision arithmetic, and returns xRPeaks in single precision.

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

collapse all

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₁ and its two surrounding points y₀ and y₂ from y along Dimension. These three samples form the vector Y such that

$Y = [\begin{matrix} y_{0} \\ y_{1} \\ y_{2} \end{matrix}] = [\begin{matrix} y [x_{PI} - 1] \\ y [x_{PI}] \\ y [x_{PI} + 1] \end{matrix}] .$

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.

Signal Amplitude Data Preprocessing

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₀, y₁, and y₂ have these values:

Logarithmic Processing	Exponential Processing
y₀ = log₁₀\|y[x_PI – 1]\|	y₀ = \|y[x_PI – 1]\|^k
y₁ = log₁₀\|y[x_PI]\|	y₁ = \|y[x_PI]\|^k
y₂ = log₁₀\|y[x_PI + 1]\|	y₂ = \|y[x_PI + 1]\|^k

For logarithmic processing, specify the exponentk using KernelPower.

Peak Location Index Refinement

The function refines the peak location index by performing quadratic least square estimation using y₀, y₁, and y₂ and differentiates the least-square quadratic polynomial to obtain a subsample estimate of the peak location, x_RPI, where

$x_{RPI} = x_{PI} + Δ x_{PI} = x_{PI} - \frac{y_{2} - y_{0}}{2 y_{0} - 4 y_{1} + 2 y_{2}} .$

The function reverts the values of y₀, y₁, and y₂ to their original values before preprocessing.

Refined Peak Scaled Location

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 of xRPeaks from x_RP, such that

$x_{RPeak} = x_{1} + Δ x_{PI} (\frac{x_{2} - x_{0}}{2}) = x_{1} - \frac{(y_{2} - y_{0}) (x_{2} - x_{0})}{4 y_{0} - 8 y_{1} + 4 y_{2}},$
where x₀ = x[x_RPI − 1], x₁ = x[x_RPI], and x₂ = x[x_RPI + 1].
If you do not specify x, then the function sets each element of xRPeaks such that x_RPeak = x_RPI.

Refined Peak Amplitude

The function calculates the refined peak amplitude y_RPeak from the least-square quadratic polynomial such that

$y_{RPeak} = a {(x_{RPeak})}^{2} + b (x_{RPeak}) + c,$

where a = (y₀ – 2y₁ + y₂)/2, b = (y₂ – y₀)/2, and c = y₁.

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 ≅ f (X; λ) = λ_{0} sinc((X - λ_{1}) λ_{2})s$

where X = [–1 0 1]^T and λ = [λ₀ λ₁ λ₂]^T.

Peak Iterative Refinement

The refinepeaks function estimates the parameter vector λ to refine the peak location and amplitude such that λ minimizes the residual error S,

$S = \sum_{i = 0}^{2} {(y_{i} - f (x_{i}, λ))}^{2} .$

The refinepeaks function uses Gauss–Newton optimization, for which it finds the gradients with respect to λ₀, λ₁, and λ₂ such that

$\begin{array}{l} \frac{\partial f (X; λ)}{\partial λ_{0}} = sinc(λ_{2} (X - λ_{1})) \\ \frac{\partial f (X; λ)}{\partial λ_{1}} = \frac{λ_{0} [sinc(λ_{2} (X - λ_{1})) - cos(π λ_{2} (X - λ_{1}))]}{X - λ_{1}} \\ \frac{\partial f (X; λ)}{\partial λ_{2}} = \frac{λ_{0} [cos(π λ_{2} (X - λ_{1})) - sinc(λ_{2} (X - λ_{1}))]}{λ_{2}} . \end{array}$

Then, refinepeaks performs the optimization as an iterative process, with initial parameter values λ₀ = |y[x_PI]|, λ₁ = eps, and λ₂ = 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 = [\begin{matrix} \frac{\partial f (X; λ_{m})}{\partial λ_{0}} & \frac{\partial f (X; λ_{m})}{\partial λ_{1}} & \frac{\partial f (X; λ_{m})}{\partial λ_{2}} \end{matrix}], \\ Δ Y = Y - f (X; λ_{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} Δ λ = (^{J^{T} J) - 1} J^{T} Δ Y, \\ λ_{m + 1} = λ_{m} + Δ λ . \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.

Refined Peak Scaled Location

The refinepeaks function first calculates the subsample estimate of the peak location x_RPI using the final value of λ₁ from the Peak Iterative Refinement process, where

$x_{RPI} = x_{PI} + Δ x_{PI} = x_{PI} + λ_{1} .$

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 of xRPeaks from x_RP such that

$x_{RPeak} = x_{1} + Δ x_{PI} (\frac{x_{2} - x_{0}}{2}) = x_{1} + λ_{1} (\frac{x_{2} - x_{0}}{2}),$
where x₀ = x[x_RPI − 1], x₁ = x[x_RPI], and x₂ = x[x_RPI + 1].
If you do not specify x, then the function sets each element of xRPeaks such that x_RPeak = x_RPI.

Refined Peak Amplitude

The refinepeaks function calculates the refined peak amplitude y_RPeak from the fitted sinc function such that

$y_{RPeak} = λ_{0} sinc((x_{RPeak} - λ_{1}) λ_{2}),$

where the coefficients λ₀, λ₁, and λ₂ 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

Logarithmic Processing	Exponential Processing
y₀ = log₁₀\|y[x_PI – 1]\|	y₀ = \|y[x_PI – 1]\|^k
y₁ = log₁₀\|y[x_PI]\|	y₁ = \|y[x_PI]\|^k
y₂ = log₁₀\|y[x_PI + 1]\|	y₂ = \|y[x_PI + 1]\|^k

refinepeaks

Syntax

Description

Examples

Refine Known Peak from Signal

Quadratic Least Squares Peak Refinement

Nonlinear Least Squares Peak Refinement

Find and Refine Peaks from Three-Dimensional Data

Input Arguments

y — Signal amplitude data vector | matrix | N-D array | timetable

xPeaksIdx — Location indices of initially estimated peaks scalar, vector, or matrix of positive integers | N-D array of logical values

x — Scaled locations of signal amplitude data real-valued vector | datetime vector | duration vector

Name-Value Arguments

Method — Peak refinement method "QLS" (default) | "NLS"

Dimension — Dimension along which to refine peaks positive integer scalar

QLSKernel — Amplitude preprocessing kernel "none" (default) | "log" | "exp"

Dependencies

KernelPower — Exponent of exponential QLS kernel 2 (default) | positive scalar up to 100

Dependencies

LobeWidth — Width of interpolating sinc curves on peaks 2 (default) | scalar | column vector

Dependencies

MaxIterations — Maximum number of iterations 8 (default) | positive integer scalar greater than or equal to 3

Dependencies

SquaredErrorThreshold — Maximum squared error threshold eps (default) | positive scalar

Dependencies

Output Arguments

yRPeaks — Amplitudes of refined peaks column vector

xRPeaks — Scaled locations of refined peaks column vector

Tips

Algorithms

QLS Algorithm

NLS Algorithm

References

Extended Capabilities

C/C++ Code Generation Generate C and C++ code using MATLAB® Coder™.

Version History

See Also

`y` — Signal amplitude data
vector | matrix | N-D array | timetable

`xPeaksIdx` — Location indices of initially estimated peaks
scalar, vector, or matrix of positive integers | N-D array of logical values

`x` — Scaled locations of signal amplitude data
real-valued vector | `datetime` vector | `duration` vector

`Method` — Peak refinement method
`"QLS"` (default) | `"NLS"`

`Dimension` — Dimension along which to refine peaks
positive integer scalar

`QLSKernel` — Amplitude preprocessing kernel
`"none"` (default) | `"log"` | `"exp"`

`KernelPower` — Exponent of exponential QLS kernel
`2` (default) | positive scalar up to 100

`LobeWidth` — Width of interpolating sinc curves on peaks
`2` (default) | scalar | column vector

`MaxIterations` — Maximum number of iterations
`8` (default) | positive integer scalar greater than or equal to 3

`SquaredErrorThreshold` — Maximum squared error threshold
`eps` (default) | positive scalar

`yRPeaks` — Amplitudes of refined peaks
column vector

`xRPeaks` — Scaled locations of refined peaks
column vector

C/C++ Code Generation
Generate C and C++ code using MATLAB® Coder™.