КИХ фильтр на Verilog

title img

В этой статье я хочу рассказать о своих экспериментах по созданию простого параметрического цифрового КИХ фильтра на Verilog HDL. До сегодняшнего дня я старался избегать тем цифровой обработки сигналов на нашем сайте: все таки это довольно сложно. Ну а когда сложно, то приходится много времени тратить просто на самостоятельное изучение вопроса. Вот сейчас только неделю читал сайт http://www.dsplib.ru/content/filters/ch10/ch10.html - очень интересный ресурс, где все довольно лаконично и понятно. После прочтения нескольких статей с сайта dsplib.ru я понял, что теорию мне лучше не писать. А вот свое понимание и свою реализацию фильтра на Verilog я пожалуй смогу предложить.

Итак, цифровые фильтры - это устройства цифровой обработки сигналов.

Цифровой сигнал - это последовательность чисел, которые представляют собой результат измерения какой-то физической величины через одинаковые интервалы времени. Можно измерять напряжение, ток, освещенность, температуру, угол наклона, скорость и прочее.

Однако, есть нюансы.

Во-первых, говорят, что измерения нужно делать с частотой более чем в 2 раза выше, чем самая высокая частота присутствующая в спектре сигнала. Об этом говорит теорема Котельникова. Практическое следствие из этого - перед измерением сигнала нужно убедиться, что в его спектре нет гармонических составляющих с частотами выше половины частоты измерений. Значит перед измерителем требуется установить соответствующий аналоговый фильтр. Если этого не сделать, то в спектре оцифрованного сигнала могут появиться частоты, которых там раньше не было, исходный сигнал будет искажен. Это вредное явление называется "алиасинг".

Во-вторых, важна разрядность измерительного прибора, АЦП, Аналогово-Цифрового Преобразователя. Вот тут, к сожалению, я не смогу дать никаких рекомендаций. Почему-то во многих статьях, которые я прочитал, говорится о выборе частоты дискретизации по теореме Котельникова, но не говорится о разрядности АЦП. Это немного странно, так как сама теорема Котельникова подразумевает взятые отсчеты с абсолютной точностью, чего в реальной жизни быть не может. Еще нужно посмотреть на уровень шумов во входном сигнале, может вообще оказаться, что младшие разряды АЦП "дрожат". Вероятно, нет смысла брать АЦП с разрядностью допускающей "дрожание" нескольких младших битов из-за шумов во входном сигнале.

Так или иначе, после того, как входной сигнал оцифрован и представляет из себя поток чисел, его можно обрабатывать численными методами в микроконтроллерах или в ПЛИС.

Одна из задач цифровой обработки сигналов - фильтрация. Фильтрация обозначает подавление одних частот в спектре и пропускание/усиление других. Про спектр я писал.

Существуют фильтры низких частот (low-pass filter), которые пропускают только частоты ниже заданной частоты среза. Есть фильтры высоких частот (high-pass), которые наоборот, подавляют частоты ниже частоты среза. Бывают полосовые фильтры, которые пропускают или подавляют только частоты в заданной полосе (band-pass & band-stop filters).

По внутреннему устройству, как мне кажется, самые простые - это КИХ фильтры. КИХ фильтр - это фильтр с конечной импульсной характеристикой. В англоязычной литературе КИХ - это FIR, то есть Finite Impulse Response filter.

Структура этого фильтра вот такая:

fir filter

Здесь блоки Z-1 - это линии задержки. На практике цепочка блоков Z-1 - это просто последовательность регистров для хранения данных. По каждому фронту тактовой частоты, которая и является частотой дискретизации, входные данные движутся от первого регистра к последующим. Фильтр может содержать цепочку регистров произвольной длины. Дополнительно, от каждого блока Z-1 идет передача данных на умножитель. Каждый из умножителей умножает свои собственные данные на некоторую константу. Потом, все результаты умножения складываются и вот получается отфильтрованный сигнал. Особая наука - это рассчитать длину фильтра и коэффициенты для умножителей. От коэффициентов зависит характеристика фильтра.

КИХ фильтр имеет несколько важных особенностей:

  • его поведение устойчиво и предсказуемо;
  • при желании можно реализовать очень быстродействующий фильтр, так как в его алгоритме данные движутся только вперед. Довольно легко сделать конвейерную обработку данных (pipeline).

Поясню, как я понимаю оба этих пункта.
Про устойчивость и предсказуемость. Если посмотреть внимательно на структуру фильтра, то становится понятно движение данных в нем. Если входные данные 16-ти битные и коэффициенты для умножителей 16-ти битные, то после умножения результат становится 32-х битный. Далее идет суммирование. Сложение двух 32-х битных чисел дает 33-х битное число. Зная длину фильтра можно точно сказать, какой разрядности будет результат. Фильтр длиной 4 будет на выходе иметь 34-х битное число и точка. Математические операции внутри фильтра не накапливают ошибку, не приводят к переполнению результата. Некоторые другие фильтры, такие как БИХ (фильтр с Бесконечно Импульсной Характеристикой, IIR - Infinite Impulse Response filter) содержат цепи обратной связи. например, вот такие:

iir filter

Вот тут все иначе. Здесь есть сумматор, который складывает входной сигнал с задержанным и умноженным на некоторый коэффициент выходной сигнал. У сумматора фиксированная разрядность. Понятно, что такой вычислитель может легко переполниться. Вот подать на такую цепь сигнал с постоянной составляющей, то через некоторое время гарантированно произойдет переполнение аккумулятора результата. Вот уже получается непредсказуемость и неустойчивость. Применять БИХ фильтры нужно с осторожностью, зараннее зная характер входного сигнала.

Теперь о быстродействии КИХ фильтра. Поскольку данные движутся только вперед, то довольно просто использовать преимущества конвейерной обработки. Например, пока умножители умножают очередные данные из задерживающих регистров, сумматор может вычислять сумму от предыдущих результатов умножителей. И сами умножители могут использовать внутренний pipeline и сумматор может использовать внутренний pipeline. Да, результат из фильтра получится задержанным на несколько тактов, но только и всего. Если же посмотреть на звено БИХ фильтра, как на рисунке выше, то понятно, что из-за обратной связи складывать нельзя, пока умножитель не вычислит свой результат.

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

Теперь расскажу, как рассчитать коэффициенты КИХ фильтра и его длину. Точнее, расскажу, как не рассчитывать - для ленивых, как я. Есть несколько online ресурсов, которые помогут рассчитать фильтры. Вот некоторые из них:

  1. https://www-users.cs.york.ac.uk/~fisher/mkfilter
  2. http://www.micromodeler.com/dsp/
  3. http://t-filter.engineerjs.com/

Воспользуемся инструментом TFilter.
Например, я хочу построить цифровой фильтр нижних частот. Частота дискретизации 20МГц. Хочу на частотах выше 2х МГц значительно подавить все частоты, не менее -40Дб.

TFilter lowpass

Ввожу требуемые мне данные в полях ввода внизу страницы и нажимаю кнопку "Design FIlter". Через пару секунд мой фильтр рассчитан и результат справа в окошке, как список коэффициентов для умножителей: -510, -520, -625, -575, -287, 306, 1232, 2467, 3927, 5477, 6948, 8162, 8962, 9241, 8962, 8162, 6948, 5477, 3927, 2467, 1232, 306, -287, -575, -625, -520, -510. Всего коэффициентов 27 - значит вот такая длина фильтра и будет.

Например, я хочу построить полосовой фильтр. Частота дискретизации 20МГц. Пусть фильтром нужно подавлять частоты ниже 1МГц и выше 4х МГц.

TFilter bandpass

Ввожу данные, делаю рассчет, получаю коэффициенты фильтра: -801, -1026, -210, 1914, 4029, 3905, 330, -5174, -8760, -7040, -152, 7700, 11130, 7700, -152, -7040, -8760, -5174, 330, 3905, 4029, 1914, -210, -1026, -801. Здесь длина фильтра получилась 25 тапов.

Вот коэффициенты фильтров есть. Что с ними делать? Хорошо бы сделать параметрический модуль на Verilog HDL, чтобы его можно было легко настраивать. По логике, параметры модуля фильтра - это его длина, то есть количество линий задержки, разрядность входных данных и разрядность коэффициентов для умножителей.

Если разрядность входных данных - это IWIDTH, а число регистров в цепочке линии задержки - это TAPS, то довольно просто описать саму цепочку регистров:


module fir ( clk, coefs, in, out );

parameter IWIDTH = 16; //input data (signal) width
parameter TAPS = 2; //number of filter taps

input wire clk;
input wire [IWIDTH-1:0]in;

genvar i;
generate
  for( i=0; i<TAPS; i=i+1 )
  begin:tap
    //make tap register chain
    reg [IWIDTH-1:0]r=0;
    if(i==0)
    begin
      //1st tap takes signal from input
      always @(posedge clk)
        r <= in;
    end
    else
    begin
      //tap reg takes signal from prev tap reg
      always @(posedge clk)
        tap[i].r <= tap[i-1].r;
    end
end
endgenerate
endmodule


Конструкция generate-endgenerate языка Verilog позволяет динамически создать нужное число регистров и посоединять их как надо в цепочку. Получится вот так:

regs

Таким образом, длина фильтра может быть задана через параметр модуля Verilog. Аналогичным образом, внутри цикла for() внутри generate-endgenerate создаются умножители:


genvar i;
generate
  for( i=0; i<TAPS; i=i+1 )
  begin:tap
    ........
    //get tap multiplication constant coef
    wire [CWIDTH-1:0]c;
    assign c = coefs[((TAPS-1-i)*32+CWIDTH-1):(TAPS-1-i)*32];

    //calculate multiplication and fix result in register
    reg [MWIDTH-1:0]m;
    always @(posedge clk)
      m <= $signed(r) * $signed( c );
    ............


Таким образом, дополнительно к регистрам линии задержки будут сформированы и подключены умножители со своими коэффициентами:

schema

Я решил, что коэффициенты для умножителей можно передавать в модуль объединенные в шину {...} как входной сигнал, вот так:


fir #( .TAPS(27) ) fir_lp_inst(
  .clk(tb_clk),
  .coefs( {
    -32'd510,
    -32'd520,
    .......
     32'd575,
     32'd625,
    -32'd520,
    -32'd510
    } ),
  .in(),
  .out()
);


Здесь каждый коэффициент фильтра занимает в шине фиксированное число бит 32 (решил, что этого должно хватить). Но из этих 32х бит модуль выберет нужное объявленное число бит для коэффициентов. Например, parameter CWIDTH = 16; Напомню, что от количества бит в коэффициенте для умножения зависит в конечном итоге ширина результата.

Для простоты реализации я пока не стал делать pipeline для сумматоров. Их можно было бы сделать побыстрее, например, складывать значения из умножителей попарно, потом складывать попарно результаты предыдущих парных сложений и так далее.. Пока я сделал один большой сумматор, как одну большую и сложную комбинаторную функцию. Наверняка в ПЛИС это не сможет работать на высоких частотах, но пока я на это закрываю глаза. Сделать бы чтобы просто работало.

Полный текст модуля будет выглядеть вот так:

module fir ( clk, coefs, in, out );

parameter IWIDTH = 16;	//input data (signal) width
parameter CWIDTH = 16;	//tap coef data width (should be less then 32 bit)
parameter TAPS   = 2;	//number of filter taps
localparam MWIDTH = (IWIDTH+CWIDTH); //multiplied width
localparam RWIDTH = (MWIDTH+TAPS-1); //filter result width

input  wire clk;
input  wire [IWIDTH-1:0]in;
input  wire [TAPS*32-1:0]coefs; //all input coefficient concatineted
output wire [RWIDTH-1:0]out; //output takes only top bits part of result

genvar i;
generate
	for( i=0; i<TAPS; i=i+1 )
	begin:tap
		//make tap register chain
		reg [IWIDTH-1:0]r=0;
		if(i==0)
		begin
			//1st tap takes signal from input
			always @(posedge clk)
				r <= in;
		end
		else
		begin
			//tap reg takes signal from prev tap reg
			always @(posedge clk)
				tap[i].r <= tap[i-1].r;
		end

		//get tap multiplication constant coef
		wire [CWIDTH-1:0]c;
		assign c = coefs[((TAPS-1-i)*32+CWIDTH-1):(TAPS-1-i)*32];

		//calculate multiplication and fix result in register
		reg [MWIDTH-1:0]m;
		always @(posedge clk)
			m <= $signed(r) * $signed( c );
			
		//make combinatorial adders
		reg [MWIDTH-1+i:0]a;
		if(i==0)
		begin
			always @*
				tap[i].a = $signed(tap[i].m);
		end
		else
		begin
			always @*
				tap[i].a = $signed(tap[i].m)+$signed(tap[i-1].a);
		end
	end
endgenerate

//fix calculated taps summa in register
reg [RWIDTH-1:0]result;
always @(posedge clk)
	result <= tap[TAPS-1].a;

//deliver output
assign out = result;

endmodule

Если будет время, то постараюсь улучшить этот модуль, добавить pipeline сумматоров.

Теперь главный вопрос - а как убедиться, что это вообще работает? Ну как? да просто! - есть симулятор Verilog. Напишу тестбенч, который будет иммитировать входной сигнал на разных частотах и посмотрю отклик на выходе фильтра. Даже более того. Если тестбенч будет синтезировать входной сигнал как синусоиду и если плавно менять ее частоту в некотором диапазоне, от низкой частоты к более высокой, то наблюдая отклик на выходе фильтра я прямо увижу его амплитудно-частотную характеристику!

Структура программы testbench.v:

verilog testbench

Синусоиду для тестбенча я уже когда-то показывал. Возьму тот старый код, немного подправлю его и сделаю плавное изменение частоты от 100кГц до 4МГц с шагом в 1000Гц.

Внутрь тестбенча установлю два экземпляра модуля своего fir() фильтра, но с разными коэффициентами, для фильтра нижних частот и для полосового фильтра.

Весь код тестбенча вот здесь:

`timescale 1ns / 1ns

module testbench();

reg tb_clk;
initial tb_clk=0;
always
	#25 tb_clk = ~tb_clk;

real PI=3.14159265358979323846;
real last_time=0; //Sec
real current_time=0; //Sec
real angle=0;	//Rad
real frequency=100; //Hz
integer freq_x100kHz=0; //*100kHz
reg signed [15:0]sin16;

//function which calculates Sinus(x)
function real sin;
input x;
real x;
real x1,y,y2,y3,y5,y7,sum,sign;
 begin
  sign = 1.0;
  x1 = x;
  if (x1<0)
  begin
   x1 = -x1;
   sign = -1.0;
  end
  while (x1 > PI/2.0)
  begin
   x1 = x1 - PI;
   sign = -1.0*sign;
  end  
  y = x1*2/PI;
  y2 = y*y;
  y3 = y*y2;
  y5 = y3*y2;
  y7 = y5*y2;
  sum = 1.570794*y - 0.645962*y3 +
      0.079692*y5 - 0.004681712*y7;
  sin = sign*sum;
 end
endfunction

task set_freq;
input f;
real f;
begin
	frequency = f;
	freq_x100kHz = f/100000.0;
end
endtask

always @(posedge tb_clk)
begin
	current_time = $realtime;
	angle = angle+(current_time-last_time)*2*PI*frequency/1000000000.0;
	//$display("%f %f",current_time,angle);
	while ( angle > PI*2.0 )
	begin
		angle = angle-PI*2.0;
	end 
	sin16 = 32000*sin(angle);
	last_time = current_time;
end

//low-pass filter
wire [57:0]out_lowpass;
fir #( .TAPS(27) ) fir_lp_inst(
	.clk(tb_clk),
	.coefs( { 
		-32'd510,
		-32'd520,
		-32'd625,
		-32'd575,
		-32'd287,
		 32'd306,
		 32'd1232,
		 32'd2467,
		 32'd3927,
		 32'd5477,
		 32'd6948,
		 32'd8162,
		 32'd8962,
		 32'd9241,
		 32'd8962,
		 32'd8162,
		 32'd6948,
		 32'd5477,
		 32'd3927,
		 32'd2467,
		 32'd1232,
		 32'd306,
		-32'd287,
		-32'd575,
		-32'd625,
		-32'd520,
		-32'd510
		} ),
	.in(sin16),
	.out(out_lowpass)
	);

//band-pass
wire [55:0]out_bandpass;
fir #( .TAPS(25) ) fir_bp_inst(
	.clk(tb_clk),
	.coefs( { 
		-32'd801,
		-32'd1026,
		-32'd210,
		 32'd1914,
		 32'd4029,
		 32'd3905,
		 32'd330,
		-32'd5174,
		-32'd8760,
		-32'd7040,
		-32'd152,
		 32'd7700,
		 32'd11130,
		 32'd7700,
		-32'd152,
		-32'd7040,
		-32'd8760,
		-32'd5174,
		 32'd330,
		 32'd3905,
		 32'd4029,
		 32'd1914,
		-32'd210,
		-32'd1026,
		-32'd801
		} ),
	.in(sin16),
	.out(out_bandpass)
	);

integer i;
real f;

initial
begin
	$dumpfile("out.vcd");
	$dumpvars(0,testbench);
	f=100000;
	for(i=0; i<4000; i=i+1)
	begin
		set_freq(f);
		#1000;
		f=f+1000;
	end
	$finish;
end

endmodule

Симуляцию провожу в icarus verilog.
Компилирую:
 > iverilog -o qqq testbench.v fir.v
Симулирую:
 > vvp qqq
Смотрю временные диаграммы в GtkWave:
 > gtkwave out.vcd

В результате симуляции программой vvp получается выходной файл с временными диаграммами сигналов out.vcd. Их можно посмотреть в программе GtkWave. Вот что видно: входной сигнал sin16 для обоих фильтров меняет частоту, она растет:

wave freq grow

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

gtkwave FIR filter

Входной сигнал для sin16 равномерен по амплитуде во всей полосе частот. "Рваные" сигналы в области низких частот - это просто GtkWave так отрисовывает, при увеличении (zoom) там все нормально выглядит.

А вот выходные сигналы, такие какие и положены быть на выходе фильтров. Фильтр низких частот от 500кГц начинается плавный спад амлитуды синусоиды out_lowpass и после 2х МГц фильтр уже практически не пропускает. Полосовой фильтр хорошо пропускает где-то в полосе 2-3МГц, а по краям пологие спады (сигнал out_bandpass).

Чтобы внимательно рассмотреть получившуюся картину нужно использовать специальные возможности просмотра в программа GtkWave. Вот несколько полезных советов:

  • для знаковых сигналов, как sin16 или out_lowpass, out_bandpass кликнуть правой кнопкой мыши по сигналу и выбрать сперва Data Format => Sogned Decimal, а потом Data Format => Analog => Step;
  • чтобы отобразить аналоговый сигнал на графике в большую высоту используйте правый клик мыши на сигнале и потом пункт popup меню Insert Analog Height Extension.

Еще совет: GtkWave может отображать аналоговый сигнал отмасштабированный по высоте с учетом всех выборок симуляции по времени или с учетом только выборок, которые помещаются в окне. Если поставить Data Format => Analog => Resizing => Screen Data, то есть если масштабирование по вертикали будет учитывать только помещающиеся в окне выборки, то можно увидеть боковые лепестки амплитудно-частотной характеристики фильтра:

gtkwave2

Таким образом, думаю мне удалось создать параметрический КИХ (FIR) фильтр и получилось посмотреть, как он работает.

Справедливости ради нужно сказать, что среда проектирования Altera Quartus II и Quartus Prime уже имеют встроенные мегафункции для FIR фильтров. Создать фильтр в среде Quartus по имеющимся коэффициентам довольно просто. Однако, имеющиеся мегафункции квартуса доступны в версиях Web Edition и Lite только для ознакомления. Они работают ограниченное время. Требуется приобретение полной версии САПР, чтобы применять FIR в своих проектах без ограничений.

 

Добавить комментарий