Основы цифровой обработки сигналов
03. Анализ и обработка во временной области
Анализируем и обрабатываем сигнал с датчика акустических волн. Стандартные задачи временной обработки - выделить из сигнала информацию, изменить его характеристики и статистические показатели, масштабировать его по уровню или сгладить его форму.
В данном примере выполним основные этапы технических вычислений - чтение данных из файла стандартного формата, визуализация, предобработка, оценка метрик, поиск локальных экстремумов, изменение формы сигнала и запись результатов в таблицу.
Импорт и визуализация сигнала
Отсчёты выхода АЦП с пропусками хранятся в файле data_missing.txt. Считаем их в виде матрицы с двумя столбцами в переменную datamat:
using DelimitedFiles
datamat = readdlm("data_missing.txt")
Выделим из неё отдельные вектора для времени и исходных значений с пропусками, отобразим на графике во временной области:
t = datamat[:,1];
original = datamat[:,2];
plot(t, original, xguide = "Время (с)",
title = "Исходный сигнал",
legend = false)
Определим программно индексы пропущенных отсчётов (NaN). Заменим их в исходном векторе нулями:
indices = isequal.(original, NaN);
zerosig = datamat[:,2];
zerosig[indices] .= 0;
plot(t, zerosig, xguide = "Время (с)",
title = "Cигнал с нулями в пропусках",
legend = false)
Необходимо "смоделировать" пропущенные значения тем или иным способом. Рассмотрим интерполяцию и авторегрессию.
Заполнение пропусков
Используем функцию fillgaps библиотеки EngeeDSP:
using EngeeDSP.Functions
# filled = fillgaps(original,3,1);
filled = fillgaps(original,40,40);
plot(t, filled, xguide = "Время (с)",
title = "Заполнение пропусков",
legend = false)
Удаление постоянной составляющей
Избавимся от постоянной составляющей при помощи функции detrend из набора функций EngeeDSP. Под капотом функции применяется метод приближения значений сигнала полиномом произвольного порядка, а затем рассчитанный полином вычитается из сигнала. Используем полином 12-го порядка:
notrend = detrend(filled, 12);
trend = filled - notrend;
plot(t, filled, xguide = "Время (с)",
title = "Постоянная составляющая",
label = "Исходный сигнал")
plot!(t, trend, lw=:5, label = "Тренд")
plot(t, notrend, xguide = "Время (с)",
title = "Сигнал без постоянной составляющей",
leg = false)
Статистика сигнала
Посчитаем самые "популярные" статистические метрики сигнала - экстремумы, средние значения, размах и прочее. Но сначала определим из вектора отсчётов времени период дискретизации сигнала:
dt = EngeeDSP.Functions.mean(diff(t))
И частоту дискретизации сигнала:
fs = 1/dt
Определим длительность дискретного сигнала - количество отсчётов в векторе:
nsamples = length(notrend)
Найдём максимальное значение и порядковый номер максимального отсчёта в численном векторе:
maxval, maxidx = EngeeDSP.Functions.max(notrend)
То же самое для минимального значения. Это - абсолютные экстремумы сигнала:
minval, minidx = EngeeDSP.Functions.min(notrend)
Зная их, можно определить размах сигнала:
range = abs(maxval - minval)
Масштабирование сигнала
Зная максимальное значение, можно поделить каждый отсчёт вектора на это значение, тем самым обеспечив колебания сиганала в пределах от -1 до +1:
sig = notrend ./ maxval;
plot(t, sig, xguide = "Время, (с)",
title = "Масштабированный сигнал",
leg = false)
Подсчитаем среднее арифметическое (обратите внимание, оно по значению близко к нулю):
EngeeDSP.Functions.mean(sig)
Найдём среднеквадратическое значение (RMS):
rmsval = EngeeDSP.Functions.rms(sig)
Наконец, подсчитаем отношение величин двух самых больших выбросов в сигнале:
peak2peak(sig)
И отношение максимального значения сигнала к среднеквадратическому:
peak2rms(sig)
Поиск пиков сигнала
Зачастую, полезно уметь программно определять локальные экстремумы сигнала, или так называемые пики. Для этого воспользуемся функцией findpeaks. Функция возвращает величины и индексы пиков, а также (опционально) их ширину и выраженность. В качестве дополнительных входных аргументов укажем минимальную высоту пика и минимальное расстояние между соседними пиками:
pks, locs, w, p = EngeeDSP.Functions.findpeaks(sig, out=:data,
MinPeakHeight = 0.3,
MinPeakDistance = 200);
plot(t, sig, xguide = "Время (с)", label = false)
scatter!(t[locs], pks, label = "Пики")
Выведем амплитуды и соответствующие им моменты времени программно. Также выведем ширину и выраженность:
hcat(t[locs], pks, w, p)
Фильтрация скользящим окном
Избавимся от одиночных выбросов при помощи нелинейного фильтра, а именно - бегущей медианы с размером окна в три отсчёта. Для этого применим функцию movmedian:
nospikes = movmedian(sig, 3);
plot(t,sig, xguide = "Время (с)", label = "До фильтра")
plot!(t, nospikes, linewidth = 1.5, label = "Медианный фильтр")
И немного сгладим форму сигнала фильтром бегущего среднего, который мы вызовем из функции smoothdata:
smooth, metric = smoothdata(nospikes,"movmean",5)
plot(t, nospikes, xguide = "Время (с)", label = "До фильтра")
plot!(t, smooth, linewidth = 2, label = "Сглаживающий фильтр")
Оценка огибающей сигнала
Есть различные способы оценки огибающей сигнала при помощи функции envelope: фильтр Гильберта, оценка СКО скользящим окном, и приближение сплайновой интерполяцией по локальным экстремумам. Рассмотрим последний способ, и укажем, что точки, по которым будет проводиться интерполяция, должны располагаться минимум через 14 отсчётов сигнала друг от друга:
plot(t, smooth, label = "Сигнал")
yupper, ylower = envelope(smooth, 14, "peak", out=:data);
plot!(t, yupper, lw=2, label = "Верх")
plot!(t, ylower, lw=2, label = "Низ")
Сохранение результата
Воспользуемся библиотекой DataFrames для записи результатов в таблицу и библиотекой XLSX для сохранения данных в файл Excel:
using DataFrames, XLSX
mytable = DataFrame(
Время = t,
Сигнал = smooth,
Верх = yupper,
Низ = ylower)
XLSX.writetable("output.xlsx", mytable)
В завершении послушаем обработанный сигнал как аудио-фрагмент:
include("audioplayer.jl");
longsig = repeat(smooth,4);
audioplayer(longsig*2,4800)