Main Content

峰值分析

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

查找最大值或峰值

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

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

Figure contains an axes object. The axes object with title Find All Peaks contains 2 objects of type line.

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

测量峰值之间的距离

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

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

Figure contains an axes object. The axes object with title Find Prominent Peaks contains 2 objects of type line.

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

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

Figure contains an axes object. The axes object with title Histogram of Peak Intervals (years) contains an object of type histogram.

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

Figure contains 2 axes objects. Axes object 1 with title Detecting Saturated Peaks contains 2 objects of type line. Axes object 2 with title Filtering Out Saturated Peaks contains 2 objects of type line.

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

测量峰的幅值

此示例显示 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

Figure contains an axes object. The axes object with title Signal with a Trend contains an object of type line. This object represents Noisy ECG Signal.

数据去趋势

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

[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')

Figure contains an axes object. The axes object with title Detrended ECG Signal contains an object of type line. This object represents 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')

Figure contains an axes object. The axes object with title R wave and S wave in Noisy ECG Signal contains 3 objects of type line. These objects represent ECG Signal, R waves, S waves.

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

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

Figure contains an axes object. The axes object with title Filtering Noisy ECG Signal contains 2 objects of type line. These objects represent Noisy ECG Signal, Filtered 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')

Figure contains an axes object. The axes object with title Thresholding Peaks in Signal contains 4 objects of type line. These objects represent 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)

Figure contains an axes object. The axes object with title QRS-complex in an ECG signal contains 11 objects of type line, text. These objects represent QRS-Complex, Peak, Minima.

另请参阅