Пиковый анализ

Этот пример показывает, как выполнить базовый пиковый анализ. Это поможет вам ответить на такие вопросы, как: Как я могу найти peaks в своем сигнале? Как измерить расстояние между 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')

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

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

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

Peaks в сигнале появляются, по-видимому, с регулярными интервалами. Однако некоторые из peaks очень близки друг к другу. The MinPeakProminence свойство может использоваться, чтобы отсеять эти пики. Рассмотрим peaks, которые падают с обеих сторон как минимум на 40 относительных чисел солнечных пятен, прежде чем встретить большее значение.

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

Figure contains an axes. The axes 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. The axes 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 лет между 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])

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

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

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

Этот пример показывает пиковый анализ в сигнале ЭКГ (электро-кардиограмма). ЭКГ является мерой электрической активности сердца с течением времени. Сигнал измеряется электродами, прикрепленными к коже, и чувствителен к нарушениям порядка, как интерференция источника степени и шум от программных продуктов движения.

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

После удаления тренда, найдите комплекс QRS, который является самым заметным повторяющимся пиком в сигнале ECG. Комплекс QRS соответствует деполяризации правого и левого желудочков человеческого сердца. Его можно использовать, чтобы определить частоту сердечных сокращений пациента или предсказать нарушения в работе сердца. Следующий рисунок показывает форму комплекса QRS в сигнале ECG.

Пороговое значение для поиска интересующих вас пиков

Комплекс QRS состоит из трех основных компонентов: Q-волна, R-волна, S-волна. Волны R могут быть обнаружены с помощью пороговых достигать максимума выше 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')

Figure contains an axes. The axes 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-волн. Установка порогового peaks для определения местоположения Q-волн приводит к обнаружению нежелательного peaks, когда Q-волны зарыты в шум. Сначала фильтруем сигнал, а затем находим peaks. Фильтрация Савицкого-Голея используется для удаления шума в сигнале.

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

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

См. также

Для просмотра документации необходимо авторизоваться на сайте