Фильтрация в Engee

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

Подключение библиотек

using Plots; # Билиотека отрисовки графиков
using Noise; # Библиотека добавления шумов
using DigitalComm; # Библиотека систем связи
using DSP; # Библиотека цифровой обработки сигналов
using FFTW; # Библиотека преобразований Фурье
plotlyjs(); # Подключения метода отрисовки

Задание синусоиды

Над данной синусоидой будут проведены все преобразования в примере.

fs = 5000; # Шагов в секунду
t = [0:1/fs:1;]; # Время

y = sin.(2*pi*10*t); # Задание синусоиды

plot(t,y)

interactive-scripts/images/dsp_Filtration/6cae79378dfd677717f98e600a9dc878853165b0

Проведём оценку спектральной плотности мощности сигнала через periodogram — это функция построения периодограммы, основанная на вычислении квадрата модуля преобразования Фурье для последовательности данных.

Спектральная плотность мощности (СПМ) в физике и обработке сигналов — функция, описывающая распределение мощности сигнала в зависимости от частоты, то есть мощность, приходящаяся на единичный интервал частоты.

y_p = periodogram(y, onesided=true, nfft=length(y), fs=fs);
plot(freq(y_p)[1:20], power(y_p)[1:20], xlabel="частота, Гц", ylabel="спектральная плотность мощности", legend=false)

interactive-scripts/images/dsp_Filtration/12b466a10e4a1e9958d53d3028914a3838830fe1

Квантование

Следующим этапом преобразования нашего сигнала станет его квантование — разбиение диапазона отсчётных значений сигнала на конечное число уровней и округление этих значений до одного из двух ближайших к ним уровней.

y_quant = quantization(y, 80, minv=-1, maxv=1); # Квантование исходной синусоиды

# Вывод результата
plot(y[1:250], label="исходный сигнал")
plot!(y_quant[1:250], label="квантованный сигнал по 80 уровням")

interactive-scripts/images/dsp_Filtration/e77caedbbcb8f38b25028f6451e6dba3a76e3284

y_p = periodogram(y_quant, onesided=true, nfft=length(y_quant), fs=fs);
plot(freq(y_p)[1:20], power(y_p)[1:20], xlabel="частота, Гц", ylabel="спектральная плотность мощности", legend=false)

interactive-scripts/images/dsp_Filtration/ffc205518e00eab3f8eb0ae2fc89b4ffd9c54ddd

Как мы видим, данное преобразование не влияет на спектральную плотность мощности сигнала и не смещает его частоту.

Прореживание данных

Следующее действие, которое мы проведём с исходной синусоидой, — это прореживание данных. Мы возьмём каждый второй отсчёт из исходного сигнала.

y_d = y_quant[1:2:end]
t_new = t[1:2:end]
plot(t_new,y_d)

interactive-scripts/images/dsp_Filtration/1251afdb86475738428918a019a221cba66913d5

y_p = periodogram(y_d, onesided=true, nfft=length(y_d), fs=fs);
plot(freq(y_p)[1:20], power(y_p)[1:20], xlabel="частота, Гц", ylabel="спектральная плотность мощности", legend=false)

interactive-scripts/images/dsp_Filtration/dc981cb4cc3f7626398fd96efff06d0956106b73

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

Наложение шумов

Шум — это нежелательные явления, мешающие нам получать информацию из полезного сигнала.

Шум случаен по своей природе.

В цифровой обработке сигналов чаще всего используется аддитивный белый гауссовский шум (АБГШ).

Его характеристики:

  • равномерная спектральная плотность мощности (поэтому он называется белым);

  • нормальное распределение временных значений (поэтому он называется гауссовским);

  • суммируется с полезным сигналом (поэтому он называется аддитивным);

  • АБГШ статистически независим от полезного сигнала.

signal = addNoise(y_d, 20); # Добавление шумов

plot(t_new, signal, label="синусоида с шумом")
plot!(t_new, y_d, label="синусоида без шума")

interactive-scripts/images/dsp_Filtration/304548ed163831cc4236ba1c7879c1b56688294e

Важно отметить, что функция addNoise задаёт кортеж, из которого нам нужен только первый вектор. Далее в коде показаны типы данных входа addNoise, его выхода и результата извлечения из кортежа первого вектора.

typeof(y_d)
Vector{Float64} (alias for Array{Float64, 1})
typeof(signal)
Tuple{Vector{Float64}, Vector{Float64}}
signal=signal[1];
typeof(signal)
Vector{Float64} (alias for Array{Float64, 1})

Фильтрация сигнала в Engee

В данном примере мы разберём три фильтра, но в документации Engge вы можете найти любой необходимый вам фильтр. Первый способ, который мы рассмотрим, — это задание фильтра через коэффициенты фильтра, в виде двух векторов a и b, затем реализуем представление фильтра через коэффициенты числителя b и знаменатель a передаточной функции с помощью функции PolynomialRatio, содержащейся в библиотеке DSP.

Примечание: b и a могут быть указаны как объекты Polynomial, или векторы, упорядоченные от наибольшей степени к наименьшей.

b = [0.75, 0.5];
a = [1, 0.2, 0.1];
f = PolynomialRatio(b,a);

Out_sig = filt(f,signal);

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

size(Out_sig,1)==size(signal,1)
true
plot(signal, label="исходный сигнал")
plot!(Out_sig, label="отфильтрованный сигнал")

interactive-scripts/images/dsp_Filtration/108a121de19d3eba356a7a16f63339f1b6e3e779

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

h = impresp(f);
Out_sig2 = conv(signal,h);

plot(signal, label="исходный сигнал")
plot!(Out_sig2, label="отфильтрованный сигнал")

interactive-scripts/images/dsp_Filtration/6735166d78243d23c0099675e11335e93849d596

Сравним выход функции filt и выход функции conv.

plot(Out_sig[1:2500]-Out_sig2[1:2500])

interactive-scripts/images/dsp_Filtration/2f7b7f8fa6210a934c23f50e41c531f78e4569e9

Как мы видим из графика, выходы функции filt и выход функции conv практически полностью совпадают.

Третий вариант фильтрации, который мы разберём, — это задание оконного КИХ-фильтра. Он будет отличаться от предыдущих, так как имеет задержку относительно входного сигнала, равную размеру окна.

responsetype = Lowpass(5; fs=30) # Определение полосы пропускания
designmethod = FIRWindow(hanning(128)) # Определение размеров окна

Out_fir = filt(digitalfilter(responsetype, designmethod), Out_sig)

plot(Out_sig, label="исходный сигнал")
plot!(Out_fir, label="исходный сигнал")

interactive-scripts/images/dsp_Filtration/6c05cb2675b4b894333ba8796a4edea130cd519b

Для сравнения выходов фильтров воспользуемся периодограммой Уэлча.

n=div(length(Out_sig),8);
y_out = welch_pgram(Out_sig, n, div(n,2), onesided=true, nfft=length(Out_sig), fs=fs);

n=div(length(Out_sig2),8);
y_in = welch_pgram(Out_sig2, n, div(n,2), onesided=true, nfft=length(Out_sig2), fs=fs);

n=div(length(Out_fir),8);
y_fir = welch_pgram(Out_fir, n, div(n,2), onesided=true, nfft=length(Out_fir), fs=fs);

plot(freq(y_in)[1:25], power(y_in)[1:25], xlabel="частота, Гц", ylabel="спектральная плотность мощности", label="filt")
plot!(freq(y_out)[1:25], power(y_out)[1:25], xlabel="частота, Гц", ylabel="спектральная плотность мощности", label="conv")
plot!(freq(y_fir)[1:25], power(y_fir)[1:25], xlabel="частота, Гц", ylabel="спектральная плотность мощности", label="fir")

interactive-scripts/images/dsp_Filtration/7124ccf3aecc75feb3894e3a0ed8fe11145d8ca8

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

Изменение частоты сигнала

Для того, чтобы увеличить или уменьшить количество отсчётов, в Engee используется команда resample, которая позволяет выполнить upsampling и downsampling сигнала.

Примечание: upsampling и downsampling задаются через отношение n//k, где:

n: количество точек, которые необходимо создать;

k: количество ближайших соседей, которых следует учитывать для каждой точки.

# Downsampling
rs = resample(Out_fir, 1//2)

n=div(length(rs),8);
y_out = welch_pgram(rs, n, div(n,2), onesided=true, nfft=length(rs), fs=fs);
plot(freq(y_out)[1:25], power(y_out)[1:25], xlabel="частота, Гц", ylabel="спектральная плотность мощности", legend=false)

interactive-scripts/images/dsp_Filtration/5c3be998787aacb2913f80f3556b042521042b81

Как мы видим, частота отфильтрованного сигнала увеличилась вдвое.

# Upsampling
rs2 = resample(Out_fir, 2//1)

n=div(length(rs2),8);
y_out = welch_pgram(rs2, n, div(n,2), onesided=true, nfft=length(rs2), fs=fs);
plot(freq(y_out)[1:25], power(y_out)[1:25], xlabel="частота, Гц", ylabel="спектральная плотность мощности", legend=false)

interactive-scripts/images/dsp_Filtration/ae8b12f8265bbc419e603fad3485fcd16439ea72

Как мы видим, частота отфильтрованного сигнала уменьшилась вдвое.

println("Размер входа: ",size(Out_fir,1))
println("Размер при downsampling: ",size(rs,1))
println("Размер при upsampling: ",size(rs2,1))
Размер вохда: 2501
Размер при downsampling: 1251
Размер при upsampling: 5001
plot(rs)
plot!(rs2)

interactive-scripts/images/dsp_Filtration/07f47b8307f4816605031af5119a412cd22f8e6c

Те же эффекты мы можем наблюдать при сравнении сигналов и их размеров. При downsampling уменьшается количество отсчётов сигнала, а при upsampling, наоборот, увеличивается относительно исходного сигнала.

Вывод

В данном примере мы продемонстрировали возможности взаимодействия с сигналами и показали, что Engee может быть в полной мере применим для цифровой обработки сигналов.

Загрузить пример в Engee

Чтобы скачать пример себе в файловый браузер, скопируйте и выполните следующие команды. Пример появится по адресу start/examples в соответствующей папке.

# Проверка наличия директории для примера
if isdir("/user/start/examples/dsp/filtration")
	rm("/user/start/examples/dsp/filtration"; force = true, recursive = true)
end

# Копирование примера на диск Engee
run(`git clone https://git.engee.com/learn-engee/examples/dsp/filtration.git /user/start/examples/dsp/filtration`)