以下は、SCILABの eqiir または iir
を使って、IIRフィルターの係数を求めるもであるが、次数が大きくなると動作が不安定で上手く動かないことがある。 また、直列の2次のIIRフィル
ターで分解した表示のcellsは、計算精度を考慮した最適な極と零の組み
合わせにはなっていない。
//-------------------------------------------------------------------------------------------
// Analog to IIR digital filter design
//
//
//-------------------------------------------------------------------------------------------
// IIR digital filter design
// by scilab eqiir
//
//-------------------------------------------------------------------------------------------
//
// IIR digital de-emphasis filter design
// as almost same as analog de-emphasis filter 50/15uS
//
// for window scilab-4.1.2
//
//-------------------------------------------------------------------------------------------
// PLEASE PAY ATTENTION that this program may have some bugs and also
you may modify program
// to fit your circumstances.
// If you use this program, everything should be done at your own risk.
// Thanks.
//===========================================================================================
// global variables No. 1
c50=50.0; // constant 1
c15=15.0; // constant 2
fs=44100.0; // sampling frequency
iiro=4; // oder of iir filter
start_freq=20.0;
end_freq=20000.0;
step_freq=1000.0; // many is good
scale_mode=0; // 0 linear, 1 log
trans_mode=1; // 0 least-square filter design, 1 bilinear
transform (= SOU SENKEI HENKAN in jappanese ) filter design
c15_change=0.0; // c15 = c15 * (1.0 + c15_change/100) This is available
only when trans_mode=1
//
z=poly(0,'z');
H1=z/(z-.5);
H0=syslin('d',H1); // dummy set, digital discrete time system
H1=syslin('d',H1); // dummy set, digital discrete time system
s=poly(0,'s');
S1=(1+2*s)/s^2;
S1A=syslin('c',S1); // dummy set, analog continuous time
system
//
zisuu=2;
ftype='bp';
approx='butt';
omf=zeros(1,4);
omf(1)= 1000.0;
omf(2)= 2000.0;
omf(3)= 3000.0;
omf(4)= 4000.0;
deltap=-0.5;
deltas=-30.0;
fact=1.; // dummy set,
fact1=1.;
zzeros=zeros(1,1); // dummy set,
zpoles=zeros(1,1); // dummy set
cells=[]; // fact
cells1=[]; // including fact
//-------------------------------------------------------------------------------------------
//
DE_EMPHASIS=0; // when DE_EMPHASIS is 1, de_emphasis menu
appears.
ANA2IIR=0;
EQIIRR=0;
NORM1=1;
SEL_CODE=x_choose(['de_emphasis design'; 'eqiir design'; 'iir design';
'analog2iir design'],['Please select one'],'default');
if SEL_CODE == 1 then
DE_EMPHASIS=1; // when DE_EMPHASIS is 1, de_emphasis
menu appears.
elseif SEL_CODE == 3 then
ANA2IIR=1;
NORM1=3;
elseif SEL_CODE == 4 then
ANA2IIR=1;
NORM1=1;
else
EQIIRR=1;
end
//
//-------------------------------------------------------------------------------------------
// choose trans_mode
//
if DE_EMPHASIS == 1 then
SEL_CODE=x_choose(['bilinear transform filter design'; 'SOU SENKEI
HENKAN + HOSEI filter design'; 'least-square filter desin'],['Please
select one'],'default');
if SEL_CODE == 1 then
trans_mode=1;
elseif SEL_CODE == 2 then
trans_mode=1;
c15_change=12.5;
disp('message: constant 2 correction factor is set 12.5');
elseif SEL_CODE == 3 then
trans_mode=0;
else
trans_mode=1;
end
if trans_mode == 0 then
disp('least-square filter design');
elseif trans_mode == 1 then
disp('bilinear transform (= SOU SENKEI HENKAN in jappanese )
filter design');
end
else
SEL_CODE = 2;
trans_mode=1;
c15_change=12.5;
end
//-------------------------------------------------------------------------------------------
function [wapprox]=selectapprox()
SEL_CODE=x_choose(['Butterworth filter' ; 'Chebyshev type 1
filter'],['Please select one'],'default');
if SEL_CODE == 1 then
wapprox='butt';
elseif SEL_CODE == 2 then
wapprox='cheb1';
else
wapprox='butt';
end
endfunction
//-------------------------------------------------------------------------------------------
function [wftype]=selectftype()
SEL_CODE=x_choose(['LPF'; 'HPF'; 'BPF'; 'SBP'],['Please select
one'],'default');
if SEL_CODE == 1 then
wftype='lp';
elseif SEL_CODE == 2 then
wftype='hp';
elseif SEL_CODE == 3 then
wftype='bp';
elseif SEL_CODE == 4 then
wftype='sb';
else
wftype='lp';
end
endfunction
//-------------------------------------------------------------------------------------------
function [wfs,wom,wdeltap,wdeltas]= set_constant2()
txt1=['sampling frequency [Hz]' ; 'cutoff frequency1 [Hz]';
'cutoff frequency2 [Hz]'; 'cutoff frequency3 [Hz]';'cutoff frequency4
[Hz]';'gain at edge or ripple in the passband [dB]' ; ' gain in the
stopband [dB]'];
wstr1=sprintf('%f',fs);
wstr2=sprintf('%f',omf(1));
wstr3=sprintf('%f',omf(2));
wstr4=sprintf('%f',omf(3));
wstr5=sprintf('%f',omf(4));
wstr6=sprintf('%f',deltap);
wstr7=sprintf('%f',deltas);
sig1=x_mdialog('Input constant ',txt1, [wstr1 ; wstr2 ; wstr3;
wstr4; wstr5; wstr6; wstr7 ]);
if sig1==[] then
arg1=evstr(wstr1);
arg2=evstr(wstr2);
arg3=evstr(wstr3);
arg4=evstr(wstr4);
arg5=evstr(wstr5);
arg6=evstr(wstr6);
arg7=evstr(wstr7);
else
arg1=evstr(sig1(1));
arg2=evstr(sig1(2));
arg3=evstr(sig1(3));
arg4=evstr(sig1(4));
arg5=evstr(sig1(5));
arg6=evstr(sig1(6));
arg7=evstr(sig1(7));
end
wom=zeros(1,4);
wfs=arg1;
wom(1)=arg2;
wom(2)=arg3;
wom(3)=arg4;
wom(4)=arg5;
wdeltap=arg6;
wdeltas=arg7;
endfunction
//-------------------------------------------------------------------------------------------
function [wzisuu,wfs,wom,wdeltap]= set_constant3()
txt1=['filter order (integer): even when BPF and SBP'; 'sampling
frequency [Hz]' ; 'cutoff frequency1 [Hz]'; 'cutoff frequency2 [Hz]';
'ripple in the passband [dB]' ];
wstr1=sprintf('%d',zisuu);
wstr2=sprintf('%f',fs);
wstr3=sprintf('%f',omf(1));
wstr4=sprintf('%f',omf(2));
wstr5=sprintf('%f',deltap);
sig1=x_mdialog('Input constant ',txt1, [wstr1 ; wstr2 ; wstr3;
wstr4; wstr5]);
if sig1==[] then
arg1=evstr(wstr1);
arg2=evstr(wstr2);
arg3=evstr(wstr3);
arg4=evstr(wstr4);
arg5=evstr(wstr5);
else
arg1=evstr(sig1(1));
arg2=evstr(sig1(2));
arg3=evstr(sig1(3));
arg4=evstr(sig1(4));
arg5=evstr(sig1(5));
end
wom=zeros(1,4);
wzisuu=arg1;
wfs=arg2;
wom(1)=arg3;
wom(2)=arg4;
wdeltap=arg5;
endfunction
//-------------------------------------------------------------------------------------------
function [wH0A, wS1A, wiiro2 ]=filterdesign3()
//
//-----------------------------
if ftype== 'bp' then
zisuu=zisuu /2;
elseif ftype == 'sb' then
zisuu=zisuu /2 ;
else
zisuu=zisuu;
end
if zisuu < 1 then
disp('ERROR!: zisuu is less than 1. Force to
set as 1');
zisuu=1;
end
//+++++++++++++++++++++++++++
wdeltap=0.0;
if approx=='cheb1' then
//
a=(abs(deltap) / 20);
wdeltap= 10^a;
wdeltap=wdeltap-1.0;
[poles,gain]=zpch1(zisuu,wdeltap,1.0);
wS0=gain/real(poly(poles,'s'));
//
// if wapprox=='cheb1' then
elseif approx=='butt' then
//
[pols,gain]=zpbutt(zisuu,1.0);
wS0=gain/real(poly(pols,'s'))
end // elseif wapprox=='butt' then
//++++++++++++++++++++++++++++
wS1=syslin('c',wS0);
wH0A=[];
if ftype== 'lp' then
w0= 2.0 * %pi * omf(1);
wS1B=horner( wS1,s/w0);
wS1A=syslin('c',wS1B);
w0= 2.0 * tan( (omf(1) / fs * 2. * %pi) / 2.0);
w0= fs * (w0 / (2.0 * %pi));
w0= 2.0 * %pi * w0;
wS1B=horner( wS1,s/w0);
wS1A2=syslin('c',wS1B);
wH0= horner(wS1A2,2*fs*(z-1)/(z+1));
wH0A= syslin('d', wH0);
elseif ftype== 'hp' then
w0= 2.0 * %pi * omf(1);
wS1B=horner( wS1,w0/s);
wS1A=syslin('c',wS1B);
w0= 2.0 * tan( (omf(1) / fs * 2. * %pi) / 2.0);
w0= fs * (w0 / (2.0 * %pi));
w0= 2.0 * %pi * w0;
wS1B=horner( wS1,w0/s);
wS1A2=syslin('c',wS1B);
wH0= horner(wS1A2,2*fs*(z-1)/(z+1));
wH0A= syslin('d', wH0);
elseif ftype== 'bp' then
b= 2.0 * %pi * (omf(2) - omf(1));
w0= 2.0 * %pi * ( omf(1) * omf(2))^0.5;
wS1B=horner( wS1,(s^2 + w0 * w0) / (b * s));
wS1A=syslin('c',wS1B);
w1= 2.0 * tan( (omf(1) / fs * 2. * %pi) / 2.0);
w1= fs * (w1 / (2.0 * %pi));
w2= 2.0 * tan( (omf(2) / fs * 2. * %pi) / 2.0);
w2= fs * (w2 / (2.0 * %pi));
w0= 2.0 * %pi * ( w1 * w2)^0.5;
b= 2.0 * %pi * (w2 - w1);
wS1B=horner( wS1,(s^2 + w0 * w0) / (b * s));
wS1A2=syslin('c',wS1B);
wH0= horner(wS1A2,2*fs*(z-1)/(z+1));
wH0A= syslin('d', wH0);
elseif ftype== 'sb' then
b= 2.0 * %pi * (omf(2) - omf(1));
w0= 2.0 * %pi * ( omf(1) * omf(2))^0.5;
wS1B=horner( wS1, (b * s)/ (s^2 + w0 * w0) );
wS1A=syslin('c',wS1B);
w1= 2.0 * tan( (omf(1) / fs * 2. * %pi) / 2.0);
w1= fs * (w1 / (2.0 * %pi));
w2= 2.0 * tan( (omf(2) / fs * 2. * %pi) / 2.0);
w2= fs * (w2 / (2.0 * %pi));
w0= 2.0 * %pi * ( w1 * w2)^0.5;
b= 2.0 * %pi * (w2 - w1);
wS1B=horner( wS1,(b * s)/ (s^2 + w0 * w0) );
wS1A2=syslin('c',wS1B);
wH0= horner(wS1A2,2*fs*(z-1)/(z+1));
wH0A= syslin('d', wH0);
end
//+++++++++++++++++++++++++++
if NORM1 == 1 then
[H0r,H0nx]=syssize( wH0A ); // get size
wiiro2=H0nx;
H0bunbo=denom( wH0A ); // get BUNBO
H0bunsi=numer( wH0A ); // get BUNSHI
[wcoeff_bunbo]=get_coeff( H0bunbo, wiiro2 );
[wcoeff_bunsi]=get_coeff( H0bunsi, wiiro2 );
H0bunbo=denom( wH0A )/wcoeff_bunbo(1,1);
H0bunsi=numer( wH0A )/wcoeff_bunbo(1,1);
z=poly(0,'z');
H1=H0bunsi/H0bunbo;
wH0A=syslin('d',H1);
elseif NORM1 == 3 then // call iir
disp('message: call iir of scilab');
fx=zeros(1,2);
fx(1,1)=omf(1)/fs;
fx(1,2)=omf(2)/fs;
dx=zeros(1,2);
dx(1,1)=wdeltap;
[wH0A]=iir(zisuu,ftype,approx,fx,dx);
end
[H0r,H0nx]=syssize( wH0A ); // get size
wiiro2=H0nx;
H0bunbo=denom( wH0A ); // get BUNBO
H0bunsi=numer( wH0A ); // get BUNSHI
//H0poles=roots( H0bunbo ); // get poles
//H0zeros=roots( H0bunsi ); // get zeros
[wcoeff_bunbo]=get_coeff( H0bunbo, wiiro2 );
[wcoeff_bunsi]=get_coeff( H0bunsi, wiiro2 );
//plzr( wH0 ); // plot zeros and poles
//
global ac0
global bcb0
global bc0
//
//disp('digital discrete time system');
//disp( H0nx );
//disp( '*bunbo/poles/x*');
//disp( H0bunbo);
//disp( wcoeff_bunbo(1), 'a[0] is always 1' );
for v=1:wiiro2
ac0(v)= - wcoeff_bunbo(v+1);
wstr1= 'a[' + string(v) + ']';
//disp( ac0(v), wstr1 );
end
//disp('poles are');
//disp( H0poles);
//disp( '*bunsi/zeros/o*')
//disp( H0bunsi);
bcb0=wcoeff_bunsi(1);
//disp( wcoeff_bunsi(1), 'b[0]' );
for v=1:wiiro2
bc0(v)= wcoeff_bunsi(v+1);
wstr1= 'b[' + string(v) + ']';
//disp( bc0(v), wstr1 );
end
//disp('zeros are');
//disp( H0zeros);
endfunction
//-------------------------------------------------------------------------------------------
function [wout]=xsort(xin)
// wout=gsort(xin);
sz=size(xin);
sz2=sz(1);
wout1=zeros(sz2);
wout=zeros(sz2);
for v=1:sz2
wout1(v)=abs(xin(v));
if real(xin(v)) < 0. then
wout1(v)= -1.0 * wout1(v);
end
end
[wout3,k]=gsort(wout1);
for v=1:sz2
wout(v)= xin( k(v));
end
endfunction
//------
function [wout]=xsort2(xin)
// this is increase order
sz=size(xin);
sz2=sz(1);
wout1=zeros(sz2);
wout=zeros(sz2);
for v=1:sz2
wout1(v)=abs(xin(v));
if real(xin(v)) < 0. then
wout1(v)= -1.0 * wout1(v);
end
end
[wout3,k]=gsort(wout1,'g','i');
for v=1:sz2
wout(v)= xin( k(v));
end
endfunction
//------
function [wout]=xsort3(xin)
// this is imaginary order
sz=size(xin);
sz2=sz(1);
wout1=zeros(sz2);
wout=zeros(sz2);
for v=1:sz2
wout1(v)=abs(imag(xin(v)));
if real(xin(v)) < 0. then
wout1(v)= -1.0 * wout1(v);
end
end
[wout3,k]=gsort(wout1);
for v=1:sz2
wout(v)= xin( k(v));
end
endfunction
//-------------------------------------------------------------------------------------------
function [wcells,wcells1,fact,fact1,H0poles,H0zeros]=get_cell(wH0)
//DEGWH0=1;
DEGWH0=0;
if DEGWH0 == 1 then
wH0=H0;
end
[H0r,H0nx]=syssize( wH0 ); // get size
wiiro2=H0nx;
wiiro3=int(wiiro2/2);
wiiro4=wiiro2 - 2 * wiiro3;
H0bunbo=denom( wH0 ); // get BUNBO
H0bunsi=numer( wH0 ); // get BUNSHI
H0poles=roots( H0bunbo ); // get poles
sz=size(H0poles);
H0poles2=zeros(sz(1),1);
flag=zeros(3,sz(1));
if (ftype=='bp') | (ftype=='sb') then
// like almost frequency shift order ???
H0poles2=xsort3(H0poles);
else
// ordering 1, 0.9, ..., 0, -0.1 ,...-0.9, -1
H0poles2=xsort(H0poles);
end
// check pair 1
for v=1:sz(1)
// check real ?
if imag(H0poles2(v)) == 0. then
flag(1,v)=1; // real
else
if (v < sz(1)) & ( flag(1,v) == 0) then
if ( real(H0poles2(v)) ==
real(H0poles2(v+1)) ) & ( imag(H0poles2(v)) == -1.0 *
imag(H0poles2(v+1)) ) then
flag(1,v)=2; // complex
flag(2,v)=v+1;
flag(1,v+1)=2;
flag(2,v+1)=v;
end // found pair
end
end // is real ?
end // chech pair 1
// check pair 2
for v=1:sz(1)
amari=modulo(v,2);
if (flag(1,v)==1) & (flag(2,v) == 0) then // is there a pair
?
flag(2,v)=-1; // there is No pair, single
for v2=(v+1):sz(1)
if flag(1,v2) == 1 then
flag(2,v) = v2;
flag(2,v2)= v;
break;
end
end
end
end // check pair 2
// check flag(*,v)==-1
moved_data_position=-1;
c0=0;
for v=1:sz(1)
if flag(3,v) == 0 then //not yet marked
if flag(2,v)== -1 then // when there is No pair, single,
it will move end of data
H0poles(sz(1))=H0poles2(v);
flag(3,v)=9; // marked
moved_data_position=c0+1;
disp(moved_data_position,'moved_data_position=');
else
c0=c0+1;
H0poles(c0)=H0poles2(v);
flag(3,v)=9; // marked
c0=c0+1;
H0poles(c0)=H0poles2( flag(2,v)); // set pair
flag(3,flag(2,v))=9; // marked
end
end
end
H0zeros=roots( H0bunsi ); // get zeros
sz=size(H0zeros);
H0zeros2=zeros(sz(1),1);
flag=zeros(3,sz(1));
// ordering -1,-0.9, ..., -0, 0.1 ,...0.9, 1
// from lpf to bpf to hpf
H0zeros2=xsort2(H0zeros);
// check pair 1
for v=1:sz(1)
// check real ?
if imag(H0zeros2(v)) == 0. then
flag(1,v)=1; // real
else
if (v < sz(1)) & ( flag(1,v) == 0) then
if ( real(H0zeros2(v)) ==
real(H0zeros2(v+1)) ) & ( imag(H0zeros2(v)) == -1.0 *
imag(H0zeros2(v+1)) ) then
flag(1,v)=2; // complex
flag(2,v)=v+1;
flag(1,v+1)=2;
flag(2,v+1)=v;
end // found pair
end
end // is real ?
end // chech pair 1
// check pair 2
for v=1:sz(1)
amari=modulo(v,2);
if (flag(1,v)==1) & (flag(2,v) == 0) then // is there a pair
?
flag(2,v)=-1; // there is No pair, single
for v2=(v+1):sz(1)
if flag(1,v2) == 1 then
flag(2,v) = v2;
flag(2,v2)= v;
break;
end
end
end
end // check pair 2
// check flag(*,v)==-1
moved_data_position2=-1;
c0=0;
for v=1:sz(1)
if flag(3,v) == 0 then //not yet marked
if flag(2,v)== -1 then // when there is No pair, single,
it will move end of data
H0zeros(sz(1))=H0zeros2(v);
flag(3,v)=9; // marked
moved_data_position2=c0+1;
disp(moved_data_position2,'moved_data_position2=');
else
c0=c0+1;
H0zeros(c0)=H0zeros2(v);
flag(3,v)=9; // marked
c0=c0+1;
H0zeros(c0)=H0zeros2( flag(2,v)); // set pair
flag(3,flag(2,v))=9; // marked
end
end
end
//fact=horner( H0bunsi1, 0.);
[wcoeff_bunsi]=get_coeff( H0bunsi, wiiro2 );
fact=wcoeff_bunsi(1);
global ac2
global bcb2
global bc2
ac2=zeros(mc2,2); // coefficient matrix of a[*][1,2]
bcb2=zeros(mc2); // coefficient matrix of b[*][0]
bc2=zeros(mc2,2); // coefficient matrix of b[*][1,2]
a=wiiro3+wiiro4;
a=1/a;
fact1=fact^a;
mcell=zeros(4,(wiiro3+wiiro4));
for v=1:wiiro3
w2=2*v-1;
H0nx2=2;
wcoeff_bunbo=zeros(1,3);
wcoeff_bunsi=zeros(1,3);
wcoeff_bunbo(1,1)=1.0;
wcoeff_bunbo(1,2)= 1.0 * (H0poles(w2) + H0poles(w2+1));
wcoeff_bunbo(1,3)= -1.0 * H0poles(w2) * H0poles(w2+1);
mcell(4,v)=-1.0 * wcoeff_bunbo(1,2);
mcell(3,v)=-1.0 * wcoeff_bunbo(1,3);
wcoeff_bunsi(1,1)=1.0 ; // * fact1;
wcoeff_bunsi(1,2)=-1.0 * (H0zeros(w2) + H0zeros(w2+1)); // *
fact1;
wcoeff_bunsi(1,3)= H0zeros(w2) * H0zeros(w2+1); // * fact1;
mcell(2,v)=-1.0 * (H0zeros(w2) + H0zeros(w2+1));
mcell(1,v)= H0zeros(w2) * H0zeros(w2+1);
for w=1:H0nx2
ac2(v,w)= 1.0 * wcoeff_bunbo(w+1);
wstr1= 'a[' + string(v) +','+ string(w) + ']';
//disp( ac2(v,w), wstr1 );
end
bcb2(v)=wcoeff_bunsi(1) * fact1;
//disp( bcb2(v), 'bcb2[0]' );
for w=1:H0nx2
bc2(v,w)= wcoeff_bunsi(w+1) * fact1;
wstr1= 'b[' + string(v) +',' + string(v) + ']';
//disp( bc2(v,w), wstr1 );
end
end // for v=1:sz(2)
if wiiro4 == 1 then
w2=wiiro2;
w3=wiiro3+1;
H0nx2=2;
wcoeff_bunbo=zeros(1,3);
wcoeff_bunsi=zeros(1,3);
wcoeff_bunbo(1,1)=1.0;
wcoeff_bunbo(1,2)= 1.0 * H0poles(w2);
wcoeff_bunbo(1,3)=0.0;
mcell(4,w3)=-1.0 * wcoeff_bunbo(1,2);
wcoeff_bunsi(1,1)=1.0 ; // * fact1;
wcoeff_bunsi(1,2)= -1.0 * (H0zeros(w2)); // * fact1;
wcoeff_bunsi(1,3)=0.0;
mcell(2,w3)= -1.0 * H0zeros(w2);
v=wiiro3+1;
for w=1:H0nx2
ac2(v,w)= wcoeff_bunbo(w+1) ;
wstr1= 'a[' + string(v) +','+ string(w) + ']';
//disp( ac2(v,w), wstr1 );
end
bcb2(v)=wcoeff_bunsi(1) * fact1;
//disp( bcb2(v), 'bcb2[0]' );
for w=1:H0nx2
bc2(v,w)= wcoeff_bunsi(w+1) * fact1;
wstr1= 'b[' + string(v) +',' + string(v) + ']';
//disp( bc2(v,w), wstr1 );
end
end
//
//
[wcells]=casc(mcell,z);
wcells(2)=real(wcells(2));
wcells(3)=real(wcells(3));
if wiiro4 == 1 then
wcells(1,(wiiro3+1))=wcells(1,(wiiro3+1))/ (z/z);
end
wcells1=wcells * fact1;
endfunction
//-------------------------------------------------------------------------------------------
function [wH0, cells,cells1,
fact,fact1,zzeros,zpoles,wiiro2]=call_eqiir()
om=zeros(1,4);
om(1)= omf(1) / fs * ( 2.0 * %pi);
om(2)= omf(2) / fs * ( 2.0 * %pi);
om(3)= omf(3) / fs * ( 2.0 * %pi);
om(4)= omf(4) / fs * ( 2.0 * %pi);
a=(abs(deltap) / 20);
wdeltap= 10^a;
wdeltap=wdeltap-1.0;
a=(deltas / 20);
wdeltas= 10^a;
disp( wdeltas, ' wdeltas= ', wdeltap, ' wdeltap= ');
[cells,fact,zzeros,zpoles]=eqiir(ftype,approx,om,wdeltap,wdeltas);
H1=fact*poly(zzeros,'z')/poly(zpoles,'z'); // <--- unstable
????
wH0=syslin('d',H1);
sz=size(cells);
H1=cells(1,1);
for v=2:sz(2)
H1 = H1 * cells(1,v);
end
H1=H1 * fact;
wH0=syslin('d',H1);
disp( cells );
disp( fact,'fact=')
sz=size(cells);
a=sz(2);
a=1/a;
fact1=fact^a;
disp( fact1,'fact1=')
cells1=cells * fact1;
disp( cells1,'cells1=')
disp(wH0,'H0=');
[H0r,H0nx]=syssize( wH0 ); // get size
wiiro2=H0nx;
sz=size(cells1);
//disp('poles x are');
//disp( zpoles );
//disp( 'zeros o are ')
//disp( zzeros);
//plzr( wH0 ); // plot zeros and poles
//
global ac2
global bcb2
global bc2
ac2=zeros(mc2,2); // coefficient matrix of a[*][1,2]
bcb2=zeros(mc2); // coefficient matrix of b[*][0]
bc2=zeros(mc2,2); // coefficient matrix of b[*][1,2]
for v=1:sz(2)
wH1=cells1(1,v);
H0bunbo=denom( wH1 ); // get BUNBO
H0bunsi=numer( wH1 ); // get BUNSHI
[H0r2,H0nx2]=syssize( wH1 );
[wcoeff_bunbo]=get_coeff( H0bunbo, H0nx2 );
[wcoeff_bunsi]=get_coeff( H0bunsi, H0nx2 );
for w=1:H0nx2
ac2(v,w)= - wcoeff_bunbo(w+1);
wstr1= 'a[' + string(v) +','+ string(w) + ']';
//disp( ac2(v,w), wstr1 );
end
bcb2(v)=wcoeff_bunsi(1);
//disp( bcb2(v), 'bcb2[0]' );
for w=1:H0nx2
bc2(v,w)= wcoeff_bunsi(w+1);
wstr1= 'b[' + string(v) +',' + string(v) + ']';
//disp( bc2(v,w), wstr1 );
end
end // for v=1:sz(2)
//
//
//[H0r,H0nx]=syssize( wH0 ); // get size
//wiiro2=H0nx;
H0bunbo=denom( wH0 ); // get BUNBO
H0bunsi=numer( wH0 ); // get BUNSHI
//H0poles=roots( H0bunbo ); // get poles
//H0zeros=roots( H0bunsi ); // get zeros
[wcoeff_bunbo]=get_coeff( H0bunbo, wiiro2 );
[wcoeff_bunsi]=get_coeff( H0bunsi, wiiro2 );
//plzr( wH0 ); // plot zeros and poles
//
global ac0
global bcb0
global bc0
//
//disp('digital discrete time system');
//disp( H0nx );
//disp( '*bunbo/poles/x*');
//disp( H0bunbo);
//disp( wcoeff_bunbo(1), 'a[0] is always 1' );
for v=1:wiiro2
ac0(v)= - wcoeff_bunbo(v+1);
wstr1= 'a[' + string(v) + ']';
//disp( ac0(v), wstr1 );
end
//disp('poles are');
//disp( H0poles);
//disp( '*bunsi/zeros/o*')
//disp( H0bunsi);
bcb0=wcoeff_bunsi(1);
//disp( wcoeff_bunsi(1), 'b[0]' );
for v=1:wiiro2
bc0(v)= wcoeff_bunsi(v+1);
wstr1= 'b[' + string(v) + ']';
//disp( bc0(v), wstr1 );
end
//disp('zeros are');
//disp( H0zeros);
zzeros=zzeros';
zpoles=zpoles';
endfunction
//---------------------------------------------------------------------------------------
function freq_response2(disp0, H0, str )
//
wb0=xget('window'); // stack old window
xset('window',disp0); // create new windows
clf();
//
[freqa, freqd, kosuu]=get_freqs();
[frq,respf]=repfreq( H0, freqd );
[db_fft2,phi_fft2] =dbphi(respf);
gainplot(freqa,db_fft2,phi_fft2);
wstr1 = 'IIR digital filter frequency
response Hz ' + str;
xtitle( wstr1,'','dB' );
//
xset('window',wb0); // push old windows
endfunction
//
//-------------------------------------------------------------------------------------------
function [wc50, wc15, wfs, wiiro, wc15_change]= set_constant()
if trans_mode==0 then
// 0 least-square filter design,
txt1=['constant 1 (default 50us)';'constant 2 (default 15us)';
'sampling frequency (default 44.1KHz)'; 'order of IIR filter'];
wstr1=sprintf('%f',c50);
wstr2=sprintf('%f',c15);
wstr3=sprintf('%f',fs);
wstr4=sprintf('%d',iiro);
sig1=x_mdialog('Input constant of analog de-emphasis filter,
sampling frequency and IIR order',txt1, [wstr1 ; wstr2 ; wstr3 ; wstr4
]);
if sig1==[] then
arg1=evstr(wstr1);
arg2=evstr(wstr2);
arg3=evstr(wstr3);
arg4=evstr(wstr4);
else
arg1=evstr(sig1(1));
arg2=evstr(sig1(2));
arg3=evstr(sig1(3));
arg4=evstr(sig1(4));
end
wc50=arg1;
wc15=arg2;
wfs=arg3;
wiiro=arg4;
wc15_change=0.0;
end //if trans_mode==0 then
//
//
//
if trans_mode==1 then
// 1 bilinear transform filter design
txt1=['constant 1 (default 50us)';'constant 2 (default 15us)';
'sampling frequency (default 44.1KHz)'; 'constant 2 correction factor
[%] (ex: 12.5)'];
wstr1=sprintf('%f',c50);
wstr2=sprintf('%f',c15);
wstr3=sprintf('%f',fs);
wstr4=sprintf('%f',c15_change);
sig1=x_mdialog('Input constant of analog de-emphasis filter,
sampling frequency and constant 2 correction factor ',txt1, [wstr1 ;
wstr2 ; wstr3; wstr4 ]);
if sig1==[] then
arg1=evstr(wstr1);
arg2=evstr(wstr2);
arg3=evstr(wstr3);
arg4=evstr(wstr4);
else
arg1=evstr(sig1(1));
arg2=evstr(sig1(2));
arg3=evstr(sig1(3));
arg4=evstr(sig1(4));
end
wc50=arg1;
wc15=arg2;
wfs=arg3;
wiiro=0;
wc15_change=arg4;
end //if trans_mode==1 then
endfunction
//---------------------------------------------------------------------------------------------------------------
function factor0=correct_factor( wc15)
//
fc=1.0 / ( (wc15 * 0.000001) * 2.0 * %pi );
wd= (fc / (fs /2.0)) * %pi;
fcc=tan( wd / 2.0) * 2.0 / ( 1. / fs) / ( 2.0 * %pi);
factor0= (fc /fcc - 1.0 ) * 100.0; // beacause unit is [%]
endfunction
//---------------------------------------------------------------------------------------------------------------
function [freqa, freqd, kosuu]=get_freqs()
//
mode0=scale_mode;
if end_freq > (1.1 * (fs/2.0)) then
end_freq= 0.9 * (fs/2.0);
end
if mode0==0 then // linear scale
step_freq=(end_freq-start_freq)/step_freq;
freqa=(start_freq:step_freq:end_freq);
freqd=((start_freq/fs):(step_freq/fs):(end_freq/fs));
end
if mode0==1 then // log scale
step0= (log(end_freq) - log(start_freq)) / step_freq;
bairitu0= exp( step0 );
freqa=zeros(1,2);
freqd=zeros(1,2);
freqa(1)=start_freq;
freqd(1)=freqa(1)/fs;
for v=1:int(step_freq)
freqa(v+1)=freqa(v) * bairitu0;
freqd(v+1)=freqa(v)/fs;
end
end
s0=size(freqa);
kosuu=s0(1,2);
endfunction
//---------------------------------------------------------------------------------------------------------------
function [wH0, wS1A, wiiro2]=filterdesign()
// analog_response
//
// dom='c' for a continuous time system,
// dom='d' for a discrete time system,
// n for a sampled system with sampling
period n (in seconds).
//
n-1
n-2
// B(z)
b(1)z + b(2)z + .... + b(n)
// H(z)= ---- = ---------------------------------
//
n-1 n-2
// A(z) z +
a(2)z + .... + a(n)
if trans_mode==0 then
// least-square filter design
//
s=poly(0,'s');
H1=(1+(c15 * 0.000001)*s)/(1+(c50 * 0.000001)*s),
wS1A=syslin('c',H1);
disp( wS1A );
//
[freqa, freqd, kosuu]=get_freqs();
[freqs, respa]=repfreq(wS1A, freqa);
freqz=zeros(1,2);
freqz(1)=0;
freqz(2)= fs /2.0;
[freqs, respz]=repfreq(wS1A, freqz);
wrespd=zeros(1,kosuu+2);
wfreqd=zeros(1,kosuu+2);
//
for v=1:kosuu
wrespd(v+1)= abs(respa(v));
wfreqd(v+1)= freqa(v) / (fs / 2.);
end
//
wrespd(1)=abs(respz(1));
wfreqd(1)=0.0;
wrespd(kosuu+2)=abs(respz(2));
wfreqd(kosuu+2)=1.0;
//
wH0=yulewalk(iiro,wfreqd,wrespd);
//
disp( wH0 );
end // if trans_mode==0 then
if trans_mode==1 then
// bilinear transform design
//
s=poly(0,'s');
H1=(1+(c15 * 0.000001)*s)/(1+(c50 * 0.000001)*s),
wS1A=syslin('c',H1);
disp( wS1A );
// bilinear transform from including correction factor
c15=c15 * ( 1.0 + c15_change/100.0);
H1=(1+(c15 * 0.000001)*s)/(1+(c50 * 0.000001)*s),
wS1A2=syslin('c',H1);
wS1Ass=tf2ss(wS1A2); //Now in state-space form
sl1=cls2dls(wS1Ass,1/fs); //sl1= output of cls2dls
wH0=ss2tf(sl1); // Converts in transfer form
disp( wH0 );
end // if trans_mode==1 then
[H0r,H0nx]=syssize( wH0 ); // get size
wiiro2=H0nx;
H0bunbo=denom( wH0 ); // get BUNBO
H0bunsi=numer( wH0 ); // get BUNSHI
//H0poles=roots( H0bunbo ); // get poles
//H0zeros=roots( H0bunsi ); // get zeros
[wcoeff_bunbo]=get_coeff( H0bunbo, wiiro2 );
[wcoeff_bunsi]=get_coeff( H0bunsi, wiiro2 );
//plzr( wH0 ); // plot zeros and poles
//
global ac0
global bcb0
global bc0
//
//disp('digital discrete time system');
//disp( H0nx );
//disp( '*bunbo/poles/x*');
//disp( H0bunbo);
//disp( wcoeff_bunbo(1), 'a[0] is always 1' );
for v=1:wiiro2
ac0(v)= - wcoeff_bunbo(v+1);
wstr1= 'a[' + string(v) + ']';
//disp( ac0(v), wstr1 );
end
//disp('poles are');
//disp( H0poles);
//disp( '*bunsi/zeros/o*')
//disp( H0bunsi);
bcb0=wcoeff_bunsi(1);
//disp( wcoeff_bunsi(1), 'b[0]' );
for v=1:wiiro2
bc0(v)= wcoeff_bunsi(v+1);
wstr1= 'b[' + string(v) + ']';
//disp( bc0(v), wstr1 );
end
//disp('zeros are');
//disp( H0zeros);
endfunction
//---------------------------------------------------------------------------------------------------------------
function freq_response(disp0, sub0, H0, S1A)
//
wb0=xget('window'); // stack old window
xset('window',disp0); // create new windows
clf();
//
subplot(311);
[freqa, freqd, kosuu]=get_freqs();
[freqs,respa]=repfreq(S1A, freqa);
[db_fft1,phi_fft1] =dbphi(respa);
gainplot(freqa,db_fft1,phi_fft1);
wstr1 = 'analog frequency
response
[Hz]';
if ANA2IIR == 1 then
wstr1 = wstr1 + ' filter order=' +
string(zisuu) + ' type=' + ftype + ' ' + approx;
end
xtitle( wstr1 ,'','dB');
subplot(312);
[frq,respf]=repfreq(H0, freqd );
[db_fft2,phi_fft2] =dbphi(respf);
gainplot(freqa,db_fft2,phi_fft2);
wstr1 = 'IIR digital filter frequency
response [Hz]';
if ANA2IIR == 1 then
wstr1 = wstr1 + ' filter order=' +
string(zisuu) + ' type=' + ftype + ' ' + approx;
end
xtitle( wstr1,'','dB' );
subplot(313);
gainplot(freqa,(db_fft2-db_fft1),(phi_fft2-phi_fft1));
wstr1 =
'difference
[Hz]';
xtitle( wstr1,'','dB' );
//
xset('window',wb0); // push old windows
endfunction
//===============================================================================================================
// check routine
//
//
iiro2=0; // order of the IIR filter
mc0=100; // maximum order IIR filter quantity
global ac0
global bcb0
global bc0
global wk0
ac0=zeros(1,mc0); // coefficient matrix of a[1,..,mc0]
bcb0=zeros(1,1); // coefficient matrix of b[0]
bc0=zeros(1,mc0); // coefficient matrix of b[1,..,mc0]
wk0=zeros(1,mc0*2); // data stack memory
mc1=1; // maximum 1st order IIR filter quantity
global ac1
global bcb1
global bc1
global wk1
ac1=zeros(mc1,1); // coefficient matrix a[*][1]
bcb1=zeros(mc1,1); // coefficient matrix of b[*][0]
bc1=zeros(mc1,1); // coefficient matrix of b[*][1]
wk1=zeros(mc1,2); // data stack memory
mc2=100; // maximum 2nd order IIR filter cascade quantity
global ac2
global bcb2
global bc2
global wk2
ac2=zeros(mc2,2); // coefficient matrix of a[*][1,2]
bcb2=zeros(mc2); // coefficient matrix of b[*][0]
bc2=zeros(mc2,2); // coefficient matrix of b[*][1,2]
wk2=zeros(mc2,4); // data stack memory
//
freqx=2500.0; // Hz
cyclex=30;
skip_cyclex=27;
//
sin_delta=0.1; // dummy set
sin_level=1.0; // sin's level when actual IIR calculation
//
in0=zeros(1,1); // dummy set
out0=zeros(1,1); // dummy set
m=0; // dummy set // size of input
xopt0=zeros(1,2); // dummy set
xopt1=zeros(1,2); // dummy set
//
//---------------------------------------------------------------------------------------------------------------
function [wcoeff]=get_coeff( wH0bunsi, wiiro2 )
//
z=poly(0,'z');
//
H0nx=wiiro2;
wcoeff=zeros(1,H0nx+1);
wcoeff(H0nx+1)= horner( wH0bunsi, 0.);
for v=1: H0nx
wH0bunsi=(wH0bunsi-wcoeff(H0nx+1-(v-1))) / z;
wcoeff(H0nx+1-v)= horner( wH0bunsi, 0.);
end
//
endfunction
//---------------------------------------------------------------------------------------------------------------
function [wiiro2]=system_stable_check( disp0, wH0, wS1A)
//
wb0=xget('window'); // stack old window
xset('window',disp0); // create new windows
clf();
[H0r,H0nx]=syssize( wH0 ); // get size
wiiro2=H0nx;
H0bunbo=denom( wH0 ); // get BUNBO
H0bunsi=numer( wH0 ); // get BUNSHI
H0poles=roots( H0bunbo ); // get poles
H0zeros=roots( H0bunsi ); // get zeros
[wcoeff_bunbo]=get_coeff( H0bunbo, wiiro2 );
[wcoeff_bunsi]=get_coeff( H0bunsi, wiiro2 );
plzr( wH0 ); // plot zeros and poles
//
global ac0
global bcb0
global bc0
//
disp('digital discrete time system');
disp( H0nx );
disp( '*bunbo/poles/x*');
disp( H0bunbo);
disp( wcoeff_bunbo(1), 'a[0] is always 1' );
for v=1:wiiro2
ac0(v)= - wcoeff_bunbo(v+1);
wstr1= 'a[' + string(v) + ']';
disp( ac0(v), wstr1 );
end
disp('poles are');
disp( H0poles);
disp( '*bunsi/zeros/o*')
disp( H0bunsi);
bcb0=wcoeff_bunsi(1);
disp( wcoeff_bunsi(1), 'b[0]' );
for v=1:wiiro2
bc0(v)= wcoeff_bunsi(v+1);
wstr1= 'b[' + string(v) + ']';
disp( bc0(v), wstr1 );
end
disp('zeros are');
disp( H0zeros);
//
//
xset('window',wb0); // push old windows
endfunction
//---------------------------------------------------------------------------------------------------------------
function [wiiro2]=system_stable_check2( disp0, wH0)
//
wb0=xget('window'); // stack old window
xset('window',disp0); // create new windows
clf();
[H0r,H0nx]=syssize( wH0 ); // get size
wiiro2=H0nx;
sz=size(cells1);
disp('poles x are');
disp( zpoles );
disp( 'zeros o are ')
disp( zzeros);
plzr( wH0 ); // plot zeros and poles
//
global ac2
global bcb2
global bc2
ac2=zeros(mc2,2); // coefficient matrix of a[*][1,2]
bcb2=zeros(mc2); // coefficient matrix of b[*][0]
bc2=zeros(mc2,2); // coefficient matrix of b[*][1,2]
for v=1:sz(2)
wH1=cells1(1,v);
H0bunbo=denom( wH1 ); // get BUNBO
H0bunsi=numer( wH1 ); // get BUNSHI
[H0r2,H0nx2]=syssize( wH1 );
[wcoeff_bunbo]=get_coeff( H0bunbo, H0nx2 );
[wcoeff_bunsi]=get_coeff( H0bunsi, H0nx2 );
for w=1:H0nx2
ac2(v,w)= - wcoeff_bunbo(w+1);
wstr1= 'a[' + string(v) +','+ string(w) + ']';
disp( ac2(v,w), wstr1 );
end
bcb2(v)=wcoeff_bunsi(1);
disp( bcb2(v), 'bcb2[0]' );
for w=1:H0nx2
bc2(v,w)= wcoeff_bunsi(w+1);
wstr1= 'b[' + string(v) +',' + string(v) + ']';
disp( bc2(v,w), wstr1 );
end
end // for v=1:sz(2)
//
//
//[H0r,H0nx]=syssize( wH0 ); // get size
//wiiro2=H0nx;
H0bunbo=denom( wH0 ); // get BUNBO
H0bunsi=numer( wH0 ); // get BUNSHI
H0poles=roots( H0bunbo ); // get poles
H0zeros=roots( H0bunsi ); // get zeros
[wcoeff_bunbo]=get_coeff( H0bunbo, wiiro2 );
[wcoeff_bunsi]=get_coeff( H0bunsi, wiiro2 );
//plzr( wH0 ); // plot zeros and poles
//
global ac0
global bcb0
global bc0
//
//disp('digital discrete time system');
//disp( H0nx );
//disp( '*bunbo/poles/x*');
//disp( H0bunbo);
//disp( wcoeff_bunbo(1), 'a[0] is always 1' );
for v=1:wiiro2
ac0(v)= - wcoeff_bunbo(v+1);
wstr1= 'a[' + string(v) + ']';
//disp( ac0(v), wstr1 );
end
//disp('poles are');
//disp( H0poles);
//disp( '*bunsi/zeros/o*')
//disp( H0bunsi);
bcb0=wcoeff_bunsi(1);
//disp( wcoeff_bunsi(1), 'b[0]' );
for v=1:wiiro2
bc0(v)= wcoeff_bunsi(v+1);
wstr1= 'b[' + string(v) + ']';
//disp( bc0(v), wstr1 );
end
//disp('zeros are');
//disp( H0zeros);
xset('window',wb0); // push old windows
endfunction
//-------------------------------------------------------------------------------------------
function [wfreqx, wcycle, wskip_cycle]= set_synth()
//
txt1=['frequency [Hz]'; 'cycle (integral number)'; 'start
skip cycle (integral number)'];
wstr1=sprintf('%f',freqx);
wstr2=sprintf('%d',cyclex);
wstr3=sprintf('%d',skip_cyclex);
sig1=x_mdialog('do actual IIR filter calculation',txt1,
[wstr1 ; wstr2 ; wstr3 ]);
if sig1==[] then
arg1=evstr(wstr1);
arg2=evstr(wstr2);
arg3=evstr(wstr3);
else
arg1=evstr(sig1(1));
arg2=evstr(sig1(2));
arg3=evstr(sig1(3));
end
wfreqx=arg1;
wcycle=arg2;
wskip_cycle=arg3;
//
endfunction
//---------------------------------------------------------------------------------------------------------------
function [win0,wsin_delta,xm,xsm]= make_sin()
//
wsin_delta = freqx / fs * 2. * %pi;
xm= int( ((1.0/ freqx) / (1.0 /fs ) * cyclex) + 1.0);
xsm= int( ((1.0/ freqx) / (1.0 /fs ) * skip_cyclex) + 1.0);
win0=zeros(1,xm);
for v=1:xm
win0(v)= sin_level * sin( wsin_delta * (v-1));
end
//
endfunction
//---------------------------------------------------------------------------------------------------------------
function [win0,wout0,xm]=start_skip()
xm=m-sm;
win0=zeros(1,xm);
wout0=zeros(1,xm);
for v=(sm+1):m
win0( v - sm) = in0(v);
wout0( v - sm) = out0(v);
end
//
endfunction
//---------------------------------------------------------------------------------------------------------------
function [wout0]= cal_iir2(init_flag)
//
global wk0
if init_flag == 1 then
for v=1:(mc0*2)
wk0(v)=0.0;
end
end
//
//
wout0=zeros(1, m);
//
for i=1:m
w = in0(i);
for v=1:iiro2
w= w + ac0(v) * wk0(v);
wout0(i)= wout0(i) + bc0(v) * wk0(v);
end
wout0(i) = wout0(i) + bcb0 * w;
for v=1:(iiro2-1)
wk0( iiro2 - v + 1) = wk0( iiro2 - v);
end
wk0(1)=w;
end
//
endfunction
//---------------------------------------------------------------------------------------------------------------
function [wout0]= cal_eqiir2(init_flag)
//
global wk2
if init_flag == 1 then
for v=1:mc2
wk2(v,1)=0.0;
wk2(v,2)=0.0;
wk2(v,3)=0.0;
wk2(v,4)=0.0;
end
end
//
//
wout0=zeros(1, m);
sz=size(cells1);
//
for i=1:m
// w= fact * in0(i);
w= in0(i);
for v=1:sz(2)
wo= bcb2(v) * w + bc2(v,1) * wk2(v,1) +
bc2(v,2) * wk2(v,2) + ac2(v,1) * wk2(v,3) + ac2(v,2) * wk2(v,4);
wk2(v,2)=wk2(v,1);
wk2(v,1)=w;
wk2(v,4)=wk2(v,3);
wk2(v,3)=wo;
w=wo;
end
wout0(i)=wo;
end
//
endfunction
//---------------------------------------------------------------------------------------------------------------
function plotfit(disp0)
wb0=xget('window'); // stack old window
xset('window',disp0); // create new windows
clf();
//
g0= 20.0 * log10( xopt1(1) / xopt0(1));
ph0= 360 * (xopt1(2)-xopt0(2))/2./%pi
//
freqa=[freqx,freqx];
freqd=[(freqx/fs),(freqx/fs)];
[freqs,respa]=repfreq(S1A, freqa);
[db_fft1,phi_fft1] =dbphi(respa);
[frq,respf]=repfreq(H0, freqd );
[db_fft2,phi_fft2] =dbphi(respf);
//
subplot(211);
tm=[0:1:(m-1)];
tm=tm';
plot(tm,in0','+',tm,yth(tm,xopt0),'-');
wstr1 = 'input of iir of ' + string(freqx) + 'Hz /' + '
theoretical analog gain' + string(db_fft1(1)) + ' dB' ;
xtitle( wstr1 ,'','');
//
subplot(212);
plot(tm,out0','+',tm,yth(tm,xopt1),'-');
wstr1 = 'output of iir of ' + string(g0) + ' dB ' + string( ph0)
+ ' vs ' + ' theoretical digital gain' + string(db_fft2(1)) + '
dB' ;
xtitle( wstr1 ,'','');
//
xset('window',wb0);
endfunction
//---------------------------------------------------------------------------------------------------------------
function plotfit2(disp0)
wb0=xget('window'); // stack old window
xset('window',disp0); // create new windows
clf();
//
g0= 20.0 * log10( xopt1(1) / xopt0(1));
ph0= 360 * (xopt1(2)-xopt0(2))/2./%pi
//
freqa=[freqx,freqx];
freqd=[(freqx/fs),(freqx/fs)];
//[freqs,respa]=repfreq(S1A, freqa);
//[db_fft1,phi_fft1] =dbphi(respa);
[frq,respf]=repfreq(H0, freqd );
[db_fft2,phi_fft2] =dbphi(respf);
//
subplot(211);
tm=[0:1:(m-1)];
tm=tm';
plot(tm,in0','+',tm,yth(tm,xopt0),'-');
wstr1 = 'input of iir of ' + string(freqx) + 'Hz' ;
xtitle( wstr1 ,'','');
//
subplot(212);
plot(tm,out0','+',tm,yth(tm,xopt1),'-');
wstr1 = 'output of iir of ' + string(g0) + ' dB ' + string( ph0)
+ ' vs ' + ' theoretical digital gain' + string(db_fft2(1)) + '
dB' ;
xtitle( wstr1 ,'','');
//
xset('window',wb0);
endfunction
//---------------------------------------------------------------------------------------------------------------
function y = yth(t, x)
y = x(1) * sin( sin_delta * t + x(2))
//y = x(1) * sin( t + x(2))
endfunction
//--
// myfun(x0, tm, in0, wm)
// wm .* ( yth(tm, x0) - in0)
function e = myfun(x, tm, ym, wm)
e = wm .* ( yth(tm, x) - ym )
endfunction
//--
function [xopt]=sin_fit(wym)
// ym is samples
tm=[0:1:(m-1)];
tm=tm';
wm=ones(m,1);
// estimated function
// initial parameter
wz=abs(wym(1));
for v=1:m
if abs(wym(v)) > wz then
wz=abs(wym(v));
end
end
x0 = [ wz ; 0.0];
// 1- the simplest call
[f,xopt, gropt] = leastsq(list(myfun,tm,wym',wm),x0);
disp( xopt,'xopt');
endfunction
//---------------------------------------------------------------------------------------------------------------
function save_para()
//
savefilename= input(' + enter file name for save parameters (ex:
para1.txt)',["string"]);
[fd0,err0]=mopen(savefilename,'w');
wstr1='/* ' + string(iiro2) + ' ORDER IIR DIGITAL
DE-EMPHASIS FILTER
COEFFICIENT */'
mfprintf(fd0,'%s\n',wstr1);
wstr1='/* CONSTANT1 = ' + string(c50) +
'
*/'
mfprintf(fd0,'%s\n',wstr1);
wstr1='/* CONSTANT2 = ' + string(c15) +
'
*/'
mfprintf(fd0,'%s\n',wstr1);
wstr1='/* SAMPLING FREQUENCY = ' + string(fs) +
'
*/'
mfprintf(fd0,'%s\n',wstr1);
mfprintf(fd0,'\n');
for v=1:iiro2
wstr1 = 'a[' + string(v) + ']=' + string( ac0(v)) ;
//mputstr(wstr1,fd0);
mfprintf(fd0,' %s;\n',wstr1);
end
wstr1 = 'b[0]=' + string( bcb0) ;
mfprintf(fd0,' %s;\n',wstr1);
for v=1:iiro2
wstr1 = 'b[' + string(v) + ']=' + string( bc0(v)) ;
//mputstr(wstr1,fd0);
mfprintf(fd0,' %s;\n',wstr1);
end
//mfprintf(fd0,'\n');
mfprintf(fd0,'\n');
mfprintf(fd0,'\n');
mfprintf(fd0,'\n');
wstr1='/* An example of TYOKUSETU-GATA
2
*/';
mfprintf(fd0,'%s\n',wstr1);
wstr1='/* in[i] is input, out[i] is
output
*/';
mfprintf(fd0,'%s\n',wstr1);
wstr1='/* w[] are memory to keep stack old
data
*/';
mfprintf(fd0,'%s\n',wstr1);
wstr1='/* initialize all w[1]=.....=0.0, and
then use it */';
mfprintf(fd0,'%s\n',wstr1);
for v=1:iiro2
wstr1= 'w[' + string(v) + '] = 0.0';
mfprintf(fd0,' %s;\n',wstr1);
end
mfprintf(fd0,'\n');
mfprintf(fd0,'\n');
wstr1='for (i=0; i< number ; ++i) ';
mfprintf(fd0,' %s\n',wstr1);
wstr1='{';
mfprintf(fd0,' %s\n',wstr1);
wstr1='w[0] = in[i] + ';
for v=1:iiro2
wstr1= wstr1 + 'a[' + string(v) + '] * w[' + string(v) + ']';
if v < iiro2
wstr1= wstr1 + ' + ';
end
end
mfprintf(fd0,' %s;\n',wstr1);
wstr1='out[i] = ';
for v=0:iiro2
wstr1= wstr1 + 'b[' + string(v) + '] * w[' + string(v) + ']';
if v < iiro2
wstr1= wstr1 + ' + ';
end
end
mfprintf(fd0,' %s;\n',wstr1);
mfprintf(fd0,'\n');
for v=1:iiro2
wstr1= 'w[' + string(iiro2-v+1) +'] = w[' +
string(iiro2-v) + ']';
mfprintf(fd0,' %s;\n',wstr1);
end
wstr1='}';
mfprintf(fd0,' %s\n',wstr1);
//
err0=mclose([fd0]);
endfunction
//---------------------------------------------------------------------------------------------------------------
function save_para2()
//
savefilename= input(' + enter file name for save parameters (ex:
para1.txt)',["string"]);
[fd0,err0]=mopen(savefilename,'w');
wstr1='/* ' + string(iiro2) + ' ORDER IIR DIGITAL
DE-EMPHASIS FILTER
COEFFICIENT */'
mfprintf(fd0,'%s\n',wstr1);
wstr1='/* SAMPLING FREQUENCY = ' + string(fs) +
'
*/'
mfprintf(fd0,'%s\n',wstr1);
mfprintf(fd0,'\n');
//
sz=size(cells1);
//wstr1 = 'fact=' + string( fact ) ;
//mfprintf(fd0,' %s;\n',wstr1);
//wstr1 = '/* fact1=' + string( fact1 ) + ' */' ;
//mfprintf(fd0,' %s;\n',wstr1);
for v=1:sz(2)
for w=1:2
wstr1 = 'a[' + string(v-1) + ',' + string(w) + ']=' +
string( ac2(v,w)) ;
//mputstr(wstr1,fd0);
mfprintf(fd0,' %s;\n',wstr1);
end
//end
//for v=1:sz(2)
w=0;
wstr1 = 'b[' + string(v-1) + ',' + string(w) + ']=' + string(
bcb2(v)) ;
mfprintf(fd0,' %s;\n',wstr1);
for w=1:2
wstr1 = 'b[' + string(v-1) + ',' + string(w) + ']=' + string(
bc2(v,w)) ;
//mputstr(wstr1,fd0);
mfprintf(fd0,' %s;\n',wstr1);
end
mfprintf(fd0,'\n');
end
//mfprintf(fd0,'\n');
mfprintf(fd0,'\n');
mfprintf(fd0,'\n');
mfprintf(fd0,'\n');
wstr1='/* An example of TYOKUSETU-GATA
1
*/';
mfprintf(fd0,'%s\n',wstr1);
wstr1='/* in[i] is input, out[i] is
output
*/';
mfprintf(fd0,'%s\n',wstr1);
wstr1='/* w[] are memory to keep stack old
data
*/';
mfprintf(fd0,'%s\n',wstr1);
wstr1='/* initialize all w[1]=.....=0.0, and
then use it */';
mfprintf(fd0,'%s\n',wstr1);
for v=1:sz(2)
for w=1:4
wstr1= 'w[' + string(v-1) + ',' + string(w-1) + '] = 0.0';
mfprintf(fd0,' %s;\n',wstr1);
end
end
mfprintf(fd0,'\n');
mfprintf(fd0,'\n');
wstr1='for (i=0; i< number ; ++i) ';
mfprintf(fd0,' %s\n',wstr1);
wstr1='{';
mfprintf(fd0,' %s\n',wstr1);
// wstr1=' wi=fact * in[i] ;';
wstr1=' wi=in[i] ;';
mfprintf(fd0,' %s\n',wstr1);
wstr1=' for(j=0; j< ' + string(sz(2)) + '; ++j) ';
mfprintf(fd0,' %s\n',wstr1);
wstr1=' {';
mfprintf(fd0,' %s\n',wstr1);
wstr1=' wo= b[j,0] * wi + b[j,1] * w[j,0] +
b[j,2] * w[j,1] + a[j,1] * w[j,2] + a[j,2] * w[j,3];';
mfprintf(fd0,' %s\n',wstr1);
wstr1=' w[j,1]=w[j,0];';
mfprintf(fd0,' %s\n',wstr1);
wstr1=' w[j,0]=wi;';
mfprintf(fd0,' %s\n',wstr1);
wstr1=' w[j,3]=w[j,2];';
mfprintf(fd0,' %s\n',wstr1);
wstr1=' w[j,2]=wo;';
mfprintf(fd0,' %s\n',wstr1);
wstr1=' wi=wo;';
mfprintf(fd0,' %s\n',wstr1);
wstr1=' }';
mfprintf(fd0,' %s\n',wstr1);
wstr1=' out[i]=wo;';
mfprintf(fd0,' %s\n',wstr1);
wstr1='}';
mfprintf(fd0,' %s\n',wstr1);
//
mfprintf(fd0,'\n');
mfprintf(fd0,'\n');
wstr1='/************************************************************/';
mfprintf(fd0,'%s\n',wstr1);
mfprintf(fd0,'\n');
for v=1:iiro2
wstr1 = 'a[' + string(v) + ']=' + string( ac0(v)) ;
//mputstr(wstr1,fd0);
mfprintf(fd0,' %s;\n',wstr1);
end
wstr1 = 'b[0]=' + string( bcb0) ;
mfprintf(fd0,' %s;\n',wstr1);
for v=1:iiro2
wstr1 = 'b[' + string(v) + ']=' + string( bc0(v)) ;
//mputstr(wstr1,fd0);
mfprintf(fd0,' %s;\n',wstr1);
end
err0=mclose([fd0]);
endfunction
//-------------------------------------------------------------------------------------------
impulse_steps=400;
//-------------------------------------------------------------------------------------------
function [ximpulse_steps]= set_synth2()
//
txt1=['steps (integral number)'];
wstr1=sprintf('%d',impulse_steps);
sig1=x_mdialog('impulse response based on do actual IIR
filter calculation',txt1, [wstr1 ]);
if sig1==[] then
arg1=evstr(wstr1);
else
arg1=evstr(sig1(1));
end
ximpulse_steps=arg1;
//
endfunction
//---------------------------------------------------------------------------------------------------------------
function [win0,wm]= make_impulse(ximpulse_steps)
//
win0=zeros(1,ximpulse_steps);
win0(1)=sin_level;
wm=ximpulse_steps;
//
endfunction
//---------------------------------------------------------------------------------------------------------------
function plotimpulse(disp0)
wb0=xget('window'); // stack old window
xset('window',disp0); // create new windows
clf();
//
plot(out0);
wstr1 = 'impulse response ' ;
xtitle( wstr1 ,'','');
//
xset('window',wb0);
endfunction
//---
function plotimpulse2(disp0,wdat1x)
wb0=xget('window'); // stack old window
xset('window',disp0); // create new windows
clf();
//
sz=size(wdat1x);
if sz(1,1)== 2 then
subplot(311);
plot(out0);
wstr1 = 'impulse response ' ;
xtitle( wstr1 ,'','');
subplot(312);
win1=zeros(1,impulse_steps);
for v=1:impulse_steps
win1(1,v)=wdat1x(1,v);
end
plot(win1);
wstr1 = 'wdat ch1 ' ;
xtitle( wstr1 ,'','');
subplot(313);
win1=zeros(1,impulse_steps);
for v=1:impulse_steps
win1(1,v)=wdat1x(2,v);
end
plot(win1);
wstr1 = 'wdat ch2 ' ;
xtitle( wstr1 ,'','');
elseif sz(1,1) == 1 then
subplot(211);
plot(out0);
wstr1 = 'impulse response ' ;
xtitle( wstr1 ,'','');
subplot(212);
win1=zeros(1,impulse_steps);
for v=1:impulse_steps
win1(1,v)=wdat1x(1,v);
end
plot(win1);
wstr1 = 'wdat ' ;
xtitle( wstr1 ,'','');
end
//
xset('window',wb0);
endfunction
//---------------------------------------------------------
function every_freqresp( xcells ,str0)
// delete old graphic windows
for v=10:20
xdel( v);
end
//
sz=size(xcells);
n0=sz(2);
for v=1:n0
H1=xcells(1,v);
H1=syslin('d',H1);
str= str0 + '(' + string(v) + ')';
freq_response2((10+v-1), H1, str);
end
endfunction
//---------------------------------------------------------------------------------------------------------------
//
// .wav file process
//
//
// for .wav 1
chs1=0; // channels
qty1=0; // data quantity, samples
size1=0; // actual loaded data quantity
fs1=0; // sampling frequency
bit1=0; // bits
wavfile1=''; // file name, this will be overwritten.
wdat1=zeros(1,10);// data of the .wav, only one first channel,
this will be overwritten.
wdat2=zeros(1,10);// data of synthesized .wav, this will be
overwritten.
// for others
f412=0; // flag to detect scilab 4.1.2
defch1=1; // default channel 1
f_win=1; // flag if windows
//
// max stack size of input wave file
maxduration=30; // unit is second
//-----------------------------------------------------------------
if isdir('c:\\') then
f_win=1;
else
f_win=0;
end
//-------------------------------------------------------------------
function [xwdat1,wavfile1,chs1,qty1,fs1,bit1,f412]=get_wavfile_pro()
// choose wave file
wavfile1=tk_getfile(Title="Choose xxx.wav file name");
// read channels and samples
cs1=wavread(wavfile1,'size');
chs1=cs1(2); // channel
qty1=cs1(1); // samples
//...
if chs1 > 2 then // invert data for scilab-4.1.2
chs1=cs1(1);
qty1=cs1(2);
f412=1;
else
f412=0;
end
//...
sizeof0=8;
maxstacksize= 44100 * 2 * maxduration * sizeof0; // fs=44.1KHz 2
channels * maxduration seconds
if ( chs1 * qty1 * sizeof0 ) > maxstacksize then
qty1 = int((maxstacksize / (chs1 * sizeof0)));
disp('Waring: force to decrease data quantity(samples) ');
end
sz=stacksize();
if (sz(1)-sz(2)) < maxstacksize then
stacksize( (maxstacksize+sz(1)));
disp('increase stack memory ');
end
//
if f412 == 1 then // for windows scilab-4.1.2
[xwdat1,fs1,bit1]=wavread(wavfile1,[1,qty1]);
else
[xwdat1b,fs1,bit1]=wavread(wavfile1,[1,qty1]);
xwdat1=xwdat1b';
end
//
disp('read done');
endfunction
//----------------------------------------------------------------
function snd_play2(xwdat2)
// this sound doesnot work well on windows scilab-3.1.1 and linux
scilab-3.1
// but, this sound works on windows scilab-4.1.2. But, linux
scilab-4.1.2 doesnot!
//
sz=size(xwdat2);
size2=sz(2);
if f_win == 1 then
if f412 == 1 then // for windows scilab-4.1.2
sound(xwdat2 ,fs1,bit1);
else
// for windows scilab-3.1.1
disp('This function does not work.');
end
else // if f_win == 1 then
if f412 == 1 then // for scilab-4.1.2
disp('If sound setting is good for scilab, this function maybe
work.');
sound(xwdat2 ,fs1,bit1);
else
disp('This function does not work.');
end
end //if f_win == 1 then
endfunction
//
//----------------------------------------------------------------
function snd_save2( xwdat2)
wavefilename= input(' + enter file name for saved .wav file
=>',["string"]);
//
if f412 == 1 then // for windows scilab-4.1.2
wavwrite(xwdat2 ,fs1,bit1,wavefilename );
else
// for windows scilab-3.1.1
wavwrite(xwdat2',fs1,bit1, wavefilename);
end
endfunction
//
//----------------------------------------------------------------
// function display the .wav property
function disp_pro_wav()
disp(qty1,'data quantity',bit1,
'bits',chs1,'channels',fs1,'sampling frequency');
endfunction
//--------------------------------------
function [xwdat2]=do_iir_filtering(xwdat1)
sz=size(xwdat1);
xwdat2=zeros(sz(1),sz(2));
disp( sz(1),'channels ', sz(2),'data quantity ');
//
wk0=zeros(sz(1),iiro2);
//
for i=1:sz(2) // long
a=int( i / 44100);
b= i - 44100 * a;
if b == 0 then
disp( i );
end
for p=1:sz(1) // ch
w = xwdat1(p,i);
for v=1:iiro2
w= w + ac0(v) * wk0(p,v);
xwdat2(p,i)= xwdat2(p,i) + bc0(v) * wk0(p,v);
end
xwdat2(p,i) = xwdat2(p,i) + bcb0 * w;
for v=1:(iiro2-1)
wk0( p, iiro2 - v + 1) = wk0( p,iiro2 - v);
end
wk0(p,1)=w;
if xwdat2(p,i) > 1.0 then
xwdat2(p,i) =1.0;
elseif xwdat2(p,i) < -1.0 then
xwdat2(p,i)=-1.0;
end
end
end
//
disp('do iir filter done');
endfunction
//--------------------------------------
function [xwdat2]=do_eqiir_filtering(xwdat1)
sz=size(xwdat1);
xwdat2=zeros(sz(1),sz(2));
disp( sz(1),'channels ', sz(2),'data quantity ');
//
wk2=zeros(mc2,4,sz(1));
//
szc=size(cells1);
//
for i=1:sz(2) // long
a=int( i / 44100);
b= i - 44100 * a;
if b == 0 then
disp( i );
end
for p=1:sz(1) // ch
w = xwdat1(p,i);
w= fact * w;
for v=1:szc(2)
wo= bcb2(v) * w + bc2(v,1) * wk2(v,1,p) +
bc2(v,2) * wk2(v,2,p) + ac2(v,1) * wk2(v,3,p) + ac2(v,2) * wk2(v,4,p);
wk2(v,2,p)=wk2(v,1,p);
wk2(v,1,p)=w;
wk2(v,4,p)=wk2(v,3,p);
wk2(v,3,p)=wo;
w=wo;
end
xwdat2(p,i)=wo;
if xwdat2(p,i) > 1.0 then
xwdat2(p,i) =1.0;
elseif xwdat2(p,i) < -1.0 then
xwdat2(p,i)=-1.0;
end
end //for p=1:sz(1) // ch
end // for i=1:sz(2) // long
//
disp('do iir filter done');
endfunction
//-------------------------------------------------------------------------------------
function [xwdat2]=do_iir_filtering2(wH0, xwdat1)
sz=size(xwdat1);
disp( sz(1),'channels ', sz(2),'data quantity ');
//
//
if sz(1) == 2 then
sl=tf2ss([ wH0, 0; 0, wH0]);
disp('from transfer to state-space');
xwdat2=dsimul(sl,xwdat1);
elseif sz(1) == 1 then
sl=tf2ss(wH0);
disp('from transfer to state-space');
xwdat2=dsimul(sl,xwdat1);
else
disp('error: not define yet. ');
end
//
disp('do iir filter done');
endfunction
//
//
//-------------------------------------------------------------------------------------
function [xwdat2]=do_iir_filtering3(wH0, xwdat1)
sz=size(xwdat1);
disp( sz(1),'channels ', sz(2),'data quantity ');
//
//
if sz(1) == 2 then
[H0r,H0nx]=syssize( wH0 ); // get size
viiro2=H0nx;
H0bunbo=denom( wH0 ); // get BUNBO
H0bunsi=numer( wH0 ); // get BUNSHI
vnum=[ H0bunsi,0; 0,H0bunsi];
vdenom=[ H0bunbo,0; 0,H0bunbo];
up=zeros(2,viiro2);
yp=zeros(2,viiro2);
xwdat2=rtitr(vnum,vdenom,xwdat1 ,up,yp);
elseif sz(1) == 1 then
[H0r,H0nx]=syssize( wH0 ); // get size
viiro2=H0nx;
H0bunbo=denom( wH0 ); // get BUNBO
H0bunsi=numer( wH0 ); // get BUNSHI
vnum= H0bunsi;
vdenom=H0bunbo;
up=zeros(1,viiro2);
yp=zeros(1,viiro2);
xwdat2=rtitr(vnum,vdenom,xwdat1 ,up,yp);
else
disp('error: not define yet. ');
end
//
disp('do iir filter done');
endfunction
//
//
//----------------------------------------------------------------
//===================================================================
//
// ++ MENU MENU MENU =====================
//
// Add menu buttons of which name are 'a_plotwav' and ''
//
if DE_EMPHASIS == 1 then
ADD_MENU_PLOT=1; // set ADD_MENU 1
to add menu button of some functions
else
ADD_MENU_PLOT=0;
end
//
//
if ADD_MENU_PLOT == 1 then
delmenu('step1_de_emphasis');
addmenu('step1_de_emphasis',[ '(1)de-emphasis set_constant ';
'(2)design IIR filter' ]);
step1_de_emphasis=[ '[c50, c15, fs, iiro, c15_change]=set_constant()'
;' [H0, S1A, iiro2]=filterdesign(); freq_response(0, 3, H0,
S1A);'] ;
end //if ADD_MENU_PLOT == 1 then
//
//
//
if EQIIRR==1 then
ADD_MENU_EQIIR=1; // set ADD_MENU 1 to add menu button of
some functions
else
ADD_MENU_EQIIR=0;
end
//
if ADD_MENU_EQIIR == 1 then
delmenu('step1_eqiir');
addmenu('step1_eqiir',[ '(1)select LPF,HPF,BPF,SBF'; '(2)select
Butterworth filter or Chebyshev type 1 filter'; '(3)set constant';
'(4)design IIR filter' ]);
step1_eqiir=[ '[ftype]=selectftype();' ; '[approx]=selectapprox();';
'[fs,omf,deltap,deltas]=set_constant2(); ' ;'[H0,
cells,cells1,fact,fact1,zzeros,zpoles,iiro2]=call_eqiir();
freq_response2(0, H0,''overall''); every_freqresp(cells1,''cells1'') ;'
] ;
end //if ADD_MENU_EQIIR == 1 then
//
//
//
if ANA2IIR==1 then
ADD_MENU_ANA2IIR=1; // set ADD_MENU 1 to add menu button of
some functions
else
ADD_MENU_ANA2IIR=0;
end
//
if ADD_MENU_ANA2IIR == 1 then
delmenu('step1_ana2iir');
addmenu('step1_ana2iir',[ '(1)select LPF,HPF,BPF,SBF'; '(2)select
Butterworth filter or Chebyshev type 1 filter'; '(3)set constant';
'(4)design IIR filter' ]);
step1_ana2iir=[ '[ftype]=selectftype();' ; '[approx]=selectapprox();';
'[zisuu,fs,omf,deltap]= set_constant3(); ' ;'[H0, S1A,
iiro2]=filterdesign3(); freq_response(0, 3, H0, S1A);
[cells,cells1,fact,fact1,zpoles,zzeros]=get_cell(H0);
every_freqresp(cells1,''cells1'');' ] ;
end //if ADD_MENU_ANA2IIR == 1 then
//
//
if DE_EMPHASIS == 1 then
ADD_MENU_KAISEI=1; // set ADD_MENU 1 to add menu button of
some functions
else
ADD_MENU_KAISEI=0;
end
//
if ADD_MENU_KAISEI == 1 then
delmenu('step2_de_emphasis');
addmenu('step2_de_emphasis',[ '(1)system stable check '; '(2)gain by
actual IIR calculation' ]);
step2_de_emphasis=[ '[iiro2]=system_stable_check(1, H0,S1A);' ;
'[freqx, cyclex, skip_cyclex]= set_synth(); [in0,sin_delta,m,sm]=
make_sin(); [out0]= cal_iir2(1); [in0,out0,m]=start_skip();
[xopt0]=sin_fit(in0); [xopt1]=sin_fit(out0); plotfit(3);'] ;
end //if ADD_MENU_KAISEI == 1 then
//
//
if EQIIRR==1 then
ADD_MENU_KAISEI2=1; // set ADD_MENU 1 to add menu button of
some functions
else
ADD_MENU_KAISEI2=0;
end
//
//
if ADD_MENU_KAISEI2 == 1 then
delmenu('step2_eqiir');
addmenu('step2_eqiir',[ '(1)system stable check '; '(2)gain by actual
IIR calculation' ]);
step2_eqiir=[ '[iiro2]=system_stable_check2(1, H0);' ; '[freqx, cyclex,
skip_cyclex]= set_synth(); [in0,sin_delta,m,sm]= make_sin(); [out0]=
cal_eqiir2(1); [in0,out0,m]=start_skip(); [xopt0]=sin_fit(in0);
[xopt1]=sin_fit(out0); plotfit2(3);'] ;
end //if ADD_MENU_KAISEI == 1 then
//
//
if ANA2IIR==1 then
ADD_MENU_KAISEI3=1; // set ADD_MENU 1 to add menu button of
some functions
else
ADD_MENU_KAISEI3=0;
end
//
if ADD_MENU_KAISEI3 == 1 then
delmenu('step2_ana2iir');
addmenu('step2_ana2iir',[ '(1)system stable check '; '(2)gain by actual
IIR calculation' ]);
step2_ana2iir=[ '[iiro2]=system_stable_check(1, H0,S1A);' ; '[freqx,
cyclex, skip_cyclex]= set_synth(); [in0,sin_delta,m,sm]= make_sin();
[out0]= cal_iir2(1); [in0,out0,m]=start_skip(); [xopt0]=sin_fit(in0);
[xopt1]=sin_fit(out0); plotfit(3);'] ;
end //if ADD_MENU_KAISEI3 == 1 then
//
//
if DE_EMPHASIS == 1 then
ADD_MENU_SAVE=1; // set ADD_MENU 1 to add menu button of
some functions
else
ADD_MENU_SAVE=0;
end
//
if ADD_MENU_SAVE == 1 then
delmenu('step3_de_emphasis');
addmenu('step3_de_emphasis',[ '(1)save coefficient in a text file.' ]);
step3_de_emphasis=[ 'save_para();'] ;
end //if ADD_MENU_SAVE == 1 then
//
if EQIIRR==1 then
ADD_MENU_SAVE2=1; // set ADD_MENU 1 to add menu button of
some functions
else
ADD_MENU_SAVE2=0;
end
//
if ADD_MENU_SAVE2 == 1 then
delmenu('step3_eqiir');
addmenu('step3_eqiir',[ '(1)save coefficient in a text file.' ]);
step3_eqiir=[ 'save_para2();'] ;
end //if ADD_MENU_SAVE2 == 1 then
//
if ANA2IIR==1 then
ADD_MENU_SAVE3=1; // set ADD_MENU 1 to add menu button of
some functions
else
ADD_MENU_SAVE3=0;
end
//
if ADD_MENU_SAVE3 == 1 then
delmenu('step3_ana2iir');
addmenu('step3_ana2iir',[ '(1)save coefficient in a text
file.','(2)save coefficient in a text file.' ]);
step3_ana2iir=[ 'save_para();';'save_para2();'] ;
end //if ADD_MENU_SAVE3 == 1 then
//
if DE_EMPHASIS == 1 then
ADD_MENU_IMPULSE=1; // set ADD_MENU 1 to add menu button of
some functions
else
ADD_MENU_IMPULSE=0;
end
//
if ADD_MENU_IMPULSE == 1 then
delmenu('step4_de_emphasis');
addmenu('step4_de_emphasis',[ '(1a)impulse response ', '(1b)impulse
response (state-space method)']);
step4_de_emphasis=[ '[impulse_steps]= set_synth2(); [in0,m]=
make_impulse(impulse_steps); [out0]= cal_iir2(1); plotimpulse(4)';
'[impulse_steps]= set_synth2(); [in0,m]= make_impulse(impulse_steps);
[out0]= do_iir_filtering2(H0, in0); plotimpulse(4)' ] ;
end //if ADD_MENU_IMPULSE == 1 then
////
//
if EQIIRR==1 then
ADD_MENU_IMPULSE2=1; // set ADD_MENU 1 to add menu button
of some functions
else
ADD_MENU_IMPULSE2=0;
end
//
if ADD_MENU_IMPULSE2 == 1 then
delmenu('step4_eqiir');
addmenu('step4_eqiir',[ '(1a)impulse response ', '(1b)impulse response
(state-space method)']);
step4_eqiir=[ '[impulse_steps]= set_synth2(); [in0,m]=
make_impulse(impulse_steps); [out0]= cal_eqiir2(1); plotimpulse(4)';
'[impulse_steps]= set_synth2(); [in0,m]= make_impulse(impulse_steps);
[out0]= do_iir_filtering2(H0, in0); plotimpulse(4)' ] ;
end //if ADD_MENU_IMPULSE2 == 1 then
////
//
if ANA2IIR==1 then
ADD_MENU_IMPULSE3=1; // set ADD_MENU 1 to add menu button
of some functions
else
ADD_MENU_IMPULSE3=0;
end
//
if ADD_MENU_IMPULSE3 == 1 then
delmenu('step4_ana2iir');
addmenu('step4_ana2iir',[ '(1a)impulse response ', '(1b)impulse
response (state-space method)' ,'(-)','(etc)impluse response with
wdat2','(etc)impluse response with wdat1']);
step4_ana2iir=[ '[impulse_steps]= set_synth2(); [in0,m]=
make_impulse(impulse_steps); [out0]= cal_iir2(1); plotimpulse(4)';
'[impulse_steps]= set_synth2(); [in0,m]= make_impulse(impulse_steps);
[out0]= do_iir_filtering2(H0, in0);
plotimpulse(4)';'';'[impulse_steps]= set_synth2();[in0,m]=
make_impulse(impulse_steps); [out0]= do_iir_filtering2(H0, in0);
plotimpulse2(4,wdat2)' ;'[impulse_steps]= set_synth2();[in0,m]=
make_impulse(impulse_steps); [out0]= do_iir_filtering2(H0, in0);
plotimpulse2(4,wdat1)' ] ;
end //if ADD_MENU_IMPULSE3 == 1 then
//
//-----------------------------------------------------------------
//
if DE_EMPHASIS == 1 then
ADD_MENU_WAV=1; // set ADD_MENU 1
to add menu button of some functions
else
ADD_MENU_WAV=0;
end
//
//
if ADD_MENU_WAV == 1 then
delmenu('step5_de_emphasis');
addmenu('step5_de_emphasis',[ '(1)read a .wav file (read only limited
duration)', '(2a)do IIR filtering (takes time)' , '(2b)do IIR filtering
(state-space method)' ,'(3)play the wav filtered' , '(4)save as a .wav
file' ]);
step5_de_emphasis=[ 'wdat1=zeros(1,1); wdat2=zeros(1,1);
[wdat1,wavfile1,chs1,qty1,fs1,bit1,f412]=get_wavfile_pro();
disp_pro_wav(); ','[wdat2]=do_iir_filtering(wdat1);'
,'[wdat2]=do_iir_filtering2(H0,
wdat1);','snd_play2(wdat2);','snd_save2( wdat2);'];
end //if ADD_MENU_WAV == 1 then
//
//
if EQIIRR==1 then
ADD_MENU_WAV2=1; // set ADD_MENU 1
to add menu button of some functions
else
ADD_MENU_WAV2=0;
end
//
//
if ADD_MENU_WAV2 == 1 then
delmenu('step5_eqiir');
addmenu('step5_eqiir',[ '(1)read a .wav file (read only limited
duration)', '(2a)do IIR filtering (takes time)' , '(2b)do IIR filtering
(state-space method)' ,'(2c)do IIR filtering (rtitr)' ,'(3)play the wav
filtered' , '(4)save as a .wav file' ]);
step5_eqiir=[ 'wdat1=zeros(1,1); wdat2=zeros(1,1);
[wdat1,wavfile1,chs1,qty1,fs1,bit1,f412]=get_wavfile_pro();
disp_pro_wav(); ','[wdat2]=do_eqiir_filtering(wdat1);'
,'[wdat2]=do_iir_filtering2(H0, wdat1);','[wdat2]=do_iir_filtering3(H0,
wdat1);','snd_play2(wdat2);','snd_save2( wdat2);'];
end //if ADD_MENU_WAV2 == 1 then
//
if ANA2IIR==1 then
ADD_MENU_WAV3=1; // set ADD_MENU 1
to add menu button of some functions
else
ADD_MENU_WAV3=0;
end
//
if ADD_MENU_WAV3 == 1 then
delmenu('step5_ana2iir');
addmenu('step5_ana2iir',[ '(1)read a .wav file (read only limited
duration)', '(2a)do IIR filtering (takes time)' , '(2b)do IIR filtering
(state-space method)' , '(2c)do IIR filtering (rtitr)' ,'(3)play the
wav filtered' , '(4)save as a .wav file' ]);
step5_ana2iir=[ 'wdat1=zeros(1,1); wdat2=zeros(1,1);
[wdat1,wavfile1,chs1,qty1,fs1,bit1,f412]=get_wavfile_pro();
disp_pro_wav(); ','[wdat2]=do_iir_filtering(wdat1);'
,'[wdat2]=do_iir_filtering2(H0, wdat1);','[wdat2]=do_iir_filtering3(H0,
wdat1);','snd_play2(wdat2);','snd_save2( wdat2);'];
end //if ADD_MENU_WAV3 == 1 then
//
//
ADD_MENU_TEST1=1;
//
if ADD_MENU_TEST1 == 1 then
delmenu('step_test');
addmenu('step_test',[ '(1)test','(2)disp zpole,zeros','(3)frequency
response (xcells)' ]);
step_test=[
'[wcells,wcells1,wfact,wfact1,wzpoles,wzzeros]=get_cell(H0);
every_freqresp(wcells1,''wcells1'')' ; 'disp(zpoles,''zpoles'');
disp(wzpoles,''wzploes''); disp(zzeros,''zzeros''); disp(wzzeros,
''wzzeros'');' ; ' every_freqresp( xcells, ''xcells'' )'];
end //if ADD_MENU_TEST1 == 1 then
//
//
No.5 2008年8月23日