// ------------------------------------------------------------------------------------------ // plots space of peak level and space of two signal rate, of two tube model // on window scilab-4.1.2 // // This calculation has some of incompleteness. // ------------------------------------------------------------------------------------------ // PLEASE PAY ATTENSION that this program may have some bugs and also you may modify program // to fit your circumstances. If you use this program, everything done at your own risk. // Thanks. // ------------------------------------------------------------------------------------------ // setup // r1=[-0.9:0.1:0.9]; // set section and division of calculation of r1 // r1=(A2-A1)/(A2+A1) is reflection coefficient between 1st tube and 2nd tube l1=[-0.8:0.1:0.8]; // set section and division of calculation of l1 // l1=(L2-L1)/(L2+L1) , L1 is 1st tube's length. L2 is 2nd tube's length //r1=[0.4:0.4:0.8]; //l1=[-0.8:0.4:0.8]; // ttl_Length=17; // set total length of combined tubes, Two Tubes, and etc ttl_Area=10; // set total area of combined tubes, Two Tubes, and etc for dummy display // r1s=size(r1); size_r1=r1s(2); // size of array r1 l1s=size(l1); size_l1=l1s(2); // size of array l1 // L1=zeros(1,size_l1); L2=zeros(1,size_l1); for v=1:size_l1 L2(v)= ( 1. + l1(v) ) * ttl_Length / 2.; L1(v)= ttl_Length - L2(v); end A1=zeros(1,size_r1); A2=zeros(1,size_r1); for v=1:size_r1 A2(v)= ( 1. + r1(v) ) * ttl_Area/ 2.; A1(v)= ttl_Area - A2(v); end // ydbout=zeros(size_l1,size_r1); // kekka wo kakunou suru for 1st peak dB yfrqout=zeros(size_l1,size_r1); // kekka wo kakunou suru for 1st peak freq // c0=35000; // constant. the speed of sound is round 35000 cm/second // tu1=L1./c0; // delay time in 1st tube tu2=L2./c0; // delay time in 2nd tube // tclosed=5 ; // set glottis closed duration (off time) by unit is millisecond [ms] trise=6; // set glottis opening duration (rise time) by unit is millisecond [ms] tfall=2; // set glottis closing duration (fall time) by unit is millisecond [ms] // // reflection coefficient between glottis and 1st tube rg_rise=0.9; // set some value (1. or less) when glottis is opening rg_fall=0.95; // set some value (1. or less) when glottis is closing rg_closed=1; // When glottis is closed, // reflection coefficient between glottis and 1st tube is 1 // because glottis's gate is closed and vocal tract is free from glottis's influence // Convert from Volume Velocity to Sound Pressure using High Pass Filter Model for Mouth Radiation // Parameters ( cut off frequecny ) Setup // fc=1000; // set cut off frequency of High Pass Filter by unit is [Hz] fnew=fc; // reflection coefficient between 2nd tube and mouth rl=0.9; // set adjust reduction response. rl is variable, however use one constant in this. // // Other preparation: Sampling Frequency Setup // fs=44100; // set sampling frequency by unit is [Hz] ts=1/fs; // Max_Freq=7500; // set plot maximum frequency by unit is [Hz] //-------------------------------------------------------- // +make yg // length0=0; amplitude=1.0; N1= int( tclosed / 1000 / ts ); N2= int( trise / 1000 / ts); N3= int( tfall / 1000 / ts); LL= N1+N2+N3; length0 = length0 + LL; yg= zeros(1,length0+1); tc0=1; yg(tc0)=0.; yg(tc0)=amplitude * yg(tc0); for t0=1:int(LL) tc0=tc0+1; if t0 < N1 then yg(tc0) =0; elseif t0 <= (N2 + N1) then yg(tc0)= 0.5 * ( 1 - cos( ( %pi / N2 ) * (t0 - N1) ) ); elseif t0 <= (N3+N2+N1) then yg(tc0)= cos( ( %pi / ( 2 * N3 )) * ( t0 - (N2+N1) ) ); else yg(tc0) =0; end yg(tc0)=amplitude * yg(tc0); end // t0 // -make yg //-------------------------------------------------------------------- // // //----------------------------------------------------------- // +declaration common use // frq_com=[100:25:Max_Freq]; frq_com_s=size(frq_com); resp_com=zeros(1,frq_com_s(2)); resp_yg=ones(1,frq_com_s(2)); resp_hpf=ones(1,frq_com_s(2)); // // -declaration common use // function y = two_tubes2( xw, ni, nir, rg0, rl0 ) yi = 0.5 * ( 1 + rg0 ) * ( 1 + r1(nir)) * ( 1 + rl0 ) * %e^( -1 * ( tu1(ni) + tu2(ni) ) * xw * %i); yb = 1 + r1(nir) * rg0 * %e^( -2 * tu1(ni) * xw * %i) + r1(nir) * rl0 * %e^( -2 * tu2(ni) * xw * %i) + rl0 * rg0 * %e^( -2 * (tu1(ni) + tu2(ni)) * xw * %i); y = yi / yb; endfunction function two_tubes_model_bode2(ni, nir) resp3=zeros(1,frq_com_s(2)); for v=1:frq_com_s(2) resp3(v)= two_tubes2(2 * %pi * frq_com(v), ni, nir, rg_closed ,rl); end [db,phi] =dbphi(resp3); clf(); bode(frq_com,db,phi); endfunction // //------------------------------------------------- // +yg response function [res0]=yg_response(save_load) if save_load == 1 then // load old data load('yg_resp.dat', 'res0'); else stp=0; stp=stp+N1; y00=0.; for v=1:(N2+N3+1) y00= y00 + yg(stp + v); end for w=1:frq_com_s(2) xw=2 * %pi * frq_com(w) * ts; yi=0; for v=1:(N2+N3+1) yi= yi + yg(stp + v) * %e^(-1. * xw * (v-1) * %i); end // for v res0(w)=yi / y00; //resp_yg(w)=yi / y00; end // for w // save('yg_resp.dat', res0); end //-save //[db,phi] =dbphi(res0); //clf(); //bode(frq_com,db,phi); endfunction // -yg response //------------------------------------------------- // +hpf function [res1]=hpf_response() wc=2. * %pi * fc; al1= -1. * %e^( -1. * wc * ts); bl0= wc * ts; // Then, convert them to high pass filter's coefficients function y=f_alfa( wc_org, wc_new) y= -1. * cos( (wc_org + wc_new) * ts /2. ) / cos( (wc_org - wc_new) * ts / 2. ); endfunction wnew=2. * %pi * fnew; alfa=f_alfa( wc, wnew ); // high pass filter's coefficients bh0= bl0 / (1.0 - al1 * alfa); bh1= alfa * bl0 / (1.0 - al1 * alfa); ah1= (alfa - al1) / (1.0 - al1 * alfa); // Function for disply high pass filter's frequency-phase response function y = hpf1( xw ) yi= bh0 + bh1 * cos( xw * ts) - bh1 * sin( xw * ts ) * %i; yb=1.0 + ah1 * cos( xw * ts ) - ah1 * sin( xw * ts ) * %i; y=yi ./yb; endfunction res1=hpf1(2 * %pi * frq_com ); //resp_hpf=hpf1(2 * %pi * frq_com ); //[db,phi] =dbphi(resp_hpf); //clf(); //bode(frq_com,db,phi); endfunction // -hpf //---------------------------------------------------------- // function [db0,frq0]=two_tubes_model_res(ni, nir) for v=1:frq_com_s(2) dummy0= two_tubes2(2 * %pi * frq_com(v), ni, nir, rg_closed ,rl); resp_com(v)= resp_hpf(v) * resp_yg(v) * dummy0; end [db,phi] =dbphi(resp_com); // frq0=frq_com_s(2)-v+1; db0=db(frq0); for v=(frq_com_s(2)-1):-1:1 if db(v) > db0 then frq0=v; db0=db(frq0); end end // frq0=frq_com(frq0); // transfer scale Hz endfunction // //-------------------------------------------------------------------- // function [peak_db,peak_frq, peak_count0, avr_pwr, two_db_pwr, def_frq]=two_tubes_model_res2(ni, nir, xr_flag, show_bode) for v=1:frq_com_s(2) dummy0= two_tubes2(2 * %pi * frq_com(v), ni, nir, rg_closed ,rl); resp_com(v)= resp_hpf(v) * resp_yg(v) * dummy0; end [db,phi] =dbphi(resp_com); if show_bode==1 then clf(); bode(frq_com,db,phi); end // peak_count0=0; peak_db=ones(1,frq_com_s(2)); peak_frq=ones(1,frq_com_s(2)); for v=3:(frq_com_s(2)-2) //+peak detect using 4points around the point if db(v) > db(v-2) then if db(v) > db(v-1) then if db(v) > db(v+1) then if db(v) > db(v+2) then //+detect peak of which Q is more than q0, check only higher frequency side q0=2.; dfrq0=frq_com(v) / q0; frq_w1=frq_com(v)- ( dfrq0 / 2.0 ); // lower frequency side frq_w2=frq_com(v)+ ( dfrq0 / 2.0 ); // higher frequency side w2_v=-1.; for v2=(v+1):frq_com_s(2) if (db(v) - db(v2)) > 3.0 then // diffrence is more than 3dB w2_v=v2; break; end if frq_com(v2) > frq_w2 then w2_v=-3.; break; end if db(v2) > db(v) then // not peak in Q area w2_v=-2.; break; end end //-v2 if w2_v > 0 then // peak_count0=peak_count0+1; peak_db(peak_count0)=db(v); peak_frq(peak_count0)=frq_com(v); // end //if w2_v > 0 then end end end end //-peak detect end // // if xr_flag == 1 then mess7z=' peak_count0 '; messsp0=' '; disp(peak_count0,mess7z); for v=1:peak_count0 disp(peak_frq(v),messsp0, peak_db(v),messsp0, v); end end // max two peak points under condition that they are neighborhood (next to each other) pointn=1; if peak_count0 > 1 then v=pointn; if peak_db(v) < 0. then pk1= 1.0 / 10^(-1.0 * peak_db(v) / 20.0); else pk1= 10^( peak_db(v) / 20.0); end if peak_db(v+1) < 0. then pk2= 1.0 / 10^(-1.0 * peak_db(v+1) / 20.0); else pk2= 10^( peak_db(v+1) / 20.0); end pk_ttl=pk1+pk2; peak_db_ttl=20. * log10( pk_ttl); // if peak_count0 > 2 then for v=2:(peak_count0-1) if peak_db(v) < 0. then pk1= 1.0 / 10^(-1.0 * peak_db(v) / 20.0); else pk1= 10^( peak_db(v) / 20.0); end if peak_db(v+1) < 0. then pk2= 1.0 / 10^(-1.0 * peak_db(v+1) / 20.0); else pk2= 10^( peak_db(v+1) / 20.0); end pk_ttl2=pk1+pk2; peak_db_ttl2=20. * log10( pk_ttl2); if peak_db_ttl2 > peak_db_ttl then peak_db_ttl=peak_db_ttl2; pointn=v; end // end end //-if peak_count0 > 2 then end //-if peak_count0 > 1 then if xr_flag == 1 then mess7z2=' peak_ttl of two points '; disp(pointn,mess7z2); disp(peak_db_ttl); end //+ average power pk2=0.; for v=1:frq_com_s(2) if db(v) < 0. then pk1= 1.0 / 10^(-1.0 * db(v) / 20.0); else pk1= 10^( db(v) / 20.0); end pk2=pk2+pk1; end pk2=pk2 / frq_com_s(2); avr_pwr=10. * log10( pk2); if xr_flag == 1 then mess7z3=' average power '; disp(avr_pwr,mess7z3); end //------------------------------------------------------------- //+ two db power: strength scale is dB not linear pointn2=1; if peak_count0 > 1 then two_db_pwr=( peak_db(1)+peak_db(2) ) /2.; // if peak_count0 > 2 then for v=2:(peak_count0-1) two_db_pwr2= ( peak_db(v)+peak_db(v+1)) / 2.; if two_db_pwr2 > two_db_pwr then two_db_pwr = two_db_pwr2; pointn2=v; end // end end //-if peak_count0 > 2 then end //-if peak_count0 > 1 then // //---------------------------------------------------------------------- // // define def_frq= (w2-w0)/w0 // // def_frq= (peak_frq(pointn2+1)-peak_frq(pointn2))/peak_frq(pointn2); if xr_flag == 1 then mess7z4=' peak_db_pwr of two points, def_frq[rate] '; disp(pointn2,mess7z4); disp(def_frq , two_db_pwr); end //+ average power pk2=0.; for v=1:frq_com_s(2) pk2=pk2+db(v); end pk2=pk2 / frq_com_s(2); if xr_flag == 1 then mess7z5=' average db power '; disp(pk2,mess7z3); end //- two db power //----------------------------------------------------------------- endfunction // //------------------------------------------------------------- // function input_draw() l1in=input('+input l1 '); r1in=input('+input r1 '); stack_r1=r1(1); stack_l1=l1(1); r1(1)=r1in; l1(1)=l1in; // L2(1)= ( 1. + l1(1) ) * ttl_Length / 2.; L1(1)= ttl_Length - L2(1); tu1(1)=L1(1)/c0; tu2(1)=L2(1)/c0; A2(1)= ( 1. + r1(1) ) * ttl_Area/ 2.; A1(1)= ttl_Area - A2(1); disp( L2(1),L1(1),A2(1),A1(1)); // [peak_db0,peak_frq0,peak_count0,avr_pwr0,two_db_pwr0,def_frq0]=two_tubes_model_res2(1,1,1,1); // r1(1)=stack_r1;; l1(1)=stack_l1; L2(1)= ( 1. + l1(1) ) * ttl_Length / 2.; L1(1)= ttl_Length - L2(1); tu1(1)=L1(1)/c0; tu2(1)=L2(1)/c0; A2(1)= ( 1. + r1(1) ) * ttl_Area/ 2.; A1(1)= ttl_Area - A2(1); // endfunction //+++ function input_draw_check(A1x,A2x,L1x,L2x) //l1in=(8-9)/(8+9) //r1in=(7-1)/(7+1) l1in=(L2x-L1x)/(L2x+L1x) r1in=(A2x-A1x)/(A2x+A1x) stack_r1=r1(1); stack_l1=l1(1); r1(1)=r1in; l1(1)=l1in; // L2(1)= ( 1. + l1(1) ) * ttl_Length / 2.; L1(1)= ttl_Length - L2(1); tu1(1)=L1(1)/c0; tu2(1)=L2(1)/c0; A2(1)= ( 1. + r1(1) ) * ttl_Area/ 2.; A1(1)= ttl_Area - A2(1); disp( L2(1),L1(1),A2(1),A1(1)); // [peak_db0,peak_frq0,peak_count0,avr_pwr0,two_db_pwr0,def_frq0]=two_tubes_model_res2(1,1,1,1); // r1(1)=stack_r1;; l1(1)=stack_l1; L2(1)= ( 1. + l1(1) ) * ttl_Length / 2.; L1(1)= ttl_Length - L2(1); tu1(1)=L1(1)/c0; tu2(1)=L2(1)/c0; A2(1)= ( 1. + r1(1) ) * ttl_Area/ 2.; A1(1)= ttl_Area - A2(1); // endfunction // function input_draw_area() l1in=input('+input l1 '); r1in=input('+input r1 '); stack_r1=r1(1); stack_l1=l1(1); r1(1)=r1in; l1(1)=l1in; // L2(1)= ( 1. + l1(1) ) * ttl_Length / 2.; L1(1)= ttl_Length - L2(1); tu1(1)=L1(1)/c0; tu2(1)=L2(1)/c0; A2(1)= ( 1. + r1(1) ) * ttl_Area/ 2.; A1(1)= ttl_Area - A2(1); // plot_area(1); disp( L2(1),L1(1),A2(1),A1(1)); // r1(1)=stack_r1;; l1(1)=stack_l1; L2(1)= ( 1. + l1(1) ) * ttl_Length / 2.; L1(1)= ttl_Length - L2(1); tu1(1)=L1(1)/c0; tu2(1)=L2(1)/c0; A2(1)= ( 1. + r1(1) ) * ttl_Area/ 2.; A1(1)= ttl_Area - A2(1); // endfunction // // // function input_draw_check_a() //A1,A2,L1,L2 input_draw_check(1,7,9,8); endfunction function input_draw_check_ae() //A1,A2,L1,L2 input_draw_check(1,8,4,13); endfunction // //------------------------------------------------------------- // // Function plot sqrt(area) vs length function plot_area(ni) clf(); plot2d(0,0,1,"010","",[-5,-10,20,10],[5,5,4,5]); // wake wo egaku ya1=sqrt(A1(ni) ); ya2=sqrt(A2(ni) ); xp=[0 0 L1(ni) L1(ni) L1(ni)+L2(ni) L1(ni)+L2(ni) L1(ni) L1(ni) 0 ]'; yp=[ya1 -ya1 -ya1 -ya2 -ya2 ya2 ya2 ya1 ya1 ]'; xpolys( xp, yp, [5 5 5 5 5 5 5 5] ); xar=[-2 -0.5]; yar=[0 0]; xarrows(xar, yar, 1, 5); xar=[ L1(ni)+L2(ni)+0.5 L1(ni)+L2(ni)+2]; xarrows(xar, yar, 1, 5); endfunction // function plot3d_ydbout() clf(); plot3d(l1,r1,ydbout,leg="l1@r1@dB"); legend('l1=(L2-L1) / (L2+L1)','r1=(A2-A1) / (A2+A1)',2); endfunction function plot3d_yfrqout() clf(); plot3d(l1,r1,yfrqout,leg="l1@r1@Hz"); legend('l1=(L2-L1) / (L2+L1)','r1=(A2-A1) / (A2+A1)',2); endfunction function plot3d_yfrqout2() clf(); plot3d(l1,r1,yfrqout,leg="l1@r1@rate"); legend('l1=(L2-L1) / (L2+L1)','r1=(A2-A1) / (A2+A1)',2); endfunction //------------------------------------------------------- // +Calculation // +switch 1 [resp_hpf]=hpf_response(); mess8b='+++ hpf'; mess8b // -switch1 // +switch 2 //[resp_yg]=yg_response(1); // load old data //mess9b='+++ load old data of yg'; mess9b='+++ make new data and save of yg'; mess9b [resp_yg]=yg_response(0); // make new data and save // -switch2 // old_ygfrq_load=0; // switch new make or load old data if old_ygfrq_load == 1 then messy8='+load old ydbout, yfrqout'; messy8 load('ydbfrq.dat','ydbout', 'yfrqout'); else for loopl=1:size_l1 // message loopl for loopr=1:size_r1 [db0,frq0,peak_count_c0,avr_pwr0,two_db_pwr0,def_frq0]=two_tubes_model_res2(loopl, loopr, 0, 0); ydbout(loopl, loopr)= two_db_pwr0; yfrqout(loopl, loopr)= def_frq0 ; end // end of loopr=1:size_r1 end // loopl=1:size_l1 // messy9='+save ydbout, yfrqout'; messy9 save('ydbfrq.dat',ydbout, yfrqout); end //if old_ygfrq_load == 1 then // -Calculation //------------------------------------------------------------ xset('window',1); // create new windows clf(); plot3d_yfrqout2(); xset('window',2); // create new windows clf(); plot3d_ydbout(); xset('window',0); //......................................................................................................................... // ADD_MENU=1; // if ADD_MENU == 1 then delmenu('functions'); addmenu('functions',['two_tubes_model_bode2(1,1)'; 'input_draw_area() '; 'plot3d_ydbout()'; 'plot3d_yfrqout2()'; 'input_draw()'; 'input_draw_check_a()' ; 'input_draw_check_ae()' ] ); functions=['two_tubes_model_bode2(1,1)'; 'input_draw_area() ' ; 'plot3d_ydbout()'; 'plot3d_yfrqout2()'; 'input_draw()'; 'input_draw_check_a()'; 'input_draw_check_ae()' ] ; end //if ADD_MENU == 1 then/ // //=========================================================================================== //