Блог

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

Например, предположим, что пик имеет такую форму:

 t=linspace(0,10,1e3);
w=0.25;
rw=@(t0) 2/(sqrt(3*w)*pi^0.25)*(1-((t-t0)/w).^2).*exp(-(t-t0).^2/(2*w.^2));
  

и у меня есть несколько пиков, расположенных слишком близко, поэтому их структура начинает мешать другим, и один отделен, так что вы можете видеть его структуру:

 pos=[ 2 2.5 3 8]; % positions
total_signal=0.*t;
for i=1:length(pos) 
   total_signal=total_signal rw(pos(i));
end
plot(t,total_signal);
  

Итак, цель состоит в том, чтобы найти все пики, найденные в total_signal , и проверить, совпадают ли их позиции с исходными позициями, которые использовались для их генерации в pos .

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

Я начну с переписывания вашего кода следующим образом:

 t = linspace(0,10,1e3);
w = 0.25;
rw = @(t) 2/(sqrt(3*w)*pi^0.25)*(1-((t)/w).^2).*exp(-(t).^2/(2*w.^2));

pos = [2, 2.5, 3, 8];
total_signal = zeros(size(t));
for i=1:length(pos) 
   total_signal = total_signal   rw(t-pos(i));
end
total_signal = total_signal   randn(size(t))*1e-2;
plot(t,total_signal);
  

Я изменил rw на take t-t0 , а не просто t0 . Это позволит мне позже создать чистый сигнал только с одним пиком в середине. Я также добавил шум к сигналу для более реалистичной задачи. Без шума проблему решить намного проще.

Деконволюция Винера — самый простой подход к решению этой проблемы. Короче говоря, мы предполагаем, что

 total_signal = conv(pulse_train, shape)
  

В частотной области это записывается как

 G = F .* H
  

(с .* поэлементным умножением, используя синтаксис MATLAB здесь). Деконволюция Винера — это:

 F = (conj(H) .* G) ./ (abs(H).^2   k)
  

с k некоторой константой, которую мы можем настроить для упорядочивания решения.

Мы реализуем это следующим образом:

 shape = rw(linspace(-5,5,1e3));
G = fft(total_signal);
H = fft(ifftshift(shape)); % ifftshift moves the origin to sample #0, as expected by FFT.
k = 1;
F = (conj(H) .* G) ./ (abs(H).^2   k);
pulse_train = ifft(F);
  

Теперь findpeaks (требуется набор инструментов обработки сигналов) может использоваться для поиска заметных пиков:

 findpeaks(pulse_train, 'MinPeakProminence', 0.02)
  

Обратите внимание, что четыре пика имеют примерно одинаковую высоту. Это не точно, потому что мы упорядочиваем работу с шумом. Точное решение возможно только в случае отсутствия шума. Без шума, k=0 и выражение упрощается до F = G ./ H .

Кроме того, ось x отключена на графике, созданном findpeaks , но это не должно повлиять на результаты. Местоположения, возвращаемые им, являются индексами в массиве, эти же индексы можно использовать для индексации в t и нахождения фактических местоположений пиков.