Matlab2019b信号峰值检测与提取

在信号处理的时候，有的时候希望确定信号峰值的大小以及位置，并最好对信号进行分离，可以更好地进行后续处理。本人使用了集中不同的思路和方式来进行位置确定和提取，仅供参考。

注意： 部分程序仅兼容Matlab2016b以后版本

目标

提取和检测如下图类似信号的峰值及位置

方法1

调用的findpeaks函数
参考程序，通过设定合适的峰值参数，可以输出相应的峰值

%% 找到数据中的峰值
clc;
clear;
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
Data = sum(Gauss);%产生假设信号

% NPeaks表示返回值中最大的峰值个数
[pks,locs,w,p] = findpeaks(Data,x,...pks为输出峰值，locs为峰值位置
    'NPeaks',1,...输出的最大峰值个数
    'SortStr','ascend',...输出峰值是否进行排序，升序或降序或不进行排序
    'Threshold',0,...峰值之间的最小差值阈值，如数组为[]则未找到符合参数的阈值
    'WidthReference','halfprom',...
    ...'MiniPeakWidth','0',...最小的峰值宽度
    'MaxPeakWidth',Inf,...最大的峰值宽度
    'Annotate','peaks',...在定义输出后无效
    'MinPeakHeight',50);%最小峰值高度

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

方法2

通过求取整个信号的方差和平均值，再使用逻辑与或的方式提取位置和峰值。由于方差即为随机变量和数学期望之间的偏离程度，当极度偏离方差时，信号为突变信号。首先提取位置，通过与方差比较，大于阈值的即为1，小于阈值的即为0，在于信号进行与逻辑运算，即可提取出相对准确的峰值。

clc
clear
Fs=1000;
y=xlsread('Data2','sheet1');
x=y(:,7)';%通道设置
figure('Name','原始混合信号','NumberTitle','off');
plot(x)

jun=mean(x);%求取整个信号的均值
biao=std(x);%求取整个信号的标准差
fang=var(x);%求整个信号的方差
zhong=median(x);%求整个信号的中位数
fai=1;

% x = 0:0.05:50*pi;
% x= signal+rand(1,length(signal))
L=length(x);
T=0.5*(min(x(:))+max(x(:)));

%循环判断与整个信号方差之间的关系
for i=1:L
        if(x(i)<fai*fang)%具体需要根据实际使用情况来调整系数fai
            BW1(i)=0;%未超过系数为0
        else 
            BW1(i)=1;%超过系数为1
    end
end

for i=1:L
       a(i)=x(i)*BW1(i);%与逻辑提取信号
end

for i=1:L
       b(i)=-(x(i)*(BW1(i)-1));%与逻辑
end

figure('Name','冲击信号位置','NumberTitle','off');
plot(BW1)


figure('Name','信号分离','NumberTitle','off');
subplot(2,1,1)
plot(b)
title('应变信号');
subplot(2,1,2)
plot(a)
title('冲击信号');

方法3

使用VPD方法

%一维峰值检测Matlab实现
clc;
clear

signal = 0:0.05:50*pi;
x = sin(signal);
row_acc=x;

plot(x)
% x = 0:0.05:50*pi;
% x= signal+rand(1,length(signal))
L=length(x);

% Y = fft(x);
% P2 = abs(Y/L);
% P1 = P2(1:L/2+1);% 实信号的功率谱是对称的,只用取一半即可
% P1(2:end-1) = 2*P1(2:end-1);
% f = Fs*(0:(L/2))/L;
% 
% plot(f,P1)
% title('原始信号直接进行傅里叶分析')
% grid on
% xlabel('Frequency (Hz)')
% ylabel('Amplitude')

% [c,l] = wavedec(x,2,'db2');
% approx = appcoef(c,l,'db2');
% [cd1,cd2] = detcoef(c,l,[1 2]);
% 
% subplot(3,1,1)
% plot(x)
% title('原始信号')
% subplot(3,1,2)
% plot(approx)
% title('应变信号')
% subplot(3,1,3)
% plot(cd3)
% title('冲击信号')

% subplot(4,1,3)
% plot(cd2)
% title('Level 2 Detail Coefficients')
% subplot(4,1,4)
% plot(cd1)
% title('Level 1 Detail Coefficients')
% 
% Hd=LowPass_Filter;
% output=filter(Hd,x);
% figure
% plot(output)
% 
% blo = fir1(34,0.48,chebwin(35,30));
% outlo = filter(blo,1,x);
% 
% 使用峰值检测
% % [pks,locs]=findpeaks(x,'Npeaks',1,'MinPeakHeight',250);
% % jieyue=x(:,[locs-200:L]);
% % figure
% % plot(jieyue)
% % yingbian=x(:,[1:locs-200]);
% % figure
% % plot(yingbian)
% % ylim([-50 300])
% % figure
% % plot(locs,pks)
% 
%% VPD方法

%%前三点均值滤波
row_acc = x;
m = length(row_acc);
row_acc1 = linspace(0,0,m);
row_acc1(1) = row_acc(1);
row_acc1(m) = row_acc(m);
for i=2:m-1
   row_acc1(i)=(row_acc(i-1) + row_acc(i)+row_acc(i+1))/3;
end
% figure;
% plot(row_acc1);

for i=m-1:-1:2
    row_acc(i) = (row_acc1(i-1) + row_acc1(i)+row_acc1(i+1))/3;
end

%%找到局部最小值和局部最大值及其对应的位置，波峰点、波谷点满足：
peaks = linspace(0,0,m);
valleys = linspace(0,0,m);
peakindexs = linspace(0,0,m);
valleyindexs = linspace(0,0,m);
peakindex = 1;
valleyindex = 1;

for i = 2:m-1
    if row_acc(i) >row_acc(i-1) && row_acc(i)>=row_acc(i+1)
        peaks(peakindex)=row_acc(i);
        peakindexs(peakindex)=i;
        peakindex = peakindex+1;
    end
    if row_acc(i) < row_acc(i-1) && row_acc(i)<row_acc(i+1)
        valleys(valleyindex)=row_acc(i);
        valleyindexs(valleyindex)=i;
        valleyindex=valleyindex+1;
    end
end
% 
%%计算VPD，VPD(n)表示第n个波谷点的值与第n个波峰点的值的差，VPD用来去掉那些假的波峰点
pcount = peakindex-1;
vcount = valleyindex-1;
peakindices = linspace(0,0,pcount);
for i = 1:pcount
    peakindices(i) = peakindexs(i);
end
valleyindices = linspace(0,0,vcount);
for i = 1:vcount
    valleyindices(i) = valleyindexs(i);
end
% figure;
% plot(row_acc,'-o', 'MarkerIndices',peakindices,'MarkerFaceColor','red','MarkerSize',5);
% 
% figure;
% plot(row_acc,'-s', 'MarkerIndices',valleyindices,'MarkerFaceColor','green','MarkerSize',5);
if pcount>2 && vcount>2
    if peakindexs(1) < valleyindexs(1)
        peakindex=2;
    else
        peakindex=1;
    end
    vindex=1;
end

if peakindex == 2
    for i = 1:m-1
        peaks(i)=peaks(i+1);
    end
    pcount = pcount-1;
    pindex=1;
end

vpd = linspace(0,0,m);
vpd1 = linspace(0,0,m);
for i=1:pcount
    vpd(i) = peaks(i) - valleys(i);
end

dels = linspace(0,0, pcount);
peakindexs1 = linspace(0,0,pcount);
if pcount > 2
    lastcount=pcount;
    curcount = 1;
    while lastcount ~= curcount
        lastcount = curcount;
        del_count = 0;
        for i = 2:pcount-1
            if vpd(i) <= 0.7 * (vpd(i-1) + vpd(i)+vpd(i+1)) / 3
                dels(i)=1;
            end
        end

        count = 1;
        for i = 1:pcount
          if dels(i) ~= 1
              vpd1(count) = vpd(i); 
              peakindexs1(count) = peakindexs(i);
              count = count+1;
          else
               del_count = del_count + 1;
               dels(i) = 0;
           end
        end
        pcount = pcount - del_count;
        for i = 1:pcount
            vpd(i) = vpd1(i); 
            peakindexs(i) = peakindexs1(i);
        end
        peakindexs(pcount+1) = 0;
        vpd(pcount+1) = 0;

        indices = linspace(0,0,pcount);
        for i = 1:pcount
            indices(i) = peakindexs1(i);
        end
        plot(row_acc,'-o', 'MarkerIndices',indices,'MarkerFaceColor','red','MarkerSize',10);
        curcount = pcount;
    end  
end

% [pks,locs]=findpeaks(row_acc,'Npeaks',1,'MinPeakHeight',20);
% jieyue=row_acc(:,[locs-200:L]);%从触发峰值点左边200处分离
% figure
% plot(jieyue)
% yingbian=x(:,[1:locs-200]);
% figure
% plot(yingbian)
% ylim([-50 300])
% 
% % figure
% % plot(locs,pks)
% 

Copyright © 2020 by RichardYang. All rights reserved.
仅供参考，严禁转载，感谢。