Фильтрация в 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)
Проведём оценку спектральной плотности мощности сигнала через 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)
Квантование
Следующим этапом преобразования нашего сигнала станет его квантование — разбиение диапазона отсчётных значений сигнала на конечное число уровней и округление этих значений до одного из двух ближайших к ним уровней.
y_quant = quantization(y, 80, minv=-1, maxv=1); # Квантование исходной синусоиды
# Вывод результата
plot(y[1:250], label="исходный сигнал")
plot!(y_quant[1:250], label="квантованный сигнал по 80 уровням")
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)
Как мы видим, данное преобразование не влияет на спектральную плотность мощности сигнала и не смещает его частоту.
Прореживание данных
Следующее действие, которое мы проведём с исходной синусоидой, — это прореживание данных. Мы возьмём каждый второй отсчёт из исходного сигнала.
y_d = y_quant[1:2:end]
t_new = t[1:2:end]
plot(t_new,y_d)
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)
Периодограмма показывает, что прореживание данных увеличивает частоту сигнала. Это связано с тем, что из-за уменьшения количества отсчётов в два раза частота колебаний соответственно увеличилась в два раза.
Наложение шумов
Шум — это нежелательные явления, мешающие нам получать информацию из полезного сигнала.
Шум случаен по своей природе.
В цифровой обработке сигналов чаще всего используется аддитивный белый гауссовский шум (АБГШ).
Его характеристики:
-
равномерная спектральная плотность мощности (поэтому он называется белым);
-
нормальное распределение временных значений (поэтому он называется гауссовским);
-
суммируется с полезным сигналом (поэтому он называется аддитивным);
-
АБГШ статистически независим от полезного сигнала.
signal = addNoise(y_d, 20); # Добавление шумов
plot(t_new, signal, label="синусоида с шумом")
plot!(t_new, y_d, label="синусоида без шума")
Важно отметить, что функция 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)trueplot(signal, label="исходный сигнал")
plot!(Out_sig, label="отфильтрованный сигнал")
Второй вариант задания такого же фильтра — через импульсную характеристику при помощи функции impresp, а затем вычисление выхода фильтра с помощью функции conv.
h = impresp(f);
Out_sig2 = conv(signal,h);
plot(signal, label="исходный сигнал")
plot!(Out_sig2, label="отфильтрованный сигнал")
Сравним выход функции filt и выход функции conv.
plot(Out_sig[1:2500]-Out_sig2[1:2500])
Как мы видим из графика, выходы функции 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="исходный сигнал")
Для сравнения выходов фильтров воспользуемся периодограммой Уэлча.
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")
Как мы видим, основные различия между выходами фильтров кроются в их спектральной плотности мощности.
Изменение частоты сигнала
Для того, чтобы увеличить или уменьшить количество отсчётов, в 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)
Как мы видим, частота отфильтрованного сигнала увеличилась вдвое.
# 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)
Как мы видим, частота отфильтрованного сигнала уменьшилась вдвое.
println("Размер входа: ",size(Out_fir,1))
println("Размер при downsampling: ",size(rs,1))
println("Размер при upsampling: ",size(rs2,1))Размер вохда: 2501
Размер при downsampling: 1251
Размер при upsampling: 5001plot(rs)
plot!(rs2)
Те же эффекты мы можем наблюдать при сравнении сигналов и их размеров. При downsampling уменьшается количество отсчётов сигнала, а при upsampling, наоборот, увеличивается относительно исходного сигнала.