//------------------------------------------------------------------------------------------ // // Two Tubes Model for Vocal Tract and Three Tubes Model for Vocal Tract // with Noise Source of which frequency band is limited // Continuous transformation of Tube, every pitch unit // // /h/ model with Noise Source of which frequency range is limited and Two Tubes Model for Vocal Tract // /h/ needs reflection tubes effect, // besides /s/ and /k/ should not be effect by tubes, because they origin nearer to lip // In this, in order to simulate less tube effect, uses rl value control, that when rl=small value like 0.1 then tube edge reflection is zero. // /s/ model with Noise Source output only ( only lip-around source ) // /k/ model with uncontinuous source and mixed with noise source ( only around lip-around source ) // // // /m/ /n/ model using bef(filter) to erase originally peaks and using hcf(filter) to cut off higher frequecny band signal // // /r/ model from close-mouth to open-mouth plus starting mark (with high frequency signal) // // /N/ or /M/ or /nn/ model with another two tubes // bypass vocal glottis pulse to another two Tube of which frequecny is middle and rl(reflection ratio) is lower(weak) // In this, in order to simulate, use ration Rn that is defined as bypass /( all= tubes + bypass ) // // //......................................................................................... // // PLEASE PAY ATTENTION that this program may have some bugs and also you may adjust program // to fit your circumstances. If you use this program, everything done at your own risk. // Thanks. // //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ T_DEMO=1; // set T_DEMO=1 for this version // // //-------------------------------------------------------------------- // Get Scilab Version // global scilab_version_number=4; // dummy set // ver0=getversion(); p0=strindex(ver0,'-'); p1=strindex(ver0,'.'); ver1=part(ver0,(p0(1)+1):(p1(1)-1)); scilab_version_number=str2code(ver1) //--------------------------------------------------------------------- // some functions to replace from scilab-4.1.2 to scilab-5.1.1 function bode0(a,b,c) if scilab_version_number >= 5 then bode(a',b',c'); else bode(a,b,c); end endfunction function [xmode00]=x_choose0(a,b,c) if scilab_version_number >= 5 then xmode00=x_choose(a,['Please select one and double-click it'],c); else xmode00=x_choose(a,b,c); end endfunction //function [wavfile0]=tk_getfile0(a) // if scilab_version_number >= 5 then // [wavfile0]=uigetfile(["*.wav"],"","Choose any xxx.wav file name"); // else // [wavfile0]=tk_getfile(a); // cannot this. // end //endfunction //---------------------------------------------------------------------- // Sound format has a difference among scilab version. // for sound wav file if scilab_version_number == 4 then f412=1; // flag to detect scilab 4.1.2 f511=0; elseif scilab_version_number >= 5 then f412=1; f511=1; // flag to detect scilab 5.1.1 else f412=0; f511=0; end // // xmess='This program is for windows scilab-4.1.2 . If you use other version or other os, some functions may doesnot work well. ' // // //=================================================================== // scilab global variables //=================================================================== // scilab global variables global eh_var1 global eh_var2 global eh_var3 global eh_var4 global eh_var5 // //========================================================================== // // Step -1 // //------------------------------------------------------ // global variables T_QTY=2; // default set D_QTY=7; // default set memo1=''; // dummy set memo2=''; // dummy set memo3=''; // dummy set xmode2=0; // dummy set sel_demo1=0; // dummy set // dummy set A1=zeros(1,1); // dummy set A2=zeros(1,1); // dummy set A3=zeros(1,1); // dummy set L1=zeros(1,1); // dummy set L2=zeros(1,1); // dummy set L3=zeros(1,1); // dummy set // for another two tube A1_n=zeros(1,1); // dummy set A2_n=zeros(1,1); // dummy set L1_n=zeros(1,1); // dummy set L2_n=zeros(1,1); // dummy set // tclosed=zeros(1,1); // dummy set trise=zeros(1,1); // dummy set tfall=zeros(1,1); // dummy set amplitude=zeros(1,1); // dummy set co_filter=zeros(1,1); // dummy set co_filter=[ 1 , 0, 500, 0.7 ; // dummy set 2 , 1, 1000, 0.7 ; // 3, 0, 10000, 0.7 ; // 4, 2, 5000, 2.5 ; // 5, 0, 2000,0.7 ; // 6 , 2, 1250, 10. 7 , 1, 1000, 0.7 ; // 8, 0, 10000, 0.7 ; // 9, 2, 2500, 10 ]; // dummy set // for /h/ x_coef1=zeros(1,1); // dummy set y_coef1=zeros(1,1); // dummy set // for /s/ x_coef2=zeros(1,1); // dummy set y_coef2=zeros(1,1); // dummy set x_coef3=zeros(1,1); // dummy set y_coef3=zeros(1,1); // dummy set x_coef4=zeros(1,1); // dummy set y_coef4=zeros(1,1); // dummy set // for /k/ x_coef5=zeros(1,1); // dummy set y_coef5=zeros(1,1); // dummy set x_coef6=zeros(1,1); // dummy set y_coef6=zeros(1,1); // dummy set x_coef7=zeros(1,1); // dummy set y_coef7=zeros(1,1); // dummy set x_coef8=zeros(1,1); // dummy set y_coef8=zeros(1,1); // dummy set k_gain0=zeros(1,2); // dummy set k_gain0(1)=7.0; // dummy set k_gain0(2)=0.0; // dummy set // for /r/ start mark x_coef9=zeros(1,1); // dummy set y_coef9=zeros(1,1); // dummy set x_coef10=zeros(1,1); // dummy set y_coef10=zeros(1,1); // dummy set x_coef11=zeros(1,1); // dummy set y_coef11=zeros(1,1); // dummy set r_ratio=0.3; // set ratio of start mark (high frequecny signal) vs vocal source amplitude // // for /N(nn)/ n_supend_para0=zeros(1,4); // dummy set n_supend_para0(1)=1500; // dummy set (low) bef fc n_supend_para0(2)=3; // dummy set (low) bef q n_supend_para0(3)=2500.; // dummy set (high) bef fc n_supend_para0(4)=3; // dummy set (high) bef q // //-------------------------------------------------------------------------- // global select_offset_page1=0; // for select page 1 select_offset_page2=100; // for select page 2 select_offset_page3=200; // for select page 3 // function [xmode2,sel_demo1]=set_select() if T_DEMO==1 then pagex0=1; for xloop=1:30 if pagex0==1 then xmode2=x_choose0(['/a/';'/ha/';'/sa/';'/ka/';'/ta/';'/ma/';'/na/';'/ra/';'/N(nn)/';'/mu/'; '-';'(next page)';],['Please select one'],'default') if xmode2 == 12 then // when (next page) pagex0=pagex0+1; else break; end elseif pagex0==2 then xmode2=x_choose0(['/a/';'/i/';'/u/';'/e(before a)/';'/e(after a)/';'/o/';'/wa(ua)/';'/ya(ia)/';'/ae/';'/au/'; '(previous page)'; '(next page)';],['Please select one'],'default') if xmode2 == 11 then // when (previous page) pagex0=pagex0-1;; elseif xmode2 == 12 then // when (next page) pagex0=pagex0+1;; else if (xmode2 > 0) & (xmode2 < 12) then xmode2=xmode2 + select_offset_page2; end break; end elseif pagex0==3 then xmode2=x_choose0(['/a/';'/ga/';'/da/';'/pa/';'-';'-';'-';'-';'-';'-'; '-';'(previous page)';],['Please select one'],'default') if xmode2 == 12 then // when (previous page) pagex0=pagex0-1;; else if (xmode2 > 0) & (xmode2 < 12) then xmode2=xmode2 + select_offset_page3; end break; end else pagex0=1; end xmode2=0; end // for xloop=1:20 if xmode2 <= 0 then sel_demo1=1; else sel_demo1=xmode2; end xmode2=0; else xmode2=x_choose0(['edit mode';'/a/';'/ae/';'/o/';'nasal /u/';'/i/'],['Please select one'],'default') if xmode2 <= 0 then xmode2=1; else xmode2=xmode2-1; // for compatibility to 412b.sci // xmode2==0, edit mode end end endfunction // // 1st call set_select() [xmode2,sel_demo1]=set_select(); // // // // Physical Parameters (length and area) Setup // if xmode2 == 0 then // else D_QTY=7; // set dimension of list. This defines length of generated wave // Generated wave's length will be around D_QTY * (N1+N2+N3). end //----------------------------- if xmode2 == 5 then T_QTY=5; elseif xmode2 == 4 then T_QTY=4; elseif xmode2 == 3 then T_QTY=3; elseif xmode2 == 1 then T_QTY=2; end //---------------------------- //T_QTY=5; // T_QTY defines which Tubes Model is used. // set T_QTY 2 or other than below values when Two Tubes Model // set T_QTY 3 when Three Tubes Model of serial type // set T_QTY 4 when Three Tubes Model of "T" type // set T_QTY 5 when Two Tubes Model with "rl" noise if xmode2==1 then // sample length & area for /a/ from problems 3.8 in Digital Processing of Speech Signals by L.R.Rabiner and R.W.Schafer // set T_QTY=2 and use Two Tubes Model L1=[9, 9, 9, 9, 9, 9, 9]; // set list of 1st tube's length by unit is [cm] A1=[1, 1, 1, 1, 1, 1, 1]; // set list of 1st tube's area by unit is [cm^2] L2=[8, 8, 8, 8, 8, 8, 8]; // set list of 2nd tube's length by unit is [cm] A2=[7, 7, 7, 7, 7, 7, 7]; // set list of 2nd tube's area by unit is [cm^2] end if xmode2==2 then // sample length & area for /ae/ from problems 3.8 in Digital Processing of Speech Signals by L.R.Rabiner and R.W.Schafer // set T_QTY=2 Two Tubes Model L1=[4, 4, 4, 4, 4, 4, 4]; // set list of 1st tube's length by unit is [cm] A1=[1, 1, 1, 1, 1, 1, 1]; // set list of 1st tube's area by unit is [cm^2] L2=[13, 13, 13, 13, 13, 13, 13]; // set list of 2nd tube's length by unit is [cm] A2=[8, 8, 8, 8, 8, 8, 8]; // set list of 2nd tube's area by unit is [cm^2] end if xmode2==5 then // sample length & area for /i/ from problems 3.8 in Digital Processing of Speech Signals by L.R.Rabiner and R.W.Schafer // (((set T_QTY=2 and use Two Tubes Model))) // set T_QTY=5 and use Two Tubes Model with "rl" noise L1=[9, 9, 9, 9, 9, 9, 9]; // set list of 1st tube's length by unit is [cm] A1=[8, 8, 8, 8, 8, 8, 8]; // set list of 1st tube's area by unit is [cm^2] L2=[6, 6, 6, 6, 6, 6, 6]; // set list of 2nd tube's length by unit is [cm] A2=[1, 1, 1, 1, 1, 1, 1]; // set list of 2nd tube's area by unit is [cm^2] end // sample for /u/ // set T_QTY=2 and use Two Tubes Model //L1=[10, 10, 10, 10, 10, 10, 10]; // set list of 1st tube's length by unit is [cm] //A1=[7, 7, 7, 7, 7, 7, 7]; // set list of 1st tube's area by unit is [cm^2] //L2=[7, 7, 7, 7, 7, 7, 7]; // set list of 2nd tube's length by unit is [cm] //A2=[3, 3, 3, 3, 3, 3, 3]; // set list of 2nd tube's area by unit is [cm^2] if xmode2==3 then // sample length & area for /o/ modoki // set T_QTY=3 and use Three Tubes Model of serial type L1=[3, 3, 3, 3, 3, 3, 3]; // set list of 1st tube's length by unit is [cm] A1=[6, 6, 6, 6, 6, 6, 6]; // set list of 1st tube's area by unit is [cm^2] L2=[7, 7, 7, 7, 7, 7, 7]; // set list of 2nd tube's length by unit is [cm] A2=[1, 1, 1, 1, 1, 1, 1]; // set list of 2nd tube's area by unit is [cm^2] L3=[8, 8, 8, 8, 8, 8, 8]; // set list of 3rd tube's length by unit is [cm] A3=[7, 7, 7, 7, 7, 7, 7]; // set list of 3rd tube's area by unit is [cm^2] end if xmode2==4 then // sample for nasal /u/ This is not enough for nasal sound effect. // set T_QTY=4 and use Three Tubes Model of "T" type L1=[5, 5, 5, 5, 5, 5, 5]; // set list of 1st tube's length by unit is [cm] A1=[4, 4, 4, 4, 4, 4, 4]; // set list of 1st tube's area by unit is [cm^2] L2=[9, 9, 9, 9, 9, 9, 9]; // set list of 2nd tube's length by unit is [cm] A2=[2, 2, 2, 2, 2, 2, 2]; // set list of 2nd tube's area by unit is [cm^2] L3=[6, 6, 6, 6, 6, 6, 6]; // set list of 3rd tube's length by unit is [cm] A3=[2, 2, 2, 2, 2, 2, 2]; // set list of 3rd tube's area by unit is [cm^2] end //---------------------------------------------------------------------------- // for another two tubes to purpose a nose // // l1_nose_default=-0.3; // set for l1 r1_nose_default=-0.5; // set for r1 rl_nose_default=0.5; // set for rl // global l1_nose=l1_nose_default; r1_nose=r1_nose_default; rl_nose=rl_nose_default; // //---------------------------------------------------------------------------- c0=35000; // constant. the speed of sound is round 35000 cm/second //---------------------------------------------------------------- function [tu1,tu2,tu3,r1,r2,r12,r21,r31,tu1n,tu2n,r1n]=make_tu_r() tu3=zeros(1,1); // dummy set r1=zeros(1,1); // dummy set r2=zeros(1,1); // dummy set r12=zeros(1,1); // dummy set r21=zeros(1,1); // dummy set r31=zeros(1,1); // dummy set if T_QTY == 3 then //Three Tubes Model of serial type tu1=L1./c0; // delay time in 1st tube tu2=L2./c0; // delay time in 2nd tube tu3=L3./c0; // delay time in 3rd tube r1=(A2-A1)./(A2+A1); // reflection coefficient between 1st tube and 2nd tube r2=(A3-A2)./(A3+A2); // reflection coefficient between 2nd tube and 3rd tube elseif T_QTY == 4 then //Three Tubes Model of T type tu1=L1./c0; // delay time in 1st tube tu2=L2./c0; // delay time in 2nd tube tu3=L3./c0; // delay time in 3rd tube r12=((A2+A3)-A1)./(A3+A2+A1); // reflection coefficient between 1st tube and others r21=(A2-(A1+A3))./(A3+A2+A1); // reflection coefficient between 2nd tube and others r31=(A3-(A1+A2))./(A3+A2+A1); // reflection coefficient between 3rd tube and others elseif T_QTY == 2 then //Two Tubes Model tu1=L1./c0; // delay time in 1st tube tu2=L2./c0; // delay time in 2nd tube r1=(A2-A1)./(A2+A1); // reflection coefficient between 1st tube and 2nd tube elseif T_QTY == 5 then //Two Tubes Model tu1=L1./c0; // delay time in 1st tube tu2=L2./c0; // delay time in 2nd tube r1=(A2-A1)./(A2+A1); // reflection coefficient between 1st tube and 2nd tube end // for another two tube tu1n=L1_n./c0; // delay time in 1st tube tu2n=L2_n./c0; // delay time in 2nd tube r1n=(A2_n-A1_n)./(A2_n+A1_n); // reflection coefficient between 1st tube and 2nd tube // endfunction //--------------------------------------------------------------------------- // Input Glottal Volume Velocity for Two Tubes Model for Vocal Tract // A.E.Rosenberg's formula as Glottal Volume Velocity // Physical Parameters (duration) Setup // // default tclosed, trise, tfall tclosed0=5.; trise0=6.; tfall0=2.; if xmode2 == 0 then // else //D_QTY // set dimension of list. This defines length of generated wave tclosed=[5, 5, 5, 5, 5, 5, 5] ; // set glottis closed duration (off time) by unit is millisecond [ms] trise= [6, 6, 6, 6, 6, 6, 6]; // set glottis opening duration (rise time) by unit is millisecond [ms] tfall= [2, 2, 2, 2, 2, 2, 2]; // set glottis closing duration (fall time) by unit is millisecond [ms] // below (tfall is smaller) is for /i/ that needs higher frequency signals as input //tclosed=[4, 4, 4, 4, 4, 4, 4] ; // set glottis closed duration (off time) by unit is millisecond [ms] //trise= [6, 6, 6, 6, 6, 6, 6]; // set glottis opening duration (rise time) by unit is millisecond [ms] //tfall= [0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7]; // set glottis closing duration (fall time) by unit is millisecond [ms] amplitude=[1, 1, 1, 1, 1, 1, 1]; // set value of multiplier end // 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 // Control flag, process yg bef as nose effect BEF_VOCAL_TRACT=0; // when 1, ON; now is OFF; //--------------------------------------------------------------------------- // Output Sound Pressure for Two Tubes Model for Vocal Tract // 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; //------------------------------------------------------------------------------------------------------ // // setting for nose sound // another_tube_on=0; //dummy set, when this value is 1, another tube, RN, and BEF control is available fc_nose=2000; // set for another two tubes, cut off frequency of High Pass Filter by unit is [Hz] nose_hpf_times=2; // if this value is 2, process hpf two times to reduce lower frequency elements fc_nose_lpf=2500; // this lpf nose_lpf_times=2; // if this value is 2, process lpf two times to reduce higher frequency elements magic_gain_nose=10.; // add a mix gain parameter, because another two tubes sound volume is too small than vocal tract sound. // reflection coefficient between 2nd tube and mouth rl_default=0.9; // set adjust reduction response. rl is variable, however use one constant in this. // //--------------------------------------------------------------------------- // "rl" noise is, When Volume Velocity of the flow is above certain value, // Noise occured by turblent flow at narrow space near outside (mouth) // and as new sound source, it returns through Tubes Model. // On calculation, this noise will be added to at "rl" reflection inside. // In this, one easy model is used, // Noise on/off is decided value of smoothing Volume Velocity by Low Pass Filter. // Addional Noise level is constant, beside actual depeds on value of Volume Velocity. // // Property of noise Setup threshold_occur_minus = -1.1; // set threshold of Volume Velocity of minus side when noise occurs threshold_occur_plus = 1.1; // set threshold of Volume Velocity of plus side when noise occurs // two side threshold, minus side and plus side, is kunikuno-saku. mix_amplitude=0.05; // set mix amplitude at "rl" point frq_center= 3000; // set center frequency of noise by unit is [Hz] freq_band= 2000; // set frequency band of noise by unit is [Hz] nstep_fir=129; // set step number of band pass filter of linear-phase FIR filter // to transfer random noise to noise with limited frequency band // Low Pass Filter to smooth Volume Velocity fsm=1000; // set cut off frequency of Low Pass Filter by unit is [Hz] // // // reflection coefficient of 3rd tube's terminal side. 3rd tube is supposed that one side is closed. // This parameter is used only when T_QTY is set to 4 as Three Tubes Model of "T" type rl3=-0.97; // set certain a value, beside ideal value is -1 //--------------------------------------------------------------------------- // Other preparation: Sampling Frequency Setup // fs=44100; // set sampling frequency by unit is [Hz] ts=1/fs; // Min_Freq=100; // set plot minimum frequency by unit is [Hz] Max_Freq=10000; // set plot maximum frequency by unit is [Hz] // //========================================================================== // // Step 0 Edit // //------------------------------------------------------ // global variables ttl_Length=17; // set total length of combined tubes, Two Tubes, Three Tubes, and etc ttl_Area2=10; // set total area of combined tubes, Two Tubes, for dummy display ttl_Area3=14; // set total area of combined tubes, Three Tubes, for dummy display lrpa=zeros(1,1); // rlap is matrix of n-th,l1,r1,pitch(+-%),amplitude(0-1),vocal-source(0/h=1/s=2/k=3), lip-source(0/s=1/r_start=2/rl_noise=3), nose-suppression(0/n=1) // // +addition for another two tube ratio_nose=1.0; // set ratio of nose tract vs mouth vocal tract ttl_Length_n=ttl_Length * ratio_nose; // set total length of another combined tubes, Two Tubes ttl_Area2_n=2; // set total area of another combined tubes, Two Tubes // // //------------------------------------------------- function [D_QTY,T_QTY,lrpa,memo1x,co_filter]= set_QTY() txt1=['duration, pitch cycle number'; 'number of Tubes (2 or 3)';'memo']; wstr1=sprintf('%d',D_QTY); wstr2=sprintf('%d',T_QTY); wstr3=memo1; sig1=x_mdialog('set duration and type of Tube Model',txt1, [wstr1; wstr2; wstr3]); if sig1==[] then arg1=evstr(wstr1); arg2=evstr(wstr2); memo1x=wstr3; else arg1=evstr(sig1(1)); arg2=evstr(sig1(2)); memo1x=sig1(3); end D_QTY=int(arg1); T_QTY=int(arg2); // if (T_QTY == 2) | (T_QTY == 5) then lrpa=zeros(D_QTY,8); for v=1:D_QTY lrpa(v,1)=v; lrpa(v,2)=0.; lrpa(v,3)=0.8; lrpa(v,4)=0.; lrpa(v,5)=1.; lrpa(v,6)=0; lrpa(v,7)=0; lrpa(v,8)=9; end elseif T_QTY == 3 then lrpa=zeros(D_QTY,10); for v=1:D_QTY lrpa(v,1)=v; lrpa(v,2)=0.4; lrpa(v,3)=-0.7; lrpa(v,4)=0.07; lrpa(v,5)=0.75; lrpa(v,6)=0.; lrpa(v,7)=1.; lrpa(v,8)=0; lrpa(v,9)=0; lrpa(v,10)=9; end end // endfunction //-------------------------------------------------------------------------- function [A1,A2,L1,L2]=cal_AL2(r1,l1) // for Two Tubes Model 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_Area2/ 2.; A1(v)= ttl_Area2 - A2(v); end endfunction //-------------------------------------------------------------------------- function [A1_n,A2_n,L1_n,L2_n]=cal_AL2_n(r1_n,l1_n) // for another two tube // for Two Tubes Model //r1s=size(r1); //size_r1=r1s(2); // size of array r1 //l1s=size(l1); //size_l1=l1s(2); // size of array l1 // size_r1=D_QTY; size_l1=D_QTY; L1_n=zeros(1,size_r1); L2_n=zeros(1,size_l1); for v=1:size_l1 L2_n(v)= ( 1. + l1_n ) * ttl_Length_n / 2.; L1_n(v)= ttl_Length_n - L2_n(v); end A1_n=zeros(1,size_r1); A2_n=zeros(1,size_r1); for v=1:size_r1 A2_n(v)= ( 1. + r1_n ) * ttl_Area2_n/ 2.; A1_n(v)= ttl_Area2_n - A2_n(v); end endfunction //-------------------------------------------------------------------------- function [A1,A2,A3,L1,L2,L3]=cal_AL3(r1,l1,r2,l2) // for Three Tubes Model // tube #1 #2 #3 // ___ ___ ___ // glottis=> |___|___|___| => out // L1 L2 L3 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); L3=zeros(1,size_l1); for v=1:size_l1 y0=l1(v) * l2(v) + ( l1(v) - l2(v) ) + 3.; L1(v)= ( l2(v) - 1. ) * ( l1(v) - 1. ) / y0 * ttl_Length; L2(v)= -1.0 * ( 1. + l1(v) ) * ( l2(v) - 1. ) / y0 * ttl_Length; L3(v)= ( l2(v) + 1. ) * ( l1(v) + 1. ) / y0 * ttl_Length; end A1=zeros(1,size_r1); A2=zeros(1,size_r1); A3=zeros(1,size_r1); for v=1:size_r1 y0=r1(v) * r2(v) + ( r1(v) - r2(v) ) + 3.; A1(v)= ( r2(v) - 1. ) * ( r1(v) - 1. ) / y0 * ttl_Area3; A2(v)= -1.0 * ( 1. + r1(v) ) * ( r2(v) - 1. ) / y0 * ttl_Area3; A3(v)= ( r2(v) + 1. ) * ( r1(v) + 1. ) / y0 * ttl_Area3; end endfunction //------------------------------------------------------------------ function [lrpar,memo2x]=edit_lrpa() // if (T_QTY == 2) | (T_QTY == 5) then [lrpar]=x_matrix('Edit no, l1, r1, pitch(+-%), amplitude(0-1), vocal-source(0/h=1/s=2/k=3/t=4), lip-around-source(0/s=1), rl(if=9,use rl_default)',lrpa); memo2x=' no, l1, r1, pitch(+-%), amplitude(0-1), vocal-source(0/h=1/s=2/k=3),lip-around-source(0/s=1),rl(if=9,use rl_default)'; elseif T_QTY == 3 then [lrpar]=x_matrix('Edit no, l1, r1, l2, r2 ,pitch(+-%), amplitude(0-1), vocal-source(0/h=1/s=2/k=3/t=4), lip-around-source(0/s=1), rl(if=9,use rl_default)',lrpa); memo2x=' no, l1, r1, l2, r2 ,pitch(+-%), amplitude(0-1), vocal-source(0/h=1/s=2/k=3),lip-around-source(0/s=1),rl(if=9,use rl_default)'; end // if lrpar==[] then // set old setting s0=size(lrpa); lrpar=zeros(s0(1),s0(2)); for v1=1:s0(1) for v2=1:s0(2) lrpar(v1,v2)=lrpa(v1,v2); end end else // end //... endfunction //------------------------------------------------------------------ function [co_filterx,memo3x]=edit_h_co_filter() // s0=size(co_filter); co_filter2=zeros(1,s0(2)); for v=1:s0(2) co_filter2(1,v)=co_filter(1,v); end [co_filterx2]=x_matrix('Edit no(/h/=1), type(LPF=0), f0, Q ,',co_filter2); memo3x='no, type(LPF=0,HPF=1,BPF=2), f0, Q'; // if co_filterx2==[] then // set old setting s0=size(co_filter); co_filterx=zeros(s0(1),s0(2)); for v1=1:s0(1) for v2=1:s0(2) co_filterx(v1,v2)=co_filter(v1,v2); end end else s0=size(co_filter); co_filterx=zeros(s0(1),s0(2)); for v1=1:s0(1) for v2=1:s0(2) if v1==1 then co_filterx(v1,v2)=co_filterx2(1,v2); else co_filterx(v1,v2)=co_filter(v1,v2); end end end end //... endfunction //------------------------------------------------------------------ function [co_filterx,memo3x]=edit_s_co_filter() // s0=size(co_filter); co_filter2=zeros(3,s0(2)); for w=2:4 for v=1:s0(2) co_filter2(w-1,v)=co_filter(w,v); end end [co_filterx2]=x_matrix('Edit no(/s/=2,3,4), type(LPF=0,HPF=1,BPF=2), f0, Q ,',co_filter2); memo3x='no, type(LPF=0,HPF=1,BPF=2), f0, Q'; // if co_filterx2==[] then // set old setting s0=size(co_filter); co_filterx=zeros(s0(1),s0(2)); for v1=1:s0(1) for v2=1:s0(2) co_filterx(v1,v2)=co_filter(v1,v2); end end else s0=size(co_filter); co_filterx=zeros(s0(1),s0(2)); for v1=1:s0(1) for v2=1:s0(2) if (v1>=2) & (v1<=4) then co_filterx(v1,v2)=co_filterx2(v1-1,v2); else co_filterx(v1,v2)=co_filter(v1,v2); end end end end //... endfunction //------------------------------------------------------------------- function [co_filterx,memo3x]=edit_r_co_filter() // s0=size(co_filter); co_filter2=zeros(3,s0(2)); for w=7:9 for v=1:s0(2) co_filter2(w-6,v)=co_filter(w,v); end end [co_filterx2]=x_matrix('Edit no(/r/ start mark=7,8,9), type(LPF=0,HPF=1,BPF=2), f0, Q ,',co_filter2); memo3x='no, type(LPF=0,HPF=1,BPF=2), f0, Q'; // if co_filterx2==[] then // set old setting s0=size(co_filter); co_filterx=zeros(s0(1),s0(2)); for v1=1:s0(1) for v2=1:s0(2) co_filterx(v1,v2)=co_filter(v1,v2); end end else s0=size(co_filter); co_filterx=zeros(s0(1),s0(2)); for v1=1:s0(1) for v2=1:s0(2) if (v1>=7) & (v1<=9) then co_filterx(v1,v2)=co_filterx2(v1-6,v2); else co_filterx(v1,v2)=co_filter(v1,v2); end end end end //... endfunction //------------------------------------------------------------------ function [co_filterx,memo3x,k_gain0x]=edit_k_co_filter() // s0=size(co_filter); co_filter2=zeros(1,s0(2)+1); for w=5:6 for v=1:s0(2) co_filter2(w-4,v)=co_filter(w,v); end co_filter2(w-4,s0(2)+1)=k_gain0(w-4); end [co_filterx2]=x_matrix('Edit no(/k/=5,6), type(LPF=0/BPF=2), f0, Q ,beat gain',co_filter2); memo3x='no, type(LPF=0,HPF=1,BPF=2), f0, Q'; // if co_filterx2==[] then // set old setting s0=size(co_filter); co_filterx=zeros(s0(1),s0(2)); for v1=1:s0(1) for v2=1:s0(2) co_filterx(v1,v2)=co_filter(v1,v2); end end k_gain0x=zeros(1,2); for v=1:2 k_gain0x(v)=k_gain0(v); end else s0=size(co_filter); co_filterx=zeros(s0(1),s0(2)); for v1=1:s0(1) for v2=1:s0(2) if (v1>=5) & (v1<=6) then co_filterx(v1,v2)=co_filterx2(v1-4,v2); else co_filterx(v1,v2)=co_filter(v1,v2); end end end for v=1:2 k_gain0x(v)=co_filterx2(v,(s0(2)+1)); end end //... endfunction //-------------------------------------------------------------------------- function [n_supend_para0x,memo4x]=edit_n_supend_para() // s0=size(n_supend_para0); n_supend_para0y=zeros(1,s0(2)); for w=1:1 for v=1:s0(2) n_supend_para0y(w,v)=n_supend_para0(w,v); end end [n_supend_para0z]=x_matrix('(low)bef fc, q, (high)bef fc, q',n_supend_para0y); memo4x='(low)bef fc, q, (high)bef fc, q'; // if n_supend_para0z==[] then // set old setting s0=size(n_supend_para0); n_supend_para0x=zeros(s0(1),s0(2)); for v1=1:s0(1) for v2=1:s0(2) n_supend_para0x(v1,v2)=n_supend_para0(v1,v2); end end else s0=size(n_supend_para0); n_supend_para0x=zeros(s0(1),s0(2)); for v1=1:s0(1) for v2=1:s0(2) if (v1>=1) & (v1<=1) then n_supend_para0x(v1,v2)=n_supend_para0z(v1,v2); else n_supend_para0x(v1,v2)=n_supend_para0(v1,v2); end end end end //... endfunction //-------------------------------------------------------------------------- function [A1,A2,A3,L1,L2,L3,A1_n,A2_n,L1_n,L2_n,tclosed,trise,tfall,amplitude]=make_paras() A1=zeros(1,1); // dummy set A2=zeros(1,1); // dummy set A3=zeros(1,1); // dummy set L1=zeros(1,1); // dummy set L2=zeros(1,1); // dummy set L3=zeros(1,1); // dummy set // r1=zeros(1,D_QTY); r2=zeros(1,D_QTY); l1=zeros(1,D_QTY); l2=zeros(1,D_QTY); // s0=size(lrpa); if (T_QTY == 2) | (T_QTY ==5) then vp=4; va=5; for v=1:D_QTY l1(v)=lrpa(v,2); r1(v)=lrpa(v,3); end [A1,A2,L1,L2]=cal_AL2(r1,l1); elseif T_QTY == 3 then vp=6; va=7; for v=1:D_QTY l1(v)=lrpa(v,2); r1(v)=lrpa(v,3); l2(v)=lrpa(v,4); r2(v)=lrpa(v,5); end [A1,A2,A3,L1,L2,L3]=cal_AL3(r1,l1,r2,l2); end // common tclosed=zeros(1,D_QTY); trise=zeros(1,D_QTY); tfall=zeros(1,D_QTY); amplitude=zeros(1,D_QTY); pl0=tclosed0+trise0+tfall0; for v=1:D_QTY amplitude(v)=lrpa(v,va); // adjp=lrpa(v,vp); delta= pl0 * adjp / 100.; if delta > 0 then // plus, more high speed if delta < tclosed0 then tclosed(v) = tclosed0 - delta; else tclosed(v) = 0.0; disp('waring: cannot pitch control because tclosed0 is not enough'); end else // when minus, more low speed tclosed(v)=tclosed0+delta; end trise(v)=trise0; tfall(v)=tfall0; end // // for another two tube [A1_n,A2_n,L1_n,L2_n]=cal_AL2_n(r1_nose,l1_nose); // disp(' - makeing parameters was done. '); // endfunction //-------------------------------------------------------------------------- function save_lrpa() filo=input(' + enter file name for save =>',["string"]); save( filo,D_QTY,T_QTY,lrpa,memo1,co_filter,k_gain0,n_supend_para0); disp(' - save was done. '); endfunction function [D_QTY,T_QTY,lrpa,memo1,co_filter,k_gain0,n_supend_para0]=load_lrpa() if scilab_version_number >= 5 then filo=uigetfile("","","Choose Choose load file"); else filo=tk_getfile(Title="Choose load file"); end load( filo,'D_QTY','T_QTY','lrpa','co_filter','k_gain0','n_supend_para0'); memo1=''; load( filo,'memo1'); disp(' - load was done. '); endfunction function [l1_nose,r1_nose,rl_nose,n_supend_para0]=set_nose(sel_demo1) if sel_demo1 == 6 then // /ma/ // two peaks with one dip between the peaks l1_nose=0.; // set for l1 r1_nose=0.; // set for r1 rl_nose=0.7; // set for rl n_supend_para0=[2000,1,2500,3]; elseif sel_demo1 == 7 then // /na/ l1_nose=0.; // set for l1 r1_nose=0.; // set for r1 rl_nose=0.7; // set for rl n_supend_para0=[1000,1,2500,3]; elseif sel_demo1 == 9 then // /NN/ l1_nose=0.; // set for l1 r1_nose=0.; // set for r1 rl_nose=0.7; // set for rl n_supend_para0=[1000,1,2500,3]; else l1_nose=l1_nose_default; // set for l1 r1_nose=r1_nose_default; // set for r1 rl_nose=rl_nose_default; // set for rl n_supend_para0=[1500,3,2500,3]; end endfunction //-------------------------------------------------------------------------- // // //========================================================================== // // Step 1 Input Glottal Volume Velocity Generation // //------------------------------------------------------ // global variables rg=zeros(1,1); // dummy set ygn=zeros(1,1); // dummy set for bypass like nose rl=zeros(1,1); // dummy set yg=zeros(1,1); // dummy set N1=zeros(1,1); // dummy set N2=zeros(1,1); // dummy set N3=zeros(1,1); // dummy set LL=zeros(1,1); // dummy set length0=0; // function [rg,ygn,yg,rl,N1,N2,N3,LL,length0]=step1() disp('+++enter step 1 Input Glottal Volume Velocity Generation'); // length0 is wave length length0=0; N1= zeros(1, D_QTY); N2= zeros(1, D_QTY); N3= zeros(1, D_QTY); LL= zeros(1, D_QTY); for loop=1:D_QTY N1(loop)= int( tclosed(loop) / 1000 / ts ); N2(loop)= int( trise(loop) / 1000 / ts); N3(loop)= int( tfall(loop) / 1000 / ts); LL(loop)= N1(loop)+N2(loop)+N3(loop); length0 = length0 + LL(loop); end // for bypass if (T_QTY == 2) | (T_QTY ==5) then idy0=9; elseif T_QTY == 3 then idy0=11; end // yg is Glottal Volume Velocity. rg is reflection coefficient between glottis and 1st tube yg= zeros(1,length0+1); rg= zeros(1,length0+1); ygn= zeros(1,length0+1); yg_temp= zeros(1,length0+1); ygn_temp= zeros(1,length0+1); rl= zeros(1,length0+1); tc0=1; yg(tc0)=0.; yg(tc0)=amplitude(1) * yg(tc0); for loop=1:D_QTY for t0=1:int(LL(loop)) tc0=tc0+1; if t0 < N1(loop) then yg(tc0) =0; rg(tc0)=rg_closed; elseif t0 <= (N2(loop) + N1(loop)) then yg(tc0)= 0.5 * ( 1 - cos( ( %pi / N2(loop) ) * (t0 - N1(loop)) ) ); rg(tc0)=rg_rise; elseif t0 <= (N3(loop)+N2(loop)+N1(loop)) then yg(tc0)= cos( ( %pi / ( 2 * N3(loop) )) * ( t0 - (N2(loop)+N1(loop)) ) ); rg(tc0)=rg_fall; else yg(tc0) =0; rg(tc0)=rg_closed; end if another_tube_on == 1 then RNx0=lrpa(loop,idy0); else RNx0=0.; end temp0=amplitude(loop) * yg(tc0); yg(tc0)=(1.- RNx0) * temp0; yg_temp(tc0)=yg(tc0); ygn(tc0)= RNx0 * temp0; ygn_temp(tc0)=ygn(tc0); rl(tc0)=rl_default; end // t0 end // loop // // process bef filtering on ygn as another tubes input if another_tube_on == 1 then if (T_QTY == 2) | (T_QTY ==5) then idy0=10; elseif T_QTY == 3 then idy0=12; end spp0=1; epp0=0; eppx=0; sd_qty=1; ed_qty=0; xcode0=lrpa(1,idy0); spp1=0; epp1=0; xcode1=-1; for loop=1:D_QTY // when only 1 D_QTY if D_QTY==1 then epp0=int(LL(loop)); ed_qty=1; // when last frame elseif loop==D_QTY then if lrpa(loop,idy0) == xcode0 then epp0=eppx+int(LL(loop)); ed_qty=D_QTY; else epp0=eppx; ed_qty=loop-1; xcode1= lrpa(loop,idy0); spp1=eppx+1; epp1=eppx+int(LL(loop)); end // find lrpa term end elseif lrpa(loop,idy0) <> xcode0 then epp0=eppx; ed_qty=loop-1; end // process if epp0 > spp0 then if xcode0 == 1 then // low length1= epp0 - spp0 + 1; y_ran1=zeros(1, length1); for w=1:length1 y_ran1(w)=ygn_temp(w+spp0-1); end [out0]=do_iir_filtering3(y_ran1,x_coef7,y_coef7); tc0=spp0; for v=sd_qty:ed_qty for t0=1:int(LL(v)) ygn(tc0)=out0(tc0-spp0+1); tc0=tc0+1; end end elseif xcode0 == 2 then // high length1= epp0 - spp0 + 1; y_ran1=zeros(1, length1); for w=1:length1 y_ran1(w)=ygn_temp(w+spp0-1); end [out0]=do_iir_filtering3(y_ran1,x_coef8,y_coef8); tc0=spp0; for v=sd_qty:ed_qty for t0=1:int(LL(v)) ygn(tc0)=out0(tc0-spp0+1); tc0=tc0+1; end end end //common spp0=epp0+1; epp0=0; sd_qty=loop; end // process for only last frame, as same as process above if epp1 > spp1 then if xcode1 == 1 then // low length1= epp1 - spp1 + 1; y_ran1=zeros(1, length1); for w=1:length1 y_ran1(w)=ygn_temp(w+spp1-1); end [out0]=do_iir_filtering3(y_ran1,x_coef7,y_coef7); for v=spp1:epp1 ygn(v)=out0(v-spp1+1); end elseif xcode0 == 2 then // high length1= epp1 - spp1 + 1; y_ran1=zeros(1, length1); for w=1:length1 y_ran1(w)=ygn_temp(w+spp1-1); end [out0]=do_iir_filtering3(y_ran1,x_coef8,y_coef8); for v=spp1:epp1 ygn(v)=out0(v-spp1+1); end end end // stack up xcode0=lrpa(loop,idy0); eppx=eppx+int(LL(loop)); end // loop // process bef filtering on yg to vocal tract if BEF_VOCAL_TRACT == 1 then if (T_QTY == 2) | (T_QTY ==5) then idy0=10; elseif T_QTY == 3 then idy0=12; end spp0=1; epp0=0; eppx=0; sd_qty=1; ed_qty=0; xcode0=lrpa(1,idy0); spp1=0; epp1=0; xcode1=-1; for loop=1:D_QTY // when only 1 D_QTY if D_QTY==1 then epp0=int(LL(loop)); ed_qty=1; // when last frame elseif loop==D_QTY then if lrpa(loop,idy0) == xcode0 then epp0=eppx+int(LL(loop)); ed_qty=D_QTY; else epp0=eppx; ed_qty=loop-1; xcode1= lrpa(loop,idy0); spp1=eppx+1; epp1=eppx+int(LL(loop)); end // find lrpa term end elseif lrpa(loop,idy0) <> xcode0 then epp0=eppx; ed_qty=loop-1; end // process if epp0 > spp0 then if xcode0 == 1 then // low length1= epp0 - spp0 + 1; y_ran1=zeros(1, length1); for w=1:length1 y_ran1(w)=yg_temp(w+spp0-1); end [out0]=do_iir_filtering3(y_ran1,x_coef7,y_coef7); tc0=spp0; for v=sd_qty:ed_qty for t0=1:int(LL(v)) yg(tc0)=out0(tc0-spp0+1); tc0=tc0+1; end end elseif xcode0 == 2 then // high length1= epp0 - spp0 + 1; y_ran1=zeros(1, length1); for w=1:length1 y_ran1(w)=yg_temp(w+spp0-1); end [out0]=do_iir_filtering3(y_ran1,x_coef8,y_coef8); tc0=spp0; for v=sd_qty:ed_qty for t0=1:int(LL(v)) yg(tc0)=out0(tc0-spp0+1); tc0=tc0+1; end end end //common spp0=epp0+1; epp0=0; sd_qty=loop; end // process for only last frame, as same as process above if epp1 > spp1 then if xcode1 == 1 then // low length1= epp1 - spp1 + 1; y_ran1=zeros(1, length1); for w=1:length1 y_ran1(w)=yg_temp(w+spp1-1); end [out0]=do_iir_filtering3(y_ran1,x_coef7,y_coef7); for v=spp1:epp1 yg(v)=out0(v-spp1+1); end elseif xcode0 == 2 then // high length1= epp1 - spp1 + 1; y_ran1=zeros(1, length1); for w=1:length1 y_ran1(w)=yg_temp(w+spp1-1); end [out0]=do_iir_filtering3(y_ran1,x_coef8,y_coef8); for v=spp1:epp1 yg(v)=out0(v-spp1+1); end end end // stack up xcode0=lrpa(loop,idy0); eppx=eppx+int(LL(loop)); end // loop end //if BEF_VOCAL_TRACT == 1 then end // if another_tube_on == 1 then // /k/ org_mix_rate0=1.0; // /h/ model with Noise Source of which frequency range is limited if (T_QTY == 2) | (T_QTY ==5) then idy0=6; elseif T_QTY == 3 then idy0=8; end spp0=1; epp0=0; eppx=0; sd_qty=1; ed_qty=0; xcode0=lrpa(1,idy0); spp1=0; epp1=0; xcode1=-1; for loop=1:D_QTY // when only 1 D_QTY if D_QTY==1 then epp0=int(LL(loop)); ed_qty=1; // when last frame elseif loop==D_QTY then if lrpa(loop,idy0) == xcode0 then epp0=eppx+int(LL(loop)); ed_qty=D_QTY; else epp0=eppx; ed_qty=loop-1; xcode1= lrpa(loop,idy0); spp1=eppx+1; epp1=eppx+int(LL(loop)); end // find lrpa term end elseif lrpa(loop,idy0) <> xcode0 then epp0=eppx; ed_qty=loop-1; end // process if epp0 > spp0 then if xcode0 == 1 then // /ha/ rand('normal'); // random generator is set to a Gaussian length1= epp0 - spp0 + 1; y_ran1=rand(1, length1); [out0]=do_iir_filtering3(y_ran1,x_coef1,y_coef1); tc0=spp0; for v=sd_qty:ed_qty for t0=1:int(LL(v)) yg(tc0)=amplitude(v) * out0(tc0-spp0+1); tc0=tc0+1; end end elseif xcode0 == 2 then // /s/ rand('normal'); // random generator is set to a Gaussian length1= epp0 - spp0 + 1; y_ran1=rand(1, length1); [out0x]=do_iir_filtering3(y_ran1,x_coef2,y_coef2); [out0]=do_iir_filtering3(out0x,x_coef3,y_coef3); tc0=spp0; for v=sd_qty:ed_qty for t0=1:int(LL(v)) yg(tc0)=amplitude(v) * out0(tc0-spp0+1); tc0=tc0+1; end end elseif (xcode0 == 3) | (xcode0 == 4) then // /k/ et /t/ rand('normal'); // random generator is set to a Gaussian length1= epp0 - spp0 + 1; y_ran1=rand(1, length1); [out0]=do_iir_filtering3(y_ran1,x_coef5,y_coef5); [out1]=do_iir_filtering3(out0,x_coef6,y_coef6); tc0=spp0; if xcode0 == 3 then temp0=k_gain0(1); // only /k/ addtion else temp0=0.; end for v=sd_qty:ed_qty for t0=1:int(LL(v)) yg(tc0)=amplitude(v) * (org_mix_rate0 * out0(tc0-spp0+1) + temp0 * out1(tc0-spp0+1) ); tc0=tc0+1; end end end //common spp0=epp0+1; epp0=0; sd_qty=loop; end // process for only last frame, as same as process above if epp1 > spp1 then if xcode1 == 1 then // /ha/ rand('normal'); // random generator is set to a Gaussian length1= epp1 - spp1 + 1; y_ran1=rand(1, length1); [out0]=do_iir_filtering3(y_ran1,x_coef1,y_coef1); for v=spp1:epp1 yg(v)=amplitude(loop) * out0(v-spp1+1); end elseif xcode0 == 2 then // /s/ rand('normal'); // random generator is set to a Gaussian length1= epp1 - spp1 + 1; y_ran1=rand(1, length1); [out0x]=do_iir_filtering3(y_ran1,x_coef2,y_coef2); [out0]=do_iir_filtering3(out0x,x_coef3,y_coef3); for v=spp1:epp1 yg(v)=amplitude(loop) * out0(v-spp1+1); end elseif (xcode0 == 3) | (xcode0 == 4) then // /k/ et /t/ rand('normal'); // random generator is set to a Gaussian length1= epp1 - spp1 + 1; y_ran1=rand(1, length1); [out0]=do_iir_filtering3(y_ran1,x_coef5,y_coef5); [out1]=do_iir_filtering3(out0,x_coef6,y_coef6); if xcode0 == 3 then temp0=k_gain0(1); // only /k/ addtion else temp0=0.; end for v=spp1:epp1 yg(v)=amplitude(loop) * (org_mix_rate0 * out0(v-spp1+1) + temp0 * out1(v-spp1+1)); end end end // stack up xcode0=lrpa(loop,idy0); eppx=eppx+int(LL(loop)); end // loop if (T_QTY == 2) | (T_QTY ==5) then idy0=8; elseif T_QTY == 3 then idy0=10; end spp0=1; epp0=0; eppx=0; sd_qty=1; ed_qty=0; xcode0=lrpa(1,idy0); spp1=0; epp1=0; xcode1=-1; for loop=1:D_QTY // when only 1 D_QTY if D_QTY==1 then epp0=int(LL(loop)); ed_qty=1; // when last frame elseif loop==D_QTY then if lrpa(loop,idy0) == xcode0 then epp0=eppx+int(LL(loop)); ed_qty=D_QTY; else epp0=eppx; ed_qty=loop-1; xcode1= lrpa(loop,idy0); spp1=eppx+1; epp1=eppx+int(LL(loop)); end // find lrpa term end elseif lrpa(loop,idy0) <> xcode0 then epp0=eppx; ed_qty=loop-1; end // process if epp0 > spp0 then if xcode0 <> 9 then // other than xcode0 = 9 tc0=spp0; for v=sd_qty:ed_qty for t0=1:int(LL(v)) rl(tc0)=xcode0; tc0=tc0+1; end end end //common spp0=epp0+1; epp0=0; sd_qty=loop; end // process for only last frame, as same as process above if epp1 > spp1 then if xcode1 <> 9 then // other than xcode0 = 9 for v=spp1:epp1 rl(v)=xcode1; end end end // stack up xcode0=lrpa(loop,idy0); eppx=eppx+int(LL(loop)); end // loop endfunction // step1 // function save yg (Glottal Volume Velocity) as one sound wave file function save_as_sound_file_yg( filename) max01=0.; for loop=1:length0 if abs( yg(loop) ) > max01 then max01 = abs( yg(loop)); end end y3out2=zeros(1,length0+1); if max01 > 0. then for loop=1:length0 y3out2(loop)= yg(loop) / max01; end end wavwrite(y3out2, fs, 16, filename ); endfunction function save_as_sound_file_ygn( filename) max01=0.; for loop=1:length0 if abs( y3outn(loop) ) > max01 then max01 = abs( y3outn(loop)); end end y3out2=zeros(1,length0+1); if max01 > 0. then for loop=1:length0 y3out2(loop)= y3outn(loop) / max01; end end wavwrite(y3out2, fs, 16, filename ); endfunction // Function for yg (Glottal Volume Velocity) frequency-phase response function yg_bode(ni) frq=[100:25:Max_Freq]; frqs=size(frq); resp=zeros(1,frqs(2)); stp=0; if ni > 1 then for v=1:(ni-1) stp=stp + LL(v); end end stp=stp+N1(ni); y00=0.; for v=1:(N2(ni)+N3(ni)+1) y00= y00 + yg(stp + v); end for w=1:frqs(2) xw=2 * %pi * frq(w) * ts; yi=0; for v=1:(N2(ni)+N3(ni)+1) yi= yi + yg(stp + v) * %e^(-1. * xw * (v-1) * %i); end // for v resp(w)=yi / y00; end // for w [db,phi] =dbphi(resp); clf(); bode0(frq,db,phi); endfunction // // End of Step 1 Input Glottal Volume Velocity Generation // //========================================================================== // // Step 2 Three Tubes/Two Tubes Model for Vocal Tract Calculation // //------------------------------------------------------ // global variables M1= zeros(1, 1); // dummy set M2= zeros(1, 1); // dummy set M3= zeros(1, 1); // dummy set y2tm= zeros(1,1); // dummy set y2tmn= zeros(1,1); // dummy set // function [y2tm,M1,M2,M3]=step2() if T_QTY == 3 then //Three Tubes Model of serial type disp('+++enter step 2 Three Tubes Model of serial type for Vocal Tract Calculation'); // yt is output of Two Tubes Model for Vocal Tract M1= zeros(1, D_QTY); M2= zeros(1, D_QTY); M3= zeros(1, D_QTY); y2tm= zeros(1,length0+1); for loop=1:D_QTY M1(loop)= int( tu1(loop) / ts ) + 1; M2(loop)= int( tu2(loop) / ts ) + 1; M3(loop)= int( tu3(loop) / ts ) + 1; end MAX_M1=M1(1); MAX_M2=M2(1); MAX_M3=M3(1); for loop=1:D_QTY if M1(loop) > MAX_M1 then MAX_M1= M1(loop) end if M2(loop) > MAX_M2 then MAX_M2= M2(loop) end if M3(loop) > MAX_M3 then MAX_M3= M3(loop) end end ya1= zeros(1,MAX_M1); ya2= zeros(1,MAX_M1); yb1= zeros(1,MAX_M2); yb2= zeros(1,MAX_M2); yc1= zeros(1,MAX_M3); yc2= zeros(1,MAX_M3); tc0=1; ya1(tc0)=0.; ya2(tc0)=0.; yb1(tc0)=0.; yb2(tc0)=0.; yc1(tc0)=0.; yc2(tc0)=0.; y2tm(tc0)=0.; for loop=1:D_QTY for t0=1:int(LL(loop)) tc0=tc0+1; for jl=MAX_M1:-1:2 ya1(jl)=ya1(jl-1); ya2(jl)=ya2(jl-1); end for jl=MAX_M2:-1:2 yb1(jl)=yb1(jl-1); yb2(jl)=yb2(jl-1); end for jl=MAX_M3:-1:2 yc1(jl)=yc1(jl-1); yc2(jl)=yc2(jl-1); end ya1(1)= ((1. + rg(tc0) ) / 2.) * yg(tc0) + rg(tc0) * ya2(M1(loop)); ya2(1)= -1. * r1(loop) * ya1(M1(loop)) + ( 1. - r1(loop) ) * yb2(M2(loop)); yb1(1)= ( 1 + r1(loop) ) * ya1(M1(loop)) + r1(loop) * yb2(M2(loop)); yb2(1)= -1. * r2(loop) * yb1(M2(loop)) + ( 1. - r2(loop) ) * yc2(M3(loop)); yc1(1)= ( 1 + r2(loop) ) * yb1(M2(loop)) + r2(loop) * yc2(M3(loop)); yc2(1)= -1. * rl(tc0) * yc1(M3(loop)); y2tm(tc0)= (1 + rl(tc0)) * yc1(M3(loop)); end // t0 end //loop elseif T_QTY == 4 then //Three Tubes Model of T type disp('+++enter step 2 Three Tubes Model of T type for Vocal Tract Calculation'); // yt is output of Two Tubes Model for Vocal Tract M1= zeros(1, D_QTY); M2= zeros(1, D_QTY); M3= zeros(1, D_QTY); y2tm= zeros(1,length0+1); for loop=1:D_QTY M1(loop)= int( tu1(loop) / ts ) + 1; M2(loop)= int( tu2(loop) / ts ) + 1; M3(loop)= int( tu3(loop) / ts ) + 1; end MAX_M1=M1(1); MAX_M2=M2(1); MAX_M3=M3(1); for loop=1:D_QTY if M1(loop) > MAX_M1 then MAX_M1= M1(loop) end if M2(loop) > MAX_M2 then MAX_M2= M2(loop) end if M3(loop) > MAX_M3 then MAX_M3= M3(loop) end end ya1= zeros(1,MAX_M1); ya2= zeros(1,MAX_M1); yb1= zeros(1,MAX_M2); yb2= zeros(1,MAX_M2); yc1= zeros(1,MAX_M3); yc2= zeros(1,MAX_M3); tc0=1; ya1(tc0)=0.; ya2(tc0)=0.; yb1(tc0)=0.; yb2(tc0)=0.; yc1(tc0)=0.; yc2(tc0)=0.; y2tm(tc0)=0.; for loop=1:D_QTY for t0=1:int(LL(loop)) tc0=tc0+1; for jl=MAX_M1:-1:2 ya1(jl)=ya1(jl-1); ya2(jl)=ya2(jl-1); end for jl=MAX_M2:-1:2 yb1(jl)=yb1(jl-1); yb2(jl)=yb2(jl-1); end for jl=MAX_M3:-1:2 yc1(jl)=yc1(jl-1); yc2(jl)=yc2(jl-1); end ya1(1)= ((1. + rg(tc0) ) / 2.) * yg(tc0) + rg(tc0) * ya2(M1(loop)); ya2(1)= -1. * r12(loop) * ya1(M1(loop)) + ( 1. - r12(loop) ) * yb2(M2(loop)) + ( 1. - r12(loop)) * yc2(M3(loop)); yb1(1)= ( 1 + r21(loop) ) * ya1(M1(loop)) + r21(loop) * yb2(M2(loop)) + ( 1. + r21(loop)) * yc2(M3(loop)); yb2(1)= -1. * rl(tc0) * yb1(M2(loop)); y2tm(tc0)= (1 + rl(tc0)) * yb1(M2(loop)); yc1(1)= ( 1. + r31(loop) ) * ya1(M1(loop)) + r31(loop) * yc2(M3(loop)) + ( 1. + r31(loop) ) * yb2(M2(loop)) ; yc2(1)= -1. * rl3 * yc1(M3(loop)); end // t0 end //loop else //Two Tubes Model disp('+++enter step 2 Two Tubes Model for Vocal Tract Calculation'); // yt is output of Two Tubes Model for Vocal Tract M1= zeros(1, D_QTY); M2= zeros(1, D_QTY); y2tm= zeros(1,length0+1); for loop=1:D_QTY M1(loop)= int( tu1(loop) / ts ) + 1; M2(loop)= int( tu2(loop) / ts ) + 1; end MAX_M1=M1(1); MAX_M2=M2(1); for loop=1:D_QTY if M1(loop) > MAX_M1 then MAX_M1= M1(loop) end if M2(loop) > MAX_M2 then MAX_M2= M2(loop) end end ya1= zeros(1,MAX_M1); ya2= zeros(1,MAX_M1); yb1= zeros(1,MAX_M2); yb2= zeros(1,MAX_M2); tc0=1; ya1(tc0)=0.; ya2(tc0)=0.; yb1(tc0)=0.; yb2(tc0)=0.; y2tm(tc0)=0.; for loop=1:D_QTY for t0=1:int(LL(loop)) tc0=tc0+1; for jl=MAX_M1:-1:2 ya1(jl)=ya1(jl-1); ya2(jl)=ya2(jl-1); end for jl=MAX_M2:-1:2 yb1(jl)=yb1(jl-1); yb2(jl)=yb2(jl-1); end ya1(1)= ((1. + rg(tc0) ) / 2.) * yg(tc0) + rg(tc0) * ya2(M1(loop)); ya2(1)= -1. * r1(loop) * ya1(M1(loop)) + ( 1. - r1(loop) ) * yb2(M2(loop)); yb1(1)= ( 1 + r1(loop) ) * ya1(M1(loop)) + r1(loop) * yb2(M2(loop)); yb2(1)= -1. * rl(tc0) * yb1(M2(loop)); y2tm(tc0)= (1 + rl(tc0)) * yb1(M2(loop)); end // t0 end //loop end //End of Two Tubes model endfunction // step2() //------------------------------------------------------------------- // Function for disply two tubes model's frequency-phase response when glottis is closed function y = two_tubes( xw, ni, rg0, rl0 ) yi = 0.5 * ( 1 + rg0 ) * ( 1 + r1(ni)) * ( 1 + rl0 ) * %e^( -1 * ( tu1(ni) + tu2(ni) ) .* xw * %i); yb = 1 + r1(ni) * rg0 * %e^( -2 * tu1(ni) .* xw * %i) + r1(ni) * rl0 * %e^( -2 * tu2(ni) .* xw * %i) + rl0 * rg0 * %e^( -2 * (tu1(ni) + tu2(ni)) .* xw * %i); y = yi ./ yb; endfunction function y = two_tubes2( xw, ni, rg0, rl0 ) yi = 0.5 * ( 1 + rg0 ) * ( 1 + r1(ni)) * ( 1 + rl0 ) * %e^( -1 * ( tu1(ni) + tu2(ni) ) .* xw * %i); yb = 1 + r1(ni) * rg0 * %e^( -2 * tu1(ni) * xw * %i) + r1(ni) * rl0 * %e^( -2 * tu2(ni) * xw * %i) + rl0 * rg0 * %e^( -2 * (tu1(ni) + tu2(ni)) * xw * %i); y = yi / yb; endfunction function y = two_tubes2n( xw, ni, rg0, rl0 ) yi = 0.5 * ( 1 + rg0 ) * ( 1 + r1n(ni)) * ( 1 + rl0 ) * %e^( -1 * ( tu1n(ni) + tu2n(ni) ) .* xw * %i); yb = 1 + r1n(ni) * rg0 * %e^( -2 * tu1n(ni) * xw * %i) + r1n(ni) * rl0 * %e^( -2 * tu2n(ni) * xw * %i) + rl0 * rg0 * %e^( -2 * (tu1n(ni) + tu2n(ni)) * xw * %i); y = yi / yb; endfunction function two_tubes_model_bode(ni) frq=[100:25:Max_Freq]'; [db,phi] =dbphi(two_tubes(2 * %pi * frq, ni,rg_closed ,rl_default)); clf(); bode0(frq,db,phi); endfunction // Function for disply three tubes model's frequency-phase response when glottis is closed function y = three_tubes( xw, ni, rg0, rl0 ) if T_QTY == 4 then // "T" type yi = 0.5 * ( 1. + rg0 ) * ( 1. + r21(ni)) * ( 1. + rl0 ) * %e^( -1. * ( tu1(ni) + tu2(ni)) * xw * %i); yb = 1.0; yb = yb + r12(ni) * rg0 * %e^( -2 * tu1(ni) * xw * %i) + r21(ni) * rl0 * %e^( -2 * tu2(ni) * xw * %i); yb = yb + rg0 * rl0 * ( 1. - r12(ni) + r21(ni) ) * %e^( -2 * (tu1(ni) + tu2(ni)) * xw * %i); yc = rl3 * ( 1. + r31(ni)); yc = yc * (1. + rg0 * %e^( -2 * tu1(ni) * xw * %i)); yc = yc * (1. - rl0 * %e^( -2 * tu2(ni) * xw * %i)); yc = yc * %e^( -1 * tu3(ni) * xw * %i); yd = 1. - rl3 * %e^( -2 * tu3(ni) * xw * %i); if yd == 0. then // There is room for reconsideration in this limit value. y=0.; // because yc/yd is infinity when yd = 0 else yc = yc / yd; y = yi / ( yb + yc); end else // serial type yi = 0.5 * ( 1. + rg0 ) * ( 1. + r1(ni)) * ( 1. + r2(ni)) * ( 1. + rl0 ) * %e^( -1. * ( tu1(ni) + tu2(ni) + tu3(ni)) * xw * %i); yb = 1 + rg0 * r1(ni) * %e^( -2 * tu1(ni) * xw * %i) + r1(ni) * r2(ni) * %e^( -2 * tu2(ni) * xw * %i) + r2(ni) * rl0 * %e^( -2 * tu3(ni) * xw * %i); yb = yb + rg0 * r2(ni) * %e^( -2 * (tu1(ni) + tu2(ni)) * xw * %i) + r1(ni) * rl0 * %e^( -2 * (tu2(ni) + tu3(ni)) * xw * %i); yb = yb + rg0 * r1(ni) * r2(ni) * rl0 * %e^( -2 * (tu1(ni) + tu3(ni)) * xw * %i); yb = yb + rg0 * rl0 * %e^( -2 * (tu1(ni) + tu2(ni) + tu3(ni)) * xw * %i); y = yi / yb; end endfunction function three_tubes_model_bode0(ni) frq3=[100:25:Max_Freq]; frq3s=size(frq3); resp3=zeros(1,frq3s(2)); for v=1:frq3s(2) resp3(v)= three_tubes(2 * %pi * frq3(v), ni,rg_closed ,rl_default); end [db,phi] =dbphi(resp3); clf(); bode(frq3,db,phi); 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 if T_QTY == 4 then //Three Tubes Model of T type ya1=sqrt(A1(ni) ); ya2=sqrt(A2(ni) ); ya3=sqrt(A3(ni) ); if ya2 > ya1 then yy0=ya2; else yy0=ya1; end yy0=10- (yy0 * 2); // Area scale is double than other Tube Models. yy1=ya1 + yy0; yy2=ya2 + yy0; yy3=yy0 - L3(ni); xx1=L1(ni) - (ya3 / 2); xx2=L1(ni) + (ya3 / 2); xp=[0 0 xx1 xx1 xx2 xx2 L1(ni)+L2(ni) L1(ni)+L2(ni) L1(ni) L1(ni) 0 ]'; yp=[yy1 yy0 yy0 yy3 yy3 yy0 yy0 yy2 yy2 yy1 yy1]'; xpolys( xp, yp, [5 5 5 5 5 5 5 5 5 5 5] ); xar=[-2 -0.5]; yar=[yy0+ ( ya1 / 2 ) yy0+ ( ya1 / 2 ) ]; if scilab_version_number <= 3 then xarrows(xar, yar, 1, 5); else xarrows(xar, yar, 20, 5); end yar=[yy0 + ( ya2 / 2 ) yy0+ ( ya2 / 2 ) ]; xar=[ L1(ni)+L2(ni)+0.5 L1(ni)+L2(ni)+2]; if scilab_version_number <= 3 then xarrows(xar, yar, 1, 5); else xarrows(xar, yar, 20, 5); end xstring(xx1, yy3 - 1, "CLOSED"); xstring(0,-9,"Attension: Area scale is double than others Tube Models"); // elseif T_QTY == 3 then //Three Tubes Model of serial type ya1=sqrt(A1(ni) ); ya2=sqrt(A2(ni) ); ya3=sqrt(A3(ni) ); xp=[0 0 L1(ni) L1(ni) L1(ni)+L2(ni) L1(ni)+L2(ni) L1(ni)+L2(ni)+L3(ni) L1(ni)+L2(ni)+L3(ni) L1(ni)+L2(ni) L1(ni)+L2(ni) L1(ni) L1(ni) 0 ]'; yp=[ya1 -ya1 -ya1 -ya2 -ya2 -ya3 -ya3 ya3 ya3 ya2 ya2 ya1 ya1 ]'; xpolys( xp, yp, [5 5 5 5 5 5 5 5 5 5 5 5] ); xar=[-2 -0.5]; yar=[0 0]; if scilab_version_number <= 3 then xarrows(xar, yar, 1, 5); else xarrows(xar, yar, 20, 5); end xar=[ L1(ni)+L2(ni)+L3(ni)+0.5 L1(ni)+L2(ni)+L3(ni)+2]; if scilab_version_number <= 3 then xarrows(xar, yar, 1, 5); else xarrows(xar, yar, 20, 5); end else //Two Tubes Model 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]; if scilab_version_number <= 3 then xarrows(xar, yar, 1, 5); else xarrows(xar, yar, 20, 5); end xar=[ L1(ni)+L2(ni)+0.5 L1(ni)+L2(ni)+2]; if scilab_version_number <= 3 then xarrows(xar, yar, 1, 5); else xarrows(xar, yar, 20, 5); end if T_QTY==5 then xar=[ L1(ni)+L2(ni)-0.5 L1(ni)+L2(ni)-2]; if scilab_version_number <= 3 then xarrows(xar, yar, 1, 3); else xarrows(xar, yar, 20, 3); end end end endfunction // // End of Step 2 Two Tubes Model for Vocal Tract Calculation // //========================================================================== // // Sub-Step 2-2 "rl" noise Calculation // //------------------------------------------------------ // global variables y_ran2=zeros(1,1); // dummy set y2tm_lpf= zeros(1,1); // dummy set al1sm= -1.; // dummy set bl0sm= 1.; // dummy set bl1sm= 0.; // dummy set // function [y2tm,y2tm_lpf,y_ran2,al1sm,bl0sm,bl1sm]=step2_2() if T_QTY == 5 then disp('+++enter sub-step 2-2 rl noise Calculation') ; // band pass filter of linear-phase FIR filter // 0 sc0 then sc0=int(loop / sizex); disp( sc0 * sizex ); // counter display end ydz = 0.; for loop2=1:nstep_fir ydz = ydz + wft(loop2) * y_ran1(1 + nstep_fir - loop2 + loop); end y_ran2(loop)=ydz; end // disp('+++enter sub-step 2-3 Re-Calculation of Two Tubes Model for Vocal Tract ') ; // // low pass filter's coefficients wcsm=2. * %pi * fsm; al1sm= -1. * %e^( -1. * wcsm * ts); bl0sm= wcsm * ts; bl1sm= 0.; // IIR type low pass filter LI1sm=2; LI2sm=2; ya1sm= zeros(1,LI1sm); yb1sm= zeros(1,LI2sm); // // yt is output of Two Tubes Model for Vocal Tract M1= zeros(1, D_QTY); M2= zeros(1, D_QTY); y2tm= zeros(1,length0+1); y2tm_lpf= zeros(1,length0+1); // lip-source(0/s=1/r_start=2/rl_noise=3) if (T_QTY == 2) | (T_QTY ==5) then idy0=7; elseif T_QTY == 3 then idy0=9; end for loop=1:D_QTY M1(loop)= int( tu1(loop) / ts ) + 1; M2(loop)= int( tu2(loop) / ts ) + 1; end MAX_M1=M1(1); MAX_M2=M2(1); for loop=1:D_QTY if M1(loop) > MAX_M1 then MAX_M1= M1(loop) end if M2(loop) > MAX_M2 then MAX_M2= M2(loop) end end ya1= zeros(1,MAX_M1); ya2= zeros(1,MAX_M1); yb1= zeros(1,MAX_M2); yb2= zeros(1,MAX_M2); y_ran3=zeros(1, (length0 + 1)); tc0=1; ya1(tc0)=0.; ya2(tc0)=0.; yb1(tc0)=0.; yb2(tc0)=0.; y2tm(tc0)=0.; y2tm_lpf(tc0)=0.; for loop=1:D_QTY //---lip-source(0/s=1/r_start=2/rl_noise=3) xcode0=lrpa(loop,idy0); xamp0 =lrpa(loop,(idy0-2)); //--- for t0=1:int(LL(loop)) tc0=tc0+1; for jl=MAX_M1:-1:2 ya1(jl)=ya1(jl-1); ya2(jl)=ya2(jl-1); end for jl=MAX_M2:-1:2 yb1(jl)=yb1(jl-1); yb2(jl)=yb2(jl-1); end ya1(1)= ((1. + rg(tc0) ) / 2.) * yg(tc0) + rg(tc0) * ya2(M1(loop)); ya2(1)= -1. * r1(loop) * ya1(M1(loop)) + ( 1. - r1(loop) ) * yb2(M2(loop)); yb1(1)= ( 1 + r1(loop) ) * ya1(M1(loop)) + r1(loop) * yb2(M2(loop)); // y2tm(tc0)= (1 + rl(tc0)) * yb1(M2(loop)); // smoothing y2tm // low pass filter ya1sm(LI1sm)=ya1sm(1); yb1sm(LI2sm)=yb1sm(1); ya1sm(1)=y2tm(tc0); yb1sm(1)= bl0sm * ya1sm(1) + bl1sm * ya1sm(LI1sm) - al1sm * yb1sm(LI2sm); y2tm_lpf(tc0)=yb1sm(1); if xcode0 == 3 then // ---lip-source(0/s=1/r_start=2/rl_noise=3) if y2tm_lpf(tc0) > (threshold_occur_plus * xamp0) then ynoise = y_ran2(tc0); elseif y2tm_lpf(tc0) < (threshold_occur_minus * xamp0) then ynoise = y_ran2(tc0); else ynoise = 0.; end else if y2tm_lpf(tc0) > threshold_occur_plus then ynoise = y_ran2(tc0); elseif y2tm_lpf(tc0) < threshold_occur_minus then ynoise = y_ran2(tc0); else ynoise = 0.; end end y_ran3(tc0)=ynoise; if xcode0 == 3 then // ---lip-source(0/s=1/r_start=2/rl_noise=3) yb2(1)= -1. * rl(tc0) * yb1(M2(loop)) + mix_amplitude * ynoise * xamp0; // mix with noise else yb2(1)= -1. * rl(tc0) * yb1(M2(loop)) + mix_amplitude * ynoise; // mix with noise end end // t0 end //loop // end // if T_QTY == 5 then endfunction // step2_2() // Function for disply two tubes model with "rl" noise frequency-phase response when glottis is closed // function y = two_tubes_rl_noise( xw, ni, rg0, rl0 ) //yi = 0.5 * ( 1. + rg0 ) * ( 1. + r1(ni)) * ( 1. + rl0 ) * %e^( -1. * ( tu1(ni) + tu2(ni)) * xw * %i); // this is for ug yi = r1(ni) * r1(ni) * rg0 * %e^( -2. * tu1(ni) * xw * %i) + r1(ni) * %e^( -2. * tu2(ni) * xw * %i) + (1. - r1(ni) * r1(ni)) * rg0 * %e^( -2 * (tu1(ni) + tu2(ni)) * xw * %i); yi = (1. + rl0 ) * yi; yb = 1 + r1(ni) * rg0 * %e^( -2 * tu1(ni) * xw * %i) + r1(ni) * rl0 * %e^( -2 * tu2(ni) * xw * %i) + rl0 * rg0 * %e^( -2 * (tu1(ni) + tu2(ni)) * xw * %i); y = yi / yb; endfunction function two_tubes_rl_noise_bode(ni) frq3=[100:25:Max_Freq]; frq3s=size(frq3); resp3=zeros(1,frq3s(2)); for v=1:frq3s(2) resp3(v)= two_tubes_rl_noise(2 * %pi * frq3(v), ni,rg_closed ,rl_default); end [db,phi] =dbphi(resp3); clf(); bode0(frq3,db,phi); endfunction // Function FFT y_ran2(noise with limited frequency band) and plot frequency-phase // input argument lng is fft points 1024 or others // points 1024, when nzyou=10 // points 2048, when nzyou=11 // points 4096, when nzyou=12 function ntz=fft_plot_ran2( lng ) nzyou= int(log2( lng )); ntz=2^nzyou; // check if data length is enough for fft ct0=0; for loop=2:D_QTY ct0=ct0+LL(loop); end if ct0 > ntz then // when actual data > ntz xin1=zeros(1,ntz); for loop=1:ntz xin1(loop) = y_ran2( loop); end win_hn=window('hn',ntz); // make hanning windows data xin2= xin1 .* win_hn; fftout=fft(xin2,-1); dfs= fs / ntz; p100= int(100. / dfs); p10000= int(Max_Freq / dfs); frq=[(dfs * p100):dfs:(dfs * p10000)]; frqs=size(frq); resp=zeros(1,frqs(2)); ct0=1; for loop=p100:p10000 resp(ct0)=fftout( loop + 1) / fftout(1); ct0=ct0+1; end [db,phi] =dbphi(resp); clf(); bode0(frq,db,phi); end // when actual data > ntz endfunction // Function for disply low pass filter's frequency-phase response function y = lpf1( xw ) yi= bl0sm + bl1sm * cos( xw * ts) - bl1sm * sin( xw * ts ) * %i; yb=1.0 + al1sm * cos( xw * ts ) - al1sm * sin( xw * ts ) * %i; y=yi ./yb; endfunction function y = lpf1n( xw ) yi= bl0smb + bl1smb * cos( xw * ts) - bl1smb * sin( xw * ts ) * %i; yb=1.0 + al1smb * cos( xw * ts ) - al1smb * sin( xw * ts ) * %i; y=yi ./yb; endfunction function lpf_bode() frq=[100:25:Max_Freq]'; [db,phi] =dbphi(lpf1(2 * %pi * frq )); clf(); bode0(frq,db,phi); endfunction //-------------------------------------------------------------------------- function [y2tm0]=after_step2() s0=size(y2tm); y2tm0=zeros(s0(1),s0(2)); for v=1:s0(2) y2tm0(v)=y2tm(v); end s1=size(LL); l0=0; for v=1:s1(2) l0=l0+LL(v); end if s0(2) xcode0 then epp0=eppx; ed_qty=loop-1; end // process if epp0 > spp0 then if xcode0 == 1 then // /sa/ length1= epp0 - spp0 + 1; y2tm0d = zeros(1,length1); for v=spp0:epp0 y2tm0d(v-spp0+1)=y2tm(v); end [out0]=do_iir_filtering3(y2tm0d,x_coef4,y_coef4); for v=spp0:epp0 y2tm0(v)=out0(v-spp0+1); end //.. elseif xcode0 == 2 then // addtion /r/ start mark length1= epp0 - spp0 + 1; y_ran1=rand(1, length1); [out0x]=do_iir_filtering3(y_ran1,x_coef9,y_coef9); [out0x2]=do_iir_filtering3(out0x,x_coef10,y_coef10); [out0]=do_iir_filtering3(out0x2,x_coef11,y_coef11); tc0=spp0; for v=sd_qty:ed_qty for t0=1:int(LL(v)) y2tm0(tc0)= y2tm(tc0) + r_ratio * amplitude(v) * out0(tc0-spp0+1); tc0=tc0+1; end end //.. end //common spp0=epp0+1; epp0=0; sd_qty=loop; end // process for only last frame, as same as process above if epp1 > spp1 then if xcode1 == 1 then // /sa/ length1= epp1 - spp1 + 1; y2tm0d = zeros(1,length1); for v=spp1:epp1 y2tm0d(v-spp0+1)=y2tm(v); end [out0]=do_iir_filtering3(y2tm0d,x_coef4,y_coef4); for v=spp1:epp1 y2tm0(v)=out0(v-spp0+1); end //.. elseif xcode0 == 2 then // addtion /r/ start mark length1= epp1 - spp1 + 1; y_ran1=rand(1, length1); [out0x]=do_iir_filtering3(y_ran1,x_coef9,y_coef9); [out0x2]=do_iir_filtering3(out0x,x_coef10,y_coef10); [out0]=do_iir_filtering3(out0x2,x_coef11,y_coef11); tc0=spp1; for v=spp1:epp1 y2tm0(v)= y2tm(v) + r_ratio * amplitude(loop) * out0(v-spp1+1); end //.. end end // stack up xcode0=lrpa(loop,idy0); eppx=eppx+int(LL(loop)); end // loop endfunction // after_step2() // //-------------------------------------------------------------------------- // global al1smb= -1.; // dummy set bl0smb= 1.; // dummy set bl1smb= 0.; // dummy set // function [y2tm,al1smb,bl0smb,bl1smb]=step2_n() // //+++++++++++++++++++++ //Two Tubes Model disp('+++enter step 2n Two Tubes Model for another two Tube'); M1= zeros(1, D_QTY); M2= zeros(1, D_QTY); y2tm_temp= zeros(1,length0+1); y2tm= zeros(1,length0+1); rl=rl_nose; for loop=1:D_QTY M1(loop)= int( tu1n(loop) / ts ) + 1; M2(loop)= int( tu2n(loop) / ts ) + 1; end MAX_M1=M1(1); MAX_M2=M2(1); for loop=1:D_QTY if M1(loop) > MAX_M1 then MAX_M1= M1(loop) end if M2(loop) > MAX_M2 then MAX_M2= M2(loop) end end ya1= zeros(1,MAX_M1); ya2= zeros(1,MAX_M1); yb1= zeros(1,MAX_M2); yb2= zeros(1,MAX_M2); tc0=1; ya1(tc0)=0.; ya2(tc0)=0.; yb1(tc0)=0.; yb2(tc0)=0.; y2tm_temp(tc0)=0.; // for loop=1:D_QTY for t0=1:int(LL(loop)) tc0=tc0+1; for jl=MAX_M1:-1:2 ya1(jl)=ya1(jl-1); ya2(jl)=ya2(jl-1); end for jl=MAX_M2:-1:2 yb1(jl)=yb1(jl-1); yb2(jl)=yb2(jl-1); end ya1(1)= ((1. + rg(tc0) ) / 2.) * ygn(tc0) + rg(tc0) * ya2(M1(loop)); ya2(1)= -1. * r1n(loop) * ya1(M1(loop)) + ( 1. - r1n(loop)) * yb2(M2(loop)); yb1(1)= ( 1 + r1n(loop) ) * ya1(M1(loop)) + r1n(loop) * yb2(M2(loop)); yb2(1)= -1. * rl * yb1(M2(loop)); y2tm_temp(tc0)= (1 + rl) * yb1(M2(loop)); end // t0 end //loop // next is low pass filtering // low pass filter's coefficients wcsm=2. * %pi * fc_nose_lpf; al1smb= -1. * %e^( -1. * wcsm * ts); bl0smb= wcsm * ts; bl1smb= 0.; // IIR type low pass filter LI1sm=2; LI2sm=2; ya1sm= zeros(1,LI1sm); yb1sm= zeros(1,LI2sm); ya1sm2= zeros(1,LI1sm); yb1sm2= zeros(1,LI2sm); // tc0=1; for loop=1:D_QTY for t0=1:int(LL(loop)) tc0=tc0+1; // low pass filter ya1sm(LI1sm)=ya1sm(1); yb1sm(LI2sm)=yb1sm(1); ya1sm(1)=y2tm_temp(tc0); yb1sm(1)= bl0smb * ya1sm(1) + bl1smb * ya1sm(LI1sm) - al1smb * yb1sm(LI2sm); if nose_lpf_times == 2 then ya1sm2(LI1sm)=ya1sm2(1); yb1sm2(LI2sm)=yb1sm2(1); ya1sm2(1)= yb1sm(1); //y2tm_temp(tc0); yb1sm2(1)= bl0smb * ya1sm2(1) + bl1smb * ya1sm2(LI1sm) - al1smb * yb1sm2(LI2sm); y2tm(tc0)=yb1sm2(1); else y2tm(tc0)=yb1sm(1); end end // t0 end //loop //+++++++++++++++++++++ endfunction // step2_n() //-------------------------------------------------------------------------- // global // variable_filter_on=1; // when this is 1, then, do variable filtering variable_filter_on_hcf=1; // when this is 1, then, do variable hcf filtering hcf_order_on=4; // set variable hcf iir filter order ripple_db_on=1; // set variable hcf iir Chebyshev filter pass band ripple [dB] xi_var=zeros(1,6); // dummy set yi_var=zeros(1,6); // dummy set resp_vari=zeros(1,1); // dummy set xi_var_hcf=zeros(1,(hcf_order_on+1)); // dummy set yi_var_hcf=zeros(1,(hcf_order_on+1)); // dummy set resp_vari_hcf=zeros(1,1); // dummy set wakeme=10000; // when double_bef, fc1= int(fc/wakeme) and fc2=fc-wakeme*fc1 function [xi_var0,yi_var0,resp_vari0,xi_var1,yi_var1,resp_vari1]=set_variable_coeff() // init xi_var0=zeros(length0+1,6); // change from 3 to 6 for double yi_var0=zeros(length0+1,6); for v=1:(length0+1) xi_var0(v,1)=1.; yi_var0(v,1)=1.; xi_var0(v,4)=1.; // for double yi_var0(v,4)=1.; end frq=[Min_Freq:25:Max_Freq]; sfrq=size(frq); resp_vari0=ones(D_QTY,sfrq(2)); if variable_filter_on==1 then disp(' ...calculating variable coefficients'); //// variable filter type(off=0/lpf=1/hpf=2/bpf=3/bef=4/double_bef=5),fc,q,gain if (T_QTY == 2) | (T_QTY ==5) then idy0=12; elseif T_QTY == 3 then idy0=14; end // when only 1 frame xloop=1; xcode0=lrpa(xloop,idy0); xspp0=1; xepp0=int(LL(xloop)); xoffset0=0; if D_QTY==1 then if xcode0 > 0 then for v=xspp0:xepp0 w0 = lrpa(xloop, idy0+1) * 2. * %pi; q0 = lrpa(xloop, idy0+2); g0= lrpa(xloop, idy0+3); if xcode0 == 1 then // lpf [x,y]=set_lpf_dig( fs, g0 , w0, q0, 2); elseif xcode0 == 2 then // hpf [x,y]=set_hpf_dig( fs, g0 , w0, q0, 2); elseif xcode0 == 3 then // bpf [x,y]=set_bpf_dig( fs, g0 , w0, q0); elseif xcode0 == 4 then // bef [x,y]=set_bef_dig2( fs, 1.0 , w0, q0, g0); elseif xcode0 == 5 then // double_bef w0a=int(lrpa(xloop, idy0+1) / wakeme); w0b=lrpa(xloop, idy0+1) - w0a * wakeme; w0a=w0a * 2. * %pi; w0b=w0b * 2. * %pi; [x,y]=set_bef_dig2( fs, 1.0 , w0a, q0, g0); [x2,y2]=set_bef_dig2( fs, 1.0 , w0b, q0, g0); end xi_var0(v,1)=x(1); xi_var0(v,2)=x(2); xi_var0(v,3)=x(3); yi_var0(v,1)=y(1); yi_var0(v,2)=y(2); yi_var0(v,3)=y(3); if v==(xspp0+xoffset0) then for v2=1:sfrq(2) resp_vari0(xloop,v2)=gz0(x,y, 2 * %pi * frq(v2) ) end end if xcode0 >= 5 then // for double xi_var0(v,4)=x2(1); xi_var0(v,5)=x2(2); xi_var0(v,6)=x2(3); yi_var0(v,4)=y2(1); yi_var0(v,5)=y2(2); yi_var0(v,6)=y2(3); if v==(xspp0+xoffset0) then for v2=1:sfrq(2) resp_vari0(xloop,v2)=resp_vari0(xloop,v2) * gz0(x2,y2, 2 * %pi * frq(v2) ) end end end end end // if xcode0> 0 then else // end of when only 1 frame // if D_QTY > = 2 then xx=zeros(1,6); yy=zeros(1,6); spp0=0; epp0=0; spp1=-1; epp1=-1; sframe=-1; eframe=-2; for loop=1:D_QTY xcode0=lrpa(loop,idy0); spp0=epp0+1; epp0=epp0+int(LL(loop)); offset0=spp0 + int((epp0-spp0)/2); if xcode0== 0 then // nothing done else // check portion start? start_flag=1; if loop== 1 sframe=loop; spp1=spp0; elseif lrpa(loop-1,idy0) == 0 then sframe=loop; spp1=spp0; elseif xcode0 <> lrpa(loop-1,idy0) then sframe=loop; spp1=spp0; else start_flag=0; end // check portion end ? end_flag=1; if loop== D_QTY eframe=loop; epp1=epp0; elseif lrpa(loop+1,idy0) == 0 then eframe=loop; epp1=epp0; elseif xcode0 <> lrpa(loop+1,idy0) then eframe=loop; epp1=epp0; else end_flag=0; end // when end, process all stack frames if end_flag == 1 then if eframe == sframe then xloop=sframe; xspp0=spp1; xepp0=epp1; xoffset0=0; if xcode0> 0 then for v=xspp0:xepp0 w0 = lrpa(xloop, idy0+1) * 2. * %pi; q0 = lrpa(xloop, idy0+2); g0 = lrpa(xloop, idy0+3); if xcode0 == 1 then // lpf [x,y]=set_lpf_dig( fs, g0 , w0, q0, 2); elseif xcode0 == 2 then // hpf [x,y]=set_hpf_dig( fs, g0 , w0, q0, 2); elseif xcode0 == 3 then // bpf [x,y]=set_bpf_dig( fs, g0 , w0, q0); elseif xcode0 == 4 then // bef [x,y]=set_bef_dig2( fs, 1.0 , w0, q0, g0); elseif xcode0 == 5 then // double_bef w0a=int(lrpa(xloop, idy0+1) / wakeme); w0b=lrpa(xloop, idy0+1) - w0a * wakeme; w0a=w0a * 2. * %pi; w0b=w0b * 2. * %pi; [x,y]=set_bef_dig2( fs, 1.0 , w0a, q0, g0); [x2,y2]=set_bef_dig2( fs, 1.0 , w0b, q0, g0); end xi_var0(v,1)=x(1); xi_var0(v,2)=x(2); xi_var0(v,3)=x(3); yi_var0(v,1)=y(1); yi_var0(v,2)=y(2); yi_var0(v,3)=y(3); if v==(xspp0+xoffset0) then for v2=1:sfrq(2) resp_vari0(xloop,v2)=gz0(x,y, 2 * %pi * frq(v2) ) end end if xcode0 >= 5 then // for double xi_var0(v,4)=x2(1); xi_var0(v,5)=x2(2); xi_var0(v,6)=x2(3); yi_var0(v,4)=y2(1); yi_var0(v,5)=y2(2); yi_var0(v,6)=y2(3); if v==(xspp0+xoffset0) then for v2=1:sfrq(2) resp_vari0(xloop,v2)=resp_vari0(xloop,v2) * gz0(x2,y2, 2 * %pi * frq(v2) ) end end end end end // if xcode0> 0 then elseif eframe == (sframe+1) then xloop=sframe; xspp0=spp1; xepp0=epp1; xoffset0=0; if xcode0> 0 then for v=xspp0:xepp0 w0 = ((v-xspp0) * (lrpa(xloop+1, idy0+1) - lrpa(xloop, idy0+1))/(xepp0-xspp0) + lrpa(xloop, idy0+1)) * 2. * %pi; q0 = ((v-xspp0) * (lrpa(xloop+1, idy0+2) - lrpa(xloop, idy0+2))/(xepp0-xspp0) + lrpa(xloop, idy0+2)); g0 = ((v-xspp0) * (lrpa(xloop+1, idy0+3) - lrpa(xloop, idy0+3))/(xepp0-xspp0) + lrpa(xloop, idy0+3)); if xcode0 == 1 then // lpf [x,y]=set_lpf_dig( fs, g0 , w0, q0, 2); elseif xcode0 == 2 then // hpf [x,y]=set_hpf_dig( fs, g0 , w0, q0, 2); elseif xcode0 == 3 then // bpf [x,y]=set_bpf_dig( fs, g0 , w0, q0); elseif xcode0 == 4 then // bef [x,y]=set_bef_dig2( fs, 1.0 , w0, q0, g0); elseif xcode0 == 5 then // double_bef w0a1=int(lrpa(xloop+1, idy0+1) / wakeme); w0b1=lrpa(xloop+1, idy0+1) - w0a1 * wakeme; w0a2=int(lrpa(xloop, idy0+1) / wakeme); w0b2=lrpa(xloop, idy0+1) - w0a2 * wakeme; w0a = ((v-xspp0) * (w0a1 - w0a2)/(xepp0-xspp0) + w0a2) * 2. * %pi; w0b = ((v-xspp0) * (w0b1 - w0b2)/(xepp0-xspp0) + w0b2) * 2. * %pi; [x,y]=set_bef_dig2( fs, 1.0 , w0a, q0, g0); [x2,y2]=set_bef_dig2( fs, 1.0 , w0b, q0, g0); end xi_var0(v,1)=x(1); xi_var0(v,2)=x(2); xi_var0(v,3)=x(3); yi_var0(v,1)=y(1); yi_var0(v,2)=y(2); yi_var0(v,3)=y(3); if xcode0 >= 5 then // double_bef xi_var0(v,4)=x2(1); xi_var0(v,5)=x2(2); xi_var0(v,6)=x2(3); yi_var0(v,4)=y2(1); yi_var0(v,5)=y2(2); yi_var0(v,6)=y2(3); end end spp2=spp1; for v=sframe:eframe xx(1)=xi_var0(spp2,1); xx(2)=xi_var0(spp2,2); xx(3)=xi_var0(spp2,3); yy(1)=yi_var0(spp2,1); yy(2)=yi_var0(spp2,2); yy(3)=yi_var0(spp2,3); for v2=1:sfrq(2) resp_vari0(v,v2)=gz0(xx,yy, 2 * %pi * frq(v2) ) end if xcode0 >= 5 then // double_bef xx(1)=xi_var0(spp2,4); xx(2)=xi_var0(spp2,5); xx(3)=xi_var0(spp2,6); yy(1)=yi_var0(spp2,4); yy(2)=yi_var0(spp2,5); yy(3)=yi_var0(spp2,6); for v2=1:sfrq(2) resp_vari0(v,v2)=resp_vari0(v,v2) * gz0(xx,yy, 2 * %pi * frq(v2) ) end end spp2=spp2 + int(LL(v)); end end // if xcode0> 0 then elseif eframe > (sframe+1) then xtpp0=0; xp0=zeros(1,(eframe-sframe+1)); for vf0=sframe:(eframe-1) if vf0==sframe then xspp0=spp1; xtpp0=spp1+int(LL(vf0)); xepp0=xtpp0+int(LL(vf0+1))/2; xp0(vf0-sframe+1)=spp1; elseif vf0==(eframe-1) then xp0(vf0-sframe+1)=xtpp0; xp0(vf0-sframe+2)=xtpp0+int(LL(vf0)); xspp0=xepp0+1; xepp0=epp1; else xspp0=xepp0+1; xp0(vf0-sframe+1)=xtpp0; xtpp0=xtpp0+int(LL(vf0)); xepp0=xtpp0+int(LL(vf0+1))/2; end ////if vf0 == sframe then //// disp ('vf0, LL(vf0), xspp0, xepp0, xp0(vf0)'); ////end ////disp( xp0(vf0), xepp0, xspp0, LL(vf0), vf0); xloop=vf0; if xcode0> 0 then for v=xspp0:xepp0 w0 = ((v-xspp0) * (lrpa(xloop+1, idy0+1) - lrpa(xloop, idy0+1))/(xepp0-xspp0) + lrpa(xloop, idy0+1)) * 2. * %pi; q0 = ((v-xspp0) * (lrpa(xloop+1, idy0+2) - lrpa(xloop, idy0+2))/(xepp0-xspp0) + lrpa(xloop, idy0+2)); g0 = ((v-xspp0) * (lrpa(xloop+1, idy0+3) - lrpa(xloop, idy0+3))/(xepp0-xspp0) + lrpa(xloop, idy0+3)); if xcode0 == 1 then // lpf [x,y]=set_lpf_dig( fs, g0 , w0, q0, 2); elseif xcode0 == 2 then // hpf [x,y]=set_hpf_dig( fs, g0 , w0, q0, 2); elseif xcode0 == 3 then // bpf [x,y]=set_bpf_dig( fs, g0 , w0, q0); elseif xcode0 == 4 then // bef [x,y]=set_bef_dig2( fs, 1.0 , w0, q0, g0); elseif xcode0 == 5 then // double_bef w0a1=int(lrpa(xloop+1, idy0+1) / wakeme); w0b1=lrpa(xloop+1, idy0+1) - w0a1 * wakeme; w0a2=int(lrpa(xloop, idy0+1) / wakeme); w0b2=lrpa(xloop, idy0+1) - w0a2 * wakeme; w0a = ((v-xspp0) * (w0a1 - w0a2)/(xepp0-xspp0) + w0a2) * 2. * %pi; w0b = ((v-xspp0) * (w0b1 - w0b2)/(xepp0-xspp0) + w0b2) * 2. * %pi; [x,y]=set_bef_dig2( fs, 1.0 , w0a, q0, g0); [x2,y2]=set_bef_dig2( fs, 1.0 , w0b, q0, g0); end xi_var0(v,1)=x(1); xi_var0(v,2)=x(2); xi_var0(v,3)=x(3); yi_var0(v,1)=y(1); yi_var0(v,2)=y(2); yi_var0(v,3)=y(3); if xcode0 >= 5 then // double_bef xi_var0(v,4)=x2(1); xi_var0(v,5)=x2(2); xi_var0(v,6)=x2(3); yi_var0(v,4)=y2(1); yi_var0(v,5)=y2(2); yi_var0(v,6)=y2(3); end end end // if xcode0> 0 then end //for vf0=sframe:(eframe-1) if xcode0> 0 then for v=sframe:eframe spp2=xp0(v-sframe+1); xx(1)=xi_var0(spp2,1); xx(2)=xi_var0(spp2,2); xx(3)=xi_var0(spp2,3); yy(1)=yi_var0(spp2,1); yy(2)=yi_var0(spp2,2); yy(3)=yi_var0(spp2,3); for v2=1:sfrq(2) resp_vari0(v,v2)=gz0(xx,yy, 2 * %pi * frq(v2) ) end if xcode0 >= 5 then // double_bef xx(1)=xi_var0(spp2,4); xx(2)=xi_var0(spp2,5); xx(3)=xi_var0(spp2,6); yy(1)=yi_var0(spp2,4); yy(2)=yi_var0(spp2,5); yy(3)=yi_var0(spp2,6); for v2=1:sfrq(2) resp_vari0(v,v2)=resp_vari0(v,v2) * gz0(xx,yy, 2 * %pi * frq(v2) ) end end end end // if xcode0> 0 then end // eframe > (sframe+1) then // reset value after done spp1=-1; epp1=-1; sframe=-1; eframe=-2; end //if end_flag == 1 then end // if xcode0 <> 0 then end // loop end //if D_QTY > = 2 then end //if variable_filter_on==1 then //................. // for hcf // init xi_var1=zeros(length0+1,(hcf_order_on+1)); yi_var1=zeros(length0+1,(hcf_order_on+1)); for v=1:(length0+1) xi_var1(v,1)=1.; yi_var1(v,1)=1.; end frq=[Min_Freq:25:Max_Freq]; sfrq=size(frq); resp_vari1=ones(D_QTY,sfrq(2)); if variable_filter_on_hcf==1 then disp(' ...calculating variable hcf coefficients. this work may take some long time.'); //// hcf(off=0/on=1),fc,gain if (T_QTY == 2) | (T_QTY ==5) then idy0=16; elseif T_QTY == 3 then idy0=18; end // when only 1 frame xloop=1; xcode0=lrpa(xloop,idy0); xspp0=1; xepp0=int(LL(xloop)); xoffset0=0; if D_QTY==1 then if xcode > 0 then for v=xspp0:xepp0 w0 = lrpa(xloop, idy0+1) * 2. * %pi; g0 = lrpa(xloop, idy0+2); if xcode0 == 1 then // lpf [x,y]=set_lpf_cheb_dig( fs, hcf_order_on, w0, ripple_db_on, g0); end for w=1:(hcf_order_on+1) xi_var1(v,w)=x(w); yi_var1(v,w)=y(w); end if v==(xspp0+xoffset0) then for v2=1:sfrq(2) resp_vari1(xloop,v2)=gzx(x,y, 2 * %pi * frq(v2) ) end end end end // if xcode0> 0 then else // end of when only 1 frame // if D_QTY > = 2 then xx=zeros(1,(hcf_order_on+1)); yy=zeros(1,(hcf_order_on+1)); spp0=0; epp0=0; spp1=-1; epp1=-1; sframe=-1; eframe=-2; for loop=1:D_QTY //disp(loop); xcode0=lrpa(loop,idy0); spp0=epp0+1; epp0=epp0+int(LL(loop)); offset0=spp0 + int((epp0-spp0)/2); if xcode0== 0 then // nothing done else // check portion start? start_flag=1; if loop== 1 sframe=loop; spp1=spp0; elseif lrpa(loop-1,idy0) == 0 then sframe=loop; spp1=spp0; elseif xcode0 <> lrpa(loop-1,idy0) then sframe=loop; spp1=spp0; else start_flag=0; end // check portion end ? end_flag=1; if loop== D_QTY eframe=loop; epp1=epp0; elseif lrpa(loop+1,idy0) == 0 then eframe=loop; epp1=epp0; elseif xcode0 <> lrpa(loop+1,idy0) then eframe=loop; epp1=epp0; else end_flag=0; end // when end, process all stack frames if end_flag == 1 then if eframe == sframe then xloop=sframe; xspp0=spp1; xepp0=epp1; xoffset0=0; if xcode0> 0 then wstr1=sprintf(' enter for v=xspp0:xepp0, %d -> %d', xspp0, xepp0); disp(wstr1); for v=xspp0:xepp0 w0 = lrpa(xloop, idy0+1) * 2. * %pi; g0 = lrpa(xloop, idy0+2); if xcode0 == 1 then // lpf [x,y]=set_lpf_cheb_dig( fs, hcf_order_on, w0, ripple_db_on, g0); end for w=1:(hcf_order_on+1) xi_var1(v,w)=x(w); yi_var1(v,w)=y(w); end if v==(xspp0+xoffset0) then for v2=1:sfrq(2) resp_vari1(xloop,v2)=gzx(x,y, 2 * %pi * frq(v2) ) end end end end // if xcode0> 0 then elseif eframe == (sframe+1) then xloop=sframe; xspp0=spp1; xepp0=epp1; xoffset0=0; if xcode0> 0 then wstr1=sprintf(' enter for v=xspp0:xepp0, %d -> %d', xspp0, xepp0); disp(wstr1); for v=xspp0:xepp0 w0 = ((v-xspp0) * (lrpa(xloop+1, idy0+1) - lrpa(xloop, idy0+1))/(xepp0-xspp0) + lrpa(xloop, idy0+1)) * 2. * %pi; g0 = ((v-xspp0) * (lrpa(xloop+1, idy0+2) - lrpa(xloop, idy0+2))/(xepp0-xspp0) + lrpa(xloop, idy0+2)); if xcode0 == 1 then // lpf [x,y]=set_lpf_cheb_dig( fs, hcf_order_on, w0, ripple_db_on, g0); end for w=1:(hcf_order_on+1) xi_var1(v,w)=x(w); yi_var1(v,w)=y(w); end end spp2=spp1; for v=sframe:eframe for w=1:(hcf_order_on+1) xx(w)=xi_var1(spp2,w); yy(w)=yi_var1(spp2,w); end for v2=1:sfrq(2) resp_vari1(v,v2)=gzx(xx,yy, 2 * %pi * frq(v2) ) end spp2=spp2 + int(LL(v)); end end // if xcode0> 0 then elseif eframe > (sframe+1) then xtpp0=0; xp0=zeros(1,(eframe-sframe+1)); for vf0=sframe:(eframe-1) if vf0==sframe then xspp0=spp1; xtpp0=spp1+int(LL(vf0)); xepp0=xtpp0+int(LL(vf0+1))/2; xp0(vf0-sframe+1)=spp1; elseif vf0==(eframe-1) then xp0(vf0-sframe+1)=xtpp0; xp0(vf0-sframe+2)=xtpp0+int(LL(vf0)); xspp0=xepp0+1; xepp0=epp1; else xspp0=xepp0+1; xp0(vf0-sframe+1)=xtpp0; xtpp0=xtpp0+int(LL(vf0)); xepp0=xtpp0+int(LL(vf0+1))/2; end ////if vf0 == sframe then //// disp ('vf0, LL(vf0), xspp0, xepp0, xp0(vf0)'); ////end ////disp( xp0(vf0), xepp0, xspp0, LL(vf0), vf0); xloop=vf0; if xcode0> 0 then wstr1=sprintf(' enter for v=xspp0:xepp0, %d -> %d', xspp0, xepp0); disp(wstr1); for v=xspp0:xepp0 w0 = ((v-xspp0) * (lrpa(xloop+1, idy0+1) - lrpa(xloop, idy0+1))/(xepp0-xspp0) + lrpa(xloop, idy0+1)) * 2. * %pi; g0 = ((v-xspp0) * (lrpa(xloop+1, idy0+2) - lrpa(xloop, idy0+2))/(xepp0-xspp0) + lrpa(xloop, idy0+2)); if xcode0 == 1 then // lpf [x,y]=set_lpf_cheb_dig( fs, hcf_order_on, w0, ripple_db_on, g0); end for w=1:(hcf_order_on+1) xi_var1(v,w)=x(w); yi_var1(v,w)=y(w); end end end // if xcode0> 0 then end //for vf0=sframe:(eframe-1) if xcode0> 0 then for v=sframe:eframe spp2=xp0(v-sframe+1); for w=1:(hcf_order_on+1) xx(w)=xi_var1(spp2,w); yy(w)=yi_var1(spp2,w); end for v2=1:sfrq(2) resp_vari1(v,v2)=gzx(xx,yy, 2 * %pi * frq(v2) ) end end end // if xcode0> 0 then end // eframe > (sframe+1) then // reset value after done spp1=-1; epp1=-1; sframe=-1; eframe=-2; end //if end_flag == 1 then end // if xcode0 <> 0 then end // loop end //if D_QTY > = 2 then end //if variable_filter_on_hcf==1 then endfunction //---------------------------------------------------- function [y2tm0]=step2_variable_filter() s0=size(y2tm); y2tm0=zeros(1,s0(2)); for v=1:s0(2) y2tm0(v)=y2tm(v); end if variable_filter_on == 1 then disp('+++enter step 2v variable filtering'); //// variable filter type(off=0/lpf=1/hpf=2/bpf=3/bef=4/double_bef=5),fc,q,gain if (T_QTY == 2) | (T_QTY ==5) then idy0=12; elseif T_QTY == 3 then idy0=14; end yx=zeros(1,6); yx2=zeros(1,6); idx=1; x0=0.; k0=0; for loop=1:D_QTY xcode0=lrpa(loop,idy0); for t0=1:int(LL(loop)) k0=k0+1; yx(3+idx)= y2tm(k0); x0= yi_var(k0, 0+idx) * yx(3+idx) + yi_var(k0, 1+idx) * yx(4+idx) + yi_var(k0,2+idx) * yx(5+idx); yx(0+idx) = x0 + yx(1+idx) * xi_var(k0,1+idx) + yx(2+idx) * xi_var(k0,2+idx); //* shift for next time */ yx(2+idx)=yx(1+idx); yx(1+idx)=yx(0+idx); yx(5+idx)=yx(4+idx); yx(4+idx)=yx(3+idx); y2tm0(k0)=yx(0+idx); ////if xcode0 >= 5 then // for double //// Set yi_var(4,*) and xi_var(4,*) =0 others =0 if xcode < 5 yx2(3+idx)= y2tm0(k0); x0= yi_var(k0, 3+idx) * yx2(3+idx) + yi_var(k0, 4+idx) * yx2(4+idx) + yi_var(k0,5+idx) * yx2(5+idx); yx2(0+idx) = x0 + yx2(1+idx) * xi_var(k0,4+idx) + yx2(2+idx) * xi_var(k0,5+idx); yx2(2+idx)=yx2(1+idx); yx2(1+idx)=yx2(0+idx); yx2(5+idx)=yx2(4+idx); yx2(4+idx)=yx2(3+idx); y2tm0(k0)=yx2(0+idx); ////else //// yx2(1)=0.; yx2(2)=0.; yx2(3)=0.; yx2(4)=0.; yx2(5)=0.; yx2(6)=0.; ////end // if xcode0 >= 5 then // end //for loop=1:D_QTY end //for t0=1:int(LL(loop)) // ++++++++++++++++++++++++++++ end //if variable_filter_on == 1 then if variable_filter_on_hcf == 1 then s0=size(y2tm); y2tmx=zeros(1,s0(2)); for v=1:s0(2) y2tmx(v)=y2tm0(v); end disp('+++enter step 2vb variable hcf filtering'); n0=hcf_order_on+1; yx=zeros(1,(n0*2)); idx=1; x0=0.; k0=0; for loop=1:D_QTY for t0=1:int(LL(loop)) k0=k0+1; yx(n0+1)= y2tmx(k0); x0=0.; for v=1:n0 x0 = x0 + yi_var_hcf(k0, v) * yx(n0+v); end yx(1)=x0; for v=2:n0 yx(1) = yx(1) + yx(v) * xi_var_hcf(k0,v); end //* shift for next time */ for v=0:(n0-2) yx(n0-v)=yx(n0-v-1); yx((n0*2)-v)=yx((n0*2)-v-1); end // y2tm0(k0)=yx(1); end //for loop=1:D_QTY end //for t0=1:int(LL(loop)) // ++++++++++++++++++++++++++++ end //if variable_filter_on_hcf == 1 then endfunction // step2_variable_filter() //-------------------------------------------------------------------------- // for debug to test filter as an input make_sin_freq=500; make_sin_amp=0.3; switch_on=0; function [y2tm0]=make_sin_y2tm() s0=size(y2tm); y2tm0=zeros(1,s0(2)); for v=1:s0(2) y2tm0(v)=y2tm(v); end if switch_on == 1 then disp(' * y2tm becomes sin wave for debug to test filter '); for v=1:s0(2) w= make_sin_freq * (2.0 * %pi / fs); y2tm0(v)= make_sin_amp * sin( w * v ); end end endfunction //--------------------------------------------------------------------------- // // // End of Sub-Step 2-2 "rl" noise Calculation // //========================================================================== // // Step 3 Output Sound Pressure for Two Tubes Model for Vocal Tract Calculation // //------------------------------------------------------ // global variables y3out= zeros(1,1); // dummy set y3outn= zeros(1,1); // dummy set ah1=1.; // dummy set bh0=1.; // dummy set bh1=0.; // dummy set ah1b=1.; // dummy set bh0b=1.; // dummy set bh1b=0.; // dummy set // function [y3out,y3outn,ah1,bh0,bh1,ah1b,bh0b,bh1b]=step3() disp('+++enter step 3 Output Sound Pressure for Two Tubes Model for Vocal Tract Calculation'); // At first, low pass filter's coefficients 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 ); alfa= -1. * cos( (wc + wnew) * ts /2. ) / cos( (wc - wnew) * ts / 2. ); // 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); // At first, low pass filter's coefficients for another two tubes wc=2. * %pi * fc_nose; 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 * fc_nose; //alfa=f_alfa( wc, wnew ); alfa= -1. * cos( (wc + wnew) * ts /2. ) / cos( (wc - wnew) * ts / 2. ); // high pass filter's coefficients bh0b= bl0 / (1.0 - al1 * alfa); bh1b= alfa * bl0 / (1.0 - al1 * alfa); ah1b= (alfa - al1) / (1.0 - al1 * alfa); // IIR type high pass filter y3out= zeros(1,length0+1); y3outn= zeros(1,length0+1); LI1=2; LI2=2; ya1= zeros(1,LI1); yb1= zeros(1,LI2); ya1b= zeros(1,LI1); yb1b= zeros(1,LI2); ya1b2= zeros(1,LI1); yb1b2= zeros(1,LI2); //loop for loop=1:length0 ya1(LI1)=ya1(1); yb1(LI2)=yb1(1); ya1(1)=y2tm(loop); yb1(1)= bh0 * ya1(1) + bh1 * ya1(LI1) - ah1 * yb1(LI2); y3out(loop)=yb1(1); //--- nose ya1b(LI1)=ya1b(1); yb1b(LI2)=yb1b(1); ya1b(1)=y2tmn(loop); yb1b(1)= bh0b * ya1b(1) + bh1b * ya1b(LI1) - ah1b * yb1b(LI2); ya1b2(LI1)=ya1b2(1); yb1b2(LI2)=yb1b2(1); ya1b2(1)=yb1b(1); yb1b2(1)= bh0b * ya1b2(1) + bh1b * ya1b2(LI1) - ah1b * yb1b2(LI2); if nose_hpf_times == 2 then y3outn(loop)=yb1b2(1); else y3outn(loop)=yb1b(1); end // last sum y3out(loop)=y3out(loop)+ magic_gain_nose * y3outn(loop); // bypass addtion based on Rn // // end //loop endfunction // step3() // // 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 function y = hpf1n( xw ) yi= bh0b + bh1b * cos( xw * ts) - bh1b * sin( xw * ts ) * %i; yb=1.0 + ah1b * cos( xw * ts ) - ah1b * sin( xw * ts ) * %i; y=yi ./yb; endfunction function hpf_bode() frq=[100:25:Max_Freq]'; [db,phi] =dbphi(hpf1(2 * %pi * frq )); clf(); bode0(frq,db,phi); endfunction // function save y3out (Sound Pressure) as one sound wave file function save_as_sound_file( filename) max01=0.; for loop=1:length0 if abs( y3out(loop) ) > max01 then max01 = abs( y3out(loop)); end end y3out2=zeros(1,length0+1); if max01 > 0. then for loop=1:length0 y3out2(loop)= y3out(loop) / max01; end end wavwrite(y3out2, fs, 16, filename ); endfunction //------------------------ function [use_flag]=check_ygn_use() s1=size(ygn); use_flag=0; for w=1:s1(2) if ygn(w) <> 0. then use_flag=1; break; end end endfunction // // // End of Step 3 Output Sound Pressure for Two Tubes Model for Vocal Tract // //================================================================================ // // Step 4 Plot Glottal Volume Velocity and Sound Pressure // function step4() disp('+++enter step 4 Plot generated sound wave'); plot_yg_yout(); endfunction // function plot_yg_yout() clf(); plot(yg,'b'); [use_flag]=check_ygn_use(); if use_flag == 1 then plot(ygn,'g'); end plot(y3out,'r'); wstr1='generated wave' + ' : ' + memo1; xtitle(wstr1,'Samples (Time)', 'Amplitude'); if use_flag == 1 then legend('Glottal Volume Velocity','another','Sound Pressure',2); else legend('Glottal Volume Velocity','Sound Pressure',2); end endfunction // function plot_ygn_y3outn() clf(); plot(ygn,'b'); plot(y3outn,'r'); wstr1='ygn y3outn from another two tube like nose'; xtitle(wstr1,'Samples (Time)', 'Amplitude'); legend('Glottal Volume Velocity','Sound Pressure',2); endfunction // function plot_yg_y2tm_y2tm_lpf_yo() clf(); plot(yg,'b'); plot(y2tm,'y'); plot(y2tm_lpf,'g'); plot(y3out,'r'); //plot(y2tm,'r'); //plot(y_ran3,'g'); // noise plot xtitle('generated wave','Samples (Time)', 'Amplitude'); legend('Glottal Volume Velocity','2nd Tube Edge Volume Velocity','Smoothed Volume Velocity','Sound Pressure',2); endfunction // Function FFT y3out(Sound Pressure) and plot frequency-phase // FFT starts from 2nd section // input argument lng is fft points 1024 or others // points 1024, when nzyou=10 // points 2048, when nzyou=11 // points 4096, when nzyou=12 function ntz=fft_plot( lng ) nzyou= int(log2( lng )); ntz=2^nzyou; // check if data length is enough for fft ct0=0; for loop=2:D_QTY ct0=ct0+LL(loop); end if ct0 > ntz then // when actual data > ntz xin1=zeros(1,ntz); ofs=LL(1); // FFT starts from 2nd section for loop=1:ntz xin1(loop) = y3out( loop + ofs); end win_hn=window('hn',ntz); // make hanning windows data xin2= xin1 .* win_hn; fftout=fft(xin2,-1); dfs= fs / ntz; p100= int(100. / dfs); p10000= int(Max_Freq / dfs); frq=[(dfs * p100):dfs:(dfs * p10000)]; frqs=size(frq); resp=zeros(1,frqs(2)); ct0=1; for loop=p100:p10000 resp(ct0)=fftout( loop + 1) / fftout(1); ct0=ct0+1; end [db,phi] =dbphi(resp); clf(); bode(frq,db,phi); end // when actual data > ntz endfunction //------------------------------------------------------------------------- // + others function library // //------------------------------------------------------ // global variables sp1=1; // start point for actual display of section 1 ep1=1; // end point for actual display of section 1 fft_point1=512; // default set db_fft1=zeros(1,1); // dummy set phi_fft1=zeros(1,1); // dummy set db_all=zeros(1,1); // dummy set phi_all=zeros(1,1); // dummy set //---------------------------------------------------------------- function [fft_point1]=set_fft_points() l1=list('points',3,['128','256','512','1024']); wrep=x_choices('Select FFT points',list(l1)); fft_point1=512; // default if wrep== 1 then fft_point1=128; elseif wrep== 2 then fft_point1=256; elseif wrep== 3 then fft_point1=512; elseif wrep== 4 then fft_point1=1024; end endfunction //------------------------------------------------- function [frq1,sfrq1,is1,ie1]=set_frq( spec1024 ) // spec1024, tube model response no syuhasuu bunkainou wo ageru tame fs1=fs; // to keep old source compatibility dfsw= fs1 / fft_point1; is1= ceil(Min_Freq / dfsw); ie1= int(Max_Freq / dfsw); if spec1024 == 1024 then dfsw2= fs1 / 1024.0; else dfsw2=dfsw; end frq1=[(dfsw * is1):dfsw2:(dfsw * ie1)]; frqs=size(frq1); sfrq1=frqs(2); endfunction //----------------------------------------------- function [db_fft1,phi_fft1,rt_code]=do_fft_wav(wdat1) // size1=length0; // to keep old source compatibility [size1]=rt_size(wdat1); rt_code=-1; if size1 >= ( sp1 + fft_point1 - 1 ) then win_hn=window('hn',fft_point1); // make hanning windows data wdatw=zeros(1, fft_point1); for v=1:fft_point1 wdatw(v)=wdat1(sp1 + v - 1); end wdatw2= wdatw .* win_hn; fftwout=fft( wdatw2, -1 ); // [frq1,sfrq1,is1,ie1]=set_frq(0); // respw=zeros(1,sfrq1); ct0=1; for loop=is1:ie1 respw(ct0)=fftwout(loop+1); ct0=ct0+1; end [db_fft1,phi_fft1] =dbphi(respw); rt_code=0; end // -if size1 <= ( sp1 + fft_point1 - 1 ) then endfunction //----------------------------------------------- function [wsp1, wep1]= set_sp1_ep1(wdat1) // wsize1=length0; // to keep old source compatibility [wsize1]=rt_size(wdat1); txt1=['start point(sec.1)';'end point(sec.1)']; wstr1=sprintf('%d',sp1); if ep1 > wsize1 then ep2 = wsize1; else ep2 = ep1; end if ep2 <= sp1 then ep2 = wsize1; end wstr2=sprintf('%d',ep2); sig1=x_mdialog('Input start point and end point for portion display.',txt1, [wstr1 ; wstr2 ]); if sig1==[] then arg1=evstr(wstr1); arg2=evstr(wstr2); else arg1=evstr(sig1(1)); arg2=evstr(sig1(2)); end wsp1=sp1; wep1=ep2; if arg1 >= 1 then if arg2 <= wsize1 then if arg1 < arg2 then wsp1=arg1; wep1=arg2; end end end endfunction //========================================================================== // filter design // //* w0= 2 * pai * f0 */ function [x,y]=set_bpf_ana( g0, w0, q0) y=zeros(1,3); x=zeros(1,3); idx=1; y(0+idx)=0.0; y(1+idx)= g0 * w0 / q0; y(2+idx)=0.0; x(0+idx)=1.0; x(1+idx)= w0 / q0; x(2+idx)= w0 * w0; endfunction //----- // some adjustment bef filter bottom level drift0=0.01; drift02=0.001; function [x,y]=set_bef_ana( g0, w0, q0) y=zeros(1,3); x=zeros(1,3); idx=1; y(0+idx)=g0; y(1+idx)=0.; y(2+idx)=g0 * w0 * w0 * (1. + drift0); x(0+idx)=1.0; x(1+idx)= w0 / q0; x(2+idx)= w0 * w0; endfunction function [x,y]=set_bef_ana2( g0, w0, q0) y=zeros(1,3); x=zeros(1,3); idx=1; y(0+idx)=g0; y(1+idx)=0.; y(2+idx)=g0 * w0 * w0 * (1. + drift02); x(0+idx)=1.0; x(1+idx)= w0 / q0; x(2+idx)= w0 * w0; endfunction //--- function [x,y]=set_hpf_ana(g0, w0, q0, xorder) y=zeros(1,3); x=zeros(1,3); idx=1; if xorder == 1 then // 1-st order hpf y(0+idx)= g0; y(1+idx)= 0.0; y(2+idx)= 0.0; x(0+idx)= 1.0; x(1+idx)= w0; x(2+idx)= 0.0 ; else // 2-nd order hpf y(0+idx)= g0; y(1+idx)= 0.0; y(2+idx)=0.0; x(0+idx)=1.0; x(1+idx)= w0 / q0; x(2+idx)= w0 * w0; end endfunction //--- function [x,y]=set_lpf_ana(g0, w0, q0, xorder) y=zeros(1,3); x=zeros(1,3); idx=1; if xorder == 1 then // 1-st order hpf y(0+idx)= 0.0; y(1+idx)= g0 * w0; y(2+idx)= 0.0; x(0+idx)= 1.0; x(1+idx)= w0; x(2+idx)= 0.0 ; else // 2-nd order hpf y(0+idx)= 0.0; y(1+idx)= 0.0; y(2+idx)= w0 * w0 * g0; x(0+idx)=1.0; x(1+idx)= w0 / q0; x(2+idx)= w0 * w0; end endfunction //-------------------------------------------------- function [wout]=souitizi_trans(w0,t0) //* w0 is target f0 * 2 * PAI */ //* t0 is sampling time */ wout = ((w0/(2.0 * %pi)) / (1.0 / t0)) * 2.0 * %pi; //* w0-> digital w he */ wout = tan( wout /2.0); wout = (wout * 2.0 / t0); endfunction //-------------------------------------------------- function [x,y]=set_bpf_dig( fs, g0, w0, q0) //* cal new _w0, and t0 */ t0x= 1.0/fs; [new_w0]= souitizi_trans( w0, t0x); //* cal analog x,y */ [x0,y0]=set_bpf_ana(g0,new_w0,q0); x=zeros(1,3) y=zeros(1,3); idx=1; //* cal X */ y(0+idx)= (2.0 * y0(1+idx) / t0x) + y0(2+idx); y(1+idx)= 2.0 * y0(2+idx); y(2+idx)= (-2.0 * y0(1+idx) / t0x) + y0(2+idx); //* cal Y */ x(0+idx)= (4.0/ (t0x * t0x) ) + (2.0 * x0(1+idx) / t0x) + x0(2+idx); x(1+idx)= (-8.0 / t0x /t0x) + (2.0 * x0(2+idx)); x(2+idx)= (4.0/ (t0x * t0x) ) + (-2.0 * x0(1+idx) / t0x) + x0(2+idx); //* hyokusetu gata sono 1 he henkei */ y(0+idx)=y(0+idx)/x(0+idx); y(1+idx)=y(1+idx)/x(0+idx); y(2+idx)=y(2+idx)/x(0+idx); x(1+idx)= -1.0 * x(1+idx)/x(0+idx); x(2+idx)= -1.0 * x(2+idx)/x(0+idx); x(0+idx)=x(0+idx)/x(0+idx); endfunction //-------------------------------------------------- function [x,y]=bpf_peak_0db(x0,y0,w0) //* cal new _w0, and t0 */ t0x= 1.0/fs; [new_w0]= souitizi_trans( w0, t0x); a0=gz0( x0 , y0 , new_w0 ); a1=abs(a0); s0=size(y); x=zeros(1,s0(2)); y=zeros(1,s0(2)); for v=1:s0(2) y(v)=y0(v) / a1; x(v)=x0(v); end endfunction //-------------------------------------------------- function [x,y]=set_hpf_dig( fs, g0, w0, q0, xorder) //* cal new _w0, and t0 */ t0x= 1.0/fs; [new_w0]= souitizi_trans( w0, t0x); //* cal analog x,y */ [x0,y0]=set_hpf_ana(g0,new_w0,q0, xorder); x=zeros(1,3) y=zeros(1,3); idx=1; if xorder == 1 then y(0+idx)=y0(0+idx); y(1+idx)= -1.0 * y0(0+idx); y(2+idx)= 0.0; x(0+idx)= (t0x * x0(1+idx) / 2.0) + 1.0; x(1+idx)= (t0x * x0(1+idx) / 2.0 ) - 1.0; x(2+idx)=0.0; else //* cal X */ y(0+idx)= (4.0 * y0(0+idx) / ( t0x * t0x)); y(1+idx)= ( -8.0 * y0(0+idx) / ( t0x * t0x)); y(2+idx)= (4.0 * y0(0+idx) / ( t0x * t0x)); //* cal Y */ x(0+idx)= (4.0/ (t0x * t0x) ) + (2.0 * x0(1+idx) / t0x) + x0(2+idx); x(1+idx)= (-8.0 / t0x /t0x) + (2.0 * x0(2+idx)); x(2+idx)= (4.0/ (t0x * t0x) ) + (-2.0 * x0(1+idx) / t0x) + x0(2+idx); end //* hyokusetu gata sono 1 he henkei */ y(0+idx)=y(0+idx)/x(0+idx); y(1+idx)=y(1+idx)/x(0+idx); y(2+idx)=y(2+idx)/x(0+idx); x(1+idx)= -1.0 * x(1+idx)/x(0+idx); x(2+idx)= -1.0 * x(2+idx)/x(0+idx); x(0+idx)=x(0+idx)/x(0+idx); endfunction //-------------------------------------------------- function [x,y]=set_lpf_dig( fs, g0, w0, q0, xorder) //* cal new _w0, and t0 */ t0x= 1.0/fs; [new_w0]= souitizi_trans( w0, t0x); //* cal analog x,y */ [x0,y0]=set_lpf_ana(g0,new_w0,q0, xorder); x=zeros(1,3) y=zeros(1,3); idx=1; if xorder == 1 then y(0+idx)=y0(1+idx); y(1+idx)=y0(1+idx); y(2+idx)= 0.0; x(0+idx)= (2.0 / t0x) + x0(1+idx); x(1+idx)= x0(1+idx) - (2.0 / t0x); x(2+idx)=0.0; else //* cal X */ y(0+idx)= (2.0 * y0(1+idx) / t0x) + y0(2+idx); y(1+idx)= 2.0 * y0(2+idx); y(2+idx)= (-2.0 * y0(1+idx) / t0x) + y0(2+idx); //* cal Y */ x(0+idx)= (4.0/ (t0x * t0x) ) + (2.0 * x0(1+idx) / t0x) + x0(2+idx); x(1+idx)= (-8.0 / t0x /t0x) + (2.0 * x0(2+idx)); x(2+idx)= (4.0/ (t0x * t0x) ) + (-2.0 * x0(1+idx) / t0x) + x0(2+idx); end //* hyokusetu gata sono 1 he henkei */ y(0+idx)=y(0+idx)/x(0+idx); y(1+idx)=y(1+idx)/x(0+idx); y(2+idx)=y(2+idx)/x(0+idx); x(1+idx)= -1.0 * x(1+idx)/x(0+idx); x(2+idx)= -1.0 * x(2+idx)/x(0+idx); x(0+idx)=x(0+idx)/x(0+idx); endfunction //-------------------------------------------------------------- function [x,y]=set_bef_dig( fs, g0, w0, q0) //* cal new _w0, and t0 */ t0x= 1.0/fs; [new_w0]= souitizi_trans( w0, t0x); //* cal analog x,y */ [x0,y0]=set_bef_ana(g0,new_w0,q0); x=zeros(1,3) y=zeros(1,3); idx=1; //* cal X */ y(0+idx)= (2.0 * y0(1+idx) / t0x) + y0(2+idx) + (4.0 * y0(0+idx) / ( t0x * t0x)); y(1+idx)= (2.0 * y0(2+idx)) + ( -8.0 * y0(0+idx) / ( t0x * t0x)); y(2+idx)= (-2.0 * y0(1+idx) / t0x) + y0(2+idx) + (4.0 * y0(0+idx) / ( t0x * t0x)); //* cal Y */ x(0+idx)= (4.0/ (t0x * t0x) ) + (2.0 * x0(1+idx) / t0x) + x0(2+idx); x(1+idx)= (-8.0 / t0x /t0x) + (2.0 * x0(2+idx)); x(2+idx)= (4.0/ (t0x * t0x) ) + (-2.0 * x0(1+idx) / t0x) + x0(2+idx); //* hyokusetu gata sono 1 he henkei */ y(0+idx)=y(0+idx)/x(0+idx); y(1+idx)=y(1+idx)/x(0+idx); y(2+idx)=y(2+idx)/x(0+idx); x(1+idx)= -1.0 * x(1+idx)/x(0+idx); x(2+idx)= -1.0 * x(2+idx)/x(0+idx); x(0+idx)=x(0+idx)/x(0+idx); endfunction //---------------------------------------------------- function [x,y]=set_bef_dig2( fs, g0, w0, q0, mix0db) // mix0= 10^( mix0db / 20. ); // out= (1 - mix0) * bef + mix0 * original //* cal new _w0, and t0 */ t0x= 1.0/fs; [new_w0]= souitizi_trans( w0, t0x); //* cal analog x,y */ [x0,y0]=set_bef_ana2(g0,new_w0,q0); x=zeros(1,3) y=zeros(1,3); idx=1; //* cal X */ y(0+idx)= (2.0 * y0(1+idx) / t0x) + y0(2+idx) + (4.0 * y0(0+idx) / ( t0x * t0x)); y(1+idx)= (2.0 * y0(2+idx)) + ( -8.0 * y0(0+idx) / ( t0x * t0x)); y(2+idx)= (-2.0 * y0(1+idx) / t0x) + y0(2+idx) + (4.0 * y0(0+idx) / ( t0x * t0x)); //* cal Y */ x(0+idx)= (4.0/ (t0x * t0x) ) + (2.0 * x0(1+idx) / t0x) + x0(2+idx); x(1+idx)= (-8.0 / t0x /t0x) + (2.0 * x0(2+idx)); x(2+idx)= (4.0/ (t0x * t0x) ) + (-2.0 * x0(1+idx) / t0x) + x0(2+idx); // y(0+idx)=mix0 * x(0+idx) + (1.0 - mix0) * y(0+idx); y(1+idx)=mix0 * x(1+idx) + (1.0 - mix0) * y(1+idx); y(2+idx)=mix0 * x(2+idx) + (1.0 - mix0) * y(2+idx); //* hyokusetu gata sono 1 he henkei */ y(0+idx)=y(0+idx)/x(0+idx); y(1+idx)=y(1+idx)/x(0+idx); y(2+idx)=y(2+idx)/x(0+idx); x(1+idx)= -1.0 * x(1+idx)/x(0+idx); x(2+idx)= -1.0 * x(2+idx)/x(0+idx); x(0+idx)=x(0+idx)/x(0+idx); endfunction //=================================================================== // scilab global variables //=================================================================== // global stack_infor1 stack_infor_size=30; stack_infor_data_xstart=10; stack_infor_data_ystart=20; stack_infor1=zeros(1,stack_infor_size); check_warning1=0; //---------------------------------------------------- function [x,y]=set_lpf_cheb_dig( fs, n_order0, w0, ripple_db, mix0db) global stack_infor1 if stack_infor1(1) == 0 then cal_flag=1; // this is call 1st time else cal_flag=0; if (( stack_infor1(2) <> fs) | ( stack_infor1(3) <> n_order0) | (stack_infor1(4) <> w0 ) | (stack_infor1(5) <> ripple_db)) then cal_flag =1; // constant has been changed, re-cal end end // Chebyshev type 1 filter design mix0= 10^( mix0db / 20. ); // out= (1 -mix0) * lpf + mix0 * original // this is not good. stop this work. // because phase rotate around fc and make dip due to cancel, if lpf and original mix is same level // mix0= 0.; // force to this work off if check_warning1 == 1 then if mix0db > -20 then disp('warning: mix0db value is not good, phase rotate will disturb frequency response. In set_lpf_cheb_dig'); end end // out= (mix0 -1) * bef + mix0 * original //* cal new _w0, and t0 */ //t0x= 1.0/fs; //[new_w0]= souitizi_trans( w0, t0x); new_w0=w0; // step 1, call iir // cheb1 lpf // order nth0= n_order0; //"""""""""""""""""""""""""""""""""""" if cal_flag == 1 then frqch0=zeros(1,2); frqch0(1)= ( new_w0 / ( 2. * %pi )) / fs; frqch0(2)=0.; // ripple_db= -1.; deltac0=zeros(1,2); deltac0(1)= 10^( ripple_db / 20. ) - 1.0; deltac0(2)=0.; [wH0A]=iir(nth0,'lp','cheb1', frqch0, deltac0); // for debug // disp( wH0A ); // [hzm,fr]=frmag(wH0A,256); // wb0=xget('window'); // stack old window // xset('window',2); // clf(); // plot2d(fr',hzm'); // xset('window',wb0); // end of for debug [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 [x]=get_coeff( H0bunbo, wiiro2 ); [y]=get_coeff( H0bunsi, wiiro2 ); // stack data for future use stack_infor1(2)=fs; stack_infor1(3)=n_order0; stack_infor1(4)=w0; stack_infor1(5)=ripple_db; for v=1:(n_order0+1) stack_infor1(stack_infor_data_xstart+v-1)=x(v); stack_infor1(stack_infor_data_ystart+v-1)=y(v); end stack_infor1(1)=stack_infor1(1)+1; else // when cal_flag == 0 x=zeros(1,(n_order0+1)); y=zeros(1,(n_order0+1)); for v=1:(n_order0+1) x(v)=stack_infor1(stack_infor_data_xstart+v-1); y(v)=stack_infor1(stack_infor_data_ystart+v-1); end stack_infor1(1)=stack_infor1(1)+1; end //if cal_flag == 1 then //"""""""""""""""""""""""""""""""""""""""""""" for v=1:(nth0+1) y(v)=mix0 * x(v) + (1.0 - mix0) * y(v); end x00=x(1); for v=1:(nth0+1) if v==1 then x(v)=x(v)/x00; else x(v)= -1.0 * x(v)/x00; end y(v)=y(v)/x00; end endfunction //-------------------------------------------------- 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 [out0]=do_iir_filtering3(sp,xi,yi) // sp is wave-input yx=zeros(1,6); idx=1; x0=0.; s0=size(sp); l0=s0(2); out0=zeros(1,l0); for k0=1:l0 yx(3+idx)= sp(k0); x0= yi(0+idx) * yx(3+idx) + yi(1+idx) * yx(4+idx) + yi(2+idx) * yx(5+idx); yx(0+idx) = x0 + yx(1+idx) * xi(1+idx) + yx(2+idx) * xi(2+idx); //* shift for next time */ yx(2+idx)=yx(1+idx); yx(1+idx)=yx(0+idx); yx(5+idx)=yx(4+idx); yx(4+idx)=yx(3+idx); out0(k0)=yx(0+idx); end endfunction //********************************************************************/ //* G(z)= y[0]z0 + y[1]z-1 + y[2]z-2/ 1.0 - (x[1]z-1 + x[2]z-2) */ //* z0 = tannni kakeru dake x[0]=1.0 function y = gz0( x ,y ,xw ) ibx=1; yi= y(0+ibx) + y(1+ibx) * cos( xw * ts) - y(1+ibx) * sin( xw * ts ) * %i + y(2+ibx) * cos( xw * ts * 2.0) - y(2+ibx) * sin( xw * ts * 2.0) * %i; yb=1.0 - ( x(1+ibx) * cos( xw * ts ) - x(1+ibx) * sin( xw * ts ) * %i + x(2+ibx) * cos( xw * ts * 2.0) - x(2+ibx) * sin( xw * ts * 2.0) * %i) ; y=yi ./yb; endfunction //--- function [x1,y1]=set_h_iir_et_plot_bode(disp0) frq=[100:25:Max_Freq]'; if disp0 >= 1 then wb0=xget('window'); // stack old window xset('window',disp0); // create new windows clf(); //(1) //subplot(211); else disp('+set_h_iir_et_plot_bode'); end idx=1; if co_filter(idx,2) == 0 then // lpf w0 = co_filter(idx,3) * 2. * %pi; q0 = co_filter(idx,4); [x,y]=set_lpf_dig( fs, 1.0 , w0, q0, 2); elseif co_filter(idx,2) == 2 then // bpf w0 = co_filter(idx,3) * 2. * %pi; q0 = co_filter(idx,4); [x,y]=set_bpf_dig( fs, 1.0 , w0, q0); elseif co_filter(idx,2) == 1 then // hpf w0 = co_filter(idx,3) * 2. * %pi; q0 = co_filter(idx,4); [x,y]=set_hpf_dig( fs, 1.0 , w0, q0, 2); end if disp0 >= 1 then [db,phi] =dbphi( gz0(x,y, 2 * %pi * frq )); bode0(frq,db,phi); wstr1='noise vocal source filter for /h/'; xtitle( wstr1 ); end s0=size(y); x1=zeros(1,s0(2)); y1=zeros(1,s0(2)); for v=1:s0(2) y1(v)=y(v); x1(v)=x(v); end //(2) //subplot(212); //idx=2; //if co_filter(idx,2) == 0 then // w0 = co_filter(idx,3) * 2. * %pi; // q0 = co_filter(idx,4); // [x,y]=set_lpf_dig( fs, 1.0 , w0, q0, 2); //elseif co_filter(idx,2) == 2 then // w0 = co_filter(idx,3) * 2. * %pi; // q0 = co_filter(idx,4); // [x,y]=set_bpf_dig( fs, 1.0 , w0, q0); //elseif co_filter(idx,2) == 1 then // hpf // w0 = co_filter(idx,3) * 2. * %pi; // q0 = co_filter(idx,4); // [x,y]=set_hpf_dig( fs, 1.0 , w0, q0, 2); //end //[db,phi] =dbphi( gz0(x,y, 2 * %pi * frq )); //bode0(frq,db,phi); //wstr1='noise source filter for /s/'; //xtitle( wstr1 ); //s0=size(y); //x2=zeros(1,s0(2)); //y2=zeros(1,s0(2)); //for v=1:s0(2) // y2(v)=y(v); // x2(v)=x(v); //end if disp0 >= 1 then xset('window',wb0); // push old windows end // if disp0 >= 1 then endfunction //--------------------------------------------------------- function [x1,y1,x2,y2,x3,y3]=set_s_iir_et_plot_bode(disp0) frq=[100:25:Max_Freq]'; if disp0 >= 1 then wb0=xget('window'); // stack old window xset('window',disp0); // create new windows clf(); //(1) subplot(311); else disp('+set_s_iir_et_plot_bode'); end idx=2; if co_filter(idx,2) == 0 then // lpf w0 = co_filter(idx,3) * 2. * %pi; q0 = co_filter(idx,4); [x,y]=set_lpf_dig( fs, 1.0 , w0, q0, 2); elseif co_filter(idx,2) == 2 then // bpf w0 = co_filter(idx,3) * 2. * %pi; q0 = co_filter(idx,4); [x,y]=set_bpf_dig( fs, 1.0 , w0, q0); elseif co_filter(idx,2) == 1 then // hpf w0 = co_filter(idx,3) * 2. * %pi; q0 = co_filter(idx,4); [x,y]=set_hpf_dig( fs, 1.0 , w0, q0, 2); end if disp0 >= 1 then [db,phi] =dbphi( gz0(x,y, 2 * %pi * frq )); bode0(frq,db,phi); wstr1='noise source filter for /s/'; xtitle( wstr1 ); end s0=size(y); x1=zeros(1,s0(2)); y1=zeros(1,s0(2)); for v=1:s0(2) y1(v)=y(v); x1(v)=x(v); end if disp0 >= 1 then //(2) subplot(312); end idx=3; if co_filter(idx,2) == 0 then w0 = co_filter(idx,3) * 2. * %pi; q0 = co_filter(idx,4); [x,y]=set_lpf_dig( fs, 1.0 , w0, q0, 2); elseif co_filter(idx,2) == 2 then w0 = co_filter(idx,3) * 2. * %pi; q0 = co_filter(idx,4); [x,y]=set_bpf_dig( fs, 1.0 , w0, q0); elseif co_filter(idx,2) == 1 then // hpf w0 = co_filter(idx,3) * 2. * %pi; q0 = co_filter(idx,4); [x,y]=set_hpf_dig( fs, 1.0 , w0, q0, 2); end if disp0 >= 1 then [db,phi] =dbphi( gz0(x,y, 2 * %pi * frq )); bode0(frq,db,phi); wstr1='noise source filter for /s/'; xtitle( wstr1 ); end s0=size(y); x2=zeros(1,s0(2)); y2=zeros(1,s0(2)); for v=1:s0(2) y2(v)=y(v); x2(v)=x(v); end if disp0 >= 1 then //(3) subplot(313); end idx=4; if co_filter(idx,2) == 0 then w0 = co_filter(idx,3) * 2. * %pi; q0 = co_filter(idx,4); [x,y]=set_lpf_dig( fs, 1.0 , w0, q0, 2); elseif co_filter(idx,2) == 2 then w0 = co_filter(idx,3) * 2. * %pi; q0 = co_filter(idx,4); [x,y]=set_bpf_dig( fs, 1.0 , w0, q0); elseif co_filter(idx,2) == 1 then // hpf w0 = co_filter(idx,3) * 2. * %pi; q0 = co_filter(idx,4); [x,y]=set_hpf_dig( fs, 1.0 , w0, q0, 2); end if disp0 >= 1 then [db,phi] =dbphi( gz0(x,y, 2 * %pi * frq )); bode0(frq,db,phi); wstr1='noise lip-around source filter for /s/'; xtitle( wstr1 ); end s0=size(y); x3=zeros(1,s0(2)); y3=zeros(1,s0(2)); for v=1:s0(2) y3(v)=y(v); x3(v)=x(v); end if disp0 >= 1 then xset('window',wb0); // push old windows end // if disp0 >= 1 then endfunction //--------------------------------------------------------- function [x1,y1,x2,y2,x3,y3]=set_r_iir_et_plot_bode(disp0) frq=[100:25:Max_Freq]'; if disp0 >= 1 then wb0=xget('window'); // stack old window xset('window',disp0); // create new windows clf(); //(1) subplot(311); else disp('+set_r_start_mark_iir et_plot_bode'); end idx=7; if co_filter(idx,2) == 0 then // lpf w0 = co_filter(idx,3) * 2. * %pi; q0 = co_filter(idx,4); [x,y]=set_lpf_dig( fs, 1.0 , w0, q0, 2); elseif co_filter(idx,2) == 2 then // bpf w0 = co_filter(idx,3) * 2. * %pi; q0 = co_filter(idx,4); [x,y]=set_bpf_dig( fs, 1.0 , w0, q0); elseif co_filter(idx,2) == 1 then // hpf w0 = co_filter(idx,3) * 2. * %pi; q0 = co_filter(idx,4); [x,y]=set_hpf_dig( fs, 1.0 , w0, q0, 2); end if disp0 >= 1 then [db,phi] =dbphi( gz0(x,y, 2 * %pi * frq )); bode0(frq,db,phi); wstr1='noise source filter for /r/ start mark'; xtitle( wstr1 ); end s0=size(y); x1=zeros(1,s0(2)); y1=zeros(1,s0(2)); for v=1:s0(2) y1(v)=y(v); x1(v)=x(v); end if disp0 >= 1 then //(2) subplot(312); end idx=8; if co_filter(idx,2) == 0 then w0 = co_filter(idx,3) * 2. * %pi; q0 = co_filter(idx,4); [x,y]=set_lpf_dig( fs, 1.0 , w0, q0, 2); elseif co_filter(idx,2) == 2 then w0 = co_filter(idx,3) * 2. * %pi; q0 = co_filter(idx,4); [x,y]=set_bpf_dig( fs, 1.0 , w0, q0); elseif co_filter(idx,2) == 1 then // hpf w0 = co_filter(idx,3) * 2. * %pi; q0 = co_filter(idx,4); [x,y]=set_hpf_dig( fs, 1.0 , w0, q0, 2); end if disp0 >= 1 then [db,phi] =dbphi( gz0(x,y, 2 * %pi * frq )); bode0(frq,db,phi); wstr1='noise source filter for /r/ start mark'; xtitle( wstr1 ); end s0=size(y); x2=zeros(1,s0(2)); y2=zeros(1,s0(2)); for v=1:s0(2) y2(v)=y(v); x2(v)=x(v); end if disp0 >= 1 then //(3) subplot(313); end idx=9; if co_filter(idx,2) == 0 then w0 = co_filter(idx,3) * 2. * %pi; q0 = co_filter(idx,4); [x,y]=set_lpf_dig( fs, 1.0 , w0, q0, 2); elseif co_filter(idx,2) == 2 then w0 = co_filter(idx,3) * 2. * %pi; q0 = co_filter(idx,4); [x,y]=set_bpf_dig( fs, 1.0 , w0, q0); elseif co_filter(idx,2) == 1 then // hpf w0 = co_filter(idx,3) * 2. * %pi; q0 = co_filter(idx,4); [x,y]=set_hpf_dig( fs, 1.0 , w0, q0, 2); end if disp0 >= 1 then [db,phi] =dbphi( gz0(x,y, 2 * %pi * frq )); bode0(frq,db,phi); wstr1='noise lip-around source filter for /r/ start mark'; xtitle( wstr1 ); end s0=size(y); x3=zeros(1,s0(2)); y3=zeros(1,s0(2)); for v=1:s0(2) y3(v)=y(v); x3(v)=x(v); end if disp0 >= 1 then xset('window',wb0); // push old windows end // if disp0 >= 1 then endfunction //-------------------------------------------------- function [x1,y1,x2,y2]=set_k_iir_et_plot_bode(disp0) frq=[100:25:Max_Freq]'; if disp0 >= 1 then wb0=xget('window'); // stack old window xset('window',disp0); // create new windows clf(); //(1) subplot(211); else disp('+set_k_iir_et_plot_bode'); end idx=5; if co_filter(idx,2) == 0 then // lpf w0 = co_filter(idx,3) * 2. * %pi; q0 = co_filter(idx,4); [x,y]=set_lpf_dig( fs, 1.0 , w0, q0, 2); elseif co_filter(idx,2) == 2 then // bpf w0 = co_filter(idx,3) * 2. * %pi; q0 = co_filter(idx,4); [x,y]=set_bpf_dig( fs, 1.0 , w0, q0); elseif co_filter(idx,2) == 1 then // hpf w0 = co_filter(idx,3) * 2. * %pi; q0 = co_filter(idx,4); [x,y]=set_hpf_dig( fs, 1.0 , w0, q0, 2); end if disp0 >= 1 then [db,phi] =dbphi( gz0(x,y, 2 * %pi * frq )); bode0(frq,db,phi); wstr1='noise source filter for /k/ et /t/'; xtitle( wstr1 ); end s0=size(y); x1=zeros(1,s0(2)); y1=zeros(1,s0(2)); for v=1:s0(2) y1(v)=y(v); x1(v)=x(v); end if disp0 >= 1 then //(2) subplot(212); end idx=6; if co_filter(idx,2) == 0 then w0 = co_filter(idx,3) * 2. * %pi; q0 = co_filter(idx,4); [x,y]=set_lpf_dig( fs, 1.0 , w0, q0, 2); elseif co_filter(idx,2) == 2 then w0 = co_filter(idx,3) * 2. * %pi; q0 = co_filter(idx,4); [x,y]=set_bpf_dig( fs, 1.0 , w0, q0); elseif co_filter(idx,2) == 1 then // hpf w0 = co_filter(idx,3) * 2. * %pi; q0 = co_filter(idx,4); [x,y]=set_hpf_dig( fs, 1.0 , w0, q0, 2); end if disp0 >= 1 then [db,phi] =dbphi( gz0(x,y, 2 * %pi * frq )); bode0(frq,db,phi); wstr1='noise beat source filter for /k/'; xtitle( wstr1 ); end s0=size(y); x2=zeros(1,s0(2)); y2=zeros(1,s0(2)); for v=1:s0(2) y2(v)=y(v); x2(v)=x(v); end if disp0 >= 1 then xset('window',wb0); // push old windows end // if disp0 >= 1 then endfunction //-------------------------------------------------- function [x1,y1,x2,y2]=set_N_iir_et_plot_bode(disp0) frq=[100:25:Max_Freq]'; if disp0 >= 1 then wb0=xget('window'); // stack old window xset('window',disp0); // create new windows clf(); //(1) subplot(211); else disp('+set_N_iir_et_plot_bode'); end idx=1; w0 = n_supend_para0(idx) * 2. * %pi; q0 = n_supend_para0(idx+1); [x,y]=set_bef_dig( fs, 1.0 , w0, q0); if disp0 >= 1 then [db,phi] =dbphi( gz0(x,y, 2 * %pi * frq )); bode0(frq,db,phi); wstr1='(low) bef to absorb'; xtitle( wstr1 ); end s0=size(y); x1=zeros(1,s0(2)); y1=zeros(1,s0(2)); for v=1:s0(2) y1(v)=y(v); x1(v)=x(v); end if disp0 >= 1 then //(2) subplot(212); end idx=3; w0 = n_supend_para0(idx) * 2. * %pi; q0 = n_supend_para0(idx+1); [x,y]=set_bef_dig( fs, 1.0 , w0, q0); if disp0 >= 1 then [db,phi] =dbphi( gz0(x,y, 2 * %pi * frq )); bode0(frq,db,phi); wstr1='(high) bef to absorb'; xtitle( wstr1 ); end s0=size(y); x2=zeros(1,s0(2)); y2=zeros(1,s0(2)); for v=1:s0(2) y2(v)=y(v); x2(v)=x(v); end if disp0 >= 1 then xset('window',wb0); // push old windows end // if disp0 >= 1 then endfunction //-------------------------------------------------- function [x1,y1,x2,y2,x3,y3]=set_N_iir_et_plot_bode2(disp0) frq=[100:25:Max_Freq]'; if disp0 >= 1 then wb0=xget('window'); // stack old window xset('window',disp0); // create new windows clf(); //(1) subplot(311); else disp('+set_N_iir_et_plot_bode2'); end idx=1; w0 = n_supend_para0(idx) * 2. * %pi; q0 = n_supend_para0(idx+1); [x,y]=set_bef_dig2( fs, 1.0 , w0, q0, 0.); if disp0 >= 1 then [db,phi] =dbphi( gz0(x,y, 2 * %pi * frq )); bode0(frq,db,phi); wstr1='bef to absorb, when mix0=0dB'; xtitle( wstr1 ); end s0=size(y); x1=zeros(1,s0(2)); y1=zeros(1,s0(2)); for v=1:s0(2) y1(v)=y(v); x1(v)=x(v); end if disp0 >= 1 then //(2) subplot(312); end idx=1; w0 = n_supend_para0(idx) * 2. * %pi; q0 = n_supend_para0(idx+1); [x,y]=set_bef_dig2( fs, 1.0 , w0, q0, -6); if disp0 >= 1 then [db,phi] =dbphi( gz0(x,y, 2 * %pi * frq )); bode0(frq,db,phi); wstr1='bef to absorb, when mix0= -6dB'; xtitle( wstr1 ); end s0=size(y); x2=zeros(1,s0(2)); y2=zeros(1,s0(2)); for v=1:s0(2) y2(v)=y(v); x2(v)=x(v); end if disp0 >= 1 then //(3) subplot(313); end idx=1; w0 = n_supend_para0(idx) * 2. * %pi; q0 = n_supend_para0(idx+1); [x,y]=set_bef_dig2( fs, 1.0 , w0, q0, -40); if disp0 >= 1 then [db,phi] =dbphi( gz0(x,y, 2 * %pi * frq )); bode0(frq,db,phi); wstr1='bef to absorb, when mix0=-40dB'; xtitle( wstr1 ); end s0=size(y); x3=zeros(1,s0(2)); y3=zeros(1,s0(2)); for v=1:s0(2) y3(v)=y(v); x3(v)=x(v); end if disp0 >= 1 then xset('window',wb0); // push old windows end // if disp0 >= 1 then endfunction //-------------------------------------------------- //********************************************************************/ //* G(z)= y[0]z0 + y[1]z-1 + y[2]z-2/ 1.0 - (x[1]z-1 + x[2]z-2) */ //* z0 = tannni kakeru dake x[0]=1.0 function [yout] = gzx( x ,y ,xwin0 ) s0=size(x); l0=s0(2); s1=size(xwin0); l1=s1(1); yout=zeros(l1,1); for v=1:l1 xw=xwin0(v); yi=y(1); yb=x(1); for ibx=1:(l0-1) yi= yi + y(1+ibx) * cos( xw * ts * ibx) - y(1+ibx) * sin( xw * ts * ibx ) * %i ; yb= yb - ( x(1+ibx) * cos( xw * ts * ibx) - x(1+ibx) * sin( xw * ts * ibx ) * %i) ; end yout(v,1)=yi / yb; end endfunction //-------------------------------------------------- function [x1,y1]=set_lpf_cheb_plot_bode(disp0,n_order0, w0x, ripple_db, mix0db1,mix0db2,mix0db3) frq=[100:25:Max_Freq]'; if disp0 >= 1 then wb0=xget('window'); // stack old window xset('window',disp0); // create new windows clf(); //(1) subplot(311); else disp('+set_lpf_cheb_et_plot_bod'); end mix0db=mix0db1; w0 = w0x * 2. * %pi; [x,y]=set_lpf_cheb_dig( fs, n_order0, w0, ripple_db, mix0db) if disp0 >= 1 then [db,phi] =dbphi( gzx(x,y, 2 * %pi * frq )); bode0(frq,db,phi); wstr1=sprintf('lpf_cheb2 mix0db=%f',mix0db); xtitle( wstr1 ); end s0=size(y); x1=zeros(1,s0(2)); y1=zeros(1,s0(2)); for v=1:s0(2) y1(v)=y(v); x1(v)=x(v); end //(2) if disp0 >= 1 then subplot(312); end mix0db=mix0db2; w0 = w0x * 2. * %pi; [x,y]=set_lpf_cheb_dig( fs, n_order0, w0, ripple_db, mix0db) if disp0 >= 1 then [db,phi] =dbphi( gzx(x,y, 2 * %pi * frq )); bode0(frq,db,phi); wstr1=sprintf('lpf_cheb2 mix0db=%f',mix0db); xtitle( wstr1 ); end s0=size(y); x2=zeros(1,s0(2)); y2=zeros(1,s0(2)); for v=1:s0(2) y2(v)=y(v); x2(v)=x(v); end //(3) if disp0 >= 1 then subplot(313); end mix0db=mix0db3; w0 = w0x * 2. * %pi; [x,y]=set_lpf_cheb_dig( fs, n_order0, w0, ripple_db, mix0db) if disp0 >= 1 then [db,phi] =dbphi( gzx(x,y, 2 * %pi * frq )); bode0(frq,db,phi); wstr1=sprintf('lpf_cheb2 mix0db=%f',mix0db); xtitle( wstr1 ); end s0=size(y); x3=zeros(1,s0(2)); y3=zeros(1,s0(2)); for v=1:s0(2) y3(v)=y(v); x3(v)=x(v); end if disp0 >= 1 then xset('window',wb0); // push old windows end // if disp0 >= 1 then endfunction //================================================== // Display Option phase_display=0; // When 1, show phase, Other case only show magnitudes: (default 0) phase_type=0; // When 1, phase representation between -360 and 0 degrees, Other case phase continuous representation: (default 0) //-------------------------------------------------- function [length1]=rt_size(wdat1) s0=size(wdat1); if s0(1) > s0(2) then length1 = s0(1); else length1 = s0(2); end endfunction //-------------------------------------------------- function display_portion(disp0,wdat0) wb0=xget('window'); // stack old window [length1]=rt_size(wdat0); if disp0 >= 1 then xset('window',disp0); // create new windows clf(); //(1) plot wave all subplot(311); xziku0=zeros(1,length1); for v=1:length1 xziku0(v)=v; end plot(xziku0,wdat0,'b'); wstr1='waveform all : ' + memo1; xtitle( wstr1 ); //(2) plot wave portion subplot(312); length12=ep1-sp1+1; xziku1=zeros(1,length12); wdat1=zeros(1,length12); for v=1:length12 xziku1(v)=v+sp1-1; wdat1(v)=wdat0(v+sp1-1); end plot(xziku1,wdat1,'r'); wstr1='waveform portion : ' ; xtitle( wstr1 ); subplot(311); plot(xziku1,wdat1,'r'); //(3) plot fft of the portion subplot(313); [db_fft1,phi_fft1,rt_code]=do_fft_wav(wdat0); if rt_code >= 0 then [frq1,sfrq1,is1,ie1]=set_frq(0); gainplot(frq1,db_fft1,phi_fft1); wstr1 = 'frequency response of waveform part of which '; wstr1 = wstr1 + ' start=' + string(sp1) + ' points=' + string( fft_point1); xtitle( wstr1 ); end end // disp0 >= 1 then xset('window',wb0); // push old windows endfunction //-------------------------------------------------- // + event handler function my_eventhandler(win,x,y,ibut) global eh_var1 global eh_var2 global eh_var3 global eh_var4 global eh_var5 global eh_fm_disp0 if ibut==-1000 then return,end [x,y]=xchange(x,y,'i2f') xinfo(msprintf('Event code %d at mouse position is (%f,%f)',ibut,x,y)) if ibut== 100 // d if (eh_var2 + eh_var3 - 1) < D_QTY then eh_var2=eh_var2+1 end plot_areas(eh_var1, eh_var2, eh_var3) elseif ibut== 117 // u if eh_var2 > 1 then eh_var2=eh_var2-1 end plot_areas(eh_var1, eh_var2, eh_var3) elseif ibut==3 // Left mouse button has been clicked eh_var4=x if (eh_var2 + eh_var3 - 1) < D_QTY then // d eh_var2=eh_var2+1 end plot_areas(eh_var1, eh_var2, eh_var3) elseif ibut==5 // Right mouse button has been clicked eh_var5=x if eh_var2 > 1 then // u eh_var2=eh_var2-1 end plot_areas(eh_var1, eh_var2, eh_var3) end // seteventhandler('my_eventhandler') // seteventhandler('') //suppress the event handler endfunction //-------------------------------------------------- function plot_areas(disp_code, nstartline, nlines) w0=xget('window'); // stack old window if disp_code >= 1 then xset('window',disp_code); // create new windows clf(); // if (T_QTY == 2) | (T_QTY ==5) then idy0=8; idy1=9; elseif T_QTY == 3 then idy0=10; idy1=11; end // frq=[Min_Freq:25:Max_Freq]; frqs=size(frq); db_fft=zeros(1,frqs(2)); phi_fft=zeros(1,frqs(2)); for vloop=1:nlines subplot(nlines ,2, 2*(vloop-1)+1); // if another_tube_on==1 then RNx0=lrpa(vloop+nstartline-1,idy1); else RNx0=0.; end if RNx0 >= 1 then // when only Rn=1 plot_area2n(vloop+nstartline-1,vloop); else plot_area2(vloop+nstartline-1,vloop); end //........................... subplot(nlines ,2, 2*(vloop-1)+2); for w=1:frqs(2) db_fft(w)= db_all(w,vloop+nstartline-1); phi_fft(w)=phi_all(w,vloop+nstartline-1); end if phase_display == 1 then bode(frq,db_fft,phi_fft); else gainplot(frq,db_fft,phi_fft); end wstr1 = 'frequency response theorycal: '; xtitle( wstr1 ); end // for vloop=1:nlines end // if disp_code >= 1 then xset('window',w0); // push old windows endfunction //---------------------------------------------------- function plot_area2(ni,nth) plot2d(0,0,1,"010","",[-5,-10,20,10],[5,5,4,5]); // wake wo egaku if nth == 1 then wstr1 = ' No. ' + string(ni) + ' :memo(' + memo1 + ')'; else wstr1 = ' No. ' + string(ni) end xstring(-4,6, wstr1); if T_QTY == 4 then //Three Tubes Model of T type ya1=sqrt(A1(ni) ); ya2=sqrt(A2(ni) ); ya3=sqrt(A3(ni) ); if ya2 > ya1 then yy0=ya2; else yy0=ya1; end yy0=10- (yy0 * 2); // Area scale is double than other Tube Models. yy1=ya1 + yy0; yy2=ya2 + yy0; yy3=yy0 - L3(ni); xx1=L1(ni) - (ya3 / 2); xx2=L1(ni) + (ya3 / 2); xp=[0 0 xx1 xx1 xx2 xx2 L1(ni)+L2(ni) L1(ni)+L2(ni) L1(ni) L1(ni) 0 ]'; yp=[yy1 yy0 yy0 yy3 yy3 yy0 yy0 yy2 yy2 yy1 yy1]'; xpolys( xp, yp, [5 5 5 5 5 5 5 5 5 5 5] ); xar=[-2 -0.5]; yar=[yy0+ ( ya1 / 2 ) yy0+ ( ya1 / 2 ) ]; if scilab_version_number <= 3 then xarrows(xar, yar, 1, 5); else xarrows(xar, yar, 20, 5); end yar=[yy0 + ( ya2 / 2 ) yy0+ ( ya2 / 2 ) ]; xar=[ L1(ni)+L2(ni)+0.5 L1(ni)+L2(ni)+2]; if scilab_version_number <= 3 then xarrows(xar, yar, 1, 5); else xarrows(xar, yar, 20, 5); end xstring(xx1, yy3 - 1, "CLOSED"); xstring(0,-9,"Attension: Area scale is double than others Tube Models"); // elseif T_QTY == 3 then //Three Tubes Model of serial type ya1=sqrt(A1(ni) ); ya2=sqrt(A2(ni) ); ya3=sqrt(A3(ni) ); xp=[0 0 L1(ni) L1(ni) L1(ni)+L2(ni) L1(ni)+L2(ni) L1(ni)+L2(ni)+L3(ni) L1(ni)+L2(ni)+L3(ni) L1(ni)+L2(ni) L1(ni)+L2(ni) L1(ni) L1(ni) 0 ]'; yp=[ya1 -ya1 -ya1 -ya2 -ya2 -ya3 -ya3 ya3 ya3 ya2 ya2 ya1 ya1 ]'; xpolys( xp, yp, [5 5 5 5 5 5 5 5 5 5 5 5] ); xar=[-2 -0.5]; yar=[0 0]; if scilab_version_number <= 3 then xarrows(xar, yar, 1, 5); else xarrows(xar, yar, 20, 5); end xar=[ L1(ni)+L2(ni)+L3(ni)+0.5 L1(ni)+L2(ni)+L3(ni)+2]; if scilab_version_number <= 3 then xarrows(xar, yar, 1, 5); else xarrows(xar, yar, 20, 5); end else //Two Tubes Model 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]; if scilab_version_number <= 3 then xarrows(xar, yar, 1, 5); else xarrows(xar, yar, 20, 5); end xar=[ L1(ni)+L2(ni)+0.5 L1(ni)+L2(ni)+2]; if scilab_version_number <= 3 then xarrows(xar, yar, 1, 5); else xarrows(xar, yar, 20, 5); end if T_QTY==5 then xar=[ L1(ni)+L2(ni)-0.5 L1(ni)+L2(ni)-2]; if scilab_version_number <= 3 then xarrows(xar, yar, 1, 3); else xarrows(xar, yar, 20, 3); end end end endfunction //------------------------------------------------------------------ function plot_area2n(ni,nth) plot2d(0,0,1,"010","",[-5,-10,20,10],[5,5,4,5]); // wake wo egaku if nth == 1 then wstr1 = ' No. ' + string(ni) + ' :memo(' + memo1 + ')' + ' another'; else wstr1 = ' No. ' + string(ni) + ' another'; end xstring(-4,6, wstr1); //Two Tubes Model ya1=sqrt(A1_n(ni) ); ya2=sqrt(A2_n(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, 20, 5); xar=[ L1(ni)+L2(ni)+0.5 L1(ni)+L2(ni)+2]; xarrows(xar, yar, 20, 5); endfunction //------------------------------------------ function [db_all,phi_all]=make_bode() // ni=1; // yg is only 1st calculate frq=[Min_Freq:25:Max_Freq]; frqs=size(frq); resp=zeros(1,frqs(2)); respt=zeros(1,frqs(2)); resph=zeros(1,frqs(2)); resphn=zeros(1,frqs(2)); respln=zeros(1,frqs(2)); respa=zeros(1,frqs(2)); respbefl=zeros(1,frqs(2)); respbefh=zeros(1,frqs(2)); db_all=zeros(frqs(2),D_QTY); phi_all=zeros(frqs(2),D_QTY); // disp(' + please wait for some time to calculate bode '); disp(' + warning: these bode are when yg is 1st. '); disp(' + warning: if yg is changed, these bode may NOT be actual ones. '); // if (T_QTY == 2) | (T_QTY ==5) then idy0=8; idy1=9; idy_noise_yg=6; elseif T_QTY == 3 then idy0=10; idy1=11; idy_noise_yg=8; end rl0=rl_default; // stp=0; if ni > 1 then for v=1:(ni-1) stp=stp + LL(v); end end stp=stp+N1(ni); y00=0.; for v=1:(N2(ni)+N3(ni)+1) y00= y00 + (yg(stp + v) + ygn(stp + v)) ; end for w=1:frqs(2) xw=2 * %pi * frq(w) * ts; yi=0; for v=1:(N2(ni)+N3(ni)+1) yi= yi + (yg(stp + v) + ygn(stp + v)) * %e^(-1. * xw * (v-1) * %i); end // for v resp(w)=yi / y00; resph(w)=hpf1(2 * %pi * frq(w) ); resphn(w)=hpf1n(2 * %pi * frq(w) ); respln(w)=lpf1n(2 * %pi * frq(w) ); respbefl(w)=gz0(x_coef7,y_coef7, 2 * %pi * frq(w) ); respbefh(w)=gz0(x_coef8,y_coef8, 2 * %pi * frq(w) ); end // for w // // bef if (T_QTY == 2) | (T_QTY ==5) then idy4=10; elseif T_QTY == 3 then idy4=12; end for vloop=1:D_QTY // if vloop >= 2 then if lrpa(vloop-1, idy_noise_yg) <> lrpa(vloop, idy_noise_yg) then ni=vloop; stp=0; if ni > 1 then for v=1:(ni-1) stp=stp + LL(v); end end stp=stp+N1(ni); y00=0.; for v=1:(N2(ni)+N3(ni)+1) y00= y00 + (yg(stp + v) + ygn(stp + v)) ; end for w=1:frqs(2) xw=2 * %pi * frq(w) * ts; yi=0; for v=1:(N2(ni)+N3(ni)+1) yi= yi + (yg(stp + v) + ygn(stp + v)) * %e^(-1. * xw * (v-1) * %i); end // for v resp(w)=yi / y00; end // for w wstr1=sprintf(' + warning: re-calculate yg bode. %d', vloop); disp(wstr1); end end // if lrpa(vloop,idy0) <> 9 then // other than xcode0 = 9 rl0=lrpa(vloop,idy0); else rl0=rl_default; end // if another_tube_on == 1 then RNx0=lrpa(vloop,idy1); else RNx0=0.; end if RNx0 >= 1 then // when only Rn=1 rl0=rl_nose; for w=1:frqs(2) respt(w)= two_tubes2n(2 * %pi * frq(w), vloop,rg_closed ,rl0); if nose_lpf_times== 2 then respa(w)=resp(w) * respt(w) * respln(w) * respln(w) * resphn(w); // lpf else respa(w)=resp(w) * respt(w) * respln(w) * resphn(w); // lpf end if lrpa(vloop,idy4) == 1 then respa(w)=respa(w) * respbefl(w); elseif lrpa(vloop,idy4) == 2 then respa(w)=respa(w) * respbefh(w); end if nose_hpf_times == 2 then respa(w)=respa(w) * resphn(w); end end else // // when only Rn=1 for w=1:frqs(2) // mark p1 if T_QTY == 2 then // respt(w)= two_tubes2(2 * %pi * frq(w), vloop,rg_closed ,rl0); elseif (T_QTY == 3) | (T_QTY == 4) then respt(w)= three_tubes(2 * %pi * frq(w), vloop,rg_closed ,rl0); elseif T_QTY == 5 then respt(w)= two_tubes_rl_noise(2 * %pi * frq(w), vloop,rg_closed ,rl_default); end respa(w)=resp(w) * respt(w) * resph(w); if BEF_VOCAL_TRACT == 1 then if lrpa(vloop,idy4) == 1 then respa(w)=respa(w) * respbefl(w); elseif lrpa(vloop,idy4) == 2 then respa(w)=respa(w) * respbefh(w); end end if variable_filter_on ==1 then //respa(w)=resp_vari(vloop,w); respa(w)=respa(w) * resp_vari(vloop,w); end if variable_filter_on_hcf ==1 then //respa(w)=resp_vari(vloop,w); respa(w)=respa(w) * resp_vari_hcf(vloop,w); end end // mark p1 end // not // when only Rn=1 // if phase_type == 1 then [phi_fft1,db_fft1] =phasemag(respa,"m"); // range -360-0 else [db_fft1,phi_fft1] =dbphi(respa); end for w=1:frqs(2) db_all(w,vloop)=db_fft1(w); phi_all(w,vloop)=phi_fft1(w); end // end disp(' - calculation was done. '); endfunction //---- // - event handler //-------------------------------------------------- // - others function library //------------------------------------------------------------------------- function [D_QTY,T_QTY,lrpa,co_filter,k_gain0,memo1x,memo2x,memo3x,variable_filter,another_tube_on,variable_filter_hcf]=sub_demo1( inum ) if inum >= select_offset_page3 then [D_QTY,T_QTY,lrpa,co_filter,k_gain0,memo1x,memo2x,memo3x,variable_filter,another_tube_on,variable_filter_hcf]=sub_demo1_page3( inum ); elseif inum >= select_offset_page2 then [D_QTY,T_QTY,lrpa,co_filter,k_gain0,memo1x,memo2x,memo3x,variable_filter,another_tube_on,variable_filter_hcf]=sub_demo1_page2( inum ); else [D_QTY,T_QTY,lrpa,co_filter,k_gain0,memo1x,memo2x,memo3x,variable_filter,another_tube_on,variable_filter_hcf]=sub_demo1_page1( inum ); end endfunction //------------------------------------------------------------------------- // // for page 1 // function [D_QTY,T_QTY,lrpa,co_filter,k_gain0,memo1x,memo2x,memo3x,variable_filter,another_tube_on,variable_filter_hcf]=sub_demo1_page1( inum ) // setting for demonstration another_tube_on=0; variable_filter=0; variable_filter_hcf=0; // common 1 co_filter=[ 1 , 0, 500, 0.7 ; 2 , 1, 1000, 0.7 ; 3, 0, 10000, 0.7 ; 4, 2, 5000, 2.5 ; 5, 0, 1500,0.7 ; 6 , 2, 1250, 10.; 7 , 1, 1000, 0.7 ; 8, 0, 10000, 0.7 ; 9, 2, 2500, 10 ]; if inum == 8 then // when /ra/ disp('+ a sample of /ra/ is setting.'); memo1x='a sample of /ra/'; D_QTY=16; T_QTY=2; memo2x=' no, l1, r1, pitch(+-%), amplitude(0-1), vocal-source(0/h=1/s=2/k=3/t=4),lip-source(0/s=1/r_start=2/rl_noise=3),rl(if=9,use rl_default)'; lrpa=[1 , 0 , -0.6 , 0 , 0.3 , 0, 2, 0.7; // add start mark 2 , 0 , -0.4 , 0 , 0.3 , 0, 2, 0.7; // add start mark 3 , 0 , -0.2 , 0 , 0.3 , 0, 0, 0.7; 4 , 0 , 0.0 , 0 , 0.3 , 0, 0, 0.7; 5 , 0 , 0.2 , 0 , 0.3 , 0, 0, 0.7; 6 , 0 , 0.4 , 0 , 0.3 , 0, 0, 9; 7 , 0 , 0.6 , 0 , 0.4 , 0, 0, 9; 8 , 0 , 0.7 , 0 , 0.5 , 0, 0, 9; 9 , 0 , 0.75, 0 , 0.6 , 0, 0, 9; 10 , 0 , 0.75 , 0 , 0.7 , 0, 0, 9; 11 , 0 , 0.8 , 0 , 0.8 , 0, 0, 9; 12 , 0 , 0.8 , 0 , 0.9 , 0, 0, 9; 13 , 0 , 0.8 , 0 , 0.7 , 0, 0, 9; 14 , 0 , 0.8 , 0 , 0.7 , 0, 0, 9; 15 , 0 , 0.8 , -3 , 0.4 , 0, 0, 9; 16 , 0 , 0.8 , -5 , 0.1 , 0, 0, 9] ; elseif inum == 11 then // when // ikinari hatuon suruno ga ta disp('+ a sample of dabug is setting. '); memo1x='a sample of dabug'; D_QTY=11; T_QTY=2; memo2x=' no, l1, r1, pitch(+-%), amplitude(0-1), vocal-source(0/h=1/s=2/k=3/t=4),lip-source(0/s=1/r_start=2/rl_noise=3),rl(if=9,use rl_default)'; lrpa=[1 , 0 , 0 , 0 , 0.05 , 4, 0, 0.1; 2 , 0 , 0 , 0 , 0.02 , 4, 0, 0.1; 3 , 0 , 0.3 , 0 , 0.3 , 0, 0, 0.5; 4 , 0 , 0.4 , 0 , 0.4 , 0, 0, 0.7; 5 , 0 , 0.5 , 0 , 0.5 , 0, 0, 9; 6 , 0 , 0.6, 0 , 0.6 , 0, 0, 9; 7 , 0 , 0.7 , 0 , 0.7 , 0, 0, 9; 8 , 0 , 0.8 , 0 , 0.8 , 0, 0, 9; 9 , 0 , 0.8 , 0 , 0.9 , 0, 0, 9; 10 , 0 , 0.8 , 0 , 0.7 , 0, 0, 9; 11 , 0 , 0.8 , 0 , 0.7 , 0, 0, 9]; co_filter(5,3)=2000; // if 2000, then it may be da elseif inum == 7 then // when /na/ disp('+ a sample of /na/ is setting.'); memo1x='a sample of /na/'; D_QTY=9; T_QTY=2; another_tube_on=1; variable_filter=1; // when variable_filter is 1, another_tube_on must be 1 ! variable_filter_hcf=1; // when variable_filter_hcf is 1, another_tube_on and variable_filter must be 1 ! memo2x=' no, l1, r1, pitch(+-%), amplitude(0-1), vocal-source(0/h=1/s=2/k=3/t=4),lip-source(0/s=1/r_start=2/rl_noise=3),rl(if=9,use rl_default), Rn(0. ... 1.), bef(off=0/low=1/high=2), always 4,variable filter type(off=0/lpf=1/hpf=2/bpf=3/bef=4/double_bef=5),fc(when double fc1x10000+fc2),q,gain(when bef,this value dip gain dB,ex -20(dB), hcf(off=0/on=1),fc,gain)'; lrpa=[ 1 , 0 , 0 , 0 , 0.6 , 0, 0, 9, 0., 0, 4, 4, 530, 0.7, -60., 1, 2000, -60; // set tune bef's fc to erase 1nd peak around 530Hz in this case. And set hcf remain other higher one peak. 2 , 0 , 0 , 0 , 0.6 , 0, 0, 9, 0., 0, 4, 4, 530, 0.7, -60., 1, 2000, -60; 3 , 0 , 0 , 0 , 0.6 , 0, 0, 9, 0., 0, 4, 4, 530, 0.7, -60., 1, 2000, -60; 4 , 0 , 0, 0 , 0.6 , 0, 0, 9, 0., 0, 4, 4, 530, 0.7, -60., 1, 2000, -40; // 5 , 0 , 0.3 , 0 , 0.7 , 0, 0, 9, 0., 0, 4, 4, 530, 0.7, -20., 1, 3000, -30; // effect on area of /a/ 6 , 0 , 0.5 , 0 , 0.8 , 0, 0, 9, 0., 0, 4, 4, 530, 0.7, -20, 1, 3000, -30; 7 , 0 , 0.6 , 0 , 0.7 , 0, 0, 9, 0., 0, 4, 4, 530, 0.7, 0, 1, 3000, -30; 8 , 0 , 0.7 ,-3 , 0.4 , 0, 0, 9, 0., 0, 4, 0, 0, 0, 0, 1, 3000, -30; 9 , 0 , 0.8 ,-5 , 0.1 , 0, 0, 9, 0., 0, 4, 0, 0, 0, 0, 1, 3000, -30] ; elseif inum == 9 then // when /N(nn)/ disp('+ a sample of /N(nn)/ is setting. '); memo1x='a sample of /N(nn)/'; D_QTY=8; T_QTY=2; another_tube_on=0; variable_filter=1; // when variable_filter is 1, another_tube_on must be 1 ! variable_filter_hcf=1; // when variable_filter_hcf is 1, another_tube_on and variable_filter must be 1 ! memo2x=' no, l1, r1, pitch(+-%), amplitude(0-1), vocal-source(0/h=1/s=2/k=3/t=4),lip-source(0/s=1/r_start=2/rl_noise=3),rl(if=9,use rl_default), Rn(0. ... 1.), bef(off=0/low=1/high=2), always 4,variable filter type(off=0/lpf=1/hpf=2/bpf=3/bef=4),fc,q,gain(when bef,this value dip gain dB,ex -20(dB), hcf(off=0/on=1),fc,gain)'; lrpa=[ 1 , 0 , 0 , 0 , 0.6 , 0, 0, 9, 0., 0, 4, 4, 1600, 0.7, -60, 1, 3000, -60; // set tune bef's fc to erase 2nd peak around 1600Hz in this case, * honrai arubeki peak wo kesu. 2 , 0 , 0 , 0 , 0.6 , 0, 0, 9, 0., 0, 4, 4, 1600, 0.7, -60, 1, 3000, -60; 3 , 0 , 0 , 0 , 0.6 , 0, 0, 9, 0., 0, 4, 4, 1600, 0.7, -60, 1, 3000, -60; 4 , 0 , 0 , 0 , 0.6 , 0, 0, 9, 0., 0, 4, 4, 1600, 0.7, -60, 1, 3000, -60; 5 , 0 , 0 , 0 , 0.6 , 0, 0, 9, 0., 0, 4, 4, 1600, 0.7, -60, 1, 3000, -60; 6 , 0 , 0 , 0 , 0.6 , 0, 0, 9, 0., 0, 4, 4, 1600, 0.7, -60, 1, 3000, -60; 7 , 0 , 0 , 0 , 0.6 , 0, 0, 9, 0., 0, 4, 4, 1600, 0.7, -60, 1, 3000, -60; 8 , 0 , 0 , 0 , 0.6 , 0, 0, 9, 0., 0, 4, 4, 1600, 0.7, -60, 1, 3000, -60]; //..................................................... // older setting //D_QTY=10; //T_QTY=2; //another_tube_on=1; //memo2x=' no, l1, r1, pitch(+-%), amplitude(0-1), vocal-source(0/h=1/s=2/k=3/t=4),lip-source(0/s=1/r_start=2/rl_noise=3),rl(if=9,use rl_default), Rn(0. ... 1.), bef(off=0/low=1/high=2)'; // //lrpa=[ 1 , 0 , 0 , 0 , 0.3 , 0, 0, 9, 1., 1; // 2 , 0 , 0 , 0 , 0.3 , 0, 0, 9, 1., 1; // 3 , 0 , 0 , 0 , 0.3 , 0, 0, 9, 1., 1; // 4 , 0 , 0 , 0 , 0.3 , 0, 0, 9, 1., 1; // 5 , 0 , 0 , 30 , 0.3 , 0, 0, 9, 1., 1; // 6 , 0 , 0 , 30 , 0.3 , 0, 0 ,9, 1., 1; // 7 , 0 , 0, 30 , 0.3 , 0, 0 , 9, 1., 1; // 8 , 0 , 0 , 30 , 0.3 , 0, 0, 9, 1., 1; // 9 , 0 , 0, 30 , 0.3 , 0, 0, 9, 1., 1; // 10 , 0 , 0 , 30 , 0.3 , 0, 0, 9, 1., 1] ; //.................................................... // elseif inum == 6 then // when /ma/ disp('+ a sample of /ma/ is setting. '); memo1x='a sample of /ma/'; D_QTY=12; T_QTY=2; another_tube_on=1; variable_filter=1; // when variable_filter is 1, another_tube_on must be 1 ! variable_filter_hcf=1; // when variable_filter_hcf is 1, another_tube_on and variable_filter must be 1 ! memo2x=' no, l1, r1, pitch(+-%), amplitude(0-1), vocal-source(0/h=1/s=2/k=3/t=4),lip-source(0/s=1/r_start=2/rl_noise=3),rl(if=9,use rl_default), Rn(0. ... 1.), bef(off=0/low=1/high=2), always 4,variable filter type(off=0/lpf=1/hpf=2/bpf=3/bef=4/double_bef=5),fc(when double fc1x10000+fc2),q,gain(when bef,this value dip gain dB,ex -20(dB), hcf(off=0/on=1),fc,gain)'; lrpa=[ 1 , 0 , 0.5 , 0 , 0.6 , 0, 0, 9, 0., 0, 4, 5, 6801400, 0.7, -60., 1, 4000, -60; //// set tune bef's fc to erase 1st and 2nd peak, And set hcf remain other higher two peaks. 2 , 0 , 0.5 , 0 , 0.6 , 0, 0, 9, 0., 0, 4, 5, 6801400, 0.7, -60., 1, 4000, -60; 3 , 0 , 0.5, 0 , 0.6 , 0, 0, 9, 0., 0, 4, 5, 6801400, 0.7, -60., 1, 4000, -40; // 4 , 0 , 0.5 , 0 , 0.7 , 0, 0, 9, 0., 0, 4, 5, 6801400, 0.7, -60., 1, 4000, -60; 5 , 0 , 0.6 , 0 , 0.7 , 0, 0, 9, 0., 0, 4, 5, 7201300, 0.7, -50., 1, 4000, -60; // escape like sound /d/, by using hcf as nose effect 6 , 0 , 0.6 , 0 , 0.7 , 0, 0, 9, 0., 0, 4, 5, 7201300, 0.7, -30., 1, 4000, -60; 7 , 0 , 0.7 , 0 , 0.7 , 0, 0, 9, 0., 0, 4, 5, 7801300, 0.7, -20., 1, 4000, -60; 8 , 0 , 0.7 , 0 , 0.7 , 0, 0, 9, 0., 0, 4, 5, 7801300, 0.7, -10., 1, 4000, -60; 9 , 0 , 0.8 , 0 , 0.7 , 0, 0, 9, 0., 0, 4, 5, 7801300, 0.7, -1., 1, 4000, -60; 10 ,0 , 0.8 , 0 , 0.6 , 0, 0, 9, 0., 0, 4, 5, 7801300, 0.7, -1., 1, 4000, -60; 11 ,0 , 0.8 , 0 , 0.5 , 0, 0, 9, 0., 0, 4, 5, 7801300, 0.7, -1., 1, 4000, -60; 12 ,0 , 0.8 , 0 , 0.3 , 0, 0, 9, 0., 0, 4, 5, 7801300, 0.7, -1., 1, 4000, -60]; elseif inum == 10 then // when /mu/ disp('+ a sample of /mu/ is setting. '); memo1x='a sample of /mu/'; D_QTY=12; T_QTY=2; another_tube_on=1; variable_filter=1; // when variable_filter is 1, another_tube_on must be 1 ! variable_filter_hcf=1; // when variable_filter_hcf is 1, another_tube_on and variable_filter must be 1 ! memo2x=' no, l1, r1, pitch(+-%), amplitude(0-1), vocal-source(0/h=1/s=2/k=3/t=4),lip-source(0/s=1/r_start=2/rl_noise=3),rl(if=9,use rl_default), Rn(0. ... 1.), bef(off=0/low=1/high=2), always 4,variable filter type(off=0/lpf=1/hpf=2/bpf=3/bef=4/double_bef=5),fc(when double fc1x10000+fc2),q,gain(when bef,this value dip gain dB,ex -20(dB), hcf(off=0/on=1),fc,gain)'; lrpa=[ 1 , 0 , 0 , 0 , 0.7 , 0, 0, 9, 0., 0, 4, 5, 5301600, 0.7, -60., 1, 4000, -60; //// set tune bef's fc to erase 1st and 2nd peak around 530Hz 1600Hz in this case, And set hcf remain other higher two peaks. 2 , 0 , 0 , 0 , 0.7 , 0, 0, 9, 0., 0, 4, 5, 5301600, 0.7, -60., 1, 4000, -60; 3 , 0 , 0, 0 , 0.7 , 0, 0, 9, 0., 0, 4, 5, 5301600, 0.7, -60., 1, 4000, -40; 4 , 0 , 0 , 0 , 0.7 , 0, 0, 9, 0., 0, 4, 5, 5301600, 0.7, -60., 1, 4000, -60; 5 , 0 , 0 , 0 , 0.7 , 0, 0, 9, 0., 0, 4, 5, 5301600, 0.7, -50., 1, 4000, -60; 6 , 0 , 0 , 0 , 0.7 , 0, 0, 9, 0., 0, 4, 5, 5301600, 0.7, -30., 1, 4000, -60; 7 , 0 , 0 , 0 , 0.7 , 0, 0, 9, 0., 0, 4, 5, 5301600, 0.7, -20., 1, 4000, -60; 8 , 0 , 0 , 0 , 0.7 , 0, 0, 9, 0., 0, 4, 5, 5301600, 0.7, -10., 1, 4000, -60; 9 , 0 , 0 , 0 , 0.7 , 0, 0, 9, 0., 0, 4, 5, 5301600, 0.7, -1., 1, 4000, -60; 10 , 0 , 0 , 0 , 0.6 , 0, 0, 9, 0., 0, 4, 5, 5301600, 0.7, -1., 1, 4000, -60; 11 , 0 , 0 , 0 , 0.5 , 0, 0, 9, 0., 0, 4, 5, 5301600, 0.7, -1., 1, 4000, -60; 12 ,0 , 0 , 0 , 0.3 , 0, 0, 9, 0., 0, 4, 5, 5301600, 0.7, -1., 1, 4000, -60]; elseif inum == 5 then // when /ta/ // ikinari hatuon suruno ga ta disp('+ a sample of /ta/ is setting. '); memo1x='a sample of /ta/'; D_QTY=11; T_QTY=2; memo2x=' no, l1, r1, pitch(+-%), amplitude(0-1), vocal-source(0/h=1/s=2/k=3/t=4),lip-source(0/s=1/r_start=2/rl_noise=3),rl(if=9,use rl_default)'; lrpa=[1 , 0 , 0 , 0 , 0.0 , 4, 0, 0.1; 2 , 0 , 0 , 0 , 0.02 , 4, 0, 0.1; 3 , 0 , 0.3 , 0 , 0.3 , 0, 0, 0.5; 4 , 0 , 0.4 , 0 , 0.4 , 0, 0, 0.7; 5 , 0 , 0.5 , 0 , 0.5 , 0, 0, 9; 6 , 0 , 0.6, 0 , 0.6 , 0, 0, 9; 7 , 0 , 0.7 , 0 , 0.7 , 0, 0, 9; 8 , 0 , 0.8 , 0 , 0.8 , 0, 0, 9; 9 , 0 , 0.8 , 0 , 0.9 , 0, 0, 9; 10 , 0 , 0.8 , 0 , 0.7 , 0, 0, 9; 11 , 0 , 0.8 , 0 , 0.7 , 0, 0, 9]; // co_filter(5,3)=2000; // if 2000, then it may be da // elseif inum == 4 then // when /ka/ disp('+ a sample of /ka/ is setting. '); memo1x='a sample of /ka/'; D_QTY=10; T_QTY=2; memo2x=' no, l1, r1, pitch(+-%), amplitude(0-1), vocal-source(0/h=1/s=2/k=3/t=4),lip-source(0/s=1/r_start=2/rl_noise=3),rl(if=9,use rl_default)'; lrpa=[1 , 0 , 0 , 0 , 0.05 , 3, 0, 0.5; // needs resonance effect, rl //2 , 0 , 0 , 0 , 0.02 , 3, 0, 0.1; // erase this line cause /t/ sound 2 , 0 , 0.3 , 0 , 0.3 , 0, 0, 0.5; 3 , 0 , 0.4 , 0 , 0.4 , 0, 0, 0.7; 4 , 0 , 0.5 , 0 , 0.5 , 0, 0, 9; 5 , 0 , 0.6, 0 , 0.6 , 0, 0, 9; 6 , 0 , 0.7 , 0 , 0.7 , 0, 0, 9; 7 , 0 , 0.8 , 0 , 0.8 , 0, 0, 9; 8 , 0 , 0.8 , 0 , 0.9 , 0, 0, 9; 9 , 0 , 0.8 , 0 , 0.7 , 0, 0, 9; 10 , 0 , 0.8 , 0 , 0.7 , 0, 0, 9]; // /ga/ setting // D_QTY=10; // T_QTY=2; //lrpa=[1 , 0 , -0.3 , 0 , 0.05 , 3, 0, 0.9; // needs resonance effect, rl // //2 , 0 , 0 , 0 , 0.02 , 3, 0, 0.1; // erase this line cause /t/ sound // 2 , 0 , 0.3 , 0 , 0.3 , 0, 0, 0.5; // 3 , 0 , 0.4 , 0 , 0.4 , 0, 0, 0.7; // 4 , 0 , 0.5 , 0 , 0.5 , 0, 0, 9; // 5 , 0 , 0.6, 0 , 0.6 , 0, 0, 9; // 6 , 0 , 0.7 , 0 , 0.7 , 0, 0, 9; // 7 , 0 , 0.8 , 0 , 0.8 , 0, 0, 9; // 8 , 0 , 0.8 , 0 , 0.9 , 0, 0, 9; // 9 , 0 , 0.8 , 0 , 0.7 , 0, 0, 9; // 10 , 0 , 0.8 , 0 , 0.7 , 0, 0, 9]; // elseif inum == 3 then // when /sa/ disp('+ a sample of /sa/ is setting. '); memo1x='a sample of /sa/'; D_QTY=17; T_QTY=2; memo2x=' no, l1, r1, pitch(+-%), amplitude(0-1), vocal-source(0/h=1/s=2/k=3/t=4),lip-source(0/s=1/r_start=2/rl_noise=3),rl(if=9,use rl_default)'; lrpa=[1 , 0 , 0 , -3 , 0.05 , 2, 1, 0.01; 2 , 0 , 0 , -3 , 0.05 , 2, 1, 0.01; 3 , 0 , 0 , -3 , 0.05 , 2, 1, 0.01; 4 , 0 , 0 , -3 , 0.05 , 2, 1, 0.01; 5 , 0 , 0 , -3 , 0.05 , 2, 1, 0.01; 6 , 0 , 0 , -3 , 0.05 , 2, 1, 0.01; 7 , 0 , 0 , -3 , 0.05 , 2, 1, 0.01; 8 , 0 , 0 , -3 , 0.05 , 2, 1, 0.01; 9 , 0 , 0 , -3 , 0.05 , 2, 1, 0.01; 10 , 0 , 0.3 , 0 , 0.2 , 0, 0, 0.5; 11 , 0 , 0.4 , 0 , 0.2 , 0, 0, 0.7; 12 , 0 , 0.5 , 0 , 0.2 , 0, 0, 9; 13 , 0 , 0.6 , 0 , 0.3 , 0, 0, 9; 14 , 0 , 0.7 , 0 , 0.5 , 0, 0, 9; 15 , 0 , 0.8, 0 , 0.6 , 0, 0, 9; 16 , 0 , 0.8, 0 , 0.7 , 0, 0, 9; 17 , 0 , 0.7, 0 , 0.8 , 0, 0, 9]; elseif inum == 2 then // when /ha/ disp('+ a sample of /ha/ is setting. '); memo1x='a sample of /ha/'; D_QTY=16; T_QTY=2; memo2x=' no, l1, r1, pitch(+-%), amplitude(0-1), vocal-source(0/h=1/s=2/k=3/t=4),lip-source(0/s=1/r_start=2/rl_noise=3),rl(if=9,use rl_default)'; lrpa=[1 , 0 , 0.6 , -3 , 0.2 , 1, 0, 9; // hazime kara kanari l1,r1 wo /a/ poku suru 2 , 0 , 0.6 , -3 , 0.2 , 1, 0, 9; 3 , 0 , 0.6 , -3 , 0.2 , 1, 0, 9; 4 , 0 , 0.6 , -3 , 0.2 , 1, 0, 9; 5 , 0 , 0.6 , -1 , 0.2 , 1, 0, 9; 6 , 0 , 0.6 , 0 , 0.3 , 1, 0, 9; 7 , 0 , 0.65 , 0 , 0.4 , 1, 0, 9; 8 , 0 , 0.7 , 0 , 0.5 , 0, 0, 9; 9 , 0 , 0.75, 0 , 0.6 , 0, 0, 9; 10 , 0 , 0.75 , 0 , 0.7 , 0, 0, 9; 11 , 0 , 0.8 , 0 , 0.8 , 0, 0, 9; 12 , 0 , 0.8 , 0 , 0.9 , 0, 0, 9; 13 , 0 , 0.8 , 0 , 0.7 , 0, 0, 9; 14 , 0 , 0.8 , 0 , 0.7 , 0, 0, 9; 15 , 0 , 0.8 , -3 , 0.4 , 0, 0, 9; 16 , 0 , 0.8 , -5 , 0.1 , 0, 0, 9] ; // elseif inum == 1 then // when /a/ disp('+ a sample of /a/ is setting. '); memo1x='a sample of /a/'; D_QTY=12; T_QTY=2; memo2x=' no, l1, r1, pitch(+-%), amplitude(0-1), vocal-source(0/h=1/s=2/k=3/t=4),lip-source(0/s=1/r_start=2/rl_noise=3),rl(if=9,use rl_default)'; lrpa=[1 , 0 , 0.5 , -3 , 0.1 , 0, 0, 9; 2 , 0 , 0.55 , -1 , 0.2 , 0, 0, 9; 3 , 0 , 0.6 , 0 , 0.3 , 0, 0, 9; 4 , 0 , 0.65 , 0 , 0.4 , 0, 0, 9; 5 , 0 , 0.7 , 0 , 0.5 , 0, 0 ,9; 6 , 0 , 0.75, 0 , 0.6 , 0, 0 , 9; 7 , 0 , 0.75 , 0 , 0.7 , 0, 0, 9; 8 , 0 , 0.8 , 0 , 0.8 , 0, 0, 9; 9 , 0 , 0.8 , 0 , 0.9 , 0, 0, 9; 10 , 0 , 0.8 , 0 , 0.7 , 0, 0, 9; 11 , 0 , 0.8 , -3 , 0.4 , 0, 0, 9; 12 , 0 , 0.8 , -5 , 0.1 , 0, 0, 9] ; // end // common 2 memo3x='no, type(LPF=0,BPF=1,HPF=2), f0, Q'; k_gain0=[7.0 , 0.0]; // check s0=size(lrpa); if D_QTY > s0(1) then disp('ERROR: D_QTY > size of lrpa: force to set D_QTY as size of lrpa. At function sub_demo1'); D_QTY = s0(1); elseif D_QTY < s0(1) then disp('WARNING: D_QTY < size of lrpa. At function sub_demo1'); end endfunction // for page 1 //------------------------------------------------------------------------- // // for page 2 // function [D_QTY,T_QTY,lrpa,co_filter,k_gain0,memo1x,memo2x,memo3x,variable_filter,another_tube_on,variable_filter_hcf]=sub_demo1_page2( inum0 ) // setting for demonstration inum=inum0 - select_offset_page2; another_tube_on=0; variable_filter=0; variable_filter_hcf=0; if inum == 11 then // when debug disp('+ This is for debug setting.'); memo1x='+ This is for debug setting.'; D_QTY=9; T_QTY=2; memo2x=' no, l1, r1, pitch(+-%), amplitude(0-1), vocal-source(0/h=1/s=2/k=3/t=4),lip-source(0/s=1/r_start=2/rl_noise=3),rl(if=9,use rl_default)'; lrpa=[1 , 0 , 0.3 , -3 , 0.05 , 2, 1, 0.1; 2 , 0 , 0.3 , -3 , 0.05 , 2, 1, 0.1; 3 , 0 , 0.3 , -3 , 0.05 , 2, 1, 0.1; 4 , 0 , 0.3 , -3 , 0.05 , 2, 1, 0.1; 5 , 0 , 0.3 , -3 , 0.05 , 2, 1, 0.1; 6 , 0 , 0.3 , -3 , 0.05 , 2, 1, 0.1; 7 , 0 , 0.3 , -3 , 0.05 , 2, 1, 0.1; 8 , 0 , 0.3 , -3 , 0.05 , 2, 1, 0.1; 9 , 0 , 0.3 , -3 , 0.05 , 2, 1, 0.1]; elseif inum == 10 then // when /au/ disp('+ a sample of /au/ is setting. '); memo1x='a sample of /au/'; D_QTY=14; T_QTY=2; memo2x=' no, l1, r1, pitch(+-%), amplitude(0-1), vocal-source(0/h=1/s=2/k=3/t=4),lip-source(0/s=1/r_start=2/rl_noise=3),rl(if=9,use rl_default)'; lrpa=[1 , 0 , 0.5 , -3 , 0.1 , 0, 0, 9; 2 , 0 , 0.55 , -1 , 0.2 , 0, 0, 9; 3 , 0 , 0.6 , 0 , 0.3 , 0, 0, 9; 4 , 0 , 0.65 , 0 , 0.4 , 0, 0, 9; 5 , 0 , 0.7 , 0 , 0.5 , 0, 0 ,9; 6 , 0 , 0.6, 0 , 0.6 , 0, 0 , 9; 7 , 0 , 0.5 , 0 , 0.7 , 0, 0, 9; 8 , 0 , 0.4 , 0 , 0.8 , 0, 0, 9; 9 , 0 , 0.3 , -2 , 0.9 , 0, 0, 9; 10 , 0 , 0.2 , -3 , 0.7 , 0, 0, 9; 11 , 0 , 0.1 , -4 , 0.5 , 0, 0, 9; 12 , 0 , 0 , -4 , 0.5 , 0, 0, 9; 13 , 0 , -0.1 , -5 , 0.4 , 0, 0, 9 ; 14 , 0 , -0.2 , -5 , 0.1 , 0, 0, 9] ; elseif inum == 9 then // when /ae/ disp('+ a sample of /ae/ is setting. '); memo1x='a sample of /ae/'; D_QTY=12; T_QTY=2; memo2x=' no, l1, r1, pitch(+-%), amplitude(0-1), vocal-source(0/h=1/s=2/k=3/t=4),lip-source(0/s=1/r_start=2/rl_noise=3),rl(if=9,use rl_default)'; lrpa=[1 , 0 , 0.5 , -3 , 0.1 , 0, 0, 9; 2 , 0 , 0.6 , -1 , 0.2 , 0, 0, 9; 3 , 0 , 0.7 , 0 , 0.3 , 0, 0, 9; 4 , 0 , 0.8 , 0 , 0.4 , 0, 0, 9; 5 , 0 , 0.8 , 0 , 0.5 , 0, 0 ,9; 6 , 0 , 0.75, 0 , 0.6 , 0, 0 , 9; 7 , -0.15 , 0.75 , -1 , 0.6 , 0, 0, 9; 8 , -0.25 , 0.75 , -2 , 0.6 , 0, 0, 9; 9 , -0.35 , 0.75 , -3 , 0.6 , 0, 0, 9; 10 , -0.45 , 0.75 , -4 , 0.6 , 0, 0, 9; 11 , -0.55 , 0.75 , -5 , 0.4 , 0, 0, 9; 12 , -0.55 , 0.75 , -5 , 0.1 , 0, 0, 9] ; elseif inum == 8 then // when /ia/ disp('+ a sample of /ya(ia)/ is setting.'); memo1x='a sample of /ya(ia)/'; D_QTY=12; T_QTY=5; memo2x=' no, l1, r1, pitch(+-%), amplitude(0-1), vocal-source(0/h=1/s=2/k=3/t=4),lip-source(0/s=1/r_start=2/rl_noise=3),rl(if=9,use rl_default)'; lrpa=[1 , -0.2 , -0.75 , -3 , 0.5 , 0, 3, 9; 2 , -0.2 , -0.75 , -3 , 0.5 , 0, 3, 9; 3 , -0.2 , -0.75 , -3 , 0.5 , 0, 3, 9; 4 , -0.2 , -0.4 , -3 , 0.5 , 0, 3, 9; 5 , -0.1 , -0.2 , -3 , 0.5 , 0, 0, 9; 6 , 0 , 0.2 , -3 , 0.5 , 0, 0, 9; 7 , 0 , 0.4 , -3 , 0.5 , 0, 0, 9; 8 , 0 , 0.6 , -3 , 0.6 , 0, 0, 9; 9 , 0 , 0.7 , -3 , 0.7 , 0, 0, 9; 10, 0 , 0.7 , -3 , 0.8 , 0, 0, 9; 11, 0 , 0.8 , -3 , 0.7 , 0, 0, 9; 12 , 0 , 0.8 , -5 , 0.1 , 0, 0, 9] ; elseif inum == 7 then // when /wa(ua)/ disp('+ a sample of /wa(ua)/ is setting.'); memo1x='a sample of /wa(ua)/'; D_QTY=12; T_QTY=2; memo2x=' no, l1, r1, pitch(+-%), amplitude(0-1), vocal-source(0/h=1/s=2/k=3/t=4),lip-source(0/s=1/r_start=2/rl_noise=3),rl(if=9,use rl_default)'; lrpa=[1 , 0 , 0. , -3 , 0.1 , 0, 0, 9; 2 , 0 , 0.2 , -1 , 0.2 , 0, 0, 9; 3 , 0 , 0.3, 0 , 0.3 , 0, 0, 9; 4 , 0 , 0.4 , 0 , 0.4 , 0, 0, 9; 5 , 0 , 0.5 , 0 , 0.5 , 0, 0 ,9; 6 , 0 , 0.6, 0 , 0.6 , 0, 0 , 9; 7 , 0 , 0.7 , 0 , 0.7 , 0, 0, 9; 8 , 0 , 0.7 , 0 , 0.8 , 0, 0, 9; 9 , 0 , 0.8 , 0 , 0.9 , 0, 0, 9; 10 , 0 , 0.8 , 0 , 0.7 , 0, 0, 9; 11 , 0 , 0.8 , -3 , 0.4 , 0, 0, 9; 12 , 0 , 0.8 , -5 , 0.1 , 0, 0, 9] ; elseif inum == 6 then // when /o/ disp('+ a sample of /o/ is setting.'); memo1x='a sample of /o/'; D_QTY=12; T_QTY=3; //This is not good resolved yet another_tube_on=1; variable_filter=1; // when variable_filter is 1, another_tube_on must be 1 ! variable_filter_hcf=1; // when variable_filter_hcf is 1, another_tube_on and variable_filter must be 1 ! memo2x=' no, l1, r1, l2, r2, pitch(+-%), amplitude(0-1), vocal-source(0/h=1/s=2/k=3/t=4),lip-source(0/s=1/r_start=2/rl_noise=3),rl(if=9,use rl_default), Rn(0. ... 1.), bef(off=0/low=1/high=2), always 4,variable filter type(off=0/lpf=1/hpf=2/bpf=3/bef=4/double_bef=5),fc(when double fc1x10000+fc2),q,gain(when bef,this value dip gain dB,ex -20(dB), hcf(off=0/on=1),fc,gain)'; lrpa=[ 1 , 0.5 , -0.6 , 0, 0.75, -3 , 0.1 , 0, 0, 9, 0., 0, 4, 0, 0, 0.7, 0., 1, 4000, -60; // use hcf to reduce pass-through direct reflection effect 2 , 0.5 , -0.6 , 0, 0.75, -1 , 0.2 , 0, 0, 9, 0., 0, 4, 0, 0, 0.7, 0., 1, 4000, -60; 3 , 0.5 , -0.6 , 0, 0.75, 0 , 0.3 , 0, 0, 9, 0., 0, 4, 0, 0, 0.7, 0., 1, 4000, -60; 4 , 0.5 , -0.6 , 0, 0.75, 0 , 0.4 , 0, 0, 9, 0., 0, 4, 0, 0, 0.7, 0., 1, 4000, -60; 5 , 0.5 , -0.6 , 0, 0.75, 0 , 0.5 , 0, 0, 9, 0., 0, 4, 0, 0, 0.7, 0., 1, 4000, -60; 6 , 0.5 , -0.6 , 0, 0.75, 0 , 0.6 , 0, 0, 9, 0., 0, 4, 0, 0, 0.7, 0., 1, 4000, -60; 7 , 0.5 , -0.6 , 0, 0.75, 0 , 0.7 , 0, 0, 9, 0., 0, 4, 0, 0, 0.7, 0., 1, 4000, -60; 8 , 0.5 , -0.6 , 0, 0.75, 0 , 0.8 , 0, 0, 9, 0., 0, 4, 0, 0, 0.7, 0., 1, 4000, -60; 9 , 0.5 , -0.6 , 0, 0.75, 0 , 0.9 , 0, 0, 9, 0., 0, 4, 0, 0, 0.7, 0., 1, 4000, -60; 10 ,0.5 , -0.6 , 0, 0.75, 0 , 0.7 , 0, 0, 9, 0., 0, 4, 0, 0, 0.7, 0., 1, 4000, -60; 11 ,0.5 , -0.6 , 0, 0.75, -3 , 0.4 , 0, 0, 9, 0., 0, 4, 0, 0, 0.7, 0., 1, 4000, -60; 12 ,0.5 , -0.6 , 0, 0.75, -5 , 0.1 , 0, 0, 9, 0., 0, 4, 0, 0, 0.7, 0., 1, 4000, -60]; //D_QTY=12; //T_QTY=3; //memo2x=' no, l1, r1, l2, r2, pitch(+-%), amplitude(0-1), vocal-source(0/h=1/s=2/k=3/t=4),lip-source(0/s=1/r_start=2/rl_noise=3),rl(if=9,use rl_default)'; // lrpa=[1 , 0.4 , -0.7 , 0, 0.75, -3 , 0.1 , 0, 0, 9; // 2 , 0.4 , -0.7 , 0, 0.75 , -1 , 0.2 , 0, 0, 9; // 3 , 0.4 , -0.7 , 0, 0.75 , 0 , 0.3 , 0, 0, 9; // 4 , 0.4 , -0.7 , 0, 0.75 , 0 , 0.4 , 0, 0, 9; // 5 , 0.4 , -0.7 , 0, 0.75 , 0 , 0.5 , 0, 0 ,9; // 6 , 0.4 , -0.7 , 0, 0.75, 0 , 0.6 , 0, 0 , 9; // 7 , 0.4 , -0.7 , 0, 0.75 , 0 , 0.7 , 0, 0, 9; // 8 , 0.4 , -0.7 , 0, 0.75 , 0 , 0.8 , 0, 0, 9; // 9 , 0.4 , -0.7 , 0, 0.75 , 0 , 0.9 , 0, 0, 9; // 10 , 0.4 , -0.7 , 0, 0.75 , 0 , 0.7 , 0, 0, 9; // 11 , 0.4 , -0.7 , 0, 0.75 , -3 , 0.4 , 0, 0, 9; // 12 , 0.4 , -0.7 , 0, 0.75, -5 , 0.1 , 0, 0, 9] ; elseif inum == 5 then // when /e(after a)/ disp('+ a sample of /e(after a)/ is setting. '); memo1x='a sample of /e(after a)/'; D_QTY=12; T_QTY=2; memo2x=' no, l1, r1, pitch(+-%), amplitude(0-1), vocal-source(0/h=1/s=2/k=3/t=4),lip-source(0/s=1/r_start=2/rl_noise=3),rl(if=9,use rl_default)'; lrpa=[1 , 0.55 , 0.75 , -3 , 0.1 , 0, 0, 9; 2 , 0.55 , 0.75 , -1 , 0.2 , 0, 0, 9; 3 , 0.55 , 0.75 , 0 , 0.3 , 0, 0, 9; 4 , 0.55 , 0.75 , 0 , 0.4 , 0, 0, 9; 5 , 0.55 , 0.75 , 0 , 0.5 , 0, 0 ,9; 6 , 0.55 , 0.75, 0 , 0.6 , 0, 0 , 9; 7 , 0.55 , 0.75 , 0 , 0.7 , 0, 0, 9; 8 , 0.55 , 0.75 , 0 , 0.8 , 0, 0, 9; 9 , 0.55 , 0.75 , 0 , 0.9 , 0, 0, 9; 10 , 0.55 , 0.75 , 0 , 0.7 , 0, 0, 9; 11 , 0.55 , 0.75 , -3 , 0.4 , 0, 0, 9; 12 , 0.55 , 0.75 , -5 , 0.1 , 0, 0, 9] ; // elseif inum == 4 then // when /e(before a)/ disp('+ a sample of /e(before a)/ is setting. '); memo1x='a sample of /e(before a)/'; D_QTY=12; T_QTY=2; memo2x=' no, l1, r1, pitch(+-%), amplitude(0-1), vocal-source(0/h=1/s=2/k=3/t=4),lip-source(0/s=1/r_start=2/rl_noise=3),rl(if=9,use rl_default)'; lrpa=[1 , -0.55 , 0.75 , -3 , 0.1 , 0, 0, 9; 2 , -0.55 , 0.75 , -1 , 0.2 , 0, 0, 9; 3 , -0.55 , 0.75 , 0 , 0.3 , 0, 0, 9; 4 , -0.55 , 0.75 , 0 , 0.4 , 0, 0, 9; 5 , -0.55 , 0.75 , 0 , 0.5 , 0, 0 ,9; 6 , -0.55 , 0.75, 0 , 0.6 , 0, 0 , 9; 7 , -0.55 , 0.75 , 0 , 0.7 , 0, 0, 9; 8 , -0.55 , 0.75 , 0 , 0.8 , 0, 0, 9; 9 , -0.55 , 0.75 , 0 , 0.9 , 0, 0, 9; 10 , -0.55 , 0.75 , 0 , 0.7 , 0, 0, 9; 11 , -0.55 , 0.75 , -3 , 0.4 , 0, 0, 9; 12 , -0.55 , 0.75 , -5 , 0.1 , 0, 0, 9] ; // elseif inum == 3 then // when /u/ disp('+ a sample of /u/ is setting. '); memo1x='a sample of /u/'; D_QTY=12; T_QTY=2; memo2x=' no, l1, r1, pitch(+-%), amplitude(0-1), vocal-source(0/h=1/s=2/k=3/t=4),lip-source(0/s=1/r_start=2/rl_noise=3),rl(if=9,use rl_default)'; lrpa=[1 , 0 , 0 , -3 , 0.1 , 0, 0, 9; 2 , 0 , 0 , -1 , 0.2 , 0, 0, 9; 3 , 0 , 0 , 0 , 0.3 , 0, 0, 9; 4 , 0 , 0 , 0 , 0.4 , 0, 0, 9; 5 , 0 , 0 , 0 , 0.5 , 0, 0 ,9; 6 , 0 , 0, 0 , 0.6 , 0, 0 , 9; 7 , 0 , 0 , 0 , 0.7 , 0, 0, 9; 8 , 0 , 0 , 0 , 0.8 , 0, 0, 9; 9 , 0 , 0 , 0 , 0.9 , 0, 0, 9; 10 , 0 , 0 , 0 , 0.7 , 0, 0, 9; 11 , 0 , 0 , -3 , 0.4 , 0, 0, 9; 12 , 0 , 0 , -5 , 0.1 , 0, 0, 9] ; elseif inum == 2 then // when /i/ disp('+ a sample of /i/ is setting. '); memo1x='a sample of /i/'; D_QTY=12; T_QTY=5; memo2x=' no, l1, r1, pitch(+-%), amplitude(0-1), vocal-source(0/h=1/s=2/k=3/t=4),lip-source(0/s=1/r_start=2/rl_noise=3),rl(if=9,use rl_default)'; lrpa=[1 , -0.2 , -0.75 , -3 , 0.1 , 0, 3, 9; 2 , -0.2 , -0.75 , -3 , 0.2 , 0, 3, 9; 3 , -0.2 , -0.75 , -3 , 0.3 , 0, 3, 9; 4 , -0.2 , -0.75 , -3 , 0.4 , 0, 3, 9; 5 , -0.2 , -0.75 , -3 , 0.5 , 0, 3, 9; // as this detect model, for some pitches, amplitude is stable is better 6 , -0.2 , -0.75 , -3 , 0.5 , 0, 3, 9; 7 , -0.2 , -0.75 , -3 , 0.5 , 0, 3, 9; 8 , -0.2 , -0.75 , -3 , 0.5 , 0, 3, 9; 9 , -0.2 , -0.75 , -3 , 0.5 , 0, 3, 9; 10, -0.2 , -0.75 , -3 , 0.5 , 0, 3, 9; 11, -0.2 , -0.75 , -3 , 0.4 , 0, 3, 9; 12 ,-0.2 , -0.75 , -5 , 0.1 , 0, 3, 9] ; elseif inum == 1 then // when /a/ disp('+ a sample of /a/ is setting. '); memo1x='a sample of /a/'; D_QTY=12; T_QTY=2; memo2x=' no, l1, r1, pitch(+-%), amplitude(0-1), vocal-source(0/h=1/s=2/k=3/t=4),lip-source(0/s=1/r_start=2/rl_noise=3),rl(if=9,use rl_default)'; lrpa=[1 , 0 , 0.5 , -3 , 0.1 , 0, 0, 9; 2 , 0 , 0.55 , -1 , 0.2 , 0, 0, 9; 3 , 0 , 0.6 , 0 , 0.3 , 0, 0, 9; 4 , 0 , 0.65 , 0 , 0.4 , 0, 0, 9; 5 , 0 , 0.7 , 0 , 0.5 , 0, 0 ,9; 6 , 0 , 0.75, 0 , 0.6 , 0, 0 , 9; 7 , 0 , 0.75 , 0 , 0.7 , 0, 0, 9; 8 , 0 , 0.8 , 0 , 0.8 , 0, 0, 9; 9 , 0 , 0.8 , 0 , 0.9 , 0, 0, 9; 10 , 0 , 0.8 , 0 , 0.7 , 0, 0, 9; 11 , 0 , 0.8 , -3 , 0.4 , 0, 0, 9; 12 , 0 , 0.8 , -5 , 0.1 , 0, 0, 9] ; // end // commoon co_filter=[ 1 , 0, 500, 0.7 ; 2 , 1, 1000, 0.7 ; 3, 0, 10000, 0.7 ; 4, 2, 5000, 2.5 ; 5, 0, 1500,0.7 ; 6 , 2, 1250, 10.; 7 , 1, 1000, 0.7 ; 8, 0, 10000, 0.7 ; 9, 2, 2500, 10 ]; memo3x='no, type(LPF=0,BPF=1,HPF=2), f0, Q'; k_gain0=[7.0 , 0.0]; // check s0=size(lrpa); if D_QTY > s0(1) then disp('ERROR: D_QTY > size of lrpa: force to set D_QTY as size of lrpa. At function sub_demo1'); D_QTY = s0(1); elseif D_QTY < s0(1) then disp('WARNING: D_QTY < size of lrpa. At function sub_demo1'); end endfunction // for page 2 //------------------------------------------------------------------------- // // for page 3 // function [D_QTY,T_QTY,lrpa,co_filter,k_gain0,memo1x,memo2x,memo3x,variable_filter,another_tube_on,variable_filter_hcf]=sub_demo1_page3( inum0 ) // setting for demonstration inum=inum0 - select_offset_page3; another_tube_on=0; variable_filter=0; variable_filter_hcf=0; // common 1 co_filter=[ 1 , 0, 500, 0.7 ; 2 , 1, 1000, 0.7 ; 3, 0, 10000, 0.7 ; 4, 2, 5000, 2.5 ; 5, 0, 1500,0.7 ; 6 , 2, 1250, 10.; 7 , 1, 1000, 0.7 ; 8, 0, 10000, 0.7 ; 9, 2, 2500, 10 ]; if inum == 4 then // when /pa/ disp('+ a sample of /pa/ is setting. '); memo1x='a sample of /pa/'; D_QTY=12; T_QTY=2; memo2x=' no, l1, r1, pitch(+-%), amplitude(0-1), vocal-source(0/h=1/s=2/k=3/t=4),lip-source(0/s=1/r_start=2/rl_noise=3),rl(if=9,use rl_default)'; lrpa=[1 , 0 , 0.7 , 0 , 0.9 , 0, 0, 0.01; // strong breath input 2 , 0 , 0.7 , 0 , 0.4 , 0, 0, 0.2; // and then resonate little 3 , 0 , 0.7 , 0 , 0.4 , 0, 0, 0.5; 4 , 0 , 0.7 , 0 , 0.4 , 0, 0, 9; 5 , 0 , 0.7 , 0 , 0.5 , 0, 0 ,9; 6 , 0 , 0.75, 0 , 0.6 , 0, 0 , 9; 7 , 0 , 0.75 , 0 , 0.7 , 0, 0, 9; 8 , 0 , 0.8 , 0 , 0.8 , 0, 0, 9; 9 , 0 , 0.8 , 0 , 0.9 , 0, 0, 9; 10 , 0 , 0.8 , 0 , 0.7 , 0, 0, 9; 11 , 0 , 0.8 , -3 , 0.4 , 0, 0, 9; 12 , 0 , 0.8 , -5 , 0.1 , 0, 0, 9] ; // elseif inum == 2 then // when /ga/ disp('+ a sample of /ga/ is setting. '); memo1x='a sample of /ga/'; D_QTY=10; T_QTY=2; memo2x=' no, l1, r1, pitch(+-%), amplitude(0-1), vocal-source(0/h=1/s=2/k=3/t=4),lip-source(0/s=1/r_start=2/rl_noise=3),rl(if=9,use rl_default)'; lrpa=[1 , 0 , -0.3 , 0 , 0.05 , 3, 0, 0.9; // needs resonance effect, rl //2 , 0 , 0 , 0 , 0.02 , 3, 0, 0.1; // erase this line cause /t/ sound 2 , 0 , 0.3 , 0 , 0.3 , 0, 0, 0.5; 3 , 0 , 0.4 , 0 , 0.4 , 0, 0, 0.7; 4 , 0 , 0.5 , 0 , 0.5 , 0, 0, 9; 5 , 0 , 0.6, 0 , 0.6 , 0, 0, 9; 6 , 0 , 0.7 , 0 , 0.7 , 0, 0, 9; 7 , 0 , 0.8 , 0 , 0.8 , 0, 0, 9; 8 , 0 , 0.8 , 0 , 0.9 , 0, 0, 9; 9 , 0 , 0.8 , 0 , 0.7 , 0, 0, 9; 10 , 0 , 0.8 , 0 , 0.7 , 0, 0, 9]; elseif inum == 3 then // when /da/ // ikinari hatuon suruno ga ta disp('+ a sample of /da/ is setting. '); memo1x='a sample of /da/'; D_QTY=11; T_QTY=2; memo2x=' no, l1, r1, pitch(+-%), amplitude(0-1), vocal-source(0/h=1/s=2/k=3/t=4),lip-source(0/s=1/r_start=2/rl_noise=3),rl(if=9,use rl_default)'; lrpa=[1 , 0 , 0 , 0 , 0.05 , 4, 0, 0.1; 2 , 0 , 0 , 0 , 0.02 , 4, 0, 0.1; 3 , 0 , 0.3 , 0 , 0.3 , 0, 0, 0.5; 4 , 0 , 0.4 , 0 , 0.4 , 0, 0, 0.7; 5 , 0 , 0.5 , 0 , 0.5 , 0, 0, 9; 6 , 0 , 0.6, 0 , 0.6 , 0, 0, 9; 7 , 0 , 0.7 , 0 , 0.7 , 0, 0, 9; 8 , 0 , 0.8 , 0 , 0.8 , 0, 0, 9; 9 , 0 , 0.8 , 0 , 0.9 , 0, 0, 9; 10 , 0 , 0.8 , 0 , 0.7 , 0, 0, 9; 11 , 0 , 0.8 , 0 , 0.7 , 0, 0, 9]; co_filter(5,3)=2000; // if 2000, then it may be da //elseif inum == 1 then // when /a/ else disp('+ a sample of /a/ is setting. '); memo1x='a sample of /a/'; D_QTY=12; T_QTY=2; memo2x=' no, l1, r1, pitch(+-%), amplitude(0-1), vocal-source(0/h=1/s=2/k=3/t=4),lip-source(0/s=1/r_start=2/rl_noise=3),rl(if=9,use rl_default)'; lrpa=[1 , 0 , 0.5 , -3 , 0.1 , 0, 0, 9; 2 , 0 , 0.55 , -1 , 0.2 , 0, 0, 9; 3 , 0 , 0.6 , 0 , 0.3 , 0, 0, 9; 4 , 0 , 0.65 , 0 , 0.4 , 0, 0, 9; 5 , 0 , 0.7 , 0 , 0.5 , 0, 0 ,9; 6 , 0 , 0.75, 0 , 0.6 , 0, 0 , 9; 7 , 0 , 0.75 , 0 , 0.7 , 0, 0, 9; 8 , 0 , 0.8 , 0 , 0.8 , 0, 0, 9; 9 , 0 , 0.8 , 0 , 0.9 , 0, 0, 9; 10 , 0 , 0.8 , 0 , 0.7 , 0, 0, 9; 11 , 0 , 0.8 , -3 , 0.4 , 0, 0, 9; 12 , 0 , 0.8 , -5 , 0.1 , 0, 0, 9] ; // end // common 2 memo3x='no, type(LPF=0,BPF=1,HPF=2), f0, Q'; k_gain0=[7.0 , 0.0]; // check s0=size(lrpa); if D_QTY > s0(1) then disp('ERROR: D_QTY > size of lrpa: force to set D_QTY as size of lrpa. At function sub_demo1'); D_QTY = s0(1); elseif D_QTY < s0(1) then disp('WARNING: D_QTY < size of lrpa. At function sub_demo1'); end endfunction // for page 3 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ // proccess step0 if xmode2==0 then if T_DEMO == 1 then [D_QTY,T_QTY,lrpa,co_filter,k_gain0,memo1,memo2,memo3,variable_filter_on,another_tube_on,variable_filter_on_hcf]=sub_demo1( sel_demo1 ); [l1_nose,r1_nose,rl_nose,n_supend_para0]=set_nose(sel_demo1); [x_coef1,y_coef1]=set_h_iir_et_plot_bode(0); [x_coef2,y_coef2,x_coef3,y_coef3,x_coef4,y_coef4]=set_s_iir_et_plot_bode(0); [x_coef5,y_coef5,x_coef6,y_coef6]=set_k_iir_et_plot_bode(0); [x_coef9,y_coef9,x_coef10,y_coef10,x_coef11,y_coef11]=set_r_iir_et_plot_bode(0); [x_coef7,y_coef7,x_coef8,y_coef8]=set_N_iir_et_plot_bode(0); [A1,A2,A3,L1,L2,L3,A1_n,A2_n,L1_n,L2_n,tclosed,trise,tfall,amplitude]=make_paras(); else [D_QTY,T_QTY,lrpa,memo1,co_filter]= set_QTY(); [lrpar,memo2]=edit_lrpa(); //[co_filter,memo3]=edit_co_filter(); // old [A1,A2,A3,L1,L2,L3,A1_n,A2_n,L1_n,L2_n,tclosed,trise,tfall,amplitude]=make_paras(); end end // if xmode2==0 then // //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ // proccess step1 to step 4 STEP_DO0=1; // = 0 or 1 or2 [rg,ygn,yg,rl,N1,N2,N3,LL,length0]=step1(); [xi_var,yi_var,resp_vari,xi_var_hcf,yi_var_hcf,resp_vari_hcf]=set_variable_coeff(); if STEP_DO0 >= 1 then [tu1,tu2,tu3,r1,r2,r12,r21,r31,tu1n,tu2n,r1n]=make_tu_r(); [y2tm,M1,M2,M3]=step2(); [y2tm,y2tm_lpf,y_ran2,al1sm,bl0sm,bl1sm]=step2_2(); [y2tm]=after_step2(); [y2tmn,bl0smb,bl1smb]=step2_n(); [y2tm]=make_sin_y2tm(); [y2tm]=step2_variable_filter(); [y3out,y3outn,ah1,bh0,bh1,ah1b,bh0b,bh1b]=step3(); step4(); end // if (T_DEMO == 1 & STEP_DO0 >= 2) then [db_all,phi_all]=make_bode(); eh_var1=1; eh_var2=1; eh_var3=3; w0x=xget('window'); xset('window', eh_var1); seteventhandler('my_eventhandler'); plot_areas(eh_var1, eh_var2, eh_var3); xset('window', w0x); // disp('+usage: key in d (=down scroll) or u (=up scroll) on Scilab Graphic (1) windows. '); end // //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ //------------------------------------------------------------------------- //================================================================================ // // other functions // //yg_bode(1) // plot yg (Glottal Volume Velocity) frequency-phase response //two_tubes_model_bode(1) // plot Two Tubes Model's frequency-phase response //three_tubes_model_bode(1) // plot Three Tubes Model's frequency-phase response //hpf_bode() // plot high pass filter's frequency-phase response //save_as_sound_file_yg( 'yg.wav' ) // save yg (Glottal Volume Velocity) as one .wav file //save_as_sound_file( 'y3out.wav' ) // save y3out (Sound Pressure) as one .wav file //plot_area(1) // plot sqrt(area) vs length //plot_yg_yout() // plot yg (Glottal Volume Velocity) and y3out (Sound Pressure) //fft_plot( 1024 ) // FFT y3out(Sound Pressure) (1024 points) and plot frequency-phase // //two_tubes_rl_noise_bode(1) // plot Two Tubes Model with "rl" noise frequency-phase response //fft_plot_ran2( 1024 ) // FFT y_ran2(noise with limited frequency band) and plot frequency-phase //plot_yg_y2tm_y2tm_lpf_yo() // plot yg (Glottal Volume Velocity) y2tm(Volume Velocity) y2tm_lpf(smoothed Volume Velocity) and y3out (Sound Pressure) //lpf_bode() // plot low pass (smoothing) filter's frequency-phase response // // //----------------------------------------------------------------- // sound setting // defch1=1; // default channel 1 f_win=1; // flag if windows bit1=16; fs1=fs; //--- check platform is windows ----------------------------------- if isdir('c:\\') then f_win=1; else f_win=0; end //---------------------------------------------------------------------- // //----------------------------------------------------------------- function [wavfile1,chs1,qty1,fs1,bit1,f412]=get_wavfile_pro() // choose wave file if scilab_version_number >= 5 then wavfile1=uigetfile(["*.wav"],"","Choose any xxx.wav file name"); else wavfile1=tk_getfile(Title="Choose any xxx.wav file name"); end // 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 //... [y,fs1,bit1]=wavread(wavfile1,[1 1]); endfunction //----------------------------------------------------------------- function snd_play2(wdat2) // 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! // for wdat2 sz=size(wdat2); size2=sz(2); if f_win == 1 then if f412 == 1 then // for windows scilab-4.1.2 sound(wdat2 ,fs1,bit1); else // for windows scilab-3.1.1 if f511== 1 then // for windows scilab-5.1.1 sound(wdat2',fs1,bit1); else if fs1 == 22050 then sound(wdat2',fs1,bit1); disp('This function may work, if luckily.'); elseif fs1 == 44100 then // down-sampling from 44100 to 22050 wdatw2=zeros(1,(int(size2/2))); for v=1:(int(size2/2)) wdatw2(v)=wdat2(2 * v ); end disp('down-sampling from 44100 to 22050'); sound(wdatw2',22050,bit1); disp('This function may work, if luckily.'); else //sound(wdatw,fs1,bit1); disp('This function does not work.'); end end // for windows scilab-5.1.1 end //// for windows scilab-4.1.2 else disp('If sound setting is good for scilab, this function maybe work.'); if f412 == 1 then // for scilab-4.1.2 sound(wdat2 ,fs1,bit1); else // for scilab-3.1.1 if f511== 1 then // for windows scilab-5.1.1 sound(wdat2',fs1,bit1); else if fs1 == 22050 then sound(wdat2',fs1,bit1); disp('This function may work, if luckily.'); elseif fs1 == 44100 then // down-sampling from 44100 to 22050 wdatw2=zeros(1,(int(size2/2))); for v=1:(int(size2/2)) wdatw2(v)=wdat2(2 * v ); end disp('down-sampling from 44100 to 22050'); sound(wdatw2',22050,bit1); disp('This function may work, if luckily.'); else //sound(wdatw,fs1,bit1); disp('This function does not work.'); end end end end endfunction // //---------------------------------------------------------------- function snd_save2(wdat2) // for wdat2 // wavefilename= input(' + enter file name for saved .wav file =>',["string"]); // if f412 == 1 then // for windows scilab-4.1.2 wavwrite(wdat2 ,fs1,bit1,wavefilename ); else // for windows scilab-3.1.1 wavwrite(wdat2',fs1,bit1, wavefilename); end endfunction //---------------------------------------------------------------- function save_text1() // save lrpa as a text file u=file('open','lrpa.txt','unknown'); wstr1=date(); dt=getdate(); wstr1=wstr1 + ' ' + string(dt(7)) + ':' + string(dt(8)) + ':' +string(dt(9)); fprintf(u,'// %s',wstr1); fprintf(u,'D_QTY=%d;',D_QTY); fprintf(u,'T_QTY=%d;',T_QTY); fprintf(u,'memo1=''%s'';',memo1); fprintf(u,'// %s',memo2); s0=size(lrpa); for v=1:s0(1) if v == 1 then wstr1='lrpa=[ '; for w=1:s0(2) if w== s0(2) then wstr1 = wstr1 + string( lrpa(v,w) ) + ' ; '; else wstr1 = wstr1 + string( lrpa(v,w) ) + ' , '; end end fprintf(u,'%s',wstr1); elseif v == s0(1) then wstr1=' '; for w=1:s0(2) if w== s0(2) then wstr1 = wstr1 + string( lrpa(v,w) ) + ' ] ;'; else wstr1 = wstr1 + string( lrpa(v,w) ) + ' , '; end end fprintf(u,'%s',wstr1); else wstr1=' '; for w=1:s0(2) if w== s0(2) then wstr1 = wstr1 + string( lrpa(v,w) ) + ' ; '; else wstr1 = wstr1 + string( lrpa(v,w) ) + ' , '; end end fprintf(u,'%s',wstr1); end end // co_filter fprintf(u,'// %s',memo3); s0=size(co_filter); for v=1:s0(1) if v == 1 then wstr1='co_filter=[ '; for w=1:s0(2) if w== s0(2) then wstr1 = wstr1 + string( co_filter(v,w) ) + ' ; '; else wstr1 = wstr1 + string( co_filter(v,w) ) + ' , '; end end fprintf(u,'%s',wstr1); elseif v == s0(1) then wstr1=' '; for w=1:s0(2) if w== s0(2) then wstr1 = wstr1 + string( co_filter(v,w) ) + ' ] ;'; else wstr1 = wstr1 + string( co_filter(v,w) ) + ' , '; end end fprintf(u,'%s',wstr1); else wstr1=' '; for w=1:s0(2) if w== s0(2) then wstr1 = wstr1 + string( co_filter(v,w) ) + ' ; '; else wstr1 = wstr1 + string( co_filter(v,w) ) + ' , '; end end fprintf(u,'%s',wstr1); end end // file('close',u); endfunction // //---------------------------------------------------------------- // // ADD_MENU_EDIT=1; // set ADD_MENU 1 to add menu button of some functions // // if ADD_MENU_EDIT == 1 then delmenu('step1_edit'); addmenu('step1_edit',['A (0)select and set','(1)set duration','(2)edit lrpa','(3-1)edit noise source filter for /h/','(3-2)edit noise source filter for /s/','(3-3)edit noise source filter for /k/ and /t/','(3-4)edit noise source filter for /r/ start mark','(3-5)edit bef filter for another tubes','B (4)make parameters','(-)save lrpa as a file','(-)load a lrpa file','(--)save lrpa as a text file as lrpa.txt']); step1_edit=['[xmode2,sel_demo1]=set_select(); [D_QTY,T_QTY,lrpa,co_filter,k_gain0,memo1,memo2,memo3,variable_filter_on,another_tube_on,variable_filter_on_hcf]=sub_demo1( sel_demo1 );[l1_nose,r1_nose,rl_nose,n_supend_para0]=set_nose(sel_demo1);','[D_QTY,T_QTY,lrpa,memo1]= set_QTY()','[lrpa,memo2]=edit_lrpa()','[co_filter,memo3]=edit_h_co_filter(); [x_coef1,y_coef1]=set_h_iir_et_plot_bode(2);','[co_filter,memo3]=edit_s_co_filter(); [x_coef2,y_coef2,x_coef3,y_coef3,x_coef4,y_coef4]=set_s_iir_et_plot_bode(2);','[co_filter,memo3,k_gain0]=edit_k_co_filter(); [x_coef5,y_coef5,x_coef6,y_coef6]=set_k_iir_et_plot_bode(2);','[co_filter,memo3]=edit_r_co_filter(); [x_coef9,y_coef9,x_coef10,y_coef10,x_coef11,y_coef11]=set_r_iir_et_plot_bode(2);','[n_supend_para0,memo4]=edit_n_supend_para(); [x_coef7,y_coef7,x_coef8,y_coef8]=set_N_iir_et_plot_bode(2);','[A1,A2,A3,L1,L2,L3,A1_n,A2_n,L1_n,L2_n,tclosed,trise,tfall,amplitude]=make_paras()','save_lrpa()','[D_QTY,T_QTY,lrpa,memo1,co_filter,k_gain0,n_supend_para0]=load_lrpa()','save_text1()']; end //if ADD_MENU == 1 then/ // //---------------------------------------------------------------- // // ADD_MENU_GEN=1; // set ADD_MENU 1 to add menu button of some functions // // if ADD_MENU_GEN == 1 then delmenu('step2_wave_generate'); addmenu('step2_wave_generate',['C (1)wave generate','(-)plot portion','(--)set fft points ']); step2_wave_generate=['[rg,ygn,yg,rl,N1,N2,N3,LL,length0]=step1(); [xi_var,yi_var,resp_vari,xi_var_hcf,yi_var_hcf,resp_vari_hcf]=set_variable_coeff(); [tu1,tu2,tu3,r1,r2,r12,r21,r31,tu1n,tu2n,r1n]=make_tu_r(); [y2tm,M1,M2,M3]=step2();[y2tm,y2tm_lpf,y_ran2,al1sm,bl0sm,bl1sm]=step2_2(); [y2tm]=after_step2(); [y2tmn,bl0smb,bl1smb]=step2_n(); [y2tm]=make_sin_y2tm(); [y2tm]=step2_variable_filter(); [y3out,y3outn,ah1,bh0,bh1,ah1b,bh0b,bh1b]=step3();step4();disp('' - done. '');' ,'[sp1, ep1]= set_sp1_ep1(y3out); display_portion(2,y3out)' ,'[fft_point1]=set_fft_points()']; end //if ADD_MENU == 1 then/ //----------------------------------------------------------------- // ADD_MENU_EH=1; // set ADD_MENU 1 to add menu button of some functions // // if ADD_MENU_EH == 1 then delmenu('step3_check_area'); addmenu('step3_check_area',['(1)preparation bode','(2) set event handler window(1), Mouse left click(=down) or Mouse right click(=up) / key in d(=down) or u(=up)','(e) suppress handler']); step3_check_area=[ '[db_all,phi_all]=make_bode();', 'eh_var1=1; eh_var2=1; eh_var3=3; w0x=xget(''window''); xset(''window'', eh_var1); seteventhandler(''my_eventhandler''); plot_areas(eh_var1, eh_var2, eh_var3); xset(''window'', w0x);' , 'seteventhandler('''') ' ]; end // if ADD_MENU_EH == 1 then // //----------------------------------------------------------------- // ADD_MENU_PLAY=1; // set ADD_MENU 1 to add menu button of some functions // // if ADD_MENU_PLAY == 1 then delmenu('menu_sound'); addmenu('menu_sound',['(1)check wav version','D (2)play sound','(3)save sound']); menu_sound=['[xwavfile1,xchs1,xqty1,xfs1,xbit1,f412]=get_wavfile_pro()','snd_play2(y3out)','snd_save2(y3out)']; end //if ADD_MENU == 1 then/ // //=========================================================================================== // // Add menu button of which name is 'functions' // ADD_MENU=1; // set ADD_MENU 1 to add menu button of some functions // // if ADD_MENU == 1 then delmenu('functions'); if T_QTY == 4 then addmenu('functions',['three_tubes_model_bode(1)'; 'fft_plot(1024)'; 'save_as_sound_file(''y3out.wav'')'; 'plot_area(1) '; 'plot_yg_yout()' ; 'plot_ygn_y3outn()'] ); functions=['three_tubes_model_bode(1)'; 'fft_plot(1024)'; 'save_as_sound_file(''y3out.wav'')'; 'plot_area(1) '; 'plot_yg_yout()' ; 'plot_ygn_y3outn()'] ; elseif T_QTY == 3 then addmenu('functions',['three_tubes_model_bode(1)'; 'fft_plot(1024)'; 'save_as_sound_file(''y3out.wav'')'; 'plot_area(1) '; 'plot_yg_yout()' ; 'plot_ygn_y3outn()'] ); functions=['three_tubes_model_bode(1)'; 'fft_plot(1024)'; 'save_as_sound_file(''y3out.wav'')'; 'plot_area(1) '; 'plot_yg_yout()' ; 'plot_ygn_y3outn()'] ; elseif T_QTY == 5 then addmenu('functions',['two_tubes_model_bode(1)'; 'fft_plot(1024)'; 'save_as_sound_file(''y3out.wav'')'; 'plot_area(1) '; 'plot_yg_yout()'; 'two_tubes_rl_noise_bode(1)'; 'fft_plot_ran2(1024)'; 'plot_yg_y2tm_y2tm_lpf_yo()'; 'plot_ygn_y3outn()' ]); functions=['two_tubes_model_bode(1)'; 'fft_plot(1024)'; 'save_as_sound_file(''y3out.wav'')'; 'plot_area(1) '; 'plot_yg_yout()' ; 'two_tubes_rl_noise_bode(1)'; 'fft_plot_ran2(1024)'; 'plot_yg_y2tm_y2tm_lpf_yo()'; 'plot_ygn_y3outn()'] ; else addmenu('functions',['two_tubes_model_bode(1)'; 'fft_plot(1024)'; 'save_as_sound_file(''y3out.wav'')'; 'plot_area(1) '; 'plot_yg_yout()' ; 'save_yg(''ygout.wav'')'; 'plot_ygn_y3outn()' ; 'snd_play2(y3outn)'; 'snd_save2(y3outn)'] ); functions=['two_tubes_model_bode(1)'; 'fft_plot(1024)'; 'save_as_sound_file(''y3out.wav'')'; 'plot_area(1) '; 'plot_yg_yout()' ; 'save_as_sound_file_yg(''ygout.wav'')' ; 'plot_ygn_y3outn()' ; 'snd_play2(y3outn)'; 'snd_save2(y3outn)'] ; end end //if ADD_MENU == 1 then/ // //=========================================================================================== // ADD_MENU_STATUS=1; if ADD_MENU_STATUS == 1 then function status1() ver0= getversion(); wstr1=sprintf('scilab version is %s',ver0); disp(wstr1); endfunction function status2() sz=stacksize(); gsz=gstacksize(); wstr1=sprintf('current stack size is %d of %d ,(%f)', sz(2),sz(1), sz(2)/sz(1) ); wstr2=sprintf('current global stack size is %d of %d ,(%f)', gsz(2),gsz(1), gsz(2)/gsz(1)); disp( wstr1); disp( wstr2); endfunction delmenu('status_'); addmenu('status_',[ 'scilab version' ; 'check stack memory size' ]); status_=[ 'status1()' ; 'status2()' ] ; end //if ADD_MENU_STATUS == 1 then // //-------------------------------------------------------------------- // ADD_MENU_DEBUG = 1; if ADD_MENU_DEBUG == 1 then delmenu('for_debug'); addmenu('for_debug',[ 'test 1' ; 'test 2' ]); for_debug=[ '[x1,y1,x2,y2,x3,y3]=set_N_iir_et_plot_bode2(1);' ; '[x1,y1]=set_lpf_cheb_plot_bode(1,4, 2000, 1, -60,-20,-10); ' ] ; end //if ADD_MENU_DEBUG == 1 then // //----------------------------------------------------------------------------------------------------//