group delay, ошибка вычисления или погрешность Опции

[shf_05]

Добрый день, коллеги.
использую функцию 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');