Finding Maxima or Peaks

The Zurich sunspot relative number measures both the number and size of sunspots. Use theFindpeaksfunction to find the locations and the value of the peaks.

loadsunspot.dat年= sunspot（：，1）;relnums = sunspot（：，2）;findpeaks（relnums，年）xlabel（“年”）ylabel（'Sunspot Number') title('Find All Peaks')

The above plot shows sunspot numbers tabulated over 300 years and labels the detected peaks. The next section shows how to measure distance between these peaks.

Measuring Distance Between Peaks

Peaks in the signal seem to appear at regular intervals. However, some of the peaks are very close to each other. TheMinPeakProminenceproperty can be used filter out these peaks. Consider peaks that drop off on both sides by at least 40 relative sunspot numbers before encountering a larger value.

Findpeaks(relNums,year,'MinPeakProminence',40) xlabel(“年”）ylabel（'Sunspot Number') title('Find Prominent Peaks')

The following histogram shows the distribution of peak intervals in years:

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

AverageDistance_Peaks = mean(diff(locs))

AverageDistance_Peaks = 10.9600

The distribution shows that majority of peak intervals lie between 10 and 12 years indicating the signal has a cyclic nature. Also, the average interval of 10.96 years between the peaks matches the known cyclic sunspot activity of 11 years.

Finding Peaks in Clipped or Saturated Signals

You may want to consider flat peaks as peaks or exclude them. In the latter case, a minimum excursion which is defined as the amplitude difference between a peak and its immediate neighbors is specified using thethresholdproperty.

loadclippedpeaks.matfigure% Show all peaks in the first plotax(1) = subplot(2,1,1); findpeaks(saturatedData) xlabel('Samples'）ylabel（'Amplitude') title('Detecting Saturated Peaks')％在第二个图中指定最小游览ax(2) = subplot(2,1,2); findpeaks(saturatedData,'临界点',5) xlabel('Samples'）ylabel（'Amplitude') title('Filtering Out Saturated Peaks')% link and zoom in to show the changeslinkaxes(ax(1:2),'xy'）轴（AX，[50 70 0 250]）

The first subplot shows, that in case of a flat peak, the rising edge is detected as the peak. The second subplot shows that specifying a threshold can help to reject flat peaks.

Measuring Amplitudes of Peaks

This example shows peak analysis in an ECG (Electro-cardiogram) signal. ECG is a measure of electrical activity of the heart over time. The signal is measured by electrodes attached to the skin and is sensitive to disturbances such as power source interference and noises due to movement artifacts.

loadnoisyecg.matt = 1:length(noisyECG_withTrend); figure plot(t,noisyECG_withTrend) title(“信号有趋势”）xlabel（'Samples'); ylabel(“电压（MV）”) legend('Noisy ECG Signal') grid上

Detrending Data

上面的信号显示基线移位，因此不代表真实的幅度。为了消除趋势，将低阶多项式拟合到信号中，并使用多项式来降低它。

[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 datafigure plot(t,ECG_data) grid上ax = axis; axis([ax(1:2) -1.2 1.2]) title('Detrended ECG Signal'）xlabel（'Samples'）ylabel（“电压（MV）”) legend('Detrended ECG Signal')

After detrending, find the QRS complex, which is the most prominent repeating peak in the ECG signal. The QRS complex corresponds to the depolarization of the right and left ventricles of the human heart. It can be used to determine a patient's cardiac rate or predict abnormalities in heart function. The following figure shows the shape of the QRS complex in an ECG signal.

Thresholding to Find Peaks of Interest

QRS综合体由三个主要组成部分组成：Q wave, R wave, S wave. The R waves can be detected by thresholding peaks above 0.5 mV. Notice that the R waves are separated by more than 200 samples. Use this information to remove unwanted peaks by specifying a 'MinPeakDistance'.

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

For detection of the S-waves, find the local minima in the signal and apply thresholds appropriately.

在信号中找到本地最小值

局部最小值可以被发现豌豆ks on an inverted version of the original signal.

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

The following plot shows the R waves and S waves detected in the signal.

图上plot(t,ECG_data) plot(locs_Rwave,ECG_data(locs_Rwave),'rv','MarkerFaceColor','r') plot(locs_Swave,ECG_data(locs_Swave),'rs','MarkerFaceColor','b'）轴（[0 1850 -1.1 1.1]）网格上legend(“心电图信号”,'R waves',波浪'）xlabel（'Samples'）ylabel（“电压（MV）”) title('R wave and S wave in Noisy ECG Signal')

Next, we try and determine the locations of the Q waves. Thresholding the peaks to locate the Q waves results in detection of unwanted peaks as the Q waves are buried in noise. We filter the signal first and then find the peaks. Savitzky-Golay filtering is used to remove noise in the signal.

smoothECG = sgolayfilt(ECG_data,7,21); figure plot(t,ECG_data,'b',t,smoothECG,'r') grid上axistightxlabel('Samples'）ylabel（“电压（MV）”) legend('Noisy ECG Signal',“过滤信号”) title('Filtering Noisy ECG Signal')

我们对平稳信号和执行峰值检测use logical indexing to find the locations of the Q waves.

[~,min_locs] = findpeaks(-smoothECG,'MinPeakDistance',40);％峰值 -  0.2mV和-0.5MV之间的峰值locs_Qwave = min_locs(smoothECG(min_locs)>-0.5 & smoothECG(min_locs)<-0.2); figure hold上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上title(“信号中的阈值峰”）xlabel（'Samples'）ylabel（“电压（MV）”) ax = axis; axis([0 1850 -1.1 1.1]) legend('Smooth ECG signal','Q wave','R wave','S wave')

The above figure shows that the QRS complex successfully detected in the noisy ECG signal.

Error Between Noisy and Smooth Signal

请注意，在原始信号和脱落过滤的信号中，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 =平均值（（noisyecg_withtrend（locs_rwave）-val_rwave））

MeaneRror_rwave = 0.3476

MeaneRror_swave =平均值（（noisyecg_withtrend（locs_swave）-val_swave））

meanError_Swave = 0.1844

This demonstrates that it is essential to detrend a noisy signal for efficient peak analysis.

Peak Properties

Some important peak properties involve rise time, fall time, rise level, and fall level. These properties are computed for each of the QRS complexes in the ECG signal. The average values for these properties are displayed on the figure below.

avg_riseTime = mean(locs_Rwave-locs_Qwave);% Average Rise timeavg_fallTime = mean(locs_Swave-locs_Rwave);% Average Fall timeavg_riseLevel = mean(val_Rwave-val_Qwave);% Average Rise Levelavg_fallLevel = mean(val_Rwave-val_Swave);平均秋季水平％helperPeakAnalysisPlot(t,smoothECG,...locs_Qwave,locs_Rwave,locs_Swave,...val_Qwave,val_Rwave,val_Swave,...avg_riseTime,avg_fallTime,...avg_riselevel，avg_falllevel）