Нахождение максимумов или Peaks

Цюрихский номер родственника солнечного пятна измеряет и номер и размер солнечных пятен. Используйте findpeaks функционируйте, чтобы найти местоположения и значение peaks.

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

Вышеупомянутый график показывает, что числа солнечного пятна свели в таблицу более чем 300 лет и помечают обнаруженный peaks. Следующий раздел показывает, как измерить расстояние между этим peaks.

Измерение расстояния между Peaks

Peaks в сигнале, кажется, появляется равномерно. Однако часть peaks очень друг близко к другу. MinPeakProminence свойство может использоваться, отфильтровывают этот peaks. Рассмотрите peaks, который понижается с обеих сторон по крайней мере 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 лет между peaks совпадает с известным циклическим действием солнечного пятна 11 лет.

Нахождение Peaks в отсеченных или влажных сигналах

Можно хотеть рассмотреть плоский peaks как peaks или исключить их. В последнем случае минимальное отклонение, которое задано как амплитудное различие между пиком и его ближайшими соседями, задано с помощью 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])

Первый подграфик показывает, что в случае плоского пика, возрастающее ребро обнаруживается как пик. Второй подграфик показывает, что определение порога может помочь отклонить плоский peaks.

Измерение амплитуд Peaks

Этот пример показывает пиковый анализ в 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 соответствует деполяризации правых и левых желудочков человеческого сердца. Это может использоваться, чтобы определить сердечный уровень пациента или предсказать отклонения в сердечной функции. Следующий рисунок показывает форму комплекса QRS в сигнале ECG.

Пороговая обработка, чтобы найти Peaks интересным

Комплекс QRS состоит из трех главных компонентов: Q волна, R волна, S волна. Волны R могут быть обнаружены peaks пороговой обработки выше 0,5 мВ. Заметьте, что волны R разделяются больше чем 200 выборками. Используйте эту информацию, чтобы удалить нежелательный peaks путем определения 'MinPeakDistance'.

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

Для обнаружения S-волн найдите локальные минимумы в сигнале и примените пороги соответственно.

Нахождение локальных минимумов в сигнале

Локальные минимумы могут быть обнаружены путем нахождения peaks на инвертированной версии исходного сигнала.

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. Пороговая обработка peaks, чтобы определить местоположение результатов волн Q в обнаружении нежелательного peaks как волны Q содержится шум. Мы фильтруем сигнал сначала и затем находим peaks. Фильтрация 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')

Мы выполняем пиковое обнаружение на сглаженном сигнале и используем логическую индексацию, чтобы найти местоположения волн 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')

Вышеупомянутый рисунок показывает, что QRS объединяют успешно обнаруженный в шумном сигнале ECG.

Ошибка между шумным и сглаженным сигналом

Заметьте среднюю разницу между комплексом 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

Это демонстрирует, что важно детрендировать сигнал с шумом для эффективного пикового анализа.

Пиковые свойства

Некоторые важные пиковые свойства включают время нарастания, время спада, уровень повышения и осенний уровень. Эти свойства вычисляются для каждого из комплексов QRS в сигнале ECG. Средние значения для этих свойств отображены на рисунке ниже.

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)