Системный объект КИХ-фильтра

Системный объект - это специализированная программная конструкция в Engee со специфическими методами и свойствами для моделирования того или иного алгоритма. Системные объекты полезны для моделирования динамических систем и обработки потоковых данных. В данном примере мы рассмотрим применение системного объекта КИХ-фильтра из библиотеки EngeeDSP для моделирования потоковой цифровой обработки (ЦОС) тестового сигнала.

Анализ прототипа фильтра нижних частот (ФНЧ)

Для создания прототипа КИХ-фильтра и дальнейшего сравнения воспользуемся подключаемой библиотекой DSP.jl:

using DSP

Прототип цифрового фильтра нижних частот (ФНЧ) синтезируется оконным методом. Он будет работать на частоте дискретизации в 2 кГц, иметь границу полосы пропускания в 600 Гц и порядок 16. Рассчитаем и отобразим амплитудно-частотную характеристику (АЧХ) прототипа фильтра, воспользовавшись функцией freqresp:

fs = 2000;
b = digitalfilter(Lowpass(600/fs), FIRWindow(hanning(16)));
filter_prototype = PolynomialRatio(b,[1]);
H, w = freqresp(filter_prototype);
freq_vec = fs*w/(2*pi);
plot(freq_vec, pow2db.(abs.(H))*2, ylim = [-120, 10], legend = false, linewidth = 3, 
        xlabel = "Частота, Гц", 
        ylabel = "Усиление, дБ",
        title = "АЧХ прототипа ФНЧ")

В качестве тестового сигнала сгенерируем сумму двух синусоид на частотах 140 и 500 Гц длительностью в 0.5 секунд. Для сигнала с частотой дискретизации в 2 кГц количество отсчётов вектора long_signal будет ровно 1000. Отобразим результат прохождения тестового сигнала через прототип ФНЧ во временной области функцией plot:

t = 0:1/fs:0.5 - 1/fs;
long_signal = sum(sin.(2*pi*t*[140 500]), dims = 2);
out_signal = filt(filter_prototype, long_signal);
plot(t, long_signal)
plot!(t, out_signal, legend = false, xlim = [0,0.1], linewidth = 3, 
        xlabel = "Время, сек", 
        ylabel = "Амплитуда",
        title = "Прохождение тестового сигнала через прототип ФНЧ")

Мы наблюдаем успешно отфильтрованную синусоиду частотой в 140 Гц без разрывов фазы.

Попытка осуществить потоковую фильтрацию в цикле

Потоковые данные для системы ЦОС могут поступать по стандартному интерфейсу в виде пакетов информации, или непосредственно от буфера после аналого-цифрового преобразователя (АЦП). Промоделируем потоковые входные данные, представив входной сигнал как набор отдельных векторов меньшей длительности (по 100 отсчётов). Поместим их в столбцы матрицы для удобства индексации:

input_matrix = reshape(long_signal,(100,10))

100×10 Matrix{Float64}:
  0.0        -7.32312e-16  -1.46462e-15  …  -5.8585e-15  -7.05397e-14
  1.42578     1.42578       1.42578          1.42578      1.42578
  0.770513    0.770513      0.770513         0.770513     0.770513
 -0.0314168  -0.0314168    -0.0314168       -0.0314168   -0.0314168
  0.982287    0.982287      0.982287         0.982287     0.982287
  1.80902     1.80902       1.80902      …   1.80902      1.80902
  0.481754    0.481754      0.481754         0.481754     0.481754
 -0.937209   -0.937209     -0.937209        -0.937209    -0.937209
 -0.368125   -0.368125     -0.368125        -0.368125    -0.368125
  0.271031    0.271031      0.271031         0.271031     0.271031
 -0.951057   -0.951057     -0.951057     …  -0.951057    -0.951057
 -1.99211    -1.99211      -1.99211         -1.99211     -1.99211
 -0.844328   -0.844328     -0.844328        -0.844328    -0.844328
  ⋮                                      ⋱               
  0.844328    0.844328      0.844328         0.844328     0.844328
  1.99211     1.99211       1.99211          1.99211      1.99211
  0.951057    0.951057      0.951057     …   0.951057     0.951057
 -0.271031   -0.271031     -0.271031        -0.271031    -0.271031
  0.368125    0.368125      0.368125         0.368125     0.368125
  0.937209    0.937209      0.937209         0.937209     0.937209
 -0.481754   -0.481754     -0.481754        -0.481754    -0.481754
 -1.80902    -1.80902      -1.80902      …  -1.80902     -1.80902
 -0.982287   -0.982287     -0.982287        -0.982287    -0.982287
  0.0314168   0.0314168     0.0314168        0.0314168    0.0314168
 -0.770513   -0.770513     -0.770513        -0.770513    -0.770513
 -1.42578    -1.42578      -1.42578         -1.42578     -1.42578

Теперь посмотрим, что будет, если вызывать функцию filt с прототипом фильтра в цикле, применяя его последовательно для каждого из столбцов входной матрицы. Промежуточный результат фильтрации мы так же будем записывать в переменную output_matrix, совпадающую по размерам с входной матрицей. А визуальное сравнение будем проводить на "развёрнутом в линию" выходном векторе out_vector:

output_matrix = zeros(size(input_matrix));
for i = 1:10 
    output_matrix[:,i] = filt(filter_prototype, input_matrix[:,i])
end
out_vector = reshape(output_matrix,(1000));

plot(t, long_signal)
plot!(t, out_vector, legend = false, xlim = [0,0.1], linewidth = 3, 
        xlabel = "Время, сек", 
        ylabel = "Амплитуда",
        title = "Потоковая обработка прототипом ФНЧ")

Как мы видим, форма выходного сигнала оказалась искажена. Это связано с тем, что функция filt и сам объект фильтра filter_prototype не предусматривают хранение промежуточных состояний цифрового фильтра, что при обработке в цикле приводит к характерным нежелательным разрывам фазы.

Создание, инициализация и применение системного объекта фильтра

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

Создадим системный объект дискретного КИХ-фильтра (из состава библиотеки EngeeDSP):

fir_SO = EngeeDSP.DescretFIRFilter()

DescretFIRFilter:
    CoefSource=Dialog parameters
    FilterStructure=Direct form
    Coefficients=[0.5 0.5]
    InputProcessing=Elemets as channels
    InitialStates=0
    ShowEnablePort=false
    ExternalReset=None
    PreEmphasis_states=[0.0]

Оставим все его параметры по-умолчанию, за исключением коэффициентов фильтра. Их мы возьмём непосредственно из прототипа фильтра, то есть из вектора b:

fir_SO.Coefficients = b;
plot(fir_SO.Coefficients, line=:stem, marker=:circle, linewidth = 3, 
        xlabel = "Номер отсчёта", 
        ylabel = "Импульсная характеристика")

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

output_matrix = zeros(size(input_matrix));
EngeeDSP.setup!(fir_SO, input_matrix[1,:])
for i = 1:10 
    output_matrix[:,i] = EngeeDSP.step!(fir_SO, input_matrix[:,i])
end
out_vector = reshape(output_matrix,(1000));

plot(t, sig)
plot!(t, out_vector, legend = false, xlim = [0,0.1], linewidth = 3,
        xlabel = "Время, сек", 
        ylabel = "Амплитуда",
        title = "Потоковая обработка системным объектом")       

Сравним вектора программно и отобразим максимальное расхождение между выходом прототипа фильтра на длинном входном сигнале и выходом цикла потоковой фильтрации с применением системного объекта:

maximum(out_vector - out_signal)

3.3306690738754696e-16

Заключение

В примере мы рассмотрели основы применения системного объекта библиотеки EngeeDSP для задачи потоковой обработки входного сигнала в цикле. Также можно познакомиться с примерами применения системных объектов библиотеки EngeePhased для моделирования радиолокационной станции (РЛС), и с объектами библиотеки EngeeComms для проектирования систем беспроводной связи.