findpeaks
查找局部最大值
语法
说明
[___] = findpeaks(___, 支持上述语法中的任何输入参量，且可使用名称-值参量指定选项。Name=Value)
findpeaks(___) 不带输出参量时用于绘制信号并叠加峰值。
示例
定义一个具有三个峰值的向量。
data = [25 8 15 5 6 10 10 3 1 20 7];
找出局部最大值。findpeaks 函数按出现顺序返回峰值。该函数在峰值查找中不包括 data 的端点，因为没有左邻点或右邻点可以与端点进行比较。
pks = findpeaks(data)
pks = 1×3
15 10 20
使用不带输出参量的 findpeaks 来显示峰值。对于平峰，该函数仅突出显示索引最低的点。
findpeaks(data)
创建一个由几条高斯曲线之和组成的信号。指定每条曲线的位置、高度和宽度。
x = linspace(0,1,1000); Pos = [1 2 3 5 7 8]'/10; Hgt = [3 4 4 2 2 3]'; Wdt = [2 6 3 3 4 6]'/100; y = sum(Hgt.*(exp(-((x-Pos)./Wdt).^2)),1);
使用具有默认设置的 findpeaks 查找信号的峰值及其位置。
[pks,locs] = findpeaks(y,x);
使用 findpeaks 绘制峰值。标注每个峰值。
findpeaks(y,x) text(locs+.02,pks,num2str((1:numel(pks))'))
从最高到最低对峰值进行排序。
[psor,lsor] = findpeaks(y,x,SortStr="descend");
findpeaks(y,x)
text(lsor+.02,psor,num2str((1:numel(psor))'))
创建一个由分布在一个完整余弦周期上的几条高斯曲线之和组成的信号。指定每条曲线的位置、高度和宽度。
x = linspace(0,1,1000); base = 4*cos(2*pi*x); Pos = [1 2 3 5 7 8]'/10; Hgt = [3 7 5 5 4 5]'; Wdt = [1 3 3 4 2 3]'/100; y = sum(Hgt.*(exp(-((x-Pos)./Wdt).^2)),1) + base;
使用 findpeaks 定位并绘制相对高差至少为 4 的峰值。只有最高和最低峰值满足最小相对高差条件。
[~,locs4,~,proms4] = findpeaks(y,x,MinPeakProminence=4)
locs4 = 1×2
0.1982 0.4995
proms4 = 1×2
5.5773 4.4171
findpeaks(y,x,MinPeakProminence=4,Annotate="extents")
显示所有峰值的位置、相对高差和相对高差一半处的峰宽。第二个和第四个峰值的相对高差大于或等于 4。
[~,locs,widths,proms] = findpeaks(y,x)
locs = 1×6
0.1001 0.1982 0.2983 0.4995 0.7017 0.8018
widths = 1×6
0.0154 0.0431 0.0377 0.0625 0.0274 0.0409
proms = 1×6
2.6816 5.5773 3.1448 4.4171 2.9191 3.6363
从太阳黑子平均计数数据中提取峰值并估计平均太阳黑子周期。
加载文件 sunspot.dat，其中包含从 1700 年到 1987 年每年观测到的太阳黑子的平均数量。按年份绘制太阳黑子数量并叠加最大值。放大 1855 年到 1920 年之间的时间段。
load sunspot.dat year = sunspot(:,1); avSpots = sunspot(:,2); findpeaks(avSpots,year) xlabel("Year") xlim([1855 1920])
太阳黑子是一种周期性出现的现象。已知其数量大约每 11 年达到峰值。
通过找到峰值之间的平均间距来估计周期持续时间。首先忽略不是最大值的峰值。将可接受的峰间间距限制为大于 6 年的周期。
findpeaks(avSpots,year,MinPeakDistance=6)
xlabel("Year")
xlim([1855 1920])
使用年度数据创建一个 datetime 数组。假设太阳黑子是在每年 3 月 20 日（接近春分）统计的。找出太阳黑子的峰值年份。使用 years 函数将最小峰间距指定为 duration。
ty = datetime(year,3,20); [pk,lk] = findpeaks(avSpots,ty,MinPeakDistance=years(6)); plot(ty,avSpots,lk,pk,"o") xlabel("Year")
计算平均太阳黑子周期。估计的平均太阳黑子周期与已知的 11 年周期非常接近。
dttmCycle = years(mean(diff(lk)))
dttmCycle = 10.9600
您可以找到满足特定特征标准的峰值，其中包括最小或最大峰值距离或高度、相对高差以及峰值与其邻点差异的阈值。使用 findpeaks 约束单个或多个峰值特征。
通过音频信号的特征约束峰值
加载一段以 7418 Hz 采样的音频信号。选择 200 个采样。
load mtlb
select = mtlb(1001:1200);找到并突出显示音频信号中的所有峰值。然后，找到并突出显示间隔至少为 5 毫秒的峰值。为了应用最小峰间距约束，findpeaks 选择信号中的最高峰值，并消除在其 5 毫秒范围内的所有峰值。然后，该函数对剩余峰值中的最高峰值重复该过程并进行迭代，直到处理完要考虑的所有峰值。
figure tiledlayout("vertical") nexttile findpeaks(select,Fs) legend("No Constraints",Location="southwest") nexttile lgn = @(txt) legend(txt,FontName="FixedWidth",Location="southwest"); findpeaks(select,Fs,MinPeakDistance=0.005) lgn("MinPeakDistance")
找出并绘制分别满足以下每个约束的峰值：
找出振幅至少为 1 V 的峰值。
找出比相邻采样至少高 1 V 的峰值。
找出信号在抵达更高峰值前于任一侧下降至少 1 V 的峰值。
figure tiledlayout("vertical") nexttile findpeaks(select,Fs,MinPeakHeight=1) lgn("MinPeakHeight") nexttile findpeaks(select,Fs,Threshold=1) lgn("Threshold") nexttile findpeaks(select,Fs,MinPeakProminence=1) lgn("MinPeakProminence")
通过啁啾信号约束多个峰值特征
找到并突出显示包含 1000 个采样的啁啾信号中宽度在 2.5 到 4 毫秒之间的所有峰值。
N = 1000; x = sin(2*pi*(1:N)/N + (10*(1:N)/N).^2); figure findpeaks(x,Fs,Annotate="extents", ... MinPeakWidth=2.5e-3,MaxPeakWidth=4e-3) title("Peaks with Multiple Constraints")
如果数据大于给定的饱和点，传感器可以返回裁剪后的读数。您可以选择忽略这些没有意义的极值，也可以选择将它们纳入您的分析中。
生成一个由频率为 5 Hz 和 3 Hz 的三角函数的乘积组成的信号，嵌入方差为 0.1² 的高斯白噪声中。以 100 Hz 的速率对信号进行采样，持续时间为一秒。重置随机数生成器以获得可重现的结果。
rng("default")
Fs = 100;
t = 0:1/Fs:1;
s = sin(2*pi*5*t).*sin(2*pi*3*t) + randn(size(t))/10;通过截断绝对值大于 0.32 的指定边界的每个读数来合成经过饱和处理的测量值。然后，对信号进行归一化，使其范围在 0.1 到 1 之间。将归一化幅值转换为分贝并绘制结果。
bnd = 0.32; s(s>bnd) = bnd; s(s<-bnd) = -bnd; sNorm = 0.1+0.9*(s-min(s))/(max(s)-min(s)); sig = mag2db(sNorm); plot(t,sig) xlabel("Time (s)") ylabel("Magnitude (dB)") ylim([-25 5])
定位信号的波峰和波谷。findpeaks 仅识别每个平峰的上升沿。
[pk,lcp] = findpeaks(sig,t); % Peaks [vl,lcv] = findpeaks(-sig,t); % Valleys hold on plot(lcp,pk,"x",lcv,-vl,"*")
使用 Threshold 名称-值参量排除平坦波峰和波谷。要求波峰或波谷与其邻点之间的最小幅值差为 0.1 dB。
[pkt,lctp] = findpeaks(sig,t,Threshold=0.1); [vlt,lctv] = findpeaks(-sig,t,Threshold=0.1); plot(lctp,pkt,"o",LineWidth=1.2) plot(lctv,-vlt,"d",LineWidth=1.2)
创建一个由几条高斯曲线之和组成的信号。指定每条曲线的位置、高度和宽度。
x = linspace(0,1,1000); Pos = [1 2 3 5 7 8]'/10; Hgt = [7 6 3 2 2 3]'; Wdt = [3 8 4 3 4 6]'/100; y = sum(Hgt.*(exp(-((x-Pos)./Wdt).^2)),1);
使用相对高差的一半处和半高处作为参考来测量峰宽。
tiledlayout("flow") nexttile findpeaks(y,x,Annotate="extents") title("Half-Prominence Peak Widths") nexttile findpeaks(y,x,Annotate="extents",WidthReference="halfheight") title("Half-Height Peak Widths")
选择在 x 轴上间隔至少 0.5 个单位的最高波峰。使用相对高差的一半处和半高处作为参考来测量峰宽。
figure tiledlayout("flow") nexttile findpeaks(y,x,MinPeakDistance=0.5,Annotate="extents") title("Half-Prominence Peak Widths") nexttile findpeaks(y,x,MinPeakDistance=0.5,Annotate="extents", ... WidthReference="halfheight") title("Half-Height Peak Widths")
只有第一个和最后一个峰满足最小间距条件，因此图中显示的宽度对应于这两个峰。每个峰的范围保持不变，因此无论指定何种条件以及是否选中该峰，峰宽都保持其值不变。
使用具有 sinc 函数核的非线性最小二乘法获得一个信号中两个主要峰值的细化峰值位置估计。
生成信号
线性 FM 波形的雷达脉冲压缩产生一个正弦形状的频谱，其中峰值的频率位置与雷达和探测对象之间的距离成正比。您可以先使用findpeaks估计峰值位置和振幅，然后使用refinepeaks增强您的估计。此示例查找合成无噪声脉冲压缩信号的峰值振幅和位置，并使用 refinepeaks 改进估计。
生成一个由两个正弦波形组成的信号，峰值分别为在 4.76 kHz 和 35.8 kHz 处的 1 和 1.5。将频率间距设为 2.5 Hz。
aTg = [1 1.5]; fTg = 1e3*[4.76 35.8]; freqkHzFull = (0:0.0025:50)'; waveFull = abs(sinc([1 0.5].*(freqkHzFull-fTg/1e3)))*aTg';
以 200 为因子对信号下采样，因此采样之间的频率间距为 0.5 kHz。此示例细化下采样信号峰值的振幅和位置估计，并将改进的估计与原始信号中的值进行比较。
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")
使用非线性最小二乘细化峰值
使用 findpeaks 对信号的两个最高峰值的振幅、位置和半高宽度进行初始估计。
[PV,PL,PW] = findpeaks(wave,NPeaks=2, ... SortStr="descend",WidthReference="halfheight"); table(PV,freq(PL), RowNames="Peak estimate "+(1:numel(PV)), ... VariableNames=["Amplitude" "Frequency (kHz)"])
ans=2×2 table
Amplitude Frequency (kHz)
_________ _______________
Peak estimate 1 1.4824 36
Peak estimate 2 0.9374 5
使用 refinepeaks 增强使用非线性最小二乘 (NLS) 方法的峰值估计。指定信号的频率点和峰宽。峰值明显更接近预期值 1.5 和 1，而频率位置分别很好地逼近 35.8 kHz 和 4.76 kHz。
LW = max(PW,2); [Ypk,Xpk] = refinepeaks(wave,PL,freq,Method="NLS",LobeWidth=LW); table(Ypk,Xpk, RowNames="Refined peak "+(1:numel(Ypk)), ... VariableNames=["Amplitude" "Frequency (kHz)"])
ans=2×2 table
Amplitude Frequency (kHz)
_________ _______________
Refined peak 1 1.5063 35.8
Refined peak 2 1.0163 4.7628
在 y 轴上绘制细化峰值的振幅，在 x 轴上绘制与初始峰值估计相比的更新后的峰值位置。两个初始估计的峰值和它们周围的两个采样都相隔 0.5 kHz。由实心圆指示的细化峰值显示实际峰值位置与初始估计的峰值位置的对比以及校正后的振幅。
refinepeaks(wave,PL,freq,Method="NLS",LobeWidth=LW) yline(aTg) % Theoretical peak amplitudes errorBounds = aTg.*(1+0.03*[-1;1]); yline(errorBounds(:),":") % ±3% error bounds legend("Peak "+[1 2])
输入参数
名称-值参数
输出参量
详细信息
提示
您可以使用 findpeaks 初始估计信号峰值，然后使用 refinepeaks 增强其振幅和位置。
假设您有振幅为 y 且位置为 x 的信号。以下代码片段显示如何根据 y 和 x 估计和细化峰值。
[yPeaks,xPeaksIdx] = findpeaks(y); [yRPeaks,xRPeaks] = refinepeaks(y,xPeaksIdx,x)
