MATLAB中findpeaks函数使用

目录

语法

说明

示例

寻找向量中的峰值

查找峰值及其位置

---

        findpeaks函数的功能是寻找局部最大值。

语法

pks = findpeaks(data)
[pks,locs] = findpeaks(data)
[pks,locs,w,p] = findpeaks(data)
[___] = findpeaks(data,x)
[___] = findpeaks(data,Fs)
[___] = findpeaks(___,Name,Value)
findpeaks(___)

说明

        pks = findpeaks(data)返回带有输入信号向量data的局部最大值(峰值)的向量。局部峰值是大于其两个相邻样本或等于Inf的数据样本。峰值按出现的顺序输出。不包括非Inf信号端点。如果峰是平的，函数只返回指数最低的点。

        [pks，locs] = findpeaks(data)还会返回峰值出现的索引。

        [pks，locs，w，p] = findpeaks(data)还返回峰的宽度作为向量w，峰的突出度作为向量p。

        [___] = findpeaks(data，x)将x指定为位置向量，并返回以前语法中的任何输出参数。loc和w用x表示。

        [___] = findpeaks(data，Fs)指定数据的采样速率Fs。假设第一个数据样本是在时间零点采集的。loc和w转换为时间单位。

        [___] = find peaks(_ _ _，Name，Value)指定使用名称-值参数以及以前语法中的任何输入参数的选项。

        不带输出参数的findpeaks(___)绘制信号并覆盖峰值。

示例

寻找向量中的峰值

        定义一个有三个峰值的矢量并画出来。 

data = [25 8 15 5 6 10 10 3 1 20 7];
plot(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;
 
for n = 1:length(Pos)
    Gauss(n,:) = Hgt(n)*exp(-((x - Pos(n))/Wdt(n)).^2);
end
 
PeakSig = sum(Gauss);

        画出各条曲线及其总和。

plot(x,Gauss,'--',x,PeakSig)

        如图所示：

        使用默认设置的findpeaks查找信号的峰值及其位置。 

[pks,locs] = findpeaks(PeakSig,x);

        使用findpeaks绘制峰并标记它们。

findpeaks(PeakSig,x)
 
text(locs+.02,pks,num2str((1:numel(pks))'))

        如图所示：

        从最高到最低对山峰进行排序。 

[psor,lsor] = findpeaks(PeakSig,x,'SortStr','descend');
 
findpeaks(PeakSig,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;
 
for n = 1:length(Pos)
    Gauss(n,:) =  Hgt(n)*exp(-((x - Pos(n))/Wdt(n)).^2);
end
 
PeakSig = sum(Gauss)+base;

        画出各条曲线及其总和。

plot(x,Gauss,'--',x,PeakSig,x,base)

        如图所示：

        使用findpeaks定位并绘制显著度至少为4的峰。 

findpeaks(PeakSig,x,'MinPeakProminence',4,'Annotate','extents')

        如图所示：

        只有最高和最低的峰值满足该条件。 显示所有峰的日珥和半日珥宽度。 

[pks,locs,widths,proms] = findpeaks(PeakSig,x);
widths
widths = 1×6
 
    0.0154    0.0431    0.0377    0.0625    0.0274    0.0409
 
 
 
 
proms
proms = 1×6
 
    2.6816    5.5773    3.1448    4.4171    2.9191    3.6363

寻找分离度最小的峰

        太阳黑子是一种循环现象。已知它们的数量大约每11年达到一个峰值。 加载文件spocitizen . dat，其中包含了从1700年到1987年每年观测到的太阳黑子的平均数量。找到并画出最大值。

load sunspot.dat
 
year = sunspot(:,1);
avSpots = sunspot(:,2);
 
findpeaks(avSpots,year)

        如图所示：

        通过忽略彼此非常接近的峰值来改进对周期持续时间的估计。再次查找并绘制峰，但现在将可接受的峰-峰分离限制为大于6年的值。 

findpeaks(avSpots,year,'MinPeakDistance',6)

        如图所示：

        使用findpeaks返回的峰位置来计算最大值之间的平均间隔。 

[pks,locs] = findpeaks(avSpots,year,'MinPeakDistance',6);
 
meanCycle = mean(diff(locs))
 
 
meanCycle = 10.9600

        使用年份数据创建日期时间数组。假设太阳黑子在每年的3月20日被计数，接近春分。找出太阳黑子的高峰期。使用年数函数将最小峰间距指定为持续时间。

ty = datetime(year,3,20);
 
[pk,lk] = findpeaks(avSpots,ty,'MinPeakDistance',years(6));
 
plot(ty,avSpots,lk,pk,'o')

        如图所示：

        使用日期时间函数计算平均太阳黑子周期。 

dttmCycle = years(mean(diff(lk)))
 
dttmCycle = 10.9600

        用这些数据创建一个时间表。以年为单位指定时间变量。绘制数据。显示时间表的最后五个条目。

TT = timetable(years(year),avSpots);
plot(TT.Time,TT.Variables)

        如图所示：

entries = TT(end-4:end,:)
 
 
entries=5×1 timetable
      Time      avSpots
    ________    _______
 
    1983 yrs     66.6  
    1984 yrs     45.9  
    1985 yrs     17.9  
    1986 yrs     13.4  
    1987 yrs     29.3  

 约束峰值特征       

        加载以7418 Hz采样的音频信号。选择200个样本。

load mtlb
select = mtlb(1001:1200);

        找出间隔至少5 ms的峰。 要应用此约束，findpeaks会选择信号中最高的峰值，并消除其5毫秒内的所有峰值。然后，该函数对最高的剩余峰重复该过程，并进行迭代，直到用完要考虑的峰。

findpeaks(select,Fs,'MinPeakDistance',0.005)

        如图所示：

        找出幅度至少为1 V的峰值。 

findpeaks(select,Fs,'MinPeakHeight',1)

        如图所示：

        找出比相邻样本至少高1 V的峰值。 

findpeaks(select,Fs,'Threshold',1)

        如图所示：

        在信号达到更高值之前，找出任一侧下降至少1 V的峰值。 

findpeaks(select,Fs,'MinPeakProminence',1)

        如图所示：

饱和信号的峰值

        如果数据大于给定的饱和点，传感器可以返回削波读数。可以选择忽略这些没有意义的峰值，或者将它们纳入您的分析中。 产生一个信号，该信号由嵌入方差为0.1的高斯白噪声中的频率为5 Hz和3 Hz的三角函数的乘积组成。信号以100 Hz的速率采样一秒钟。重置随机数生成器以获得可重现的结果。

rng default
 
fs = 1e2;
t = 0:1/fs:1-1/fs;
 
s = sin(2*pi*5*t).*sin(2*pi*3*t)+randn(size(t))/10;

         通过截断每个大于指定界限0.32的读数来模拟饱和测量。绘制饱和信号。

bnd = 0.32;
s(s>bnd) = bnd;
 
plot(t,s)
xlabel('Time (s)')

        如图所示：

        找到信号的峰值。findpeaks只报告每个平峰的上升沿。 

[pk,lc] = findpeaks(s,t);
 
hold on
plot(lc,pk,'x')

        如图所示：

        使用“阈值”名称-值对来排除平峰。要求最小幅度差为10^-4 一座山峰和它的邻居之间。 

[pkt,lct] = findpeaks(s,t,'Threshold',1e-4);
 
plot(lct,pkt,'o','MarkerSize',12)

        如图所示：

确定峰宽

        创建一个由一系列钟形曲线组成的信号。指定每条曲线的位置、高度和宽度。

x = linspace(0,1,1000);
 
Pos = [1 2 3 5 7 8]/10;
Hgt = [4 4 2 2 2 3];
Wdt = [3 8 4 3 4 6]/100;
 
for n = 1:length(Pos)
    Gauss(n,:) =  Hgt(n)*exp(-((x - Pos(n))/Wdt(n)).^2);
end
 
PeakSig = sum(Gauss);

        画出各条曲线及其总和。 

plot(x,Gauss,'--',x,PeakSig)
grid

        如图所示：

        使用半日珥作为参考测量峰的宽度。 

findpeaks(PeakSig,x,'Annotate','extents')

        如图所示：

        再次测量宽度，这次使用一半高度作为参考。

findpeaks(PeakSig,x,'Annotate','extents','WidthReference','halfheight')
title('Signal Peak Widths')

        如图所示：