對(duì)信號(hào)的峰值進(jìn)行分析
Step1.
簡(jiǎn)單的尋找最大值點(diǎn)

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

圖片.png

Step2.
計(jì)算峰值點(diǎn)的距離

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

圖片.png

先用直方圖分析

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

AverageDistance_Peaks = mean(diff(locs))

圖片.png

Step3.
查看截止或者飽和信號(hào)的峰值點(diǎn)

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

圖片.png

可以通過增加門限進(jìn)行設(shè)定，否則對(duì)于一個(gè)平坦的峰，上升沿的第一個(gè)點(diǎn)會(huì)被認(rèn)為是峰值

step4.
測(cè)量峰值幅度
用ECG信號(hào)舉例說明

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

圖片.png

step5.
去除信號(hào)趨勢(shì)

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

圖片.png

下圖是一個(gè)標(biāo)準(zhǔn)的ECG信號(hào)圖

圖片.png

Step6.
感興趣的峰值點(diǎn)門限，ECG主要的三個(gè)部分，包括S波，Q波，R波
R波可以通過設(shè)定0.5mv的門限檢出

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

對(duì)于S波，要設(shè)置適當(dāng)門限

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

然后標(biāo)記出來如下

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

圖片.png

為了定位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')

圖片.png

查找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');

圖片.png

查看噪聲信號(hào)和平滑后哦信號(hào)的誤差

% 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_Rwave = mean((noisyECG_withTrend(locs_Rwave) - val_Rwave))
meanError_Swave = mean((noisyECG_withTrend(locs_Swave) - val_Swave))

圖片.png

峰值分析
主要分析上升時(shí)間，下降時(shí)間等

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)

圖片.png

總結(jié)：

1. 學(xué)會(huì)幾個(gè)函數(shù)的使用
   findpeaks， 查找峰值點(diǎn)
   linkaxes， 關(guān)聯(lián)子圖的坐標(biāo)軸
   polyfit， 擬合
   polyval，擬合后的應(yīng)用
   sgolayfilt，濾波器
   deal， 輸入到輸出的賦值

helperPeakAnalysisPlot，peak峰值點(diǎn)繪圖，內(nèi)部函數(shù)看不到介紹

2. 查找峰值點(diǎn)的一般方法
   如何設(shè)置參數(shù)范圍大小
   針對(duì)ECG如何獲取Q,R,S波峰位置