[billidean 0]

Помогите плизз.

 

Раньше с фильтрами почти не работал.

 

Необходимо написать код vhdl фильтра ФНЧ Баттерворта 6 порядка.

Частота среза = 40 кГц.

Частота дискретизации = 500 кГц.

Данные от АЦП - 14-разрядные.

 

Покопавшись в Матлабе сгенерил коэффициенты, и накидал такой код:

-- Формула фильтра Баттерворта 6-го порядка:
-- y = 1/a0 * ( b0*x + b1*x(-1) + b2*x(-2) + b3*x(-3) + b4*x(-4) + b5*x(-5) + b6*x(-6) - 
--	- a1*y(-1) - a2*y(-2) - a3*y(-3) - a4*y(-4) - a5*y(-5) - a6*y(-6) )
-- Полоса среза ФНЧ = 40 кГц
-- Порядок фильтра = 6
-- Частота дискретизации = 500 кГц
-- Из Матлаба получил следующие коэффициенты:
-- b0 = 1.0707e-004
-- b1 = 6.4241e-004
-- b2 = 1.6060e-003
-- b3 = 2.1414e-003
-- b4 = 1.6060e-003
-- b5 = 6.4241e-004
-- b6 = 1.0707e-004
-- a0 = 1.0000e+000
-- a1 = -4.0616e+000
-- a2 = 7.0995e+000
-- a3 = -6.7850e+000
-- a4 = 3.7230e+000
-- a5 = -1.1087e+000
-- a6 = 1.3966e-001
--
-- Для перехода от плавающей математики к целочисленной использую
-- масштабирующий коэффициент 2^14=16384
-- Т.к. а0=1, то 1/а0=1, и я просто убираю эту операцию из формулы
-- После умножения весовых коэффициентов на масштабирующий коэффициент и перенеся знак "-" в саму формулу получаю:
-- b0 = 1,75423488
-- b1 = 10,52524544
-- b2 = 26,312704
-- b3 = 35,0846976
-- b4 = 26,312704
-- b5 = 10,52524544
-- b6 = 1,75423488
-- a1 = -66545,2544
-- a2 = 116318,208
-- a3 = -111165,44
-- a4 = 60997,632
-- a5 = -18164,9408
-- a6 = 22881,8944

signal b0				: integer := 2;
signal b1				: integer := 11;
signal b2				: integer := 26;
signal b3				: integer := 35;
signal b4				: integer := 26;
signal b5				: integer := 11;
signal b6				: integer := 2;
--signal a0				: integer := 16384;
signal a1				: integer := 66545;
signal a2				: integer := 116318;
signal a3				: integer := 111165;
signal a4				: integer := 60998;
signal a5				: integer := 18165;
signal a6				: integer := 22882;

signal input			: integer := 0;

signal xv_0				: integer := 0;
signal xv_1				: integer := 0;
signal xv_2				: integer := 0;
signal xv_3				: integer := 0;
signal xv_4				: integer := 0;
signal xv_5				: integer := 0;
signal xv_6				: integer := 0;
signal yv_0				: integer := 0;
signal yv_1				: integer := 0;
signal yv_2				: integer := 0;
signal yv_3				: integer := 0;
signal yv_4				: integer := 0;
signal yv_5				: integer := 0;
signal yv_6				: integer := 0;

process(rst_n, clk)
begin
if rst_n = '0' then
	input			<= 0;
else
	if clk'event and clk = '1' then
		input			<= conv_integer(unsigned(din));
	end if;
end if;
end process;

process(rst_n, clk)
variable y_temp		: integer := 1;
begin
if rst_n = '0' then
	xv_0			<= 0;
	xv_1			<= 0;
	xv_2			<= 0;
	xv_3			<= 0;
	xv_4			<= 0;
	xv_5			<= 0;
	xv_6			<= 0;
	y_temp			:= 0;
	yv_0			<= 0;
	yv_1			<= 0;
	yv_2			<= 0;
	yv_3			<= 0;
	yv_4			<= 0;
	yv_5			<= 0;
	yv_6			<= 0;
else
	if clk'event and clk = '1' then
		xv_0			<= input;
		xv_1			<= xv_0;
		xv_2			<= xv_1;
		xv_3			<= xv_2;
		xv_4			<= xv_3;
		xv_5			<= xv_4;
		xv_6			<= xv_5;

		y_temp			:= b0*xv_0 + b1*xv_1 + b2*xv_2 + b3*xv_3 + b4*xv_4 + b5*xv_5 + b6*xv_6 +
						a1*yv_1 -  a2*yv_2 + a3*yv_3 - a4*yv_4 + a5*yv_5 - a6*yv_6;
		-- Делим получившийся результат на масштабирующий коэффициент
		yv_0			<= y_temp/16384;
		yv_1			<= yv_0;
		yv_2			<= yv_1;
		yv_3			<= yv_2;
		yv_4			<= yv_3;
		yv_5			<= yv_4;
		yv_6			<= yv_5;

		-- Переводим из типа integer в тип std_logic_vector на 14 разрядов
		-- для записи в ФИФО
		fifo_din		<= CONV_STD_LOGIC_VECTOR((yv_0), 14);
	end if;
end if;
end process;

 

В Моделсиме получаю такую картину:

Частота T_clk=500 кГц. Период входной пилы = 46 мкс, т.е. частота входного сигнала = 21,7 кГц. Исходя из этого я думаю, что сигнал почти весь должен пройти через фильтр без изменения. НО на картинке какая-то ересь.

 

Что мне делать, куда смотреть?

[des00 25]

Покопавшись в Матлабе сгенерил коэффициенты, и накидал такой код:

отложить моделсим и промоделировать фильтр в симулинке.

 

[billidean 0]

Я бы рад, но я с Матлабом не особо работал. А с симулинком вообще не знаком.

Как-то не было нужды до этих пор.

 

Но я чувствую, что выходная картинка должна быть другой.

 

Может быть входной сигнал какой-нибудь другой загнать, чтобы получить на выходе какой-нибудь заранее известный сигнал?

Чота я как-то вообще не могу вникнуть.

[Golikov 0]

самое простое начать все же с матлаба.

 

а сигнал то понятно какой, сумма 2 синусов, с частотой до среза и с частотой после среза, очевидно при правильном фильтре должен остаться практически один синус с частотой до среза

[billidean 0]

Решил попробовать написать на Си этот фильтр.

Вот две функции:

Вычисления с плавающей точкой:

void MainWindow::butter(float s)
{
    float b0 =  0.0001;
    float b1 =  0.0006;
    float b2 =  0.0016;
    float b3 =  0.0021;
    float b4 =  0.0016;
    float b5 =  0.0006;
    float b6 =  0.0001;
    float a0 = 1;
    float a1 = -4.0616;
    float a2 = 7.0995;
    float a3 = -6.7850;
    float a4 = 3.7230;
    float a5 = -1.1087;
    float a6 =  0.1397;

    static float s0;
    static float s1;
    static float s2;
    static float s3;
    static float s4;
    static float s5;
    static float s6;

    static float y0;
    static float y1;
    static float y2;
    static float y3;
    static float y4;
    static float y5;
    static float y6;

    float t = 0.0;

    s6 = s5;
    s5 = s4;
    s4 = s3;
    s3 = s2;
    s2 = s1;
    s1 = s0;

    y6 = y5;
    y5 = y4;
    y4 = y3;
    y3 = y2;
    y2 = y1;
    y1 = y0;

    s0 = s;
    y0 = b0*s0 + b1*s1 + b2*s2 + b3*s3 + b4*s4 + b5*s5 + b6*s6 - a1*y1 - a2*y2 - a3*y3 - a4*y4 - a5*y5 - a6*y6;
    qDebug() << y0;
}

Вычисления с фиксированной точкой:

#define MM   16384 // 2^14

void MainWindow::butter_i(int s)
{
    int b0 =  0.0001*MM;
    int b1 =  0.0006*MM;
    int b2 =  0.0016*MM;
    int b3 =  0.0021*MM;
    int b4 =  0.0016*MM;
    int b5 =  0.0006*MM;
    int b6 =  0.0001*MM;
    int a0 = 1;
    int a1 = -4.0616*MM;
    int a2 = 7.0995*MM;
    int a3 = -6.7850*MM;
    int a4 = 3.7230*MM;
    int a5 = -1.1087*MM;
    int a6 =  0.1397*MM;

    static int s0;
    static int s1;
    static int s2;
    static int s3;
    static int s4;
    static int s5;
    static int s6;

    static int y0;
    static int y1;
    static int y2;
    static int y3;
    static int y4;
    static int y5;
    static int y6;

    s6 = s5;
    s5 = s4;
    s4 = s3;
    s3 = s2;
    s2 = s1;
    s1 = s0;

    y6 = y5;
    y5 = y4;
    y4 = y3;
    y3 = y2;
    y2 = y1;
    y1 = y0;

    s0 = s*MM;
    y0 = (b0*s0 + b1*s1 + b2*s2 + b3*s3 + b4*s4 + b5*s5 + b6*s6 - a1*y1 - a2*y2 - a3*y3 - a4*y4 - a5*y5 - a6*y6)/MM;
    qDebug() << (int)(y0);
}

На вход подал постоянку (10).

Первая ф-ция дает такой результат:

Вторая:

 

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

[Golikov 0]

1. int b0 = 0.0001*MM;

 

я бы все же написал бы

int b0 = (int) (0.0001*((float)MM));

 

 

 

2.b0*s0 => b0*s0*2^14 *2^14 = b0*s0*2^28.

то есть на число у вас осталось 3 бита, если int - 32 битный...

 

 

 

 

ну и последние деление...

(int)((float)(b0*s0 + b1*s1 + b2*s2 + b3*s3 + b4*s4 + b5*s5 + b6*s6 - a1*y1 - a2*y2 - a3*y3 - a4*y4 - a5*y5 - a6*y6)/(float)MM);

 

и делить надо вроде как не на MM, а на MM*MM вы же оба числа умножили

0.1 * 0.1 == 0.1*10/10 *0.1*10/10 == 1*1 / 10 / 10 == 1*1/(10*10)

нет тут верно, чтобы сразу осталось число домноженное на MM

 

 

ну и входной параметр

void MainWindow::butter_i(int s)

не int должен быть, а флоат, если вы его умножаете в функции ...

 

 

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

[billidean 0]

Используя Си я пытаюсь преобразовать математику из плав. в фикс.точку.

 

Вместо int использовал везде long.

 

int b0 = 0.0001*MM; -- это чтобы не пересчитывать на калькуляторе каждый коэффициент, это место можно не анализировать, расчет ведется правильно, проверил отладчиком.

 

На счет:

(int)((float)(b0*s0 + b1*s1 + b2*s2 + b3*s3 + b4*s4 + b5*s5 + b6*s6 - a1*y1 - a2*y2 - a3*y3 - a4*y4 - a5*y5 - a6*y6)/(float)MM);

- когда буду реализовывать на VHDL, как мне преобразовывать в тип float? Да и на результат эти преобразования не повлияют.

 

void MainWindow::butter_i(int s)

не int должен быть, а флоат, если вы его умножаете в функции ...

Сами данные будут поступать с АЦП в виде 14-разрядов. В тип float я их не смогу преобразовать в ПЛИС. Я их сдвигаю на 14 рр и считаю, что эти 14-разрядов представляют целую часть числа в фиксированной точке.

 

Думается, что где-то кардинальная ошибка.

 

Может обратная связь при фиксированной точке дает некоторую ЗНАЧИМУЮ погрешность??

 

В-общем дошел до какого-то варианта решения:

все типы заменил на long long и получил такой график:

 

Это конечно извращение - такими типами манипулировать, но..

Теперь не знаю, как это все в VHDL буду перекидывать...

[Golikov 0]

там если один раз воткнуть как оно все делается, то потом все становиться легко и естественно... главное что становиться легко контролировать разрядность, и вылезать за 32 бита не надо будет, как у вас с лонгами.

 

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

[Maverick_ 15]

в матлабе по работе с фильтрами ничего сложного нет - попробуйте

на русском - введение

 

[billidean 0]

Спасибо всем за помощь, с этим фильтром разобрался, его реализация в ПЛИС работает достойно.

 

Теперь еще одна просьба: Помогите найти расчет коэффициентов для фильтра ФНЧ Бесселя 6-го порядка.

 

В Матлабе в fdatool для такого фильтра расчета нет. В инете тоже не могу найти.

 

Саму функцию фильтра я имею, а вот где взять коэффициенты, ума не приложу.

[Гость TSerg]

Filter Solutions

[billidean 0]

TSerg Спасибо огромное!

 

Сравнил результат этой проги по уже сделанному фильтру Баттерворта - коэффициенты совпали. => Значит можно верить программе.

[Гость TSerg]

Скорее я Матлабу не поверю

P.S.

Рад за Вас и что все получилось.

 

[Hose 0]

самое простое начать все же с матлаба.

 

а сигнал то понятно какой, сумма 2 синусов, с частотой до среза и с частотой после среза, очевидно при правильном фильтре должен остаться практически один синус с частотой до среза

 

Лучше подать много синусов) то есть дельта-фукцию, тем более, что ее значительно проще сгенерить, ну и вообще классики пишут, что фильтр на дельта-функцию откликается своей импульсной характеристикой.

[Golikov 0]

это было давно и уже не актуально. А подавать лучше 2 синуса, правду говорю;)