峰值分析

此示例说明如何执行基本峰值分析。它将帮助您回答诸如以下的问题：我如何在信号中找到峰值？我如何测量峰值之间的距离？如何测量受趋势影响的信号峰值的振幅？我如何在含噪信号中找到峰值？我如何找到局部最小值？

查找最大值或峰值

苏黎世太阳黑子相对数测量太阳黑子的数量和大小。使用 findpeaks 函数查找峰的位置和值。

load sunspot.dat
year = sunspot(:,1); 
relNums = sunspot(:,2);
findpeaks(relNums,year)
xlabel('Year')
ylabel('Sunspot Number')
title('Find All Peaks')

上图以表格形式显示 300 年来太阳黑子的数量，并标出检测到的峰值。下一节说明如何测量这些峰值之间的距离。

测量峰值之间的距离

信号的峰值似乎以固定间隔出现。然而，有些峰彼此非常接近。MinPeakProminence 属性可用于滤除这些峰。此示例只考虑两侧相对高差至少为 40 个黑子数的峰。

findpeaks(relNums,year,'MinPeakProminence',40)
xlabel('Year')
ylabel('Sunspot Number')
title('Find Prominent Peaks')

以下直方图显示以年数表示的峰值间隔分布：

figure
[pks, locs] = findpeaks(relNums,year,'MinPeakProminence',40);
peakInterval = diff(locs);
histogram(peakInterval)
grid on
xlabel('Year Intervals')
ylabel('Frequency of Occurrence')
title('Histogram of Peak Intervals (years)')

AverageDistance_Peaks = mean(diff(locs))

AverageDistance_Peaks = 10.9600

该分布显示，大多数峰值间隔在 10 至 12 年之间，表明信号具有周期性。此外，峰值之间的平均间隔为 10.96 年，这与已知的 11 年太阳黑子活动周期相匹配。

寻找裁剪或饱和信号中的峰值

您可能希望将平峰视为峰值或排除它们。如需排除平峰，请使用 threshold 属性指定最小偏移，该偏移定义为峰值与其紧邻峰值之间的振幅差。

load clippedpeaks.mat

figure

% Show all peaks in the first plot
ax(1) = subplot(2,1,1);
findpeaks(saturatedData)
xlabel('Samples')
ylabel('Amplitude')
title('Detecting Saturated Peaks')

% Specify a minimum excursion in the second plot
ax(2) = subplot(2,1,2);
findpeaks(saturatedData,'threshold',5)
xlabel('Samples')
ylabel('Amplitude')
title('Filtering Out Saturated Peaks')

% link and zoom in to show the changes
linkaxes(ax(1:2),'xy')
axis(ax,[50 70 0 250])

第一个子图显示，对于平峰，上升沿检测为峰值。第二个子图显示，指定阈值有助于排除平峰。

测量峰的振幅

此示例显示 ECG（心电图）信号的峰值分析。ECG 是对随时间变化的心跳活动电信号的测量。该信号由附着在皮肤上的电极测量，对电源干扰和由运动伪影引起的噪声等干扰十分敏感。

load noisyecg.mat
t = 1:length(noisyECG_withTrend);

figure
plot(t,noisyECG_withTrend)
title('Signal with a Trend')
xlabel('Samples');
ylabel('Voltage(mV)')
legend('Noisy ECG Signal')
grid on

数据去趋势

上述信号显示一个基线偏移，因此不表示真实振幅。为了去趋势，对信号进行某低次多项式拟合，并使用该多项式对信号去趋势。

[p,s,mu] = polyfit((1:numel(noisyECG_withTrend))',noisyECG_withTrend,6);
f_y = polyval(p,(1:numel(noisyECG_withTrend))',[],mu);

ECG_data = noisyECG_withTrend - f_y;        % Detrend data

figure
plot(t,ECG_data)
grid on
ax = axis;
axis([ax(1:2) -1.2 1.2])
title('Detrended ECG Signal')
xlabel('Samples')
ylabel('Voltage(mV)')
legend('Detrended ECG Signal')

在去趋势后，找到 QRS 复波，它是 ECG 信号中最突出的重复波峰。QRS 复波对应于人心脏左右心室的去极化。它可用于确定患者的心率或预测心脏功能异常。下图显示 ECG 信号中 QRS 复波的形状。

阈值化以找到感兴趣的峰值

QRS 复波由三个主要分量组成：Q 波、R 波、S 波。R 波可以通过设置 0.5 毫伏以上的峰值阈值来检测。注意 R 波相隔 200 多个采样。使用此信息通过指定 'MinPeakDistance' 来去除不需要的峰值。

[~,locs_Rwave] = findpeaks(ECG_data,'MinPeakHeight',0.5,...
                                    'MinPeakDistance',200);

为了检测 S 波，找到信号中的局部最小值，并适当应用阈值。

查找信号中的局部最小值

局部最小值可以通过在原始信号的反转版本上找到峰值来检测。

ECG_inverted = -ECG_data;
[~,locs_Swave] = findpeaks(ECG_inverted,'MinPeakHeight',0.5,...
                                        'MinPeakDistance',200);

下图显示在信号中检测到的 R 波和 S 波。

figure
hold on 
plot(t,ECG_data)
plot(locs_Rwave,ECG_data(locs_Rwave),'rv','MarkerFaceColor','r')
plot(locs_Swave,ECG_data(locs_Swave),'rs','MarkerFaceColor','b')
axis([0 1850 -1.1 1.1])
grid on
legend('ECG Signal','R waves','S waves')
xlabel('Samples')
ylabel('Voltage(mV)')
title('R wave and S wave in Noisy ECG Signal')

接下来，我们尝试确定 Q 波的位置。对峰值设置阈值以定位 Q 波会导致检测到不需要的峰值，因为 Q 波在噪声中被掩盖。我们先对信号进行滤波，然后找到峰值。萨维茨基-戈雷滤波用于去除信号中的噪声。

smoothECG = sgolayfilt(ECG_data,7,21);

figure
plot(t,ECG_data,'b',t,smoothECG,'r')
grid on
axis tight
xlabel('Samples')
ylabel('Voltage(mV)')
legend('Noisy ECG Signal','Filtered Signal')
title('Filtering Noisy ECG Signal')

我们对平滑信号执行峰值检测，并使用逻辑索引来找到 Q 波的位置。

[~,min_locs] = findpeaks(-smoothECG,'MinPeakDistance',40);

% Peaks between -0.2mV and -0.5mV
locs_Qwave = min_locs(smoothECG(min_locs)>-0.5 & smoothECG(min_locs)<-0.2);

figure
hold on
plot(t,smoothECG); 
plot(locs_Qwave,smoothECG(locs_Qwave),'rs','MarkerFaceColor','g')
plot(locs_Rwave,smoothECG(locs_Rwave),'rv','MarkerFaceColor','r')
plot(locs_Swave,smoothECG(locs_Swave),'rs','MarkerFaceColor','b')
grid on
title('Thresholding Peaks in Signal')
xlabel('Samples')
ylabel('Voltage(mV)')
ax = axis;
axis([0 1850 -1.1 1.1])
legend('Smooth ECG signal','Q wave','R wave','S wave')

上图显示，在含噪 ECG 信号中成功检测到 QRS 复波。

含噪信号和平滑信号之间的误差

请注意原始信号和去趋势滤波后的信号中的 QRS 复波之间的平均差异。

% Values of the Extrema
[val_Qwave, val_Rwave, val_Swave] = deal(smoothECG(locs_Qwave), smoothECG(locs_Rwave), smoothECG(locs_Swave));

meanError_Qwave = mean((noisyECG_withTrend(locs_Qwave) - val_Qwave))

meanError_Qwave = 0.2771

meanError_Rwave = mean((noisyECG_withTrend(locs_Rwave) - val_Rwave))

meanError_Rwave = 0.3476

meanError_Swave = mean((noisyECG_withTrend(locs_Swave) - val_Swave))

meanError_Swave = 0.1844

这表明，为了进行高效的峰值分析，必须对含噪信号进行去趋势。

峰值属性

一些重要的峰值属性包括上升时间、下降时间、上升水平和下降水平。为 ECG 信号中的每个 QRS 复波计算这些属性。下图显示这些属性的平均值。

avg_riseTime = mean(locs_Rwave-locs_Qwave); % Average Rise time
avg_fallTime = mean(locs_Swave-locs_Rwave); % Average Fall time
avg_riseLevel = mean(val_Rwave-val_Qwave);  % Average Rise Level
avg_fallLevel = mean(val_Rwave-val_Swave);  % Average Fall Level

helperPeakAnalysisPlot(t,smoothECG,...
                    locs_Qwave,locs_Rwave,locs_Swave,...
                    val_Qwave,val_Rwave,val_Swave,...
                    avg_riseTime,avg_fallTime,...
                    avg_riseLevel,avg_fallLevel)

另请参阅

findpeaks