Помощь - Поиск - Пользователи - Календарь
Полная версия этой страницы: Создание фильтра в MATLAB и его реальное применение
Форум разработчиков электроники ELECTRONIX.ru > Цифровая обработка сигналов - ЦОС (DSP) > Алгоритмы ЦОС (DSP)
sidy
Здравствуйте, уважаемые форумчане. В MATLAB в утилите FDATool я выбрал фильтр Баттерворта 2ого порядка, частоту среза 50 Гц, частоту дескритизации 5000 Гц и расчитал коэффициенты фильтра. После этого выполнил Convert To Single Section. Получил следующие коэффициенты:
numerator:
0.0009446918438402,
0.00188938368768,
0.0009446918438402

denominator:
1,
-1.911197067426,
0.9149758348014

Могу ли я использовать эти коэффициенты для построения фильтра в цифровом сигнальном процессоре для фильтрации входного сигнала с АЦП (входная шкала напряжения от 0 до 4096)? Или я должен сделать еше какието манипуляции с данными коэффциентами?
Т.к. используя эти коэффициенты, получаю нулевой сигнал. Посмотрев коэффициенты фильтров в различной литературе, заметил что они близки к 1 или 2. А у меня коэффициенты b получились очень маленькими.
V_G
Это распространенная проблема БИХ-фильтров - существенно отличающиеся коэффициенты а и b, не вписывающиеся в точность целочисленных DSP. Борьба с этим - в грамотном разбиении на секции, матлаб на вашей задаче вполне с этим справляется, вы зря сделали Convert To Single Section. Либо повторите расчет, либо сделайте обратно Convert to Second-Order Sections. Там вполне нормальные коэффициенты. Просто отмасштабируйте их до диапазона работы вашего конкретного DSP (умножьте все на степень двойки)
alexeyv
Если ДСП fixed-poin, то необходимо конвертировать коэффициенты в целочисленный вид
sidy
Цитата(V_G @ Feb 19 2013, 03:25) *
Либо повторите расчет, либо сделайте обратно Convert to Second-Order Sections. Там вполне нормальные коэффициенты. Просто отмасштабируйте их до диапазона работы вашего конкретного DSP (умножьте все на степень двойки)

Да спасибо очень помогло. Получил следующие коэффициенты для Low-Pass фильтра:
ACoef[3]={1, -1.911197067426073, 0.914975834801434};
BCoef[3]={1, 2, 1};
Только в Matlab fdatool я получил фазовый сдвиг выходного сигнала относительно входного 90 гр, а в моем DSP фазовый сдвиг 270 гр. Непонятно с чем это связано.

Теперь мне нужен еще один фильтр - High-Pass для фильтрации постоянной составляющей, т.е. сигналы с частотой около 0 Гц. В Matlab fdatool я рассчитал коэффициенты High-Pass фильтра Баттерворта второго порядка:
ACoef_HP[3]={1, -1.998222847291842, 0.9982244250264};
BCoef_HP[3]={1, -2, 1};
и промоделировал прохождение сигнала через него в Simulink - фазового сдвига сигнала практически не наблюдал. При занесении коэффициентов в DSP выходной сигнал сдвинут относительно входного на 180 гр.

Т.е. мой фильтр вносит сдвиг в 180 гр по-мимо расчетного в Matlab FDATool. Не понятно с чем это связано.
STAR_IK
А вы как смотрите сдвиг? Это похоже на простую инверсию, возникающую на определенном этапе обработки, а фильтр может быть и ни при чем.
sidy
Цитата(STAR_IK @ Feb 28 2013, 08:05) *
А вы как смотрите сдвиг? Это похоже на простую инверсию, возникающую на определенном этапе обработки, а фильтр может быть и ни при чем.

Да, это была моя ошибка - фильтр тут не причем.
sidy
У меня возник еще один вопрос. Теперь мне необходим фильтр Баттерворта 3ого порядка, с частотой среза 50 Гц и частотой дескритизации 5000 Гц. Я вновь рассчитал коэффициенты в Matlab Fdatool и получил не совсем понятный мне набор коэффициентов, т.е. их больше четырех для A и B (см. рис.) и не понятно как с ними работать.
Нажмите для просмотра прикрепленного файла
Если я делаю Convert to Single section, то получаю нужное количество коэффициентов A и B - по четыре на каждый, но коэффициенты B получаются очень маленькими, о чем уже обсуждалось выше. Теперь вопрос можно ли каким-то образом преобразовать коэффициенты A и B что использовать их в DSP&
Alexey K
Два фильтра 2-го порядка
sidy
Задам еще один вопрос. Столкнулся с непонятной мне вещью: есть два фильтра реализованных на C цифровых фильтра Баттерворта фторого порядка: один НЧ, другой ВЧ. Вот код реализации
Код
//*Фильтр Баттерворта 2ого порядка Low-Pass-------------------------------------------------    
    //Сдвигаем старые значения
    for(n=2; n>0; n--){
        x[n]=x[n-1];
        y[n]=y[n-1];
    }
    //Вычисляем новое выходное значение
    x[0]=NewSample;
    y[0]=BCoef[0]*x[0];
    for(n=1; n<=2; n++)
        {y[0]+=(BCoef[n]*x[n])-(ACoef[n]*y[n]);}    
//*Фильтр Баттерворта 2ого порядка High-Pass-------------------------------------------------
    //Сдвигаем старые значения
    for(n=2; n>0; n--){
        x_HP[n]=x_HP[n-1];
        y_HP[n]=y_HP[n-1];
    }
    //Вычисляем новое выходное значение
    x_HP[0]=y[0];
    y_HP[0]=BCoef_HP[0]*x_HP[0];
    for(n=1; n<=2; n++)
        {y_HP[0]+=(BCoef_HP[n]*x_HP[n])-(ACoef_HP[n]*y_HP[n]);}
    in_new=y_HP[0];

И есть две среды разработки для STM32: Keil MDK-ARM 4.53.0 и CooCox COIDE 1.7.0. В обоих средах для фильтра используется код приведенный выше.
При использовании для прошивки кода скомпилированного в Keil на выходе ЦАП после фильтрации получаю нормальную синусоиду:
Нажмите для просмотра прикрепленного файла
А при использовании для прошивки кода скомпилированного в CooCox COIDE на выходе ЦАП получаю сигнал с просечками, которые мешают мне для дальнейшей обработки.
Нажмите для просмотра прикрепленного файла
С чем все это может быть связано?
sidy
Прошу прощение, что потревожил техническую общественность. Это была моя ошибка - не имеющая отношения к алгоритму фильтрации.
MSP430F
Всем доброго времени суток!
Коллеги, помогите!
Рассчитываю в MATLAB цифровой полосовой фильтр Чебышева II рода. Сделал разбиение на секции 2-го порядка, получил в хедере коэффициенты в виде двумерных массивов NUM и DEN. Формулы секций простые, все работает. Озадачился получением стабильной версии в виде одной секции. Порядок фильтра задал 8. Применил Convert Structure... к Lattice Arma, получил стабильный фильтр. Создал сишный хедер и получил два массива.
const int KL = 16;
const real32_T K[16] = {
-0.9990231991, 0.9992837906, -0.9991765022, 0.9989774227, -0.9984907508,
0.9831924438, 0.1230626926, 0.3759146333, -0.7456214428, -0.4113018215,
-0.2647816539, 0.1375903487, 0.199342683, 0.1202145815, 0.0357035771,
0.00457067322
};
const int VL = 17;
const real32_T V[17] = {
3.237323187e-010,8.469006829e-009,-4.508386837e-007,-1.037612128e-005,0.0004177154624,
0.009835269302, 0.04770011082, -0.07836318016, -0.08220990002, 0.04402532056,
0.002590565477, -0.3493199348, -0.3774584234, 0.05220580474, 0.3641556203,
0.2648055255, 0.06757183373
};

Смутило сразу названия массивов. Вместо NUM и DEN какие-то K и V. А потом посмотрел еще структуру этого Lattice Arma и совсем запутался.
Какие формулы надо использовать для этих коэффициентов ???
Поскажите, плиз.
sidy
Назрела необходимость перейти к вычислениям с фиксированной точкой. В связи с этим возник еще один вопрос: как преобразовать мои коэффициенты
ACoef[3]={1, -1.911197067426073, 0.914975834801434};
BCoef[3]={1, 2, 1};
к целочисленным.
Как я понял я должен умножить их на 2^k, где k целые числа; я попробовал промаштабировать исходные коэффициенты умножив на 32767, получил новые коэффициенты, но работа фильтра с новыми коэффициентами существенно изменилась - исчез фазовый сдвиг (90 гр.) и фильтр потерял свои "фильтрующие" свойства. Вопрос: что я сделал неправильно? (Разрядность АЦП 12 бит, разрядность MCU 32 бита).
fontp
QUOTE (sidy @ Jul 3 2013, 15:37) *
Назрела необходимость перейти к вычислениям с фиксированной точкой. В связи с этим возник еще один вопрос: как преобразовать мои коэффициенты
ACoef[3]={1, -1.911197067426073, 0.914975834801434};
BCoef[3]={1, 2, 1};
к целочисленным.
Как я понял я должен умножить их на 2^k, где k целые числа; я попробовал промаштабировать исходные коэффициенты умножив на 32767, получил новые коэффициенты, но работа фильтра с новыми коэффициентами существенно изменилась - исчез фазовый сдвиг (90 гр.) и фильтр потерял свои "фильтрующие" свойства. Вопрос: что я сделал неправильно? (Разрядность АЦП 12 бит, разрядность MCU 32 бита).


А как Вы представите 2 с фиксированой запятой 1.15 fixed-point
В 1.15 fixed-point представимы только числа от -1 до 1
Разделите на 2 и числитель и знаменатель и тогда уже... то есть умножайте на 16384

На самом деле в Matlab Filter Design можно включить режим конечной разрядности и он сам посчитает целочисленные ( 1.N fixed-point ) коэфф.,
и покажет частотные и фазовые характеристики с эффектами разрядности
sidy
Вот что у меня получилось:
Код
//float ACoef[3]={1, -1.913780606355776, 0.917346934620046};
//float BCoef[3]={1, 2, 1};
int16_t ACoef[3]={16384, -31355, 15029};
int16_t BCoef[3]={16384, 32767, 16384};
int y[3], x[3];
int16_t NewSample;
  //*Фильтр Баттерворта 2ого порядка Low-Pass-------------------------------------------------
  for(n=2; n>0; n--) {x[n]=x[n-1]; y[n]=y[n-1];} //Сдвигаем старые значения
  //Вычисляем новое выходное значение
  x[0]=NewSample<<4; y[0]=BCoef[0]*x[0];
  y[0]=y[0]>>15;
  for(n=1; n<=2; n++) {y[0]+=(BCoef[n]*x[n])-(ACoef[n]*y[n]);  
  y[0]=y[0]>>15;}

Рассчитал коэффициенты фильтра, умножив, исходные полученные в Matlab на 16384. Также я сдвигаю полученные в 12 битном АЦП отсчеты на 4 влево (x[0]=NewSample<<4;) для приведения к формату Q1.15 и еще сдвигаю рассчитанный выходной сигнал на 15 вправо. После этого вывожу отсчеты на ЦАП. Фильтр не работает т.е. искаженная синусоида проходит без фильтрации и фазовый сдвиг вместо 90гр получился около 5гр (см. рис., входная синусоида желтая).
Нажмите для просмотра прикрепленного файла
С float арифметикой все работало, один в один как в Matlabe:
Правильная фильтрация:
Нажмите для просмотра прикрепленного файла
Что я делаю не так?
bullit
Простите что вмешиваюсь!
А чем вам float не угодил? СТМки с fpu прекрасно переваривают флоаты. Посмотрите статью http://www.compeljournal.ru/enews/2012/6/3/ особенно раздел с плавающей запятой. Но не вылезайте за границы флоата
Например расчет квадратного корня 14 тактов, проверял, правда вышло 15 тактов))))

Да и ФНЧ + ФВЧ ~ ПФ, не?


BooZe

Вот результаты моделирования Вашего фильтра с нормированными коэффициентами. Виден спад АЧХ из-за нормирования коэффициентов фильтра.
Лучше пользоваться стандартной функцией fdatool для нормировки коэффициентов.

На скриншоте виден спад из-за квантования (штриховая - необходимая Вам АЧХ).
Кроме того фильтр нестабилен, что тоже не очень хорошо.
sidy
Задам еще один вопрос. Допустим я фильтрую входной синусоидальный сигнал. В отсутствии входного сигнала вследствие шума АЦП на входе фильтра появляются отсчеты, которые не являются полезным сигналом, но также фильтруются и портят алгоритм (поскольку с помощью фильтра я определяю переход синусоиды через 0). Как в таком случае обычно поступают: перестают накладывать фильтр, или подают на него 0, или же существуют какие-либо другие варианты?
bullit
Помоему из этого:
Цитата
...также фильтруются и портят алгоритм (поскольку с помощью фильтра я определяю переход синусоиды через 0).
следует что Вы или не тот алгоритм взяли, или фильтр не тот или Вам вообще фильтрация не нужна.

И давайте лучше сначала (поскольку речь вые шла только про фильтры), Вы что хотите из сигнала "выудить"?
sidy
Выудить хочу переход через ноль, поскольку в исходном сигнале может присутсвовать достаточно сильная помеха именно в переходах, вплоть до такой как на рис.
Нажмите для просмотра прикрепленного файла
bullit
Цитата
но также фильтруются и портят алгоритм

Каким образом, покажите. И давайте больше информации, тут нет телепатов)
Верхний и нижний график это после фильтра или как?
Покажите лучше исходный сигнал, со всеми вариантами шума.
И.. попробуйте загнать данные сначала в матлаб, там его как следует обработать, ну и потом выбирайте алгоритм.

ну и Вы не до конца ответили на вопрос... зачем вам переход через ноль? Вы, надеюсь, в курсе, что выход фильтра имеет "запаздывание" по выходу?
Может Вам просто полосовой фильтр на основную частоту попробовать?
_VGA_
помогите разобраться
делаю фильтр так:
-------
BIQUADS=3;
fs=48000;
fc=450;
fn=fc/fs;
[b,a]=butter(2*BIQUADS,fn);
[sos,g]=tf2sos(b,a);
sos=sos';
d=[0,0];
fp2=fopen('low.dat','w');
for i=1:BIQUADS,
j=(i-1)*6;
sos(j+4:j+6)=-sos(j+4:j+6);
fprintf(fp2,'%2.20e,\n',sos(j+1));
fprintf(fp2,'%2.20e,\n',sos(j+2));
fprintf(fp2,'%2.20e,\n',sos(j+5));
fprintf(fp2,'%2.20e,\n',sos(j+3));
fprintf(fp2,'%2.20e,\n',sos(j+6));
fprintf(fp2,'%2.20e,\n',d(1));
fprintf(fp2,'%2.20e,\n',d(2));
end
fprintf(fp2,'%2.20e\n',g);
fclose(fp2);
----------
дальше эти коэф загружаю в adsp
хочу посмотреть график:
freqz(b,a,2048,2000);
выводится график, но частота среза раз в 10 ниже (примерно 45гц)
пока два вопроса
1 правильно ли я расчитваю коэфф
2 почему неправильно рисуется график
спс
SemperAnte
fn - это же нормализованная частота?
Она должна считаться как f/(Fs/2).

АЧХ вывожу командой:

Код
h = fvtool(b,a);
set(h, 'Fs', fs);


Вроде, совпадает с заданной.
_VGA_
fs - частота квантования (48 кгц)
fc - частота среза фильтра (450 гц)
fn =2*fc/fs - я правильно понял?
SemperAnte
Цитата(_VGA_ @ Nov 12 2013, 15:20) *
fs - частота квантования (48 кгц)
fc - частота среза фильтра (450 гц)
fn =2*fc/fs - я правильно понял?

Да, так. Отношение частоты среза к частоте Найквиста.

В Matlab примере на функцию butter это показано (по крайней мере, в моей 2012b версии).
_VGA_
там хелп есть? ну я и бестолочь sad.gif
Спасибо получилось, красиво;
теперь я хочу фнч и фвч с одинаковой частотой среза вывести на график
(результирующую)
это возможно?
Для просмотра полной версии этой страницы, пожалуйста, пройдите по ссылке.
Invision Power Board © 2001-2025 Invision Power Services, Inc.