Реализация на Verilog цифрового БИХ-фильтра

Приветствую Хабр. Не так давно здесь уже появлялись статьи на эту тему Verilog. Цифровой фильтр на RAM и Построение цифрового фильтра с конечной импульсной характеристикой. Хочу и я внести свой скромный вклад и представить вашему вниманию реализацию цифрового БИХ-фильтра на Verilog.

Не смотря на то, что занимающиеся данной темой всё из ниже описанного знают, хочу, для начала, все-таки дать совсем немного теории для широкого круга читателей.
БИХ-фильтр(IIR-фильтр) — фильтр с бесконечной импульсной характеристикой (infinite impulse response). Отличительной чертой этого типа фильтра от КИХ-фильтра (фильтр с конечной импульсной характеристикой) является наличие обратной связи, когда значение на выходе фильтра зависит не только от входных данных, но и от выходных данных, полученных фильтром на предыдущих итерациях. Структуру БИХ-фильтра можно представить следующим образом:

Структуру БИХ-фильтра
Не буду останавливаться на методиках расчета коэффициентов аналогового фильтра, с этой темой можно ознакомиться в большом количестве изданий, таких как Аналоговые и цифровые фильтры. Расчет и реализация Г.Лем или Синтез фильтров. Херреро Дж. Уиллонер Г. Так же существует много программ, позволяющих это сделать. Опишу кратко только процесс преобразования аналогового фильтра в цифровой.
Любой фильтр можно представить с помощью аналоговой передаточной функции. Например для БИХ-фильтра второго порядка она будет выглядеть так:
Аналоговая передаточная функция
Подставив в это выражение вместо нормированной комплексной переменной
нормированная комплексная переменная
получим передаточную функцию A(z), которая может быть реализована в цифровом фильтре:
Цифровая передаточная функция
Коэффициенты для фильтра 2-го порядка будут определяться по следующим формулам:

где
— тактовая частота
— частота среза фильтра
d0...2 и с0...2 — коэффициенты аналогового фильтра, предполагается что вы их посчитали заранее.
Воспользовавшись этими нехитрыми вычислениями и зная соответствующие коэффициенты для аналогового фильтра, можно рассчитать коэффициенты для цифровой передаточной функции фильтра любого порядка. (Меня, правда, в свое время хватило только до 4-го порядка, причем сложности никакой в этом нет, просто приходится много писать длинные формулы, упрощая по формулам сокращенного умножения. Зато вспомнил далекие школьные годы). Передо мной, например, стояла задача реализовать фильтр с перестраиваемыми характеристики, поэтому все коэффициенты рассчитывались на ПК и передавались в FPGA в виде 32-х разрядных чисел, преобразованные перед этим в формат с плавающей запятой.
Итак, имея коэффициенты цифровой передаточной функции, подставим их в разностное уравнение БИХ-фильтра:
разностное уравнение БИХ-фильтра
на основе которого уже не сложно написать реализацию на Verilog.

Текст модуля
module LP_FILTER
(
mhz_clk,RESET, 
D0,D1,D2,C0,C1,
X,Y,
COUNT

);

	// low pass filter INPUT32 OUTPUT32
/*
	y(i) = D2 * x(i) + D1 * x(i-1) + D0 * x(i-2) + C1 * y(i-1) + C0 * y(i-2)
	
	A(z) = (D0+D1*z+D2*z*z)/(C0+C1*z+z*z) <==> A(P) = (d0+d1*P+d2*P*P)/(c0+c1*P+c2*P*P)
											   A(P) = 1/(1+1.4*P+P*P) LPF 	

	D0 =  (  d0 - d1*l +   d2*l*l)/(c0 + c1*l + c2*l*l)
	D1 =  (2*d0        - 2*d2*l*l)/(c0 + c1*l + c2*l*l)
	D2 =  (  d0 + d1*l +   d2*l*l)/(c0 + c1*l + c2*l*l)
	C0 = -(  c0 - c1*l +   c2*l*l)/(c0 + c1*l + c2*l*l)	
	C1 = -(2*c0        - 2*c2*l*l)/(c0 + c1*l + c2*l*l)	
		  
	l = ctg ( 3.14 * f_filt / f_samp )

*/

input  mhz_clk;
input RESET;
input [17:0]	X;
input [31:0]	D0;
input [31:0]	D1;
input [31:0]	D2;
input [31:0]	C0;
input [31:0]	C1;
output [17:0]	Y;
output [5:0]	COUNT;

// adc-filter counter	

reg [5:0] COUNT;
always @( posedge mhz_clk,negedge RESET )
	if (~RESET) COUNT[5:0] = 0; //49 count
	else if (COUNT[5:0] == 6'h31) COUNT[5:0] = 0;
	else COUNT[5:0] = COUNT[5:0] + 6'b1;

//  input - COUNT[4:0] = 24:49
// output - COUNT[4:0] = 6

reg [31:0]	c0_y2;
reg [31:0]	c1_y1;
reg [31:0]	d0_x2;
reg [31:0]	d1_x1;
reg [31:0]	d2_x0;

reg [31:0]	y0;
reg [31:0]	y1;
reg [31:0]	y2;
reg [31:0]	x0;
reg [31:0]	x1;
reg [31:0]	x2;

//	y(i) = D2 * x(i) + D1 * x(i-1) + D0 * x(i-2) + C1 * y(i-1) + C0 * y(i-2)
reg [31:0]	mul_a;
reg [31:0]	mul_b;
always @(*)
	if 		( ( COUNT[5:0] >=  0 ) & ( COUNT[5:0] <=  4 ) | ( COUNT[5:0] ==  49 ) ) mul_a[31:0] = x0[31:0];
	else if ( ( COUNT[5:0] >=  5 ) & ( COUNT[5:0] <= 10 ) ) mul_a[31:0] = x1[31:0];   
	else if ( ( COUNT[5:0] >= 11 ) & ( COUNT[5:0] <= 16 ) ) mul_a[31:0] = x2[31:0];
	else if ( ( COUNT[5:0] >= 17 ) & ( COUNT[5:0] <= 22 ) ) mul_a[31:0] = y1[31:0];
	else    												mul_a[31:0] = y2[31:0]; 
	  
always @(*)
	if 		( ( COUNT[5:0] >=  0 ) & ( COUNT[5:0] <=  4 ) | ( COUNT[5:0] ==  49 ) ) mul_b[31:0] = D2[31:0];
	else if ( ( COUNT[5:0] >=  5 ) & ( COUNT[5:0] <= 10 ) ) mul_b[31:0] = D1[31:0];   
	else if ( ( COUNT[5:0] >= 11 ) & ( COUNT[5:0] <= 16 ) ) mul_b[31:0] = D0[31:0];
	else if ( ( COUNT[5:0] >= 17 ) & ( COUNT[5:0] <= 22 ) ) mul_b[31:0] = C1[31:0];
	else    												mul_b[31:0] = C0[31:0]; 

wire [31:0] mul_out;	
mul_float32 ( 1, mhz_clk, mul_a[31:0], mul_b[31:0], mul_out[31:0] );

reg [31:0] outmul;
always @(*) outmul[31:0]=mul_out[31:0]; 

always @( posedge mhz_clk,negedge RESET )
	if (~RESET)	d2_x0[31:0] = 32'h0;
	else if ( COUNT[5:0] ==   4 ) d2_x0[31:0] = mul_out[31:0];
							
always @( posedge mhz_clk,negedge RESET )	
	if (~RESET) d1_x1[31:0] = 32'h0;
	else if ( COUNT[5:0] ==  10 ) d1_x1[31:0] = mul_out[31:0]; 
							
always @( posedge mhz_clk,negedge RESET )
	if (~RESET) d0_x2[31:0] = 32'h0;
	else if ( COUNT[5:0] ==  16 ) d0_x2[31:0] = mul_out[31:0]; 
							
always @( posedge mhz_clk,negedge RESET )	
	if (~RESET) c1_y1[31:0] = 32'h0;
	else if ( COUNT[5:0] ==  22 ) c1_y1[31:0] = mul_out[31:0];
							
always @( posedge mhz_clk,negedge RESET )
	if (~RESET)	c0_y2[31:0] = 32'h0;
	else if ( COUNT[5:0] ==  28 ) c0_y2[31:0] = mul_out[31:0]; 
							
//	y(i) = D2 * x(i) + D1 * x(i-1) + D0 * x(i-2) + C1 * y(i-1) + C0 * y(i-2)
reg [31:0]	sum_a;
reg [31:0]	sum_b;

always @(*)
	if 		( ( COUNT[5:0] >=  0 ) & ( COUNT[5:0] <= 18 ) | ( COUNT[5:0] ==  49 ) ) sum_a[31:0] = d2_x0[31:0];
	else if ( ( COUNT[5:0] >= 19 ) & ( COUNT[5:0] <= 26 ) ) sum_a[31:0] = d0_x2[31:0];   
	else if ( ( COUNT[5:0] >= 27 ) & ( COUNT[5:0] <= 34 ) ) sum_a[31:0] = c1_y1[31:0];
	else    												sum_a[31:0] = c0_y2[31:0]; 

always @(*)
	if 		( ( COUNT[5:0] >=  0 ) & ( COUNT[5:0] <= 18 ) | ( COUNT[5:0] ==  49 ) ) sum_b[31:0] = d1_x1[31:0];
	else 	sum_b[31:0] = y0[31:0];  		

wire [31:0] sum_out;
sum_float32 (	1,	mhz_clk,	sum_a[31:0],	sum_b[31:0],	sum_out[31:0] );

reg [31:0] outsum;
always @(*) outsum[31:0]=sum_out[31:0]; 
			
always @( posedge mhz_clk,negedge RESET )
	if (~RESET) y0[31:0] = 32'h0;
	else if ( ( COUNT[5:0] ==  18 ) | ( COUNT[5:0] ==  26 ) | ( COUNT[5:0] ==  34 ) | ( COUNT[5:0] ==  42 ) ) y0[31:0] = sum_out[31:0];
							

reg [17:0] int_to_float_in;
always @(*) begin int_to_float_in[17:0] = X[17:0];end
wire [31:0] int_to_float_out;
int18_to_float32 (	1,	mhz_clk, int_to_float_in[17:0], int_to_float_out[31:0] );

always @( posedge mhz_clk,negedge RESET )
	if (~RESET)	x0[31:0] = 32'h0;
	else if ( COUNT[5:0] ==  49 ) x0[31:0] = int_to_float_out[31:0];
		
always @( posedge mhz_clk,negedge RESET )
	if (~RESET)	x1[31:0] = 32'h0;
	else if ( COUNT[5:0] ==  49 ) x1[31:0] = x0[31:0];
	
always @( posedge mhz_clk,negedge RESET )
	if (~RESET) x2[31:0] = 32'h0;
	else if ( COUNT[5:0] ==  49 ) x2[31:0] = x1[31:0];
	
always @( posedge mhz_clk,negedge RESET )
	if (~RESET) y1[31:0] = 32'h0;
	else if ( COUNT[5:0] ==  49 ) y1[31:0] = sum_out[31:0];	

always @( posedge mhz_clk,negedge RESET )
	if (~RESET) y2[31:0] = 32'h0;
	else if ( COUNT[5:0] ==  49 ) y2[31:0] = y1[31:0];

wire [17:0] float_to_int_out;
wire nan;
wire overflow;
wire underflow;
float32_to_int18 ( 1, mhz_clk, y0[31:0], nan, overflow, float_to_int_out[17:0], underflow ); 

reg [17:0] Y;
always @( posedge mhz_clk )	if ( COUNT[4:0] ==  6 ) Y[17:0] = float_to_int_out[17:0];	
		
endmodule



Программа представляет собой законченный модуль, готовый к включению в существующий проект. Входные данные передаются в него в виде 18-ти разрядных целых чисел. (Возможно кому-то покажется странным такая разрядность, поясню сразу, что на разработанной плате используется 18-ти разрядная АЦП AD7643). Далее они преобразовываются в 32-х разрядный формат с плавающей точкой. Это выполняется по-средствам стандартной альтеровской мегафункции int18_to_float32. Сложение и умножение реализовано так же стандартными функциями sum_float32 и mul_float32 соответственно. После выполнения необходимых вычислений результат конвертируется в целое 18-ти разрядной число функцией float32_to_int18. С помощью предварительно рассчитанных коэффициентов можно превратить этот модуль как в фильтр нижних частот, верхних частот либо полосовой.
Архив проекта в Quartus II 9.1
Помимо уже упомянутых в тексте книг и данных из Википедии при написании топика использовалось справочное руководство Полупроводниковая схемотехника. Титце У., Шенк К. и Структуры цифровых фильтров и их характеристики
Прошу прощения за некоторую сумбурность изложения, т.к. литературный талант никогда не являлся моей отличительной чертой, и спасибо за внимание.
Поделиться публикацией
AdBlock похитил этот баннер, но баннеры не зубы — отрастут

Подробнее
Реклама

Комментарии 17

    0
    2 порядок больше подходит для звена БИХ фильтра. Сколько варнингов у такого описания? Еще, что бросилось в глаза: если стоит условный оператор, то условия перечисляются через логическое и/или (&& ||). А вообще классно, было интересно читать!
      0
      Спасибо. Вылитает 16 варнингов, причем 2 относятся к тому что не используется объявленная переменная outsum, 4 к тому, что не дал внутримодульные названия мегафункциям. Это самые критичные вроде, остальные — мелочь какая-то. Насчет условного оператора, вы правы конечно, но я во всех проектах использую подобную запись, и, кажется, никаких проблем с этим не возникало. Может быть стоит начать писать правильно, мало ли где выплавит. А что вы имели ввиду по поводу звена БИХ фильтра? Нашел фразу:«БИХ-фильтры обычно реализуются с помощью звеньев второго порядка, которые называются биквадратными фильтрами». Но вроде бы эти биквадратные формы используются для упрощений расчетов фильтров больших порядков, определения нулей и полюсов фильтров. Влияет ли это на конечную реализацию устройства, например в Verilog? Всегда казалось что для повышения порядка, достаточно добавить дополнительные слагаемые в конечное разностное уравнение.
        0
        А посмотрите в проекте этот сумматор: COUNT[5:0] + 1. Просто интересно, скольки он разрядный. Просто «1» по дефолту int и преображается в 32 разряда, поэтому сумматор может получится страшным, а может и нет, если оптимизацию пройдет.
        Лично один раз спотыкался об & и && по части интерфейсов, теперь стараюсь следить.
        И да, большие БИХ фильтры можно строить звеньями, что облегчает анализ этих фильтров. Это эквивалентно последовательному включению фильтров меньшего порядка, что реализацию, в принципе, не меняет.
          0
          Хм… Да, вы правы, по этому счетчику выдавал одно предупреждение. Но, по-видимому, компилятор приводит его к 6-и разрядному виду сам. Внес изменения в код на всякий случай. Спасибо.
      0
      Не могу не запостить этот баян:
      image
        0
        Спасибо, учту…
          0
          Вы не принимайте близко к сердцу, я с благими намерениями :)
            0
            Я нормально отношусь к конструктивной критике и ничего обидного здесь не вижу. До вашего комментария я даже не задумывался об этом вопросе и как-то неудобно даже стало перед читателями за такой просчет, поэтому мое «спасибо» действительно искренне.
        +1
        хм…
        а это у Вас вообще-то работает? пробовали ли в железе или симулировали?
        Что-то меня смущает несколько моментов.
        1) похоже все операции делаются над положительными числами, это так? Не вижу signed reg или $signed… Если все числа положительны, то на выходе будет все возрастающее число — посмотрите на свою структуру: сложение многих положительных компонентов дает большее число и с каждой новой входящей выборкой число на выходе будет увеличиваться…
        2) вы как-то смело используете для всех переменных 32х битные регистры. При этом вообще-то при перемножении двух 32х битных чисел получится 64х битное число, ведь так? Почему Вы пренебрегаете этим?
        3) и как оно синтезируется для ПЛИС? для какой конкретно ПЛИС? операция умножения довольно расточительная функция и занимает много места, если только не используются встроенные умножители
        3)
          0
          Да, все хорошо работает в железе.
          1) Дело в том что при вычислениях используется формат числа с плавающей запятой. Я указал в пояснениях к модулю, что функция int18_to_float32 преобразует числа к этому виду.Вот ознакомьтесь пожалуйста Число одинарной точности.Старший бит этой 32-х битной структуры как раз и отвечает за знак числа. И числа конечно же не все положительные. Переменные С0… С2, заведомо умножены на -1, т.е. они отрицательные.
          2) Все правильно, если бы это были целые числа, то так бы оно и было. Но я использую числа с плавающей запятой, поэтому результат остается 32-х разрядный.
          3)Конкретно в этом проекте я использую Cyclone III EP3C25Q240C8.Там есть встроенные умножители. У меня там 4 таких фильтра и еще процессор Nios II, места на все хватает.
          0
          Вот пара моментов, которые прям бросились в глаза:
          1) Входы-выходы объявляются в заголовке модуля обычно. Вы объявляете регистры и цепи в процессе их появления в тексте. Мелочь, конечно, но для меня странно выглядит.
          2) Для каждого блока always используете один оператор. Почему не в операторных скобках в одном блоке always?
          3) У вас все повторяющиеся условия ( которые отрабатывают счетчик) можно свести в один case ( автомат получится). Вы это делаете, чтоб синтезировалась обработка условий через enable в триггерах вместо мультиплексора, которым реализуется автомат? Насколько это эффективнее?
          4) Я бы параметризовал многие вещи. Для гибкости. Вдруг пригодится.
          5) На какой частоте все работает? Имеется в виду оценка Тime quest'а.

          Извините, если несуразные вопросы. Стаж у самого чуть больше года. Может я что-то важное упускаю?
          Eщё нескромный вопрос. Не пробовали генерировать фильтр из матлаба?
            0
            1) Мне так удобнее ориентироваться в программе. Когда простыня на несколько экранов я ищу нужное мне место по объявлениям регистров.
            2) Ну тут много причин:
            a)Ключевое слово always указывает на процесс. Процесс — это такой кусок кода, внутри которого все операции выполняются ПОСЛЕДОВАТЕЛЬНО. Зачем мне ограничиваться, выстраивая последовательные расчеты, если я могу выполнять их ПАРАЛЛЕЛЬНО.
            b)Эти процессы зачастую тактируются по-разному. Если в этой маленькой программе этого нет (везде одна частота), то есть проекты, в которых условия срабатывания в разных процессах разные. Поэтому я привык разграничивать различные по смыслу процессы.
            3)Это такой же мультиплексор, не важно как его задавать, case-ом или if-ом. Реализация не вентильно-тригерном уровне не меняется. В case вроде бы нельзя задать выбор по диапазону, а только по конкретному значению, как в if, а раз смысл одинаковый, а функциональности больше у if, я пользуюсь обычно if. Хотя иногда, когда много условий, могу использовать и case.
            4)Да, можно было бы. А то какой-то говнокод получается:)
            5)По данным Тime quest'а максимальное значение частоты mhz_clk = 63.4MHz, но эта частота определяется частотой дискритизации АЦП. В моем проекте это 600кГц, поэтому mhz_clk я выбрал 30МГц.
            С матлабом пока дела не имел, но начинаю к нему приглядываться. Скоро опять будет проект, где нужно будет синтезировать фильтр по заданной АЧХ, причем весьма сложной. Без матлаба видимо не обойтись…
              0
              1) Вы работаете один? Я работаю в команде и стилистика, хоть и на словах, но имеется и обсуждается.
              2) а) Например блок
              Always @(posedge clk)
              begin
              a=x+y;
              b=q+g;
              c=n+m;
              end
              Все присвоения в таком случае сработают параллельно.
              5) Маловата частота конечно для такой плис, но если запас в 2 раза, то проблем нет.

              Matlab, кстати, очень рекомендую. Весь ЦОС всегда полностью моделирую, но вот генерацией HDL кода пока не занялся.
                0
                1) Нет, не один, но у нас все также пишут:)
                2)Нет, последовательно. Это можно сделать параллельно, если вместо begin...end использовать fork...join. Либо можно еще ускорить процесс присвоения регистра. Дело в том, что операция установки регистра занимает определенное время. Можно использовать неблокирующее присваивание(a<=x+y;b<=q+g), при котором перескок на следующую операцию идет сразу после того, как регистр начнет устанавливаться.
                0
                Я делал фильтры из fdatool в matlab, получааетсся очень громоздко и нечитабельно. Именно по этой причине я и написал этот пост, где брал из матлаба коэффициента фильтра, а все умножения с накоплением делал сам. Когда доделал фильтр сравнил результаты: дизайн из матлаба не влез в Cyclone2 не помню какую, мой дизайн занял 1% кристалла, аккуратно разложив сдвиговый регистр и коэффициенты в ram, которой было более чем достаточно

          Только полноправные пользователи могут оставлять комментарии. Войдите, пожалуйста.

          Самое читаемое