2014 dxdy logo

Научный форум dxdy

Математика, Физика, Computer Science, Machine Learning, LaTeX, Механика и Техника, Химия,
Биология и Медицина, Экономика и Финансовая Математика, Гуманитарные науки




Начать новую тему Ответить на тему
 
 Подсчет сегментов имульсов в сигнале MATLAB
Сообщение14.07.2022, 19:06 


17/03/20
183
Уважаемые коллеги, добрый вечер! В общем, возникла проблема следующего характера.

Имеется сигнал, достаточно большой объем точек, длительность порядка 35-40 секунд. Он представлят собой последовательность импульсов тока, длительностью примерно 5 мкс. Есть, например, такой сигнал, где последовательности сигналов (импульсов) составляют такие участки (sig1.mat):

Изображение


Суть задачи такова, что необходимо подсчитывать количество таких сегментов (выделены рамкой), и количесто импульсов в кажом таком сегменте, собственно говоря, просуммировав весь массив получим общее число рабочих импульсов. Для начала, я хотел использовать конструкцию счетчика путем применения логической индекасации, полагая, что однозначно, импульс по амплитуде имеет значение больше 200. Тогда можно для каждого импульса брать, определять начальный пик имульса, затем находить разницу до следующего пика, который больше порогового значения (200) и как только дальше идет промежуток, когда значение около нуля, т.е интервал следования между импульсами (больше нуля) использовать nnz. Т.е что-то примерно такое:

Используется синтаксис Matlab M
numberOfPulses = nnz(diff(outSig_1_curr > 1) > 0);


Однако, в сигнале могут присутствовать и такие фрагменты файл sig16.mat:

Изображение


Тогда попрбовал напрямую использовать индексацию, и вычислять интервалы вручную (к сожалению в этом случае минимальный размер окна, т.е. длительность такого фрагмента 0.1. мс и расстояние между импульсами 0.1 мс):

код: [ скачать ] [ спрятать ]
Используется синтаксис Matlab M
clc
clearvars

load('sig1.mat') % outSig_1_curr2      32800000x1
samples = numel(outSig_1_curr2);
Fs = 1000000; %  Hz
dt = 1/Fs;
t = (0:samples-1)*dt;

% извлечь сигнал между временными метками t = 21.9112 и t = 21.94 s
ind1 = round(21.9112*Fs);
ind2 = round(21.94*Fs);
t = t(ind1:ind2);
outSig_1_curr2 = outSig_1_curr2(ind1:ind2);


%% находим пики
% положительный максимум
[tf, P] = islocalmax(outSig_1_curr2,'MinProminence',200);
t_peak = t(tf);
outSig_1_peak = outSig_1_curr2(tf);

%% определить и извлечь длину корректного блока данных
blocks_distance_samples = 100; % select data blocks  that are distant above value = blocks_distance_samples (in samples)
min_contiguous_samples = 100; % select data blocks  that are contiguous with min amount of samples = min_contiguous_samples
detect_signal = zeros(size(t))+0.5; % initialisation

    % определить начало и конец сегмента выше порогового значения
        %%%%%%%%%%

        % начало и конец точек для каждого сегмента
        ind = find(tf);
        D = diff([0;ind;0]);
        ind_ends = find(D > blocks_distance_samples) - 1;
        ind_ends = ind_ends(ind_ends>0); % удаляем значения нуля
        ends = ind(ind_ends);
        begin = ind(ind_ends+1);

        % удалить некорректные значения начала и конца сегмента
        ends = ends(ends>begin(1));
        begin = begin(begin<ends(end));
        m = min(numel(begin),numel(ends));
        ends = ends(1:m);
        begin = begin(1:m);
       
       
        %%%%%%%%%%
       
    length_ind = ends - begin; % длина блока
    ind2= length_ind>min_contiguous_samples; % проверяем корректность длины (приблизительно min_contiguous_samples значение)
    begin = begin(ind2); % выбранные точки
    ends = ends(ind2); % выбранные точки
   
    t_peak2 = t(ind);
    outSig_1_peak2 = outSig_1_curr2(ind);
    % определить начало и конец импульса и значений амплитуды фрагмента
    t_peak2_begin = t(begin);
    outSig_1_curr2_begin = outSig_1_curr2(begin);
    t_peak2_ends = t(ends);
    outSig_1_curr2_ends = outSig_1_curr2(ends);
   
for ci = 1:length(begin)
    ind = (t>=t_peak2_begin(ci) & t<=t_peak2_ends(ci));
    detect_signal(ind) = 200;
    ind4 = (t_peak>=t_peak2_begin(ci) & t_peak<=t_peak2_ends(ci));
    number_of_peaks(ci) = numel(find(ind4));
   
end
number_of_peaks

% plot
figure(1)
plot(t,outSig_1_curr2,t_peak,outSig_1_peak,'dr',t(begin),outSig_1_curr2(begin),'dg',t(ends),outSig_1_curr2(ends),'dk',t,detect_signal,'k');
legend('signal',' peaks ','begin points','end points','selected signal')
 


Однако, к сожалению, в этом случае очень тяжело распознать фргаменты со второго рисунка и подсчитывать такие фрагменты. Между импульсами в таких фрагментах длительность почти нулевая, и задать размер между импульсами в 1-2 отсчета не получается, происходит ошибка индексации.

Каким образом, можно задать расстояние в 2-3 сэмпла и при этом различать такие фрагменты? Или можно же использовать конструкцию с nnz, и в полученном массиве логических единиц определять общее число импульсов и имульсов фрагмента?

К примеру:
код: [ скачать ] [ спрятать ]
Используется синтаксис Matlab M
idx1 = [0 0 0 1 1 1 0 0 1 1 1 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0 1 1 1 0 0];
diffs = diff([0 idx1 0]);

% первый индекс
loc = find(diffs == 1);

Ns = find(diffs == -1) - loc;

C = mat2cell([loc' Ns'],ones(size(loc)))

  4×1 cell array

    {[ 4 3]}
    {[ 9 3]}
    {[16 6]}
    {[29 3]}
 


Может быть есть готовые решения или более корректное использование функций islocalmax, islocalmin для подсчета ипульсов? Возможно ли задать параметр 'MinProminence как дииапазон (например, часть фргаментов, которые на втором рисунке лежит примерно выше 1500 по амлитуде, а каждый пик импульса такого фргамента приблизительно выше порогового значения max(signal)/1.03321). На данный момент я не могу предположить каким образом, можно подсчитать количесвто таких импульсов в каждом сегменте и количество таких сегментов.

Сами сигналы доступны по ссылке: https://drive.google.com/drive/folders/1NB3QQXVV023fa4fdK6MzHqLHI0WH3aLC?usp=sharing

 Профиль  
                  
 
 Re: Подсчет сегментов имульсов в сигнале MATLAB
Сообщение14.07.2022, 23:54 


10/03/16
3864
Aeroport
Для началa нужно построить огибающую, и затем находить changepoints в огибающей.

 Профиль  
                  
 
 Re: Подсчет сегментов имульсов в сигнале MATLAB
Сообщение15.07.2022, 00:14 


17/03/20
183
ozheredov
Достаточно ли построить огибающую путем преобразования Гильберта?

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

Просто проблема в том, что такие фрагменты могут быть чуть ниже и выше, и их достаточно много...

Здесь все же лучше пользоваться просто подсчётом через islocalmax, или потом пробовать findpeaks?

 Профиль  
                  
 
 Re: Подсчет сегментов имульсов в сигнале MATLAB
Сообщение15.07.2022, 02:43 


10/03/16
3864
Aeroport
Alm99 в сообщении #1560190 писал(а):
преобразования Гильберта?


Я думаю, преобразование Гильберта тут не сработает -- слишком резкие changepoint'ы.

Alm99 в сообщении #1560190 писал(а):
вот построил огибающую, дальше подсчитываю уже только число максимальных пиков только для точек, которые имеют значение огибающей большее, чем заданное пороговое?


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

Alm99 в сообщении #1560190 писал(а):
Просто проблема в том, что такие фрагменты могут быть чуть ниже и выше, и их достаточно много...


Да, это вполне вписывается.

Alm99 в сообщении #1560190 писал(а):
Здесь все же лучше пользоваться просто подсчётом через islocalmax, или потом пробовать findpeaks?


Не нужно проверять чего-то на локальный максимум или искать пики - у Вас совсем другая задача, как я понял.

 Профиль  
                  
 
 Re: Подсчет сегментов имульсов в сигнале MATLAB
Сообщение15.07.2022, 08:44 


17/03/20
183
ozheredov

Я подсчитывал с помщью пиков центр импульса и его границы (длительность)...

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

Изображение


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

 Профиль  
                  
 
 Re: Подсчет сегментов имульсов в сигнале MATLAB
Сообщение15.07.2022, 12:50 


10/03/16
3864
Aeroport
Alm99

Нужно скользящее среднее огибающей, а не исходного сигнала -- тогда картинка будет "резче".

 Профиль  
                  
 
 Re: Подсчет сегментов имульсов в сигнале MATLAB
Сообщение15.07.2022, 14:48 


17/03/20
183
ozheredov

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

Изображение


Однако, я получаю все равно ошибку индексации:
Код:
Array indices must be positive integers or logical values.

Error in movavg (line 43)
        ends = ind(ind_ends);


код: [ скачать ] [ спрятать ]
Используется синтаксис Matlab M
  1. step = 1/freq;
  2. time_sig = sig16(:,1);
  3. data = sig16(:,2);
  4. % upp1_norm - огибающая
  5. [tf, P] = islocalmin(upp1_norm,'MinProminence',min(upp1_norm));
  6. t_peak = time_sig(tf);
  7. outSig_1_peak = upp1_norm(tf);
  8.  
  9.  
  10. %% localize / extract valid data blocks
  11. blocks_distance_samples = 100; % длина блока, состоящие из импульсов, где огибающая ниже средней линии
  12. min_contiguous_samples = 100; % рамзер окна
  13. detect_signal = zeros(size(time_sig))+0.5;
  14.  
  15.     % старт и конец блока с определенным пороговым значением огибающей min(upp1_norm)
  16.    
  17.         %%%%%%%%%%
  18.         % This locates the beginning /ending points of data groups
  19.         ind = find(tf);
  20.         D = diff([0;ind;0]);
  21.         ind_ends = find(D > blocks_distance_samples) - 1;
  22.         ends = ind(ind_ends);
  23.         begin = ind(ind_ends+1);
  24.         % удалить начальные и конечне некорректные точки сегмента
  25.         ends = ends(ends>begin(1));
  26.         begin = begin(begin<ends(end));
  27.         m = min(numel(begin),numel(ends));
  28.         ends = ends(1:m);
  29.         begin = begin(1:m);
  30.        
  31.        
  32.         %%%%%%%%%%
  33.        
  34.     length_ind = ends - begin; % вычисляем длину блока
  35.     ind2= length_ind>min_contiguous_samples; % проверка корректности длины
  36.     begin = begin(ind2); % выбранные точки
  37.     ends = ends(ind2); % выбранные точки
  38.    
  39.     t_peak2 = t(ind);
  40.     outSig_1_peak2 = outSig_1_curr2(ind);
  41.  
  42.     t_peak2_begin = t(begin);
  43.     outSig_1_curr2_begin = outSig_1_curr2(begin);
  44.     t_peak2_ends = t(ends);
  45.     outSig_1_curr2_ends = outSig_1_curr2(ends);
  46.    
  47. for ci = 1:length(begin)
  48.     ind = (t>=t_peak2_begin(ci) & t<=t_peak2_ends(ci));
  49.     detect_signal(ind) = round(max(sig16(:,2))/1.03321);
  50.  
  51.     ind4 = (t_peak>=t_peak2_begin(ci) & t_peak<=t_peak2_ends(ci));
  52.     number_of_peaks(ci) = numel(find(ind4));
  53.    
  54. end
  55.  
  56. number_of_peaks;
  57.  

 Профиль  
                  
 
 Re: Подсчет сегментов имульсов в сигнале MATLAB
Сообщение15.07.2022, 17:37 


10/03/16
3864
Aeroport
Alm99
Чет картинка не хочет появляться

 Профиль  
                  
 
 Re: Подсчет сегментов имульсов в сигнале MATLAB
Сообщение15.07.2022, 18:42 


17/03/20
183
ozheredov

Попробуйте пожалуйста, еще раз, я повторно код загружу. Сигнал sig16.mat можно скачать по ссылке из начального поста
код: [ скачать ] [ спрятать ]
Используется синтаксис Matlab M
clear
clc

load('sig_no_timetable.mat', 'freq')
load('sig_no_timetable.mat', 'sig16')
step = 1/freq;
time_sig = sig16(:,1);
data = sig16(:,2);

freq = 1e6;
step = 1/freq;
time_sig = sig16(:,1);
data = sig16(:,2);
 
fl1 = 120;
[up1,lo1] = envelope(data,1000,'analytic');
upp1 = filter(ones(1,fl1)/fl1,1,up1);
upp1_norm = circshift(upp1,-95);
 
figure(1)
plot(time_sig,data)
hold on
plot(time_sig,upp1_norm);
xlim([4.435 4.445]);
hold off

% upp1_norm - огибающая
[tf, P] = islocalmin(upp1_norm,'MinProminence',min(upp1_norm));
t_peak = time_sig(tf);
outSig_1_peak = upp1_norm(tf);
 
 
%% localize / extract valid data blocks
blocks_distance_samples = 100; % длина блока, состоящие из импульсов, где огибающая ниже средней линии
min_contiguous_samples = 100; % рамзер окна
detect_signal = zeros(size(time_sig))+0.5;
 
    % старт и конец блока с определенным пороговым значением огибающей min(upp1_norm)
   
        %%%%%%%%%%
        % This locates the beginning /ending points of data groups
        ind = find(tf);
        D = diff([0;ind;0]);
        ind_ends = find(D > blocks_distance_samples) - 1;
        ends = ind(ind_ends);
        begin = ind(ind_ends+1);
        % удалить начальные и конечне некорректные точки сегмента
        ends = ends(ends>begin(1));
        begin = begin(begin<ends(end));
        m = min(numel(begin),numel(ends));
        ends = ends(1:m);
        begin = begin(1:m);
       
       
        %%%%%%%%%%
       
    length_ind = ends - begin; % вычисляем длину блока
    ind2= length_ind>min_contiguous_samples; % проверка корректности длины
    begin = begin(ind2); % выбранные точки
    ends = ends(ind2); % выбранные точки
   
    t_peak2 = t(ind);
    outSig_1_peak2 = upp1_norm(ind);
 
    t_peak2_begin = t(begin);
    outSig_1_curr2_begin = upp1_norm(begin);
    t_peak2_ends = t(ends);
    outSig_1_curr2_ends = upp1_norm(ends);
   
for ci = 1:length(begin)
    ind = (t>=t_peak2_begin(ci) & t<=t_peak2_ends(ci));
    detect_signal(ind) = min(upp1_norm);
 
    ind4 = (t_peak>=t_peak2_begin(ci) & t_peak<=t_peak2_ends(ci));
    number_of_peaks(ci) = numel(find(ind4));
   
end
 
number_of_peaks;

 

 Профиль  
                  
 
 Re: Подсчет сегментов имульсов в сигнале MATLAB
Сообщение27.10.2022, 17:33 
Аватара пользователя


26/05/12
1534
приходит весна?
А вы можете выложить отрывок ваших данных с обсуждаемым фрагментом, но весом не 200 мегов, а, скажем, не больше 5? Накладно качать ваши данные только чтобы поиграться. Хотя хотелось бы помочь.

 Профиль  
                  
Показать сообщения за:  Поле сортировки  
Начать новую тему Ответить на тему  [ Сообщений: 10 ] 

Модераторы: Karan, Toucan, PAV, maxal, Супермодераторы



Кто сейчас на конференции

Сейчас этот форум просматривают: нет зарегистрированных пользователей


Вы не можете начинать темы
Вы не можете отвечать на сообщения
Вы не можете редактировать свои сообщения
Вы не можете удалять свои сообщения
Вы не можете добавлять вложения

Найти:
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group