В этом примере показано, как выполнить основной пиковый анализ. Это поможет вам ответить на вопросы, такие как: Как я нахожу peaks в своем сигнале? Как я измеряю расстояние между peaks? Как я измеряю амплитуду peaks сигнала, который затронут трендом? Как я нахожу 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 очень друг близко к другу. 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); hist(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 или исключить их. В последнем случае минимальное отклонение, которое задано как амплитудное различие между пиком и его ближайшими соседями, задано с помощью 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.
Этот пример показывает пиковый анализ в 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.
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)