Моделирование работы КИХ фильтра в Matlab

Моделирование в Matlab часто снимает много непоняток при проектировании устройств цифровой обработки сигналов на ПЛИС. Данный пост это простейшая модель работы цифрового КИХ фильтра на ПЛИС, которая дает ответы на вопросы связанные с масштабом выходных данных  КИХ фильтра, подавлением частот вне его полосы.

Предположим что в ПЛИС необходимо создать или вернее написать на VHDL Verilog КИХ фильтр с известной полосой и коэффициентами. В этом нет ничего особенного, но моделирование работы фильтра в Matlab может прояснить некорые узкие моменты. 

Итак КИХ фильтр полоса 85-95 кГц, 51 коэффициент 16 бит со знаком. В матлаб модель такого фильтра выглядит достаточно просто:

%% инициализация

clc; 
clear all;
  
f1 = 90000;         % центральная частота Гц
f2 = 70000;         % центральная частота Гц
fs = 360000;        % частота дискретизации Гц
A = 2047;           % амплитуда сигнала  
ts = 1/fs;          % период одного дискрета
tm = 50;            % количество проходов участка моделирования  
tm_test = tm*4;			% общее количество проходов моделирования
td = 5;             % коэффициент децимации
w1=2*pi*f1;         % круговая частота 1
w2=2*pi*f2;         % круговая частота 2
tmin = 0;           % начальное время (начальная фаза)
tmax = tm*ts;       % конечное время
tmax_test = tm_test*ts;       % конечное время

% матрица коэффициентов 51 коэффициент (Equiripple 75-85-95-105 кГц)
kf = [0,-834,0,283,0,-202,0,8,0,307,0,-738,0,1260,0,-1836,0,2417, ...
        0,-2947,0,3372,0,-3647,0,3742,0,-3647,0,3372,0,-2947,0, ...   
				2417,0,-1836,0,1260,0,-738,0,307,0,8,0,-202,0,283,0,-834,0];
			
%lkf =length(kf);    % количество коэффициентов
% матрица времени 
t=linspace(tmin, tmax_test+tmin, fs*tmax_test+1);
% матрица сигнала 90 кГц
x1 = A*sin(w1*t);
x1(1:60) = 0;
x1(140:201) = 0;
% матрица сигнала 70 кГц
x2 = A*sin(w2*t);
x2(1:60) = 0;
x2(140:201) = 0;

%% моделирование

zr1=zeros(1,200);
zr2=zeros(1,200);
for tm_mod = 1:tm_test
  
  % поэлементное умножение
	if tm_mod < 150
		z1 = x1(tm_mod:(tm_mod+50)).*kf;
		z2 = x2(tm_mod:(tm_mod+50)).*kf;
	end	
  % матрица выходных данных КИХ для 70 кГц и 90 кГц 
  zr1(tm_mod) = round(sum(z1));
	zr2(tm_mod) = round(sum(z2));
  
end

% максимальное значение
fprintf('Максимальное значение 90 кГц %d\n', max(zr1));
fprintf('Максимальное значение 70 кГц %d\n', max(zr2));
fprintf('\n');

%% графика

% выборки по сигналу x1
figure(1)
plot(x1, 'r');
hold on
stem(x1);

% графика элементов массива zr1
figure(2)
plot(zr1, 'r'), grid;
hold on
stem(zr1);

% графика элементов массива zr2
figure(3)
plot(zr2, 'r');
hold on
stem(zr2);

Коротко по коду: в начале задаются частоты сигналов, первая f1 - центральная частота фильтра, вторая f2 - частота вне полосы фильтра. Далее задаются другие необходимые атрибуты сигналов, амплитуда, время существования сигналов, частота дискретизации. Коэффициенты фильтра задаются в виде матрицы. Также в виде матриц задаются сигналы.

Моделирование представляет из себя поэлементное перемножение элементов матриц сигналов и коэффицентов с суммированием результатов. Ну и простейшая графика которая визуализирует процесс работы цифрового фильтра во времени. 

На картинке выше модель сигнала, а на картинках ниже графики выходного результата работы фильтра при частоте сигнала  f1 и f2 соответственно.

Явно прослеживается полное пропускание сигнала с частотой 90 кГц и подавление сигнала с частотой 70 кГц. Наиболее существенное подавление происходит при полном 'наползании' матрицы сигнала на матрицу коэффицентов. Тоесть при данной конкретной частоте дискретизаци, коэффициентах и частоте сигнала соответствующей центральной частоте  КИХ фильтра можно сказать, что при перемножении матрицы сигнала на матрицу коэффициентов все данные имеют положительный знак. Естественно их сумма будет максимальной. При частоте сигнала вне полосы фильтра данные перемножения матриц сигнала и коэффициентов будут иметь разные знаки. При суммировании их сумма будет существенно меньше предыдущего случая.