использую функцию 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');
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');
