Цифровые фильтры, помогите разобраться Опции

[TViT]

Все вроде понятно по фильтрам. Задаешь частоты расчитываешь коэффициенты. Поэтому вопросы у меня как обычно чисто практические...

x(n)=входной массив выборок АЦП (n=0…1023)
y(n)=выходной массив

стандартная формула Баттерворта на которую натыкаешься всюду:

y(n)=( a0*x(n)+a1*x(n-1)+a2*x(n-2) - (b1*y(n)+b2*y(n-1)) )

так вот не могу понять как при начальном значении первого(нулевого) значения выборки по адресу x(0)
может быть a1*x(n-1) запрашивается значение из массива x(0-1) этот же выход из диапазона массива? Массив от 0 до 1023 –> 1024 выборки
сделал значения по порядку:

y(n)=( a0*x(n)+a1*x(n+1)+a2*x(n+2) - (b1*y(n)+b2*y(n+1)) )

все заработало для ФВЧ, а для ПФ не работает толи коэффициенты не правильно подставляю, толи формула другая для ПФ там просто коэффициентов больше 5 вроде поэтому делаю так

y(n)=( a0*x(n)+a1*x(n+1)+a2*x(n+2)+ a3*x(n+3)+ a4*x(n+4) - (b1*y(n)+b2*y(n-1)+ b3*y(n+2)+b3*y(n+3)))

видел и такую формулу где-то:

y(n)=( a0*x(n)+2*a1*x(n-1)+a2*x(n-2) - (b1*y(n)+b2*y(n-1)) )

В общем полная каша в голове, может для Чебышева иначе делается, подскажите плиз…

Сообщение отредактировал TViT - Nov 10 2010, 13:34

[Джеймс]

Алгоритм одинаковый что для ФВЧ, что для ПФ. Меняются только коэффициенты.
Посмотрите программу на стр. 12:
http://www.vlsi.ss.titech.ac.jp/~isshiki/V...stemVIII_06.pdf

[TViT]

Джеймс Да вобщем и у меня так там последовательно складываются входные значения умноженные на коэффициенты еще подробней посмотрю
расчитал в ciirf1 получается для 2-го порядка - 5 коэфф. а0,а1,а2 и в1,в2 в Матлабе всего 3 выкинул обратную связь чтоб только 3-мя коэфф. пользоваться все равно работает:

Y(n) = (a(0) * X(n) + a(1) * X(n + 1) + a(2) * X(n + 2))

Вот сравниваю тут коэффициенты может нужно для ПФ разделить коэффициенты, первую половину получится как через НЧ фильтр прошел сигнал, потом его переписать во входной массив X(n) и пропустить через вторую половину коэффициентов как бы через ФВЧ получится ПФ.

Сообщение отредактировал TViT - Nov 10 2010, 19:54

[DRUID3]

так вот не могу понять как при начальном значении первого(нулевого) значения выборки по адресу x(0)
может быть a1*x(n-1) запрашивается значение из массива x(0-1) этот же выход из диапазона массива?

Это просто показывает, что в фильтре используются предыдущие отсчеты... И разумеется в первые моменты времени там ничего нет ("0"). Надеюсь Вас не удивит, что это одна из форм т.н. "переходных процессов". Написать грамотно такой фильтр не тяжело, нужно только сеть и немного подумать.

все заработало для ФВЧ, а для ПФ не работает толи коэффициенты не правильно подставляю, толи формула другая для ПФ

...формула та же, но там возможно несколько вариантов этих формул (для IIR) и все они правильные.

Вот сравниваю тут коэффициенты может нужно для ПФ разделить коэффициенты, первую половину получится как через НЧ фильтр прошел сигнал, потом его переписать во входной массив X(n) и пропустить через вторую половину коэффициентов как бы через ФВЧ получится ПФ.

...буду краток - нет. У Вас какая-то путаница с пониманием форм построения фильтра FIR-IIR и форм АЧХ(ФНЧ,ФВЧ,полосовой и т.д.)

[TViT]

...Понимаю что гдето моя ошибка... Продумаю по внимательнее алгоритм по шагам...

Сообщение отредактировал TViT - Nov 11 2010, 08:06

[TViT]

Все разобрался просто игнорирую ошибку выход из диапазона массива в цикле перемножения и сложения выборок и тогда все отлично работает
Всем за участие спасибо!

[Guest_TSerg_*]

Занятное решение.

[TViT]

TSerg Это если быстро сказать. На самом деле привел все к соответствию формулы, взял алгоритм из DSPlib AVR32 переписал на VB, а так как там через указатели, то я просто применил возможность игнорировать ошибки на VB результаты абсолютно одинаковые что на МК передаю массив и результат обратно что через VB на компе можно конечто через смещение сделать как в примере от Джеймс.

DSPlib коэффициенты через Scilab расчитываются поэтому с этим разбираюсь от Матлаба не подходят, толи что-то не так делаю. Только из ciirf1 подходят коэффициенты.... Вот такие пироги.

Сообщение отредактировал TViT - Nov 12 2010, 07:37

[bahurin]

TSerg Это если быстро сказать. На самом деле привел все к соответствию формулы, взял алгоритм из DSPlib AVR32 переписал на VB, а так как там через указатели, то я просто применил возможность игнорировать ошибки на VB результаты абсолютно одинаковые что на МК передаю массив и результат обратно что через VB на компе можно конечто через смещение сделать как в примере от Джеймс.

DSPlib коэффициенты через Scilab расчитываются поэтому с этим разбираюсь от Матлаба не подходят, толи что-то не так делаю. Только из ciirf1 подходят коэффициенты.... Вот такие пироги.

дааа сурово все однако.

[Guest_TSerg_*]

Я - пас.
Мой мозг отказался перевести сообщение №8 в семантически значимую для меня форму.

[TViT]

Вам язык почесать негде? Или я что-то грубое сказал???

[bahurin]

Вам язык почесать негде? Или я что-то грубое сказал???

Это что еще за агрессия? Вообще то это у вас проблемы а не у нас. Если хотите получить совет то вот вам совет: Бегите от VB со всех ног в сторону С++. Это первое, второе бегите еще быстрее чем от VB от SCILab в сторону matlab или GNU Octave. Поверьте если результаты SCILab отличаются от matlab, то врет SCILab а не матлаб. Что касается dsp - matlab можно считать эталоном.

[DRUID3]

бегите еще быстрее ... в сторону matlab

...а еще никаДДа не используйте GCC... Ведь всем известно, что настоящие гуру от DSP используют исключительно кошерные решения - в виде готовых примеров поставляемых с IDE

GNU Octave.

...это жалкая пародия как на матлаб так и на скайлаб...

Что касается dsp - matlab можно считать эталоном.

да-да... Как дэлфи-кодер - эталоном программиста...

P.S.: по-ходу у аФФтАра глючит исходник, а мы тут "холиварим"...

[bahurin]

...а еще никаДДа не используйте GCC... Ведь всем известно, что настоящие гуру от DSP используют исключительно кошерные решения - в виде готовых примеров поставляемых с IDE

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

...это жалкая пародия как на матлаб так и на скайлаб...

Неправда GNU Octave практически полный аналог матлаба, совершенно бесплатный, с великим множенством тулбоксов. Оданко вы не найдете тормозных симулинков в октаве и красивого интерфейса, но что касается функций и их использования, то пожалуй октава ничем не хуже матлаба.

да-да... Как дэлфи-кодер - эталоном программиста...

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

[DRUID3]

Неправда GNU Octave практически полный аналог матлаба, совершенно бесплатный, с великим множенством тулбоксов. Оданко вы не найдете тормозных симулинков в октаве и красивого интерфейса, но что касается функций и их использования, то пожалуй октава ничем не хуже матлаба.

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

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

Вы правы - я троль... но вернемся к вопросу. Миллиарды... Им пользуются миллиарды... В столице Украины при десятках тысяч работающих программистов хорошо если есть сотня активно юзающих математические пакеты. Кстати массы добровольных тестировщиков не всегда есть критерий безошибочности. Ошибку в quickSort во FreeBSD подправили пару лет назад - а если историю жизни проекта умножить на количество пользователей как раз и получатся трилиарды гуманоидов в мирриадах звездных систем . Скайлаб в эпоху своей молодости тоже был коммерческим продуктом с радиотехническим уклоном! Для математики он юзает известные вычислительные библиотеки которыми пользуются все и далеко не только в этом проекте. Вряд ли у него есть ошибки в вычислительном ядре. Кстати, никто и нигде не дает гарантий что их нет в Матлаба . Еще с помощью Скайлаба очень удобно тестировать свои DSP исходники(вычислительную часть) не знаю как с этим в "матлабе"... Давайте говорить честно - матлаб это модно. Столько функций, а главное дофига всего в книгах и интернете - почти не надо думать. Токо копируй и в консоль, или как там оно зовеЦЦо у него, вставляй...

[bahurin]

Давайте говорить честно - матлаб это модно. Столько функций, а главное дофига всего в книгах и интернете - почти не надо думать. Токо копируй и в консоль, или как там оно зовеЦЦо у него, вставляй...

Все правильно в это и есть преймущество матлаба! Вы пользуетесь результатами других людей, которые положили много много времени чтобы дать вам этот инструмент.

[TViT]

Это что еще за агрессия? Вообще то это у вас проблемы а не у нас.

bahurin А что можно еще сказать когда ничего конкретного и конструктивного не сказано люди заходят усмехнуться и съязвить...

Бегите от VB со всех ног в сторону С++. Это первое, второе бегите еще быстрее чем от VB от SCILab в сторону matlab или GNU Octave. Поверьте если результаты SCILab отличаются от matlab, то врет SCILab а не матлаб. Что касается dsp - matlab можно считать эталоном.

Я не думаю что есть разница на чем писать. Мне если честно без разницы VB, ASM, C, Delphi просто результат нужен побыстрее, а программирую я поскольку постольку - редко, а на VB набросать быстренько инструмент для отладки пара пустяков и столько же минут. А на Си пишу для AVR, потому нет наработок и примеров для Си скачивать тоже искать нужно и разбираться на VB многое уже написано и понятно. Время крайне мало. Да и скорость выполнения не важна вот и пишу на VB.

Также не думаю что большая разница где коэффициенты считать тут больше алгоритм важен простота использования, а Матлаб платная среда можно ли результаты вычислений использовать в своих проектах если они будут платными, думаю чтоб таких проблем не возникало Atmel использует все фришное AVR32 Studio, Си(gcc), Scilab.

Еще подскажите как лучше промоделировать АЧХ у меня в программе сравнить так скажем задумка генерить с цикле синусоиды по возрастанию частоты и пропускать через фильтр и строить график

[bahurin]

А что можно еще сказать когда ничего конкретного и конструктивного не сказано люди заходят усмехнуться и съязвить...

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

Матлаб платная среда можно ли результаты вычислений использовать в своих проектах если они будут платными

Можно матлаб для того и существует. Если конечно у вас матлаб лицензионный.

Еще подскажите как лучше промоделировать АЧХ у меня в программе сравнить так скажем задумка генерить с цикле синусоиды по возрастанию частоты и пропускать через фильтр и строить график

подставить в передаточную характеристику H(z) коэффициенты которой вы посчитали z=exp(j*w) и получите комплексный к-т передачи H(exp(jw)). АЧХ при этом модуль H(exp(jw)).

[TViT]

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

вот код переписан из DSPlib AVR32:

Код

       size = 1023                                               ' //1024 выборки
For n = 0 To (size)
       sum1 = 0: sum2 = 0
       For m = 0 To (ComboNum.ListCount - 1)     ' кол-во коэф. "А" в списке 5-1 потому что в массиве от 0
              sum1 = sum1 + (num(m) * Xn(n - m))
       Next m
       For m = 0 To (ComboDen.ListCount - 1)      'коэф. "B"
              sum2 = sum2 + ((den(m) * Yn(n - m - 1)))

' //Yn(n) = (num(0) * Xn(n) + num(0) * Xn(n - 1) + num(0) * Xn(n - 2)) - (den(0) * Yn(n)+ den(1) * Yn(n-1) и т.д.) раньше так работало

       Next m
       Yn(n) = sum1 - sum2
Next n

коэффициенты из ciirf1 семплирование - 11025Гц ПФ - 800-900Гц:

Код

      секция 1
A:  0.000780  0.000000  -0.001561  0.000000  0.000780
B:  1.000000  -3.469910  4.930680  -3.332800  0.922566

Сообщение отредактировал TViT - Nov 14 2010, 14:02

[bahurin]

Код понять не могу, т.к. на VB не пишу. Но одно кидается в глаза. У фильтра принято обозначать к-ты числителя ПФ через b, а к-ты знаменателя через A. У вас обозначения наоборот, потому что A[0] должно быть равно 1. Возможно в этом ваша ошибка. Построил в матлабе АЧХ фильтра с приведенными к-тами:

Код

B=[0.000780  0.000000  -0.001561  0.000000  0.000780];
A=[1.000000  -3.469910  4.930680  -3.332800  0.922566];

h = zeros(1,1024);
h(2) = 1;
h = filter(B,A,h);

H = fft(h);
Fs = 11025;
f = Fs*(0:1023)/1024;
plot(f,20*log10(H));

фильтр как фильтр, только В с А местами поменял.

[TViT]

вот тот же код на Си из AVR32 DSPlib:

Код

{
   int n, m;
   S32 sum1, sum2;
   dsp16_t *px, *py;

   num_prediv = DSP16_QB - num_prediv; //15-3
   den_prediv = DSP16_QB - den_prediv; //15-3

   for(n=0; n<size; n++)
   {
     sum1 = 0;
     sum2 = 0;
     px = &x[n];
     for(m=0; m<num_size; m++)
       sum1 += ((S32) num[m] * (S32) px[-m]);
     py = &y[n-1];
     for(m=0; m<den_size; m++)
       sum2 += ((S32) den[m] * (S32) py[-m]);

     sum1 >>= num_prediv; //>>3
     sum2 >>= den_prediv; //>>3

     y[n] = sum1 - sum2;
   }
}

А с В поменять странно попробую, хотя все работает, немогу понять при расчете через Матлаб fdatool где там А... коэфф где В...

[bahurin]

Не видел таких конструкций.
не понятен смысл

Код

py = &y[n-1];

при n=0, не понятно на что будет указывать py. Аналогично

Код

sum1 += ((S32) num[m] * (S32) px[-m]);

в цикле по m эта конструкция должна по идее должна приводить к эксцепшену по памяти.

[TViT]

Как ни странно все работает. Думаю так для производительности может сделали чтобы небыло лишних циклов записи чтения памяти она там флеш медленная.

Вот по тем же параметрам посчитал в матлабе где тут что, где А где В.

Код

/*
* Filter Coefficients (C Source) generated by the Filter Design and Analysis Tool
*
* Generated by MATLAB(R) 7.9 and the Signal Processing Toolbox 6.12.
*
* Generated on: 15-Nov-2010 18:56:58
*
*/
/*
* Discrete-Time IIR Filter (real)
* -------------------------------
* Filter Structure    : Direct-Form II, Second-Order Sections
* Number of Sections  : 1
* Stable              : Yes
* Linear Phase        : No
*/

#define MWSPT_NSEC 3
const int NL[MWSPT_NSEC] = { 1,3,1 };
const real32_T NUM[MWSPT_NSEC][3] = {
  {
    0.02771298587,              0,              0
  },
  {
                1,              0,             -1
  },
  {
                1,              0,              0
  }
};
const int DL[MWSPT_NSEC] = { 1,3,1 };
const real32_T DEN[MWSPT_NSEC][3] = {
  {
                1,              0,              0
  },
  {
                1,   -1.721542239,   0.9445739985
  },
  {
                1,              0,              0
  }
};

Сообщение отредактировал TViT - Nov 15 2010, 16:12

[bahurin]

NUM[MWSPT_NSEC][3] это B (числитель)
DEN[MWSPT_NSEC][3] зто А (знаменатель)

[AlikM]

Но одно кидается в глаза. У фильтра принято обозначать к-ты числителя ПФ через b, а к-ты знаменателя через A. У вас обозначения наоборот, потому что A[0] должно быть равно 1.

Интересно, это где же так принято и по каким книгам Вас учат (учили)?
Вот ссылка на "библию начинающего радиоинженера" РтЦиС Гоноровского (в простонародии Гоноревич): Гоноровский
И если взлянуть на страницу 500, то увидите что там как раз наоборот.
Хотя это совсем не принципиально, а вот что принципиально- это если бы автор темы не поленился и взял и синтезировал этот фильтр 2-го порядка на бумажке собственными руками, то у него и для более сложных фильтров вопросов больше бы не возникло.
А если надо быстро расчитать фильтр и получить его программную реализацию причем в С коде, тогда вам сюда FilterSolut

[PetrovichKR]

Интересно, это где же так принято и по каким книгам Вас учат (учили)?

Во всех книгах по цифровой обработке сигналов, которые я видел, коэффициенты рекурсивной части фильтра обозначают через a, нерекурсивной части через b. Гоноровский к таким книгам не относится, поскольку в ней изложены основы радиотехники, и цифровая обработка затрагивается поверхностно.
Но, опять же, это не принципиально.

[TViT]

Ребят кто разбирается в Scilab помогите еще, нужно рассчитать коэффициенты по секциям последовательно, как в программке ciirf1. Массивом получается, а по секциям не могу понять как разбить, рассчитать в Scilab.

Вот код из AVR32 DSPlib (у меня урезанный только параметры коэффициенты и график кидаю весь код мало ли что-то по незнанию выкинул нужное):

Сообщение отредактировал TViT - Dec 18 2010, 05:56

Прикрепленные файлы

 iir.rar ( 1.55 килобайт ) Кол-во скачиваний: 25

 

[lisstret]

У фильтра принято обозначать к-ты числителя ПФ через b, а к-ты знаменателя через A.

Вопрос спорный, в России и в зарубежье они по разному обозначались, соответственно и в Матлабе они обозначаются как принято за бугром. А так как щас новые книжки ориентированны на матлаб, то и в книгах российских авторов приняты другие обозначения. В книжке Сергиенко даже есть сноска где это поясняется! Короче из за этого и различия.

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

Мне кажется что он писал так, как было принято у нас. А не из за того, что он чет не догонял.

Сообщение отредактировал lisstret - Dec 19 2010, 20:53

[TViT]

Все понятно с коэффициентами. Подскажите по существу, по Scilab...

[TViT]

Нашел функцию casc — cascade realization of filter from coefficients, что она делает и как ее юзать не пойму?

[PetrovichKR]

Мне кажется что он писал так, как было принято у нас. А не из за того, что он чет не догонял.

Наверно, так и есть. Но в современных книжках, всё же, обозначают уже одинаково. Поэтому это принято уже везде.

[bahurin]

Наверно, так и есть. Но в современных книжках, всё же, обозначают уже одинаково. Поэтому это принято уже везде.

Да ну не спорьте вы о книжках. Все они правильные. Когда я писал слово "принято" я скорее имел ввиду не литературу, а программные пакеты матлаб, октаве, скайлаб. Ибо автор вопроса на них опирается, но в его коде вроде как наоборот. Разумеется мы можем обозначить на бумаге любыми буквами, хоть иероглифами китайскими, главное, чтобы при программной реализации они не путались местами.

[TViT]

Вот описание функции:

Код

function cels=casc(x,z)
//cels=casc(x,z)
//Creates cascade realization of filter
//from a matrix of coefficients
//  x    :(4xN)-Matrix where each column is a cascade
//       :element, the first two column entries being
//       :the numerator coefficients and the second two
//       :column entries being the denominator coefficients
//  z    :Character string representing the cascade variable
//  cels :Resulting cascade representation
//
//EXAMPLE:
//  x=[ 1.     2.     3. ;
//      4.     5.     6. ;
//      7.     8.     9. ;
//      10.    11.    12. ]
//
//  cels=casc(x,'z')
//  cels      =
//
//  !             2               2               2  !
//  !   1 + 4z + z      2 + 5z + z      3 + 6z + z   !
//  !  ------------    ------------    ------------  !
//  !              2               2               2 !
//  !   7 + 10z + z     8 + 11z + z     9 + 12z + z  !
//!

cels=[];
for col=x,
      nf=[col(1:2);1];
      nd=[col(3:4);1];
      cels=[cels,syslin([],poly(nf,'z','c'),poly(nd,'z','c'))];
end,

endfunction

----------------------

Нужно как я понимаю сформировать правильно матрицу коэффициентов, а casc пересчитает их по секциям видимо 2-х порядковым, вот как и что нужно сделать, а главное что возвращается?

И что это за дроби она выдает в примере?

Код

//  !   1 + 4z + z      2 + 5z + z      3 + 6z + z   !
//  !  ------------    ------------    ------------  !
//  !              2               2               2 !
//  !   7 + 10z + z     8 + 11z + z     9 + 12z + z  !

Сообщение отредактировал TViT - Dec 23 2010, 16:52