[shf_05 0]

Добрый день, коллеги.

использую функцию group из пакета Scilab применительно к фильтру IIR, собранному из секций 2 порядка в ПО типа QED в формате 1,15.

сначала из коэффициентов фильтра формируются полиномы, потом один полином всего фильтра и из него - ГВЗ.

почему то ГВЗ строится с "выбросами", которых не должно быть и в добавок, форма ГФЗ и положение "выбросов" зависит от числа точек.

попытка привести коэф. к типу double делением коэффициентов на 32767 вместо 32768 дает еще более корявый результат.

 

подскажите, как еще можно посчитать ГВЗ фильтра, чтобы не было проблем с точностью.

 

будет ли подобное в матлабе?

 

листинг:

clc  
clear

//Y(n)= b0*x(n)+b1*x(n-1)+b2*x(n-2)-a1*y(n-1)-a2*y(n-2)
Fs= 32000;

M= [
                79 //*      4F      section  1 coefficient B0 */
        -94 //*FFFFFFA2      section  1 coefficient B1 */
         79 //*      4F      section  1 coefficient B2 */
     -14930 //*FFFFC5AE      section  1 coefficient -A2*/
      31250 //*    7A12      section  1 coefficient -A1*/       
       1696 //*     6A0      section  2 coefficient B0 */
      -3164 //*FFFFF3A4      section  2 coefficient B1 */
       1696 //*     6A0      section  2 coefficient B2 */
     -15352 //*FFFFC408      section  2 coefficient -A2*/
      31503 //*    7B0F      section  2 coefficient -A1*/         
       5664 //*    1620      section  3 coefficient B0 */
     -10924 //*FFFFD554      section  3 coefficient B1 */
       5664 //*    1620      section  3 coefficient B2 */
     -15842 //*FFFFC21E      section  3 coefficient -A2*/
      31813 //*    7C45      section  3 coefficient -A1*/          
       8845 //*    228D      section  4 coefficient B0 */
     -17189 //*FFFFBCDB      section  4 coefficient B1 */
       8845 //*    228D      section  4 coefficient B2 */
     -16222 //*FFFFC0A2      section  4 coefficient -A2*/
      32094 //*    7D5E      section  4 coefficient -A1*/
] / 32768;
G = [2 2 2 2];

SosSign= -1;

Gall= 1;    //gain for sections

Nsect= length(M)./5;
if (round(Nsect)-Nsect) ~= 0 then
    disp ('error - input matrix dimension !')
    abort;
end

Msos= matrix(M,5,Nsect);

disp(Msos)

Num= Msos(3:-1:1 , :);
Den= zeros(Num);   // [1/g1,-Msos(5),-Msos(4) , :];
Den(1,:)= 1 ./G;
Den(2:3,:)=   SosSign * Msos(5:-1:4,:);

disp(Num,'num')
disp(Den,'den')

z= poly(0,"z");

Psos= cell(1,Nsect);
Psos2= cell(1,Nsect);

//make Polynoms
for k = 1:Nsect
    Psos(k).entries = poly(Num(:,k),'z','c') ./ poly(Den(:,k),'z','c');
end
disp(Psos.entries,'Sos polynoms:')

// convert poly Z to Z^-1
for k = 1:Nsect
   Psos2(k).entries = horner(Psos(k).entries , 1/z);
end
disp(Psos2.entries,'Sos2 polynoms:')

// All polynoms in one
Pall = poly(%inf,"z")
for k = 1:Nsect
   Pall= Pall*Psos2(k).entries;
end
Pall= Gall*Pall;
disp(Pall,'All polynom:')



scf(0);
for k=1:Nsect
    [hzm,fr]=frmag(Psos2(k).entries,4096);// freq responce for each section
    plot(fr*Fs,20*log10(hzm),'--');   
end
    set(gca(),"grid",[1 1])
// freq responce for all sections
    [hzm,fr]=frmag(Pall,4096);
    plot(fr*Fs,20*log10(hzm),'r');
    title('Freq responce');

// freq magnitude responce for h = syslin
    h= syslin('d',Pall)

// time responces
scf(2);clf();
set(gca(),"grid",[1 1])  
    imprep=flts(eye(1, 1024),tf2ss(h));   //Impulse response
    plot([1:1024]/Fs,imprep./max(imprep),'b')
    
    stprep=flts(ones(1,1024),tf2ss(h));   //Step response
    plot([1:1024]/Fs,stprep,'g')
    legend(['Imp';'Step'])
    title('Time responce for all sections');

//group delay
//h= horner(poly1,z) // convert poly Z to Z^-1
scf(3);clf();
set(gca(),"grid",[1 1])
    [tg, fr] = group(512, h);
    plot(fr*Fs, tg/Fs,'.-')
    title('Group delay');

[Самурай 12]

Добрый день, коллеги.

использую функцию group из пакета Scilab применительно к фильтру IIR, собранному из секций 2 порядка в ПО типа QED в формате 1,15.

сначала из коэффициентов фильтра формируются полиномы, потом один полином всего фильтра и из него - ГВЗ.

почему то ГВЗ строится с "выбросами", которых не должно быть и в добавок, форма ГФЗ и положение "выбросов" зависит от числа точек.

попытка привести коэф. к типу double делением коэффициентов на 32767 вместо 32768 дает еще более корявый результат.

 

подскажите, как еще можно посчитать ГВЗ фильтра, чтобы не было проблем с точностью.

 

будет ли подобное в матлабе?

 

За весь Матлаб не скажу, но вот например функция grpdelay работает без нареканий. Опять же, цели провести полномасштабное тестирование вроде как не было, но на тех фильтрах, что смотрел, "выбросов" на ГВЗ не заметил.

 

Что не так в Scilab и/или у Вас в коде судить не берусь, про Scilab слышал только название, картинки с аномалиями ГВЗ тут точно бы не помешали. Попробуйте сначала на фильтре 2-ого порядка проверить.

 

А по теории можно например вот тут почитать, там же рядом исходный код функции grpdelay.

 

[shf_05 0]

в теории то все красиво, вот например на практике функция скилаба "корявит" - на первом рисунке 512 точек, на втором 128 точек - вычисленная group задержка и на третьем - производная фазы (разрывы функции фазы дают выбросы, не обращаю внимания)

ПС - проверено матлабом ГВЗ без всяких выбросов, строится красиво