Документация Engee
Notebook

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

Открыть пример в Engee

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

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

In [ ]: 
Pkg.add(["Noise", "DSP", "DigitalComm"])

   Resolving package versions...
    Updating `~/.project/Project.toml`
  [a7b11256] + DigitalComm v1.2.0
  [81d43f40] + Noise v0.3.3
    Updating `~/.project/Manifest.toml`
  [a7b11256] + DigitalComm v1.2.0
  [81d43f40] + Noise v0.3.3
  [d96e819e] + Parameters v0.12.3
  [e409e4f3] + PoissonRandom v0.4.4
  [3a884ed6] + UnPack v1.0.2
Precompiling project...
  ✓ Noise
  1 dependency successfully precompiled in 4 seconds. 23 already precompiled.
In [ ]: 
using Plots; # Билиотека отрисовки графиков
using Noise; # Библиотека добавления шумов
using DigitalComm; # Библиотека систем связи
using DSP; # Библиотека цифровой обработки сигналов
using FFTW; # Библиотека преобразований Фурье
plotlyjs(); # Подключения метода отрисовки

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

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

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

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

plot(t,y)

┌ Warning: attempting to remove probably stale pidfile
│   path = "/user/.jlassetregistry.lock"
└ @ Pidfile /usr/local/ijulia-core/packages/Pidfile/DDu3M/src/Pidfile.jl:260
Out[0]:

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

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

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

Квантование

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

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

# Вывод результата
plot(y[1:250], label="исходный сигнал")
plot!(y_quant[1:250], label="квантованный сигнал по 80 уровням")
Out[0]:
In [ ]: 
y_p = DSP.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)
Out[0]:

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

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

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

In [ ]: 
y_d = y_quant[1:2:end]
t_new = t[1:2:end]
plot(t_new,y_d)
Out[0]:
In [ ]: 
y_p = DSP.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)
Out[0]:

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

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

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

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

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

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

  • равномерная спектральная плотность мощности (поэтому он называется белым);
  • нормальное распределение временных значений (поэтому он называется гауссовским);
  • суммируется с полезным сигналом (поэтому он называется аддитивным);
  • АБГШ статистически независим от полезного сигнала.
In [ ]: 
signal = addNoise(y_d, 20); # Добавление шумов

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

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

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

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

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

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

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

Out_sig = filt(f,signal);

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

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

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

In [ ]: 
h = impresp(f);
Out_sig2 = conv(signal,h);

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

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

In [ ]: 
plot(Out_sig[1:2500]-Out_sig2[1:2500])
Out[0]:

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

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

In [ ]: 
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="исходный сигнал")
Out[0]:

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

In [ ]: 
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")
Out[0]:

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

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

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

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

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

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

In [ ]: 
# 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)
Out[0]:

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

In [ ]: 
# 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)
Out[0]:

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

In [ ]: 
println("Размер входа: ",size(Out_fir,1))
println("Размер при downsampling: ",size(rs,1))
println("Размер при upsampling: ",size(rs2,1))

Размер входа: 2501
Размер при downsampling: 1251
Размер при upsampling: 5001
In [ ]: 
plot(rs)
plot!(rs2)
Out[0]:

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

Вывод

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