Jump to content

    
Sign in to follow this  
billidean

помогите с фильтром Баттерворта

Recommended Posts

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

 

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

 

Необходимо написать код 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;

 

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

post-59925-1418655845_thumb.png

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

 

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

Share this post


Link to post
Share on other sites
Покопавшись в Матлабе сгенерил коэффициенты, и накидал такой код:

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

 

Share this post


Link to post
Share on other sites

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

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

 

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

 

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

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

Share this post


Link to post
Share on other sites

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

 

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

Share this post


Link to post
Share on other sites

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

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

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

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).

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

post-59925-1418711152_thumb.png

Вторая:

post-59925-1418711158_thumb.png

 

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

Share this post


Link to post
Share on other sites

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 должен быть, а флоат, если вы его умножаете в функции ...

 

 

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

Share this post


Link to post
Share on other sites

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

 

Вместо 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 и получил такой график:

post-59925-1418715867_thumb.png

 

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

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

Share this post


Link to post
Share on other sites

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

 

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

Share this post


Link to post
Share on other sites

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

 

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

 

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

 

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

Share this post


Link to post
Share on other sites

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

 

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

Share this post


Link to post
Share on other sites
самое простое начать все же с матлаба.

 

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

 

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

Share this post


Link to post
Share on other sites

Join the conversation

You can post now and register later. If you have an account, sign in now to post with your account.

Guest
Reply to this topic...

×   Pasted as rich text.   Paste as plain text instead

  Only 75 emoji are allowed.

×   Your link has been automatically embedded.   Display as a link instead

×   Your previous content has been restored.   Clear editor

×   You cannot paste images directly. Upload or insert images from URL.

Sign in to follow this