//----------------------------------------------------------------------------------------- // A .wav Edit and Compare with 2 tubes model's wave or 3 tubes model's wave. // // a trial of leastsq method to estimate tube model parameters including rl // Using focus weight to compare (local)part of object, evaluation by siguma |wmf(i) *(x(i)-y(i))|^2 between spectrum // And hoping into limit by addition of (1 + exp( fact0 * x0)) * (1 + exp(fact0 * x1)) // // additional noise by iir bpf and independently add it with 2 tube model output // // superpose (tyouzyou) model fitting : for example "i" 2 tubes model plus independent noise source // // "da" : tube model plus its partly emphasize // // "ka" : "k(uncontinuationly burst noise-in-Tube)+hu(continuationly noise-in-Tube)+a" // // // for window scilab-4.1.2 // // ......................................................................................... // PLEASE PAY ATTENSION that this program may have some bugs and also you may adjust program // to fit your circumstances. If you use this program, everything should be done at your own risk. // Thanks. //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ T_DEMO=1; // set T_DEMO=1 for demonstration, beside set T_DEMO=0 when normal T_SKIP1=1; // some program will be skipped when T_SKIP1=1 SEL_CODE=1; // dummy set // //-------------------------------------------------------------------- // 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 [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 //---------------------------------------------------------------------- // 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 //--------------------------------------------------------------------- // get Scilab current directory getcwd //==================================================================== // scilab's global variable global yg_res1 global tube2_res1 L1 L2 L3 A1 A2 A3 global hpf_res1 global frq_center q0 global yd_bpf xd_bpf global yd_lpf xd_lpf global yd_hpf xd_hpf global ind_noise_res1 global stable_noise_res1 global burst_noise_res1 global yd_peq xd_peq global peq_res1 global overall_res1 db_2tube phi_2tube global r1 r2 l1 l2 ttl_Length rl noise_waku_area i2nd_thd_factor nA0 global sp1 ep1 tframe fir_steps global db_fft1 phi_fft1 global yg_res9 hpf_res9 global x_init8 wm9 global xopts global r1s r2s l1s l2s ttl_Lengths rls noise_waku_areas i2nd_thd_factors nA0s global emp_space_area empA0 global emp_space_areas empA0s global stable_noise_lpf stable_noise_lpfs global WT_QTY global fact0 global fact1_2r1 fact2_2r1 fact1_2l1 fact2_2l1 fact1_rl fact2_rl limit_switch0 limit_switch2 limit_switch3 global fact1_3l1 fact2_3l1 fact1_3r1 fact2_3r1 fact1_3l2 fact2_3l2 fact1_3r2 fact2_3r2 global fact1_2r1s fact2_2r1s fact1_2l1s fact2_2l1s fact1_rls fact2_rls limit_switch0s limit_switch2s limit_switch3s global fact1_3l1s fact2_3l1s fact1_3r1s fact2_3r1s fact1_3l2s fact2_3l2s fact1_3r2s fact2_3r2s global max_nth max_value max_wm9 global mess_init0 global yg y2tm y3out y2tm_lpf y2tm_noise LL y3out_least global st_freq end_freq n_stepf end_freq2 global bpf_fb0 wft fb_out fb_am_out fb_am_out_prop global fb_out_init fb_am_out_init fb_am_out_prop_init global fb_out_work1 fb_am_out_work1 fb_am_out_prop_work1 global wdat1 global LOAD_FB_FLAG //------------------------------------------------------------------- LOAD_FB_FLAG=0; // this is optional function. to load fb data which was already calculated. //=================================================================== // global variables No. 1 // for .wav 1 chs1=0; // channels qty1=0; // data quantity, samples size1=0; // actual loaded data quantity fs1=0; // sampling frequency bit1=0; // bits wavfile1=''; // file name, this will be overwritten. wdat1=zeros(1,10);// data of the .wav, only one first channel, this will be overwritten. // for fft fft_point1=512; db_fft1=zeros(1,512); // this will be overwritten. phi_fft1=zeros(1,512); // this will be overwritten. // for many frame db_fft1s=zeros(1,512); // this will be overwritten. phi_fft1s=zeros(1,512); // this will be overwritten. Min_Freq=150; // set plot maximum frequency by unit is [Hz] Max_Freq=7500; // set plot maximum frequency by unit is [Hz] // for display width0=400; // default display width ydat0=zeros(1,width0+1); // for all data ydat1=zeros(1,width0+1); // for portion xziku0=[0:1:width0]; // for all data xziku1=[0:1:width0]; // for portion sp1=1; // start point for actual display ep1=1; // end point for actual display tframe=1; tframe0=tframe; // this is for pop and push // for many frames display PART_LIST0=[' part 1' ; ' part 2' ; ' part 3' ; ' part 4' ; ' part 5' ; ' part 6' ; ' part 7' ; 'part 8' ; 'part 9' ; 'part 10']; max_frames=10; FRAME_QTY=1; // dummy set sp1s=ones(1,max_frames); // start point for actual display ep1s=ones(1,max_frames); // end point for actual display // for others defch1=1; // default channel 1 f_win=1; // flag if windows //--------------------------------------- // function save_fb() save('fb_bak.dat',fb_out_org, fb_am_out_org, fb_am_out_prop_org); endfunction function [fb_out_org, fb_am_out_org, fb_am_out_prop_org]=load_fb() load('fb_bak.dat','fb_out_org', 'fb_am_out_org', 'fb_am_out_prop_org'); endfunction function save_fb2() save('fb_bak2.dat',fb_out_init, fb_am_out_init, fb_am_out_prop_init); endfunction function [fb_out_init, fb_am_out_init, fb_am_out_prop_init]=load_fb2() load('fb_bak2.dat','fb_out_init', 'fb_am_out_init', 'fb_am_out_prop_init'); endfunction // function save_wavprop() save('wavprop.dat',fs1,size1,bit1,wdat1,sp1,ep1,ydat0,ydat1,xziku0,xziku1); endfunction function [fs1,size1,bit1,wdat1,sp1,ep1,ydat0,ydat1,xziku0,xziku1]=load_wavprop() load('wavprop.dat','fs1','size1','bit1','wdat1','sp1','ep1','ydat0','ydat1','xziku0','xziku1'); endfunction // function save_wfft() save('wavfft.dat',fft_point1,db_fft1,phi_fft1); endfunction function [fft_point1,db_fft1,phi_fft1]=load_wfft() load('wavfft.dat','fft_point1','db_fft1','phi_fft1'); endfunction // //-------------------------------------- function [wavfile1,chs1,qty1,fs1,bit1,f412,SEL_CODE]=get_wavfile_pro() global LOAD_FB_FLAG // choose wave file if T_DEMO == 1 then //+ select //SEL_CODE=x_choose0(['ka_sample (into 5 parts)';'-'; '-'; 'manual select'],['Please select one'],'default') SEL_CODE=x_choose0(['ka_sample (into 5 parts)';'-'; '-'; '-';'(load fb)';'(save fb)'],['Please select one'],'default') if SEL_CODE == 4 then // when manual select SEL_CODE = 999; elseif SEL_CODE <= 0 then SEL_CODE=1; elseif SEL_CODE == 5 then // select load fb data and set SEL_CODE =1 LOAD_FB_FLAG=1; SEL_CODE=1; elseif SEL_CODE == 6 then // select save fb data and set SEL_CODE =1 LOAD_FB_FLAG=-1; SEL_CODE=1; end //- select if SEL_CODE == 1 then // ka_sample part 1 disp(' This is demonstration. Copy this file and ka_sample.wav to your scilab''s bin or home directory.'); [x,ierr]=fileinfo('ka_sample.wav'); xs=size(x); if xs(1)==0 then disp('Choose ka_sample.wav file'); if scilab_version_number >= 5 then wavfile1=uigetfile(["*.wav"],"","Choose ka_sample.wav file"); else wavfile1=tk_getfile(Title="Choose ka_sample.wav file"); end else wavfile1='ka_sample.wav'; end elseif SEL_CODE == 2 then // ka_sample part 2 disp(' This is demonstration. Copy this file and ka_sample.wav to your scilab''s bin or home directory.'); [x,ierr]=fileinfo('ka_sample.wav'); xs=size(x); if xs(1)==0 then disp('Choose ka_sample.wav file'); if scilab_version_number >= 5 then wavfile1=uigetfile(["*.wav"],"","Choose ka_sample.wav file"); else wavfile1=tk_getfile(Title="Choose ka_sample.wav file"); end else wavfile1='ka_sample.wav'; end elseif SEL_CODE == 3 then // ka_sample part3 disp(' This is demonstration. Copy this file and ka_sample.wav to your scilab''s bin or home directory.'); [x,ierr]=fileinfo('ka_sample.wav'); xs=size(x); if xs(1)==0 then disp('Choose ka_sample.wav file'); if scilab_version_number >= 5 then wavfile1=uigetfile(["*.wav"],"","Choose ka_sample.wav file"); else wavfile1=tk_getfile(Title="Choose ka_sample.wav file"); end else wavfile1='ka_sample.wav'; end else // include SEL_CODE == 999 if scilab_version_number >= 5 then wavfile1=uigetfile(["*.wav"],"","Choose a .wav file"); else wavfile1=tk_getfile(Title="Choose a .wav file name"); end end else if scilab_version_number >= 5 then wavfile1=uigetfile(["*.wav"],"","Choose a .wav file"); else wavfile1=tk_getfile(Title="Choose a .wav file name"); end SEL_CODE = 999; 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 display the .wav property function disp_pro_wav() disp(qty1,'data quantity',bit1, 'bits',chs1,'channels',fs1,'sampling frequency'); endfunction //-------------------------------------- function [wdat1, size1]=read_one_ch_of_wav(wavfile1) y=wavread(wavfile1); ysize=size(y); if ysize(2) == chs1 then wdat1=zeros(1,ysize(1)); if chs1 == 1 then // for v=1:ysize(1) // wdat1(v)=y(v); // end wdat1=y'; else disp('Channels are more than two, but, only 1st channel is stored.'); for v=1:ysize(1) wdat1(v)=y(v,defch1); end end size1=ysize(1); disp(size1,'stored data size is'); elseif ysize(1) == chs1 then // for scilab-4.1.2 wdat1=zeros(1,ysize(2)); if chs1 == 1 then // for v=1:ysize(2) // wdat1(v)=y(v); // end wdat1=y; else disp('Channels are more than two, but, only 1st channel is stored.'); for v=1:ysize(2) wdat1(v)=y(defch1,v); end end size1=ysize(2); disp(size1,'stored data size is'); end endfunction //----------------------------------------------- function plot_wave1(disp0) wb0=xget('window'); // stack old window xset('window',disp0); // create new windows clf(); subplot(211); //plot(xziku0,ydat0,'b'); plot2d(xziku0,ydat0,style=[color("blue")]); xtitle('waveform all'); // When linux scilab-3.1 2nd subplot does not work well. subplot(212); //plot(xziku1,ydat1,'b'); plot2d(xziku1,ydat1,style=[color("blue")]); xtitle('waveform selected'); xset('window',wb0); // push old windows endfunction //----------------------------------------------- function plot_wave1s(disp0) global sp1 ep1 wb0=xget('window'); // stack old window xset('window',disp0); // create new windows clf(); subplot( QTY_FRAME+1, 1, 1); //plot(xziku0,ydat0,'b'); plot2d(xziku0,ydat0,style=[color("blue")]); xtitle('waveform all'); // When linux scilab-3.1 2nd subplot does not work well. for v=1:QTY_FRAME subplot(QTY_FRAME+1,1, v+1); //plot(xziku1,ydat1,'b'); sp1=sp1s(v); ep1=ep1s(v); [ydat1,xziku1]=make_width_data( wdat1, size1); plot2d(xziku1,ydat1,style=[color("blue")]); wstr1='waveform selected' + PART_LIST0(v); xtitle( wstr1 ); end xset('window',wb0); // push old windows sp1=sp1s(tframe); // back to init ep1=ep1s(tframe); // back to init endfunction //-------------------------------------- function [ydat1,xziku1]=make_width_data( work1, wsize1) wsize2=wsize1 - sp1 + 1; wsize3=ep1-sp1+1; if wsize3 > (width0+1) then wstep1= wsize3 / (width0+1); for v=1:(width0+1) ydat1(v)= work1( sp1 + int(wstep1 * (v - 1))); xziku1(v)= sp1 + int(wstep1 * (v - 1)); end else for v=1:(width0+1) if v <= wsize3 then ydat1(v)=work1((sp1-1)+v); else ydat1(v)=0.; end xziku1(v)=(sp1-1)+v; end end // // plot_wave1(0); endfunction //---------------------------------------------------------------- function [wsp1, wep1]= reset_sp1_ep1(wsize1) wsp1=1; wep1=wsize1; endfunction //---------------------------------------------------------------- function [wsp1, wep1, wtframe]= set_sp1_ep1(wsize1) txt1=['start point';'end point';'no. of frame']; wstr1=sprintf('%d',sp1); wstr2=sprintf('%d',ep1); wstr3=sprintf('%d',tframe); sig1=x_mdialog('Input start point and end point for portion display.',txt1, [wstr1 ; wstr2 ; wstr3]); if sig1==[] then arg1=evstr(wstr1); arg2=evstr(wstr2); arg3=evstr(wstr3); else arg1=evstr(sig1(1)); arg2=evstr(sig1(2)); arg3=evstr(sig1(3)); end // //disp(arg2,'arg2',arg1,'arg1'); // wsp1=sp1; wep1=ep1; if arg1 >= 1 then if arg2 <= wsize1 then if arg1 < arg2 then wsp1=arg1; wep1=arg2; wtframe=arg3; disp(tframe,'tframe',wep1,'ep1',wsp1,'sp1'); end end end endfunction //---------------------------------------------------------------- function snd_play1() // 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! wdatw=zeros(ep1-sp1+1); for v=sp1:ep1 wdatw(v-sp1+1)=wdat1(v); end if f_win == 1 then if f412 == 1 then // for windows scilab-4.1.2 sound(wdatw' ,fs1,bit1); else // for windows scilab-3.1.1 if fs1 == 22050 then sound(wdatw,fs1,bit1); disp('This function may work, if luckily.'); elseif fs1 == 44100 then // down-sampling from 44100 to 22050 wdatw2=zeros((ep1-sp1)/2+1); for v=1:((ep1-sp1)/2+1) wdatw2(v)=wdat1(sp1 + 2 * (v -1) ); 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 else disp('Sorry, but, this function does not work.'); end endfunction // //---------------------------------------------------------------- function snd_save1() // wavefilename= input(' + enter file name for saved .wav file =>',["string"]); // wdatw=zeros(ep1-sp1+1); for v=sp1:ep1 wdatw(v-sp1+1)=wdat1(v); end if f412 == 1 then // for windows scilab-4.1.2 wavwrite(wdatw' ,fs1,bit1,wavefilename ); else // for windows scilab-3.1.1 wavwrite(wdatw,fs1,bit1, wavefilename); end endfunction //--- check platform is windows ----------------------------------- if isdir('c:\\') then f_win=1; else f_win=0; end // //================================================================= // // +MAIN (1)program starts // [wavfile1,chs1,qty1,fs1,bit1,f412,SEL_CODE]=get_wavfile_pro(); disp_pro_wav(); [wdat1, size1]=read_one_ch_of_wav(wavfile1); sp1=1; ep1=size1; [ydat0,xziku0]=make_width_data( wdat1, size1); QTY_FRAME=1; if T_DEMO==1 then if SEL_CODE == 1 then // ka_sample QTY_FRAME=5; // brust, un-brust, a-, a|, a= sp1s(1)=391; ep1s(1)=840; sp1s(2)=872; ep1s(2)=1345; //----------------------------------------------------------- // kato teki ni hennka siteiru tokorode tokutyou ga mada hakkiri sinai tokoro no // fittting ha umaku dekite inai. kuhuu ga hituyou de aru sp1s(3)=1870; ep1s(3)=2370; //sp1s(4)=2460; //ep1s(4)=3060; //++++++++++++++++++++++++!!!+++++++++++++++++++++++++++++++ // needs sp1-ep1 area adjusted to more feature appear sp1s(4)=2205; ep1s(4)=2805; // sp1s(5)=5000; ep1s(5)=5700; // sp1=sp1s(tframe); ep1=ep1s(tframe); elseif SEL_CODE == 2 then // ?_sample QTY_FRAME=1; sp1s(1)=3129; ep1s(1)=3700; // sp1=3129; // ?? ep1=3700; // ?? elseif SEL_CODE == 3 then // ?_sample QTY_FRAME=1; sp1s(1)=7810; ep1s(1)=8400; // sp1=7810; // ?? ep1=8400; // ?? else // include SEL_CODE == 999 [sp1, ep1, tframe]= set_sp1_ep1(size1); sp1s(tframe)=sp1; ep1s(tframe)=ep1; end end [ydat1,xziku1]=make_width_data( wdat1, size1); // ?? //plot_wave1(0); plot_wave1s(0); // // -MAIN (1)program starts // select .wav and set portion to make it fft //==============================================================++= // // // 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 agwru tame 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 [frq1,sfrq1,is1,ie1]=set_frq_fft( spec1024 ) // spec1024, tube model response no syuhasuu bunkainou wo agwru tame dfsw= fs1 / spec1024 ; is1= ceil(Min_Freq / dfsw); ie1= int(Max_Freq / dfsw); dfsw2=dfsw; frq1=[(dfsw * is1):dfsw2:(dfsw * ie1)]; frqs=size(frq1); sfrq1=frqs(2); is1=is1+1; // fft(0) is DC ie1=ie1+1; endfunction //----------------------------------------------- function plot_fft1(disp0) wb0=xget('window'); // stack old window xset('window',disp0); // create new windows [frq1,sfrq1,is1,ie1]=set_frq(0); clf(); subplot(211); gainplot(frq1,db_fft1,phi_fft1); xtitle('frequency response of waveform selected by fft analysis'); subplot(212); //plot(xziku1,ydat1,'b'); plot2d(xziku1,ydat1,style=[color("blue")]); xtitle('waveform selected'); xset('window',wb0); // push old windows endfunction // function plot_fft1s(disp0) wb0=xget('window'); // stack old window xset('window',disp0); // create new windows [frq1,sfrq1,is1,ie1]=set_frq(0); clf(); s23=size(db_fft1s); db_fft0=zeros(1,s23(2)); phi_fft0=zeros(1,s23(2)); for v=1:QTY_FRAME subplot(QTY_FRAME,1,v); for w=1:s23(2) db_fft0(1,w)=db_fft1s(v,w); phi_fft0(1,w)=phi_fft1s(v,w); end gainplot(frq1,db_fft0,phi_fft0); wstr1='frequency response of waveform selected by fft analysis' + PART_LIST0(v); xtitle(wstr1); //subplot(212); //plot(xziku1,ydat1,'b'); //plot2d(xziku1,ydat1,style=[color("blue")]); //xtitle('waveform selected'); end xset('window',wb0); // push old windows endfunction //---------------------------------------------------- function [db_fft1,phi_fft1]=do_fft_wav() 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); end // -if size1 <= ( sp1 + fft_point1 - 1 ) then endfunction //-------------------------------------------------------------- function [db_fft1s,phi_fft1s]=do_fft_wavs() global sp1 ep1 win_hn=window('hn',fft_point1); // make hanning windows data wdatw=zeros(1, fft_point1); // [frq1,sfrq1,is1,ie1]=set_frq(0); respw=zeros(1,sfrq1); // for wloop=1:QTY_FRAME sp1=sp1s(wloop); ep1=ep1s(wloop); // if size1 >= ( sp1 + fft_point1 - 1 ) then for v=1:fft_point1 wdatw(v)=wdat1(sp1 + v - 1); end wdatw2= wdatw .* win_hn; fftwout=fft( wdatw2, -1 ); ct0=1; for loop=is1:ie1 respw(ct0)=fftwout(loop+1); ct0=ct0+1; end [db_fft1,phi_fft1] =dbphi(respw); s23=size(db_fft1); if wloop == 1 then db_fft1s=zeros(QTY_FRAME,s23(2)); phi_fft1s=zeros(QTY_FRAME,s23(2)); end for v=1:s23(2) db_fft1s(wloop,v)=db_fft1(v); phi_fft1s(wloop,v)=phi_fft1(v); end end // -if size1 <= ( sp1 + fft_point1 - 1 ) then end // for v=1:QTY_FRAME endfunction //================================================================= // peak detect on FFT spectrum // // (1)2-3 smoothing on FFT data / nizi sanzi takousiki tekigou niyoru smoothing // (2)differentiation (def) / heikatsuka bibun // (3)detect peak as the point from plus to minus // sn=0; // data quantity. this will be overwritten. sm1=2; // set order of 2-3 smoothing sm2=2; // set order of 5 points differentiation (def) // if sm2=1, more peak candidate may appear // if sm2=2, sometime gentle slop (nadaraka) peak will be lost swnd= zeros(1,128); // swnd(1),swnd(2),.... this will be overwritten. sm1out=zeros(1,512); // this will be overwritten. calculate weighted average sm2out=zeros(1,512); // this will be overwritten. def of sm1out npk=0; pklist=zeros(9,512); // nth, freq, its value of peak candinates // nth, freq, its value of estimated left edge // nth, freq, its value of estimated right edge //------------------------------------------------------ function plot_fft_sm1(disp0) [frq1,sfrq1,is1,ie1]=set_frq(0); w0=xget('window'); xset('window',disp0); clf(); subplot(211); gainplot(frq1,db_fft1,phi_fft1); //plot2d(frq1,db_fft1,color("black"),logflag="ln"); xtitle('frequency response of waveform selected by fft analysis'); plot2d(frq1,sm1out,color("green"),logflag="ln"); legend(['org'; 'smoothed'],3); //xtitle('smoothed frequency response '); subplot(212); plot2d(frq1,sm2out,color("red"),logflag="ln"); xtitle('def of smoothed frequency response '); xset('window',w0); endfunction //---------------------------------------------------------- function [sn,swnd,sm1out,sm2out,npk,pklist]=smooth1(disp_code) // get db_fft1 data size ax=size(db_fft1); sn=ax(2); swnd= zeros(1,128); // preparetion smoothing 2-3 sm=sm1; snorm=(4.0 * sm * sm - 1.0) * (2.0 * sm + 3.0 ) / 3.0; sl=3.0 * sm * (sm + 1.0) -1.0; for loop=1:(sm+1) swnd(loop)= sl - 5.0 * (loop -1.0) * (loop - 1.0); end // calculate weighted average sm1out=zeros(sn); dmax0= 0.; sigma0=0.; ep0=0.; for v=1:sn if v < (sm+1) then sm1out(v)=db_fft1(v); elseif v > (sn-sm) then sm1out(v)=db_fft1(v); else sum0 = db_fft1(v) * swnd(1); for v2=1:sm sum0 = sum0 + db_fft1(v+v2) * swnd(1+v2); sum0 = sum0 + db_fft1(v-v2) * swnd(1+v2); end sm1out(v)= sum0 / snorm; // sa0=db_fft1(v) - sm1out(v); sigma0=sigma0+sa0*sa0; end ep0=ep0+db_fft1(v); if sm1out(v) > dmax0 then dmax0= sm1out(v); end end // // calculate def of smoothed sm=sm2; sm2out=zeros(sn); for v=1:sn if v<(sm+1) then sm2out(v)=0.; elseif v>(sn-sm) then sm2out(v)=0.; else sm2out(v)=0.; for v2=1:sm sm2out(v)=sm2out(v)+ ( sm1out(v+v2) - sm1out(v-v2)) * v2; end end end // find from + to - in def pklist=zeros(9,sn); [frq1,sfrq1,is1,ie1]=set_frq(0); npk=0; for v=2:sn if sm2out(v-1) > 0 then if sm2out(v) < 0 then dmax1=-9999.0; stack_v=-1; //.. .. for v2=(v-sm):(v+sm) if v2 < 1 then elseif v2 > sn then else if db_fft1(v2) > dmax1 then dmax1=db_fft1(v2); stack_v=v2; end end end //.. .. if stack_v > 0 then npk=npk+1; pklist(1,npk)=stack_v; pklist(2,npk)=frq1(stack_v); pklist(3,npk)=dmax1; end end end end // find left edge estimated based on smoothed signal for loop=1:npk if pklist(1,loop) == 1 then pklist(4,loop) = 1; pklist(5,loop) = frq1(1); pklist(6,loop) = sm1out(1); else pklist(4,loop) = 1; pklist(5,loop) = frq1(1); pklist(6,loop) = sm1out(1); for v=pklist(1,loop):-1:2 if sm1out(v-1) > sm1out(v) then pklist(4,loop) = v; pklist(5,loop) = frq1(v); pklist(6,loop) = sm1out(v); break; end end end end // for loop=1:npk //.. // find right edge estimated based on smoothed signal for loop=1:npk if pklist(1,loop) == sn then pklist(7,loop) = sn; pklist(8,loop) = frq1(sn); pklist(9,loop) = sm1out(sn); else pklist(7,loop) = sn; pklist(8,loop) = frq1(sn); pklist(9,loop) = sm1out(sn); for v=pklist(1,loop):1:(sn-1) if sm1out(v+1) > sm1out(v) then pklist(7,loop) = v; pklist(8,loop) = frq1(v); pklist(9,loop) = sm1out(v); break; end end end end // for loop=1:npk //... if disp_code == 1 then disp(''); //-- peak candinate list out disp('-- peak candinate list out --'); disp('-- peak [Hz] [dB] left edge [Hz] [dB] right edge [Hz] [dB] --'); for v=1:npk // disp([pklist(1,v) pklist(2,v) pklist(3,v)]); disp([ pklist(2,v) pklist(3,v) pklist(5,v) pklist(6,v) pklist(8,v) pklist(9,v)]); end //-- peak candinate list out disp(''); end endfunction //---- function smooths(disp0,peak_disp_code) // when peak_disp_code=1, display peak candidates, beside peak_disp_code=0, non global db_fft1 phi_fft1 s23=size(db_fft1s); db_fft1=zeros(1,s23(2)); phi_fft1=zeros(1,s23(2)); [frq1,sfrq1,is1,ie1]=set_frq(0); w0=xget('window'); xset('window',disp0); clf(); ngamen=QTY_FRAME * 2; for v=1:QTY_FRAME for w=1:s23(2) db_fft1(1,w)=db_fft1s(v,w); phi_fft1(1,w)=phi_fft1s(v,w); end disp( PART_LIST0(v)); [sn,swnd,sm1out,sm2out,npk,pklist]=smooth1(peak_disp_code); // when smooth(1), values will be displayed ! subplot(ngamen,1,2*v-1); gainplot(frq1,db_fft1,phi_fft1); wstr1=PART_LIST0(v) + ' frequency response of waveform selected by fft analysis' xtitle(wstr1); plot2d(frq1,sm1out,color("green"),logflag="ln"); legend(['org'; 'smoothed'],3); subplot(ngamen,1,2*v); plot2d(frq1,sm2out,color("red"),logflag="ln"); xtitle('def of smoothed frequency response '); end // for v=1:QTY_FRAME xset('window',w0); // endfunction // //================================================================= // // +MAIN (2)program starts // about fft if T_DEMO==1 then [db_fft1, phi_fft1]=do_fft_wav(); // ?? [db_fft1s, phi_fft1s]=do_fft_wavs(); plot_fft1s(1); ////////////////////////////////////////////////////// if T_SKIP1 <> 1 then smooths(2,0); end ////[sn,swnd,sm1out,sm2out,npk,pklist]=smooth1(5); ////plot_fft_sm1(5); ///////////////////////////////////////////////////// end // // -MAIN (2)program starts // //==============================================================++= // // For filter bank // // st_freq=500.0; end_freq=st_freq * 10; end_freq2=0.; n_stepf=10; bpf_fb0=zeros(512,3); // *,1 center frequency. this will be overwritten. // *,2 lower frequency of bpf band. this will be overwritten. // *,3 higher frequency of bpf band. this will be overwritten. fir_steps=513; // This must be odd. this will be overwritten. kwindows=1; wft=zeros(n_stepf,fir_steps); // this will be overwritten. wfm=zeros(n_stepf,fir_steps); // this will be overwritten. fr=zeros(fir_steps); // this will be overwritten. fb_out=zeros(n_stepf,1); // this will be overwritten. fb_am_out=zeros(n_stepf,1); // this will be overwritten. nmax_am_prop=5; fb_am_out_prop=zeros(n_stepf,nmax_am_prop); // *,1 average value about fb_am_out // *,2 peak value about fb_am_out cfb_list=ones(1,2); // this will be overwritten. fb_out_org=zeros(n_stepf,1); // this will be overwritten. fb_am_out_org=zeros(n_stepf,1); // this will be overwritten. fb_am_out_prop_org=zeros(n_stepf,nmax_am_prop); fb_out_init=zeros(n_stepf,1); // this will be overwritten. fb_am_out_init=zeros(n_stepf,1); // this will be overwritten. fb_am_out_prop_init=zeros(n_stepf,nmax_am_prop); fb_out_work1=zeros(n_stepf,1); // this will be overwritten. fb_am_out_work1=zeros(n_stepf,1); // this will be overwritten. fb_am_out_prop_work1=zeros(n_stepf,nmax_am_prop); //------------------------------------------------- function [st_freq, end_freq, n_stepf, fir_steps, kwindows,end_freq2]= set_fb_para() // codex defines which term.... txt1=['start freq'; 'end freq'; 'freq steps'; 'steps of fir filter (odd)'; 'priority: frq-resolution (1), drop value (2)' ; 'end freq of additional bpf']; wstr1=sprintf('%f',st_freq); wstr2=sprintf('%f',end_freq); wstr3=sprintf('%d',n_stepf); wstr4=sprintf('%d',fir_steps); wstr5=sprintf('%d',kwindows); wstr6=sprintf('%f',end_freq2); sig1=x_mdialog('set fir filter bank parameters ',txt1, [wstr1; wstr2; wstr3; wstr4; wstr5]); if sig1==[] then arg1=evstr(wstr1); arg2=evstr(wstr2); arg3=evstr(wstr3); arg4=evstr(wstr4); arg5=evstr(wstr5); arg6=evstr(wstr6); else arg1=evstr(sig1(1)); arg2=evstr(sig1(2)); arg3=evstr(sig1(3)); arg4=evstr(sig1(4)); arg5=evstr(sig1(5)); arg6=evstr(sig1(6)); end st_freq=arg1; end_freq=arg2; n_stepf=arg3; fir_steps=arg4; kwindows=arg5; end_freq2=arg6; endfunction //----------------------------------------------------------- function [bpf_fb0,wft]=make_bpf_fb0(disp0) global n_stepf bpf_fb0=zeros(n_stepf,3); //+ additional bpf if end_freq2 > end_freq then n_stepf=n_stepf-1; end //- additional bpf // log gauge divide a1=log(st_freq); b1=log(end_freq); c1=(b1-a1)/n_stepf; for v=1:n_stepf a2=exp(a1+ c1 * (v-1)); b2=exp(a1+ c1 * v); c2=(b2-a2)/2. + a2; bpf_fb0(v,1)=c2; bpf_fb0(v,2)=a2; bpf_fb0(v,3)=b2; if disp0 >= 1 then disp( [ v bpf_fb0(v,1) bpf_fb0(v,2) bpf_fb0(v,3)]); end end // for V=1:n_stepf //+ additional bpf if end_freq2 > end_freq then n_stepf=n_stepf+1; bpf_fb0(n_stepf,1)=end_freq + (end_freq2-end_freq)/2.; bpf_fb0(n_stepf,2)=end_freq; bpf_fb0(n_stepf,3)=end_freq2; if disp0 >= 1 then disp('additional bpf'); v=n_stepf; disp( [ v bpf_fb0(v,1) bpf_fb0(v,2) bpf_fb0(v,3)]); end end //- additional bpf // // ++ mk_fir_coeff // w0=xget('window'); // stack old window wft=zeros(n_stepf,fir_steps); for vloop=1:n_stepf cfreq(1)= bpf_fb0(vloop,2) / fs1; cfreq(2)= bpf_fb0(vloop,3) / fs1; fpar(1)=0.1; // dummy data fpar(2)=-0.1; // dummy data if kwindows == 2 then [wft2,wfm2,fr2]=wfir('bp',fir_steps,cfreq,'hn',fpar); // gensui riritu yuusen else [wft2,wfm2,fr2]=wfir('bp',fir_steps,cfreq,'hm',fpar); // bunkai nou yuusen end for v2=1:fir_steps wft(vloop,v2)=wft2(v2); end if disp0 >=2 then //++ only one will be done if vloop==1 then xset('window',disp0); // create new windows clf(); [frq1,sfrq1,is1,ie1]=set_frq(0); fr0=zeros(2); // freq list selected frq1f=zeros(2); // freq responce selected phf0=zeros(2); // phase selected sv0=size(fr2); end // if vloop==1 then //-- only one will be done sc0=0; for v=1:sv0(2) fx0= fr2(v) * fs1; if (fx0 >= frq1(1)) &(fx0 <= frq1(sfrq1)) then sc0=sc0+1; fr0(sc0)=fx0; frq1f(sc0)= 20. * log10( wfm2(v)); phf0(sc0)= -1.0 * ((fir_steps - 1) / 2.0 ) * ( 2. * %pi * v / fir_steps); end end // for v=1:sv0(2) gainplot(fr0,frq1f,phf0); end // vloop+1:n_stepf end // if disp0 >=2 then xset('window',w0); // push old windows // // -- mk_fir_coeff // endfunction //----------------------------------------------------------- function plot_fb(disp_code, xmessfb) w0=xget('window'); // stack old window if disp_code >= 1 then xset('window',disp_code); // create new windows clf(); wsize2=size1 - sp1 + 1; wsize3=ep1-sp1+1; zdata1=zeros(1,width0+1); zdata1p=zeros(1,width0+1); ydata1=zeros(1,width0+1); for vloop=1:n_stepf if wsize3 > (width0+1) then wstep1= (ep1-sp1) / width0; for v=1:(width0+1) zdat1(v)= fb_out(vloop,1 + int(wstep1 * (v - 1))); zdat1p(v)= fb_am_out(vloop, 1 + int(wstep1 * (v - 1))); xziku1(v)= sp1 + int(wstep1 * (v - 1)); if vloop==1 then ydat1(v)= wdat1( sp1 + int(wstep1 * (v - 1))); end end else for v=1:(width0+1) if v <= wsize3 then zdat1(v)=fb_out(vloop,v); zdat1p(v)=fb_am_out(vloop,v); if vloop==1 then ydat1(v)= wdat1((sp1-1)+v); end else zdat1(v)=0.; zdat1p(v)=0.; if vloop==1 then ydat1(v)= 0.; end end xziku1(v)=(sp1-1)+v; end end if vloop==1 then subplot((n_stepf+1),1,1); plot(xziku1,ydat1,'b'); wstr1=sprintf('waveform selected, original'); wstr1=wstr1 + ' / ' + xmessfb; xtitle(wstr1); end subplot((n_stepf+1),1,(vloop+1)); plot(xziku1,zdat1,'b'); plot(xziku1,zdat1p,'g'); wstr1=sprintf('Band Pass filtered waveform selected (%dHz- %dHz)',int(bpf_fb0(vloop,2)),int(bpf_fb0(vloop,3))); mv0=fb_am_out_prop(vloop,1); if mv0 > 0. then mv1= int(20. * log10(mv0)); else mv1=-120; end pk0=fb_am_out_prop(vloop,2); if pk0 > 0. then pk1= int( 20. * log10(pk0)); else pk1=-120; end wstr2=sprintf('%d/%d dB',mv1,pk1); // average / peak [dB] xtitle( wstr1,'',wstr2 ); end // for vloop=1:n_stepf end // if disp_code >= 1 then xset('window',w0); // push old windows endfunction //----------------------------------------------------------- function [fb_out,fb_am_out,fb_am_out_prop]=do_fb(disp0) // w0=(ep1-sp1)+1; fb_out=zeros(n_stepf,w0); dt0= int((fir_steps-1) / 2); ep1x=ep1+dt0; if ep1x > size1 then ep1x = size1; end sp1x=sp1+dt0; if sp1x < fir_steps then sp1x = fir_steps; end //+filtering+ for vloop=1:n_stepf if disp0 >= 1 then wstr1=sprintf('+entering filter bank no.%d/%d',vloop,n_stepf); disp(wstr1); end for loop=sp1x:ep1x ydz = 0.; for loop2=1:fir_steps ydz = ydz + wft(vloop,loop2) * wdat1( 1 - loop2 + loop ); end //fir_out1(loop - dt0)=ydz; fb_out(vloop,loop - dt0 -sp1 + 1)=ydz; end end // if vloop=1:n_stepf //-filtering- // //+am+ nsize_max=5000; g_work2=zeros(10,nsize_max); fb_am_out=zeros(n_stepf,w0); fb_am_out_prop=zeros(n_stepf,nmax_am_prop); fir_ou1=zeros(1,size1); fir_out1p=zeros(1,size1); for vloop=1:n_stepf if disp0 >= 1 then wstr1=sprintf('+entering AM de-modulation, peak trace no.%d/%d',vloop,n_stepf); disp(wstr1); end for v2=sp1:ep1 fir_out1(v2)= fb_out(vloop,v2-sp1+1); fir_out1p(v2)=0.; end //+from trace peak //+ spx=-1; spy=-1; epx=-1; epy=-1; c0=-1; w0=-1.; area0=-1.; if (ep1-sp1) > 4 then for v=(sp1+1):(ep1-1) if fir_out1(v) > 0. then // only postive side if (fir_out1(v) > fir_out1(v-1)) & (fir_out1(v) > fir_out1(v+1)) then spx=epx; spy=epy; epx=v; epy=fir_out1(v); c0=c0+1; // if g_work2(1,1) < nsize_max then if g_work2(1,1) <= 0 then g_work2(1,1)=1; // first start g_work2(3,1)=epx; g_work2(4,1)=epy; g_work2(5,1)=epx; g_work2(6,1)=epy; else d0=int(g_work2(1,1)); if epy >= g_work2(6,d0) then // rising or equal g_work2(2,d0)=0; //current g_work2(5,d0)=epx; g_work2(6,d0)=epy; else // falling if g_work2(2,d0)==0 then // complete data g_work2(2,d0)=1; //complete g_work2(7,d0)= g_work2(5,d0) - g_work2(3,d0); g_work2(8,d0)= g_work2(6,d0) - g_work2(4,d0); g_work2(9,d0)= g_work2(7,d0) * g_work2(8,d0) / 2.0 ; if (g_work2(7,d0) > 0.) & (g_work2(8,d0) > 0.) then if g_work2(9,d0) > area0 then area0=g_work2(9,d0); w0=g_work2(8,d0) / g_work2(7,d0); end end // atan //g_work2(10,d0)= atan( (g_work2(8,d0) / g_work2(7,d0))) / 2.0 / %pi * 360.0; // g_work2(1,1)=g_work2(1,1)+1; // set as new start point d1=int(g_work2(1,1)); g_work2(3,d1)=epx; g_work2(4,d1)=epy; g_work2(5,d1)=epx; g_work2(6,d1)=epy; elseif g_work2(2,d0) <= -1 then // restart g_work2(3,d0)=epx; g_work2(4,d0)=epy; g_work2(5,d0)=epx; g_work2(6,d0)=epy; end end end // if g_work2(1,1) <= 0 then else disp(' error: g_g_work2(1,1) >= nsize_max , in do_fb() '); end // if g_work2(1,1) < nsize_max then // // if c0 >= 1 then a0= (epy - spy) / (epx - spx); b0= epy - a0 * epx; for v2=spx:epx fir_out1p(v2)= a0 * v2 + b0; end end // if c0 >= 1 then // end // if (fir_out1(v) > fir_out1(v-1)) & (fir_out1(v) > fir_out1(v+ 1)) then end // if fir_out1(v) > 0. then end //v=(sp1+1):(ep1-1) end if c0 < 1 then epx=sp1-1; end for v=(epx+1):ep1 fir_out1p(v)=0; end pk0=0.; mv0=0.; for loop=sp1:ep1 fb_am_out(vloop,loop - sp1 + 1)= fir_out1p(loop); if fir_out1p(loop) > pk0 then pk0=fir_out1p(loop); end mv0=mv0+fir_out1p(loop); end fb_am_out_prop(vloop,1)= mv0 / (ep1-sp1+1); fb_am_out_prop(vloop,2)= pk0; //+ //-from trace peak end // for vloop=1:n_stepf //-am- endfunction //------------------------------------------------------------------ function do_fb2(wdat1x,sp1x,ep1x,st_freq1x,end_freq1x,end_freq2x,n_stepf1x,disp_no,xmessfb1x) global sp1 ep1 fir_steps global st_freq end_freq n_stepf end_freq2 global bpf_fb0 wft fb_out fb_am_out fb_am_out_prop global wdat1 sp1bak=sp1; ep1bak=ep1; fir_stepsbak=fir_steps; st_freqbak=st_freq; end_freqbak=end_freq; end_freq2bak=end_freq2; n_stepfbak=n_stepf; wdat1bak=wdat1; sp1=sp1x; ep1=ep1x; s15=size(wdat1x); if s15(2) > s15(1) then wdat1=zeros(1,s15(2)); for vloop=1:s15(2) wdat1(1,vloop)=wdat1x(1,vloop); end else wdat1=zeros(1,s15(1)); for vloop=1:s15(1) wdat1(1,vloop)=wdat1x(vloop,1); end end if fir_steps > sp1 then if (int(sp1 / 2) * 2) <> sp1 then fir_steps=sp1; else fir_steps=sp1-1; end str1=sprintf('%d',fir_steps); str1="waring: re-set fir_steps is " + str1; disp( str1 ); end st_freq=st_freq1x; end_freq=end_freq1x; end_freq2=end_freq2x; n_stepf=n_stepf1x; [bpf_fb0,wft]=make_bpf_fb0( disp_no); [fb_out,fb_am_out,fb_am_out_prop]=do_fb(1); // xmessfb=' ka_sample original (500-3000/4)'; plot_fb( disp_no+1 ,xmessfb1x) // push back sp1=sp1bak; ep1=ep1bak; fir_steps=fir_stepsbak; st_freq=st_freqbak; end_freq=end_freqbak; end_freq2=end_freq2bak; n_stepf=n_stepfbak; wdat1=wdat1bak; endfunction //------------------------------------------------------------------ function do_fb3_init(wdat1x,sp1x,ep1x,disp_no,xmessfb1x) // filter bank for init data global sp1 ep1 fir_steps global wdat1 global fb_out_init fb_am_out_init fb_am_out_prop_init sp1bak=sp1; ep1bak=ep1; fir_stepsbak=fir_steps; wdat1bak=wdat1; sp1=sp1x; ep1=ep1x; s15=size(wdat1x); if s15(2) > s15(1) then wdat1=zeros(1,s15(2)); for vloop=1:s15(2) wdat1(1,vloop)=wdat1x(1,vloop); end else wdat1=zeros(1,s15(1)); for vloop=1:s15(1) wdat1(1,vloop)=wdat1x(vloop,1); end end if fir_steps > sp1 then if (int(sp1 / 2) * 2) <> sp1 then fir_steps=sp1; else fir_steps=sp1-1; end str1=sprintf('%d',fir_steps); str1="waring: re-set fir_steps is " + str1; disp( str1 ); end [bpf_fb0,wft]=make_bpf_fb0( disp_no); [fb_out,fb_am_out,fb_am_out_prop]=do_fb(1); plot_fb( disp_no+1 ,xmessfb1x) fb_out_init=fb_out; fb_am_out_init=fb_am_out; fb_am_out_prop_init=fb_am_out_prop; // push back sp1=sp1bak; ep1=ep1bak; fir_steps=fir_stepsbak; wdat1=wdat1bak; endfunction //------------------------------------------------------------------ // experiment function [y3out_replace]=am_conv0_iso(i_start, j_end, sp1x, ep1x, disp_code, xmessfb) // // force initial one to org's AM , from fb no. i_start to fb no. j_end // global fb_out fb_am_out fb_am_out_prop global sp1 ep1 global wdat1 global fb_out_work1 fb_am_out_work1 fb_am_out_prop_work1 // force area is shorter one which itital or org sa0i=size(fb_am_out_init); sa1i=sa0i(2); sa0o=size(fb_am_out_org); sa1o=sa0o(2); if sa1i > sa1o then sa1= sa1o; disp(' org data qty is less than init one.'); else sa1= sa1i; end amp0=zeros(1,sa1); // copy init sb0=size(fb_out_init); fb_out_work1=zeros(sb0(1),sb0(2)); for vloop=1:sb0(1) for vloop2=1:sb0(2) fb_out_work1(vloop,vloop2)=fb_out_init(vloop,vloop2); end end if i_start > j_end then // when i_start is more thsn j_end, nothing done else for vloop=i_start:j_end // find max mp0=1; mv0=fb_am_out_init(vloop,1); for v=1:sa1 if fb_am_out_init(vloop,v) > mv0 then mp0=v; mv0=fb_am_out_init(vloop,v); end end // wstr1=sprintf("no. %d , max point is %d, max value is %f", vloop, mp0,mv0); disp( wstr1 ); ax0=fb_am_out_org(vloop,mp0); wstr1=sprintf("the value of kizyun is %f", ax0); disp( wstr1 ); if (( ax0 > 0.) & ( mv0 > 0.)) then ax1= mv0 / ax0; for v=1:sa1 if fb_am_out_init(vloop,v) <= 0.000001 then amp0(v)=1.; else amp0(v)=fb_am_out_org(vloop,v) * ax1 / fb_am_out_init(vloop,v); end end // for v=1:sa1 fb_out_work1(vloop,v)= amp0(v) * fb_out_work1(vloop,v); end end end // for vloop=i_start:j_end // end //if i_start >= j_end then sp1=sp1x; if (ep1x-sp1x+1) < sa1 then ep1=ep1x; else ep1=sp1x+sa1-1; end // ------------------------------- // Fade OUT and Fade IN control // fade_in_legnth = 150; fade_out_legnth = 250; offset=0; // y3out_replace=y3out; for v=(1+offset):(sa1-offset) val0=0.; for w=1:n_stepf val0=val0+ fb_out_work1(w,v); end if v < fade_in_legnth then alfa= v / fade_in_legnth; y3out_replace(sp1+v-1)= alfa * val0 + (1. - alfa) * y3out(sp1+v-1); elseif (sa1-v) < fade_out_legnth then alfa= (sa1 - v ) / fade_out_legnth; y3out_replace(sp1+v-1)= alfa * val0 + (1. - alfa) * y3out(sp1+v-1); else y3out_replace(sp1+v-1)= val0; end end if disp_code > 0 then w0=xget('window'); // stack old window if disp_code >= 1 then xset('window',disp_code); // create new windows clf(); sx=size(y3out_replace); size1=sx(2); wsize2=size1 - sp1 + 1; wsize3=ep1-sp1+1; zdata1=zeros(1,width0+1); zdata1p=zeros(1,width0+1); ydata1=zeros(1,width0+1); for vloop=1:n_stepf if wsize3 > (width0+1) then wstep1= (ep1-sp1) / width0; for v=1:(width0+1) zdat1(v)= fb_out_work1(vloop,1 + int(wstep1 * (v - 1))); //zdat1p(v)= fb_am_out(vloop, 1 + int(wstep1 * (v - 1))); xziku1(v)= sp1 + int(wstep1 * (v - 1)); if vloop==1 then ydat1(v)= y3out_replace( sp1 + int(wstep1 * (v - 1))); end end else for v=1:(width0+1) if v <= wsize3 then zdat1(v)=fb_out_work1(vloop,v); //zdat1p(v)=fb_am_out(vloop,v); if vloop==1 then ydat1(v)= y3out_replace((sp1-1)+v); end else zdat1(v)=0.; zdat1p(v)=0.; if vloop==1 then ydat1(v)= 0.; end end xziku1(v)=(sp1-1)+v; end end if vloop==1 then subplot((n_stepf+1),1,1); plot(xziku1,ydat1,'b'); wstr1=sprintf('waveform selected, '); wstr1=wstr1 + ' / ' + xmessfb; xtitle(wstr1); end subplot((n_stepf+1),1,(vloop+1)); plot(xziku1,zdat1,'b'); //plot(xziku1,zdat1p,'g'); wstr1=sprintf('Band Pass filtered waveform selected (%dHz- %dHz)',int(bpf_fb0(vloop,2)),int(bpf_fb0(vloop,3))); //mv0=fb_am_out_prop(vloop,1); //if mv0 > 0. then // mv1= int(20. * log10(mv0)); //else // mv1=-120; //end //pk0=fb_am_out_prop(vloop,2); //if pk0 > 0. then // pk1= int( 20. * log10(pk0)); //else // pk1=-120; //end //wstr2=sprintf('%d/%d dB',mv1,pk1); // average / peak [dB] //xtitle( wstr1,'',wstr2 ); end // for vloop=1:n_stepf xset('window',disp_code+1); // create new windows clf(); subplot(211); plot(y3out,'b'); xtitle('y3out' ,'Samples (Time)', 'Amplitude'); subplot(212); plot(y3out_replace,'r'); xtitle('y3out_replace' ,'Samples (Time)', 'Amplitude'); end // if disp_code >= 1 then xset('window',w0); // push old windows end endfunction //am_conv0( 1, 3) // //----------------------------------------------------------------- function [wdat2,rslt]=do_replace(disp0) rslt=ones(n_stepf,2); for v=1:n_stepf rslt(v,1)=v; end s34=size(wdat1); size1=s34(2); wdat2=zeros(1,size1); for v=1:size1 wdat2(v)=wdat1(v); end //+replacement for v=sp1:ep1 vs0=0.; for v2=1:n_stepf if rslt(v2,2)==1 then // if the fb no is Combine element vs0= vs0 + fb_out(v2, (v-sp1+1)) end end wdat2(v)=vs0; end //-replacement if disp0 >= 1 then ydat1r=zeros(1,(width0+1)); if size1 > (width0+1) then wstep1= size1 / (width0+1); for v=1:(width0+1) ydat1r(v)= wdat2( 1 + int(wstep1 * (v - 1))); // xziku1(v)= 1 + int(wstep1 * (v - 1)); end else for v=1:(width0+1) if v <= size1 then ydat1r(v)=wdat2(v); else ydat1r(v)=0.; end // xziku1(v)=(sp1-1)+v; end end // clf(); subplot(211); plot(xziku0,ydat0,'b'); //plot2d(xziku0,ydat0,style=[color("blue")]); xtitle('waveform all'); // When linux scilab-3.1 2nd subplot does not work well. subplot(212); plot(xziku0,ydat1r,'b'); //plot2d(xziku1,ydat1,style=[color("blue")]); xtitle('replacement waveform'); end // id disp0 >= 1 thne endfunction // //================================================================== // // +MAIN (fb)program starts // if T_DEMO==1 then if SEL_CODE == 1 then // ka_sample sp1bak=sp1; ep1bak=ep1; fir_stepsbak=fir_steps; // select portion of brust, un-brust of ka_sample sp1=sp1s(1); ep1=ep1s(2); if fir_steps > sp1 then if (int(sp1 / 2) * 2) <> sp1 then fir_steps=sp1; else fir_steps=sp1-1; end str1=sprintf('%d',fir_steps); str1="waring: re-set fir_steps is " + str1; disp( str1 ); end // [st_freq, end_freq, n_stepf, fir_steps, kwindows, end_freq2]= set_fb_para(); //(1) ////st_freq=500.; ////end_freq=3000.; ////end_freq2=0.; // no addtional bpf ////n_stepf=4; ////[bpf_fb0,wft]=make_bpf_fb0(30); ////[fb_out,fb_am_out,fb_am_out_prop]=do_fb(1); ////xmessfb=' ka_sample original (500-3000/4)'; ////plot_fb(31,xmessfb) //(2) st_freq=500.; end_freq=1500.; end_freq2=3000.; // addtional bpf n_stepf=6; //including addtional bpf [bpf_fb0,wft]=make_bpf_fb0(32); if LOAD_FB_FLAG == 1 then [fb_out_org, fb_am_out_org, fb_am_out_prop_org]=load_fb(); fb_out=fb_out_org; fb_am_out=fb_am_out_org; fb_am_out_prop=fb_am_out_prop_org; xmessfb=' ka_sample original (500-1500/5)'; plot_fb(33,xmessfb); else [fb_out,fb_am_out,fb_am_out_prop]=do_fb(1); xmessfb=' ka_sample original (500-1500/5)'; plot_fb(33,xmessfb); // cpy backup to *_org fb_out_org=fb_out; fb_am_out_org=fb_am_out; fb_am_out_prop_org=fb_am_out_prop; if LOAD_FB_FLAG == -1 then save_fb(); end end //(test3) //[yd_lpf, xd_lpf]=set_lpf(3000. * 2. * %pi , 0.7 ,2 ); //s15=size(wdat1); //[xy_ran2s]=make_noise_lpf( s15(2)); //wdat1bak=wdat1; //for vloop=1:s15(2) // wdat1(vloop)=xy_ran2s(vloop); //end //[fb_out,fb_am_out,fb_am_out_prop]=do_fb(1); //xmessfb=' random noise with lpf 3kHz (500-1500/5)'; //plot_fb(34,xmessfb) ////--- // set new fb coeffient.... //[fb_out,fb_am_out,fb_am_out_prop]=do_fb(1); //xmessfb=' random noise with lpf 3kHz (500-3000/4)'; //plot_fb(35,xmessfb) ////-- //wdat1=wdat1bak; // //end of (test3) // push back sp1=sp1bak; ep1=ep1bak; fir_steps=fir_stepsbak; end end // // -MAIN (fb)program starts // //==============================================================++= // // // global variables No. 2 // c0=35000; // constant. the speed of sound is round 35000 cm/second ts=1.0 / fs1; // L1=1; // this will be overwritten. L2=1; // this will be overwritten. L3=1; // this will be overwritten. A1=1; // this will be overwritten. A2=1; // this will be overwritten. A3=1; // this will be overwritten. // LL=zeros(1,10); // this will be overwritten. // ttl_Length=18.5; // set total length of combined tubes, Two Tubes, and etc ttl_Lengths=zeros(1,max_frames); ttl_Area=10; // a constant in this program. set total area of combined tubes, Two Tubes, and etc for dummy display l1=0; // this will be overwritten. l2=0; // this will be overwritten. rg=1; // a constant in this program. When glottis r1=0.8; // this will be overwritten. r2=0.8; // this will be overwritten. rl=0.9; // set adjust reduction response. rl is variable, however use one constant in this. fc=1000; // this will be overwritten. set cut off frequency of High Pass Filter by unit is [Hz] trise=6.0; // this will be overwritten. //tfall=0.7; // this will be overwritten. tfall=2.0; // this will be overwritten. tclosed=5.0; // use ONLY in waveform making // l1s=zeros(1,max_frames); l2s=zeros(1,max_frames); r1s=zeros(1,max_frames); r2s=zeros(1,max_frames); rls=zeros(1,max_frames); // frame_lengths=zeros(1,max_frames); amplitudes=zeros(1,max_frames); // yg_res1= ones(1,1024); // this will be overwritten. tube2_res1= ones(1,1024); // this will be overwritten. hpf_res1= ones(1,1024); // this will be overwritten. ind_noise_res1= zeros(1,1024); // this will be overwritten. stable_noise_res1= zeros(1,1024); // this will be overwritten. burst_noise_res1= zeros(1,1024); // this will be overwritten. peq_res1= zeros(1,1024); // this will be overwritten. overall_res1=zeros(1,1024); // this will be overwritten. db_2tube=zeros(1,1024); // this will be overwritten. phi_2tube=zeros(1,1024); // this will be overwritten. // // others WT_QTY=2; // WT_QTY=2: standard 2 tubes model WT_QTY0=WT_QTY; // this is pop and push WT_QTYs=zeros(1,max_frames); N_REPEAT=7; // how many section to generate waveform ? Wtubepropdat= 'tubeprop.dat'; // //--------------------------------------------------------------------------- // "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. // set threshold for independently output point threshold_occur_minus_id = -0.3; // set threshold of Volume Velocity of minus side when noise occurs threshold_occur_plus_id = 0.3; // set threshold of Volume Velocity of plus side when noise occurs // two side threshold, minus side and plus side, is kunikuno-saku. mix_amplitude_rl=0.05; // set mix amplitude at "rl" point mix_amplitude_1rl=0.05; // set mix amplitude at "1+rl" point mix_amplitude_indep=0.08; // set mix amplitude at independently output point. // When i_type is 4, nA0 will be used instead of mix_amplitude_indep //frq_center= 3000; // set center frequency of noise by unit is [Hz] //freq_band= 2000; // set frequency band of noise by unit is [Hz] //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ // // // Simplest independent noise parameters representation // //++++ basic definition +++++ noise_waku_area0=4.; // [cm^2] frq_center0=2500.; // the center frequecny when noise wake manseki is noise_wake_area0 i2nd_thd_factor=0.; // dai 2 koutyouha no original nitaisuru wariai // ---> changed to 2nd peak no frequecny ga nanbai ka ni henkou, supposing both base and 2nd are same high peak value // q0_0=5.; // iir bpf no Q this value is OK?? // noise_waku_area=0.; // dummy set frq_center=0.; // dummy set q0=0.; // dummy set // function [frq_center,q0]= trans_waku_area( noise_waku_area) if noise_waku_area > 0. then frq_center = frq_center0 * (noise_waku_area0 / noise_waku_area); q0=q0_0; else frq_center=0.; q0=0.; end endfunction // // //++++++++++++++++++++++++++++ //frq_center=2500.; // noise center frequecny noise_waku_area=noise_waku_area0; [frq_center,q0]= trans_waku_area( noise_waku_area); freq_band=2000.; // noise band nA0=mix_amplitude_indep; // noise amplitude // noise_waku_areas=zeros(1,max_frames); i2nd_thd_factors=zeros(1,max_frames); nA0s=zeros(1,max_frames); // ... And The range of existing independent noise low_edge_fc=1500.; // low_edge is most yokonagani kutiwo tubometa toki low_edge_band=2500.; high_edge_fc=5000.; // high_edge is most maruku kutiwo tubometa toki high_edge_band=3000.; //++++++++++++++++++++++++++++++++ function [emp_frq_center,emp_q0]= trans_space_area( emp_space_area) if emp_space_area > 0. then emp_frq_center = emp_frq_center0 * (emp_space_area0 / emp_space_area); emp_q0=emp_q0_0; else emp_frq_center=0.; emp_q0=0.; end endfunction //++++ basic definition +++++ emp_space_area0=4.; // [cm^2] emp_frq_center0=2500.; // the center frequecny of emphasised frequency band emp_q0_0=4.; emp_space_area=emp_space_area0; emp_space_areas=zeros(1,max_frames); emp_frq_center=emp_frq_center0; // dummy set emp_q0=emp_q0_0; // dummy set [emp_frq_center,emp_q0]= trans_space_area( emp_space_area); empA0=20; // emphasis amplitude empA0s=zeros(1,max_frames); // ... And The range of existing emphasis emp_low_edge_fc=1500.; // low_edge emp_high_edge_fc=5000.; // high_edge emp_low_A0=0.; emp_high_A0=30.; // max boost //+++++ stable noise ++++++ stable_noise_lpf=2000.0; // LPF fc for stable noise band stable_noise_lpfs=zeros(1,max_frames); q0_stable_noise=0.7; xorder_stable_noise=2; // stable_noise_hpf=1000.0; q0_stable_noiseh=0.7; xorder_stable_noiseh=1; //----- stable noise ------ // //++++ for burst noise ++++ ////stable_noise_lpf=2000.0; // band center frequency is stable_noise_lpf. q0_burst_noise=2.0; //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ // nstep_fir=251; // 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 // // // //===================================================================================== // xd_bpf=zeros(3,1); // this will be overwritten. yd_bpf=zeros(3,1); // this will be overwritten. xd_lpf=zeros(3,1); // this will be overwritten. yd_lpf=zeros(3,1); // this will be overwritten. xd_hpf=zeros(3,1); // this will be overwritten. yd_hpf=zeros(3,1); // this will be overwritten. iir_bpf1_res1=zeros(1,1); // this will be overwritten. iir_bpf2_res1=zeros(1,1); // this will be overwritten. iir_lpf1_res1=zeros(1,1); // this will be overwritten. iir_hpf1_res1=zeros(1,1); // this will be overwritten. db_iir1=zeros(1,1); // this will be overwritten. phi_iir1=zeros(1,1); // this will be overwritten. db_iir2=zeros(1,1); // this will be overwritten. phi_iir2=zeros(1,1); // this will be overwritten. // xd_peq=zeros(3,1); // this will be overwritten. yd_peq=zeros(3,1); // this will be overwritten. iir_peq1_res1=zeros(1,1); // this will be overwritten. iir_peq2_res1=zeros(1,1); // this will be overwritten. // function save_2tube() save(Wtubepropdat,L1,L2,L3,A1,A2,A2,ttl_Length,ttl_Area,l1,l2,rg,r1,r2,fc,trise,tfall,tclosed,yg_res1,tube2_res1,hpf_res1,overall_res1,db_2tube,phi_2tube,WT_QTY); endfunction function [L1,L2,L3,A1,A2,A3,ttl_Length,ttl_Area,l1,l2,rg,r1,r2,fc,trise,tfall,tclosed,yg_res1,tube2_res1,hpf_res1,overall_res1,db_2tube,phi_2tube,WT_QTY]=load_2tube() load(Wtubepropdat,'L1','L2','L3','A1','A2','A3','ttl_Length','ttl_Area','l1','l2','rg','r1','r2','fc','trise','tfall','tclosed','yg_res1','tube2_res1','hpf_res1','overall_res1','db_2tube','phi_2tube','WT_QTY'); endfunction // function [WT_QTY]= set_tube_model() l1=list('model',(WT_QTY-1),['2 tubes','3 tubes serial type','2 tubes plus independent noise','2 tubes plus partly emphasis','2 tubes in stable noise','2 tubes in brust noise']); //l1=list('model',1,['2 tubes']); wrep=x_choices('Select tube model',list(l1)); // if wrep== 1 then WT_QTY=2; elseif wrep== 2 then WT_QTY=3; elseif wrep== 3 then WT_QTY=21; elseif wrep== 4 then WT_QTY=41; elseif wrep== 5 then WT_QTY=51; elseif wrep== 6 then WT_QTY=52; end endfunction // // //------------------------------------------------- function [fc,trise,tfall,tclosed]= set_2tube_para() //function [r1,r2,l1,l2,ttl_Length,rl,fc,trise,tfall,tclosed, noise_waku_area,i2nd_thd_factor,nA0,emp_space_area, empA0]= set_2tube_para() global r1s r2s l1s l2s ttl_Lengths rls noise_waku_areas i2nd_thd_factors nA0s global emp_space_areas empA0s global stable_noise_lpfs wstr1=sprintf('%f',r1s(tframe)); wstr2=sprintf('%f',l1s(tframe)); wstr3=sprintf('%f',ttl_Lengths(tframe)); wstr4=sprintf('%f',fc); wstr5=sprintf('%f',trise); wstr6=sprintf('%f',tfall); wstr7=sprintf('%f',tclosed); wstr8=sprintf('%f',r2s(tframe)); wstr9=sprintf('%f',l2s(tframe)); wstr10=sprintf('%f',rls(tframe)); wstr11=sprintf('%f',noise_waku_areas(tframe)); wstr12=sprintf('%f',i2nd_thd_factors(tframe)); wstr13=sprintf('%f',nA0s(tframe)); wstr14=sprintf('%f',emp_space_areas(tframe)); wstr15=sprintf('%f',empA0s(tframe)); wstr16=sprintf('%f',stable_noise_lpfs(tframe)); wstr17='0'; // dummy data // if WT_QTY == 3 then txt1=['r1';'l1'; 'r2'; 'l2'; 'total_tube_Length'; 'rl' ; 'fc'; 'trise'; 'tfall' ; 'tclosed (use only for generation)']; sig1=x_mdialog('Input parameters',txt1, [wstr1 ; wstr2 ; wstr8 ; wstr9; wstr3 ; wstr10 ; wstr4; wstr5; wstr6; wstr7 ]); if sig1==[] then fc=evstr(wstr4); trise=evstr(wstr5); tfall=evstr(wstr6); tclosed=evstr(wstr7); else r1s(tframe)=evstr(sig1(1)); l1s(tframe)=evstr(sig1(2)); r2s(tframe)=evstr(sig1(3)); l2s(tframe)=evstr(sig1(4)); ttl_Lengths(tframe)=evstr(sig1(5)); rls(tframe)=evstr(sig1(6)); fc=evstr(sig1(7)); trise=evstr(sig1(8)); tfall=evstr(sig1(9)); tclosed=evstr(sig1(10)); end else txt1=['r1';'l1'; 'total_tube_Length'; 'rl' ; 'fc'; 'trise'; 'tfall' ; 'tclosed (use only for generation)']; sig1=x_mdialog('Input parameters',txt1, [wstr1 ; wstr2 ; wstr3 ; wstr10; wstr4; wstr5; wstr6; wstr7 ]); if sig1==[] then fc=evstr(wstr4); trise=evstr(wstr5); tfall=evstr(wstr6); tclosed=evstr(wstr7); else r1s(tframe)=evstr(sig1(1)); l1s(tframe)=evstr(sig1(2)); ttl_Lengths(tframe)=evstr(sig1(3)); rls(tframe)=evstr(sig1(4)); fc=evstr(sig1(5)); trise=evstr(sig1(6)); tfall=evstr(sig1(7)); tclosed=evstr(sig1(8)); end end if WT_QTY == 21 then txt1=['noise wake area(default=4)';'2nd/base freq factor(1-2)';'amplitude']; sig1=x_mdialog('Input independent noise parameters',txt1, [wstr11 ; wstr12 ; wstr13 ]); if sig1==[] then // else noise_waku_areas(tframe)=evstr(sig1(1)); i2nd_thd_factors(tframe)=evstr(sig1(2)); nA0s(tframe)=evstr(sig1(3)); end end if WT_QTY == 41 then txt1=['emphasis space area(default=4)'; 'amplitude']; sig1=x_mdialog('Input partly emphasis parameters',txt1, [wstr14 ; wstr15 ]); if sig1==[] then // else emp_space_areas(tframe)=evstr(sig1(1)); empA0s(tframe)=evstr(sig1(2)); end end if (WT_QTY == 51) | (WT_QTY == 52) then txt1=['stable noise band (LPF fc) or center freq for burst noise'; '-']; sig1=x_mdialog('Input parameter',txt1, [wstr16 ; wstr17 ]); if sig1==[] then // else stable_noise_lpfs(tframe)=evstr(sig1(1)); // empA0s(tframe)=evstr(sig1(2)); end end endfunction // // Function plot sqrt(area) vs length function plot_area(disp0) wb0=xget('window'); // stack old window xset('window',disp0); // create new windows // if f412 == 1 then arrowsize1=20; else arrowsize1=1; end clf(); plot2d(0,0,1,"010","",[-5,-10,30,10],[5,5,4,5]); // wake wo egaku if WT_QTY == 4 then //Three Tubes Model of T type ya1=sqrt(A1); ya2=sqrt(A2); ya3=sqrt(A3); 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; xx1=L1 - (ya3 / 2); xx2=L1 + (ya3 / 2); xp=[0 0 xx1 xx1 xx2 xx2 L1+L2 L1+L2 L1 L1 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 ) ]; xarrows(xar, yar, arrowsize1, 5); yar=[yy0 + ( ya2 / 2 ) yy0+ ( ya2 / 2 ) ]; xar=[ L1+L2+0.5 L1+L2+2]; xarrows(xar, yar, arrowsize1, 5); xstring(xx1, yy3 - 1, "CLOSED"); xstring(0,-9,"Attension: Area scale is double than others Tube Models"); // elseif WT_QTY == 3 then //Three Tubes Model of serial type ya1=sqrt(A1); ya2=sqrt(A2); ya3=sqrt(A3); xp=[0 0 L1 L1 L1+L2 L1+L2 L1+L2+L3 L1+L2+L3 L1+L2 L1+L2 L1 L1 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]; xarrows(xar, yar, arrowsize1, 5); xar=[ L1+L2+L3+0.5 L1+L2+L3+2]; xarrows(xar, yar, arrowsize1, 5); else //Two Tubes Model ya1=sqrt(A1 ); ya2=sqrt(A2 ); xp=[0 0 L1 L1 L1+L2 L1+L2 L1 L1 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, arrowsize1, 5); xar=[ L1+L2+0.5 L1+L2+2]; xarrows(xar, yar, arrowsize1, 5); if WT_QTY==5 then xar=[ L1+L2-0.5 L1+L2-2]; xarrows(xar, yar, arrowsize1, 3); end end xset('window',wb0); // push old windows endfunction // // Function plot sqrt(area) vs length function plot_areas(L1z, L2z, L3z, A1z, A2z,A3z,WT_QTY,N_REPEAT,disp0) wb0=xget('window'); // stack old window xset('window',disp0); // create new windows // if f412 == 1 then arrowsize1=20; else arrowsize1=1; end clf(); for v=1:N_REPEAT subplot(N_REPEAT,1,v); L1=L1z(v); L2=L2z(v); L3=L3z(v); A1=A1z(v); A2=A2z(v); A3=A3z(v); if v == 1 then wstr1 = ' No. ' + string(v) + ' Tube Model, Area and Length, ' + mess_init0; else wstr1 = ' No. ' + string(v); end plot2d(0,0,1,"010","",[-5,-10,30,10],[5,5,4,5]); // wake wo egaku xstring(-4,6, wstr1); if WT_QTY == 4 then //Three Tubes Model of T type ya1=sqrt(A1); ya2=sqrt(A2); ya3=sqrt(A3); 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; xx1=L1 - (ya3 / 2); xx2=L1 + (ya3 / 2); xp=[0 0 xx1 xx1 xx2 xx2 L1+L2 L1+L2 L1 L1 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 ) ]; xarrows(xar, yar, arrowsize1, 5); yar=[yy0 + ( ya2 / 2 ) yy0+ ( ya2 / 2 ) ]; xar=[ L1+L2+0.5 L1+L2+2]; xarrows(xar, yar, arrowsize1, 5); xstring(xx1, yy3 - 1, "CLOSED"); xstring(0,-9,"Attension: Area scale is double than others Tube Models"); // elseif WT_QTY == 3 then //Three Tubes Model of serial type ya1=sqrt(A1); ya2=sqrt(A2); ya3=sqrt(A3); xp=[0 0 L1 L1 L1+L2 L1+L2 L1+L2+L3 L1+L2+L3 L1+L2 L1+L2 L1 L1 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]; xarrows(xar, yar, arrowsize1, 5); xar=[ L1+L2+L3+0.5 L1+L2+L3+2]; xarrows(xar, yar, arrowsize1, 5); else //Two Tubes Model ya1=sqrt(A1 ); ya2=sqrt(A2 ); xp=[0 0 L1 L1 L1+L2 L1+L2 L1 L1 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, arrowsize1, 5); xar=[ L1+L2+0.5 L1+L2+2]; xarrows(xar, yar, arrowsize1, 5); if WT_QTY==5 then xar=[ L1+L2-0.5 L1+L2-2]; xarrows(xar, yar, arrowsize1, 3); end end end // for v=1:N_REPEAT xset('window',wb0); // push old windows endfunction // //------------------------------------------------- // +hpf function [hpf_res1]=hpf_response( point0) wc=2. * %pi * fc; al1= -1. * %e^( -1. * wc * ts); bl0= wc * ts; // Then, convert them to high pass filter's coefficients fnew=fc; wcnew=2. * %pi * fnew; alfa= -1. * cos( (wc + wcnew) * ts /2. ) / cos( (wc - wcnew) * 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); // Function for disply high pass filter's frequency-phase response [frq1,sfrq1,is1,ie1]=set_frq(point0); hpf_res1=zeros(1,sfrq1); for v=1:sfrq1 xw=2 * %pi * frq1(v); yi= bh0 + bh1 * cos( xw * ts) - bh1 * sin( xw * ts ) * %i; yb=1.0 + ah1 * cos( xw * ts ) - ah1 * sin( xw * ts ) * %i; hpf_res1(v)= yi /yb; end endfunction // -hpf //---------------------------------------------------------- // // function [ tube2_res1, L1, L2, L3, A1, A2, A3]= two_tubes2_resp() if WT_QTY==3 then // // 3 tubes model definition // // L1+L2+L3 = ttl_Length (total tubes length) // // l1= (L2-L1)/(L2+L1) // l2= (L3-L2)/(L3+L2) // // A1+A2+A3= ttl_Area // // r1= (A2-A1)/(A2+A1) // r2= (A3-A2)/(A3+A2) // bunbo= l1 * l2 + ( l1 - l2 ) + 3.0; L1= ttl_Length * ( l2 - 1.0 ) * (l1 - 1.0 ) / bunbo ; L2= ttl_Length * ( l2 - 1.0 ) * ( -1.0 * l1 - 1.0 ) / bunbo ; L3= ttl_Length * ( l2 + 1.0 ) * (l1 + 1.0 ) / bunbo ; tu1=L1/c0; tu2=L2/c0; tu3=L3/c0; bunbo= r1 * r2 + ( r1 - r2 ) + 3.0; A1= ttl_Area * ( r2 - 1.0 ) * (r1 - 1.0 ) / bunbo ; A2= ttl_Area * ( r2 - 1.0 ) * ( -1.0 * r1 - 1.0 ) / bunbo ; A3= ttl_Area * ( r2 + 1.0 ) * (r1 + 1.0 ) / bunbo ; [frq1,sfrq1,is1,ie1]=set_frq(1024); tube2_res1=zeros(1,sfrq1); for v=1:sfrq1 xw= 2 * %pi * frq1(v); //+ yi = 0.5 * ( 1. + rg ) * ( 1. + r1) * ( 1. + r2) * ( 1. + rl ) * %e^( -1. * ( tu1 + tu2 + tu3) * xw * %i); yb = 1 + rg * r1 * %e^( -2 * tu1 * xw * %i) + r1 * r2 * %e^( -2 * tu2 * xw * %i) + r2 * rl * %e^( -2 * tu3 * xw * %i); yb = yb + rg * r2 * %e^( -2 * (tu1 + tu2) * xw * %i) + r1 * rl * %e^( -2 * (tu2 + tu3) * xw * %i); yb = yb + rg * r1 * r2 * rl * %e^( -2 * (tu1 + tu3) * xw * %i); yb = yb + rg * rl * %e^( -2 * (tu1 + tu2 + tu3) * xw * %i); tube2_res1(v) = yi / yb; //- end else L2= ( 1. + l1 ) * ttl_Length / 2.; L1= ttl_Length - L2; tu1=L1/c0; tu2=L2/c0; A2= ( 1. + r1 ) * ttl_Area/ 2.; A1= ttl_Area - A2; // [frq1,sfrq1,is1,ie1]=set_frq(1024); tube2_res1=zeros(1,sfrq1); for v=1:sfrq1 xw= 2 * %pi * frq1(v); yi = 0.5 * ( 1 + rg ) * ( 1 + r1) * ( 1 + rl ) * %e^( -1 * ( tu1 + tu2 ) * xw * %i); yb = 1 + r1 * rg * %e^( -2 * tu1 * xw * %i) + r1 * rl * %e^( -2 * tu2 * xw * %i) + rl * rg * %e^( -2 * (tu1 + tu2) * xw * %i); tube2_res1(v) = yi / yb; end // +dummy data set L3=1; A3=1; // -dummy data set end endfunction // // +yg response function [yg_res1]=yg_resp(points) N2= int( trise / 1000 / ts); N3= int( tfall / 1000 / ts); sizew=max((N2+N3+1),points); yg=zeros(1,sizew); // for tc0=1:sizew if tc0 <= (N2+1) then yg(tc0)= 0.5 * ( 1 - cos( ( %pi / N2 ) * (tc0 - 1) ) ); elseif tc0<= (N2+N3+1) then yg(tc0)= cos( ( %pi / ( 2 * N3 )) * ( tc0 - (N2+1) ) ); else yg(tc0)=0.0; end end [frq1,sfrq1,is1,ie1]=set_frq(points); yg_res1=zeros(1,sfrq1); // +do points(1024) fft and get portion of data if sizew == points then fftwout=fft( yg, -1 ); v=1; is2= int(is1 * points / fft_point1); ie2= int(ie1 * points / fft_point1); for w=is2:ie2 yg_res1(v)=fftwout(is2+v) / fftwout(1); v=v+1; end else // -do 1024 fft and get portion of data // y00=0.; for v=1:(N2+N3+1) y00= y00 + yg(v); end // for w=1:sfrq1 xw=2 * %pi * frq1(w) * ts; yi=0; for v=1:(N2+N3+1) yi= yi + yg(v) * %e^(-1. * xw * (v-1) * %i); end // for v yg_res1(w)=yi / y00; end // for w end endfunction // -yg response //------------------------------------------------- // +y2tmx response function [respw,db_fft1,phi_fft1]=y2tmx_resp(y2tmx,start_frame,points) sp1=1; for v=1:start_frame sp1=sp1+int(LL(v)); end win_hn=window('hn',points); // make hanning windows data wdatw=zeros(1, points); for v=1:points wdatw(v)=y2tmx(sp1 + v - 1); end wdatw2= wdatw .* win_hn; fftwout=fft( wdatw2, -1 ); [frq1,sfrq1,is1,ie1]=set_frq_fft(points); 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); endfunction // -y2tmx response // //------------------------------------------------- // +noise function [ind_noise_res1]=noise_response( point0) [frq1,sfrq1,is1,ie1]=set_frq(point0); ind_noise_res1=zeros(1,sfrq1); if WT_QTY==21 then [iir_bpf1_res1,db_iir1,phi_iir1] = iir_bpf_resp(1,point0); [iir_bpf2_res1,db_iir2,phi_iir2] = iir_bpf_resp(2,point0); for v=1:sfrq1 ind_noise_res1(v)=iir_bpf1_res1(v) + iir_bpf2_res1(v); end plot_iir_bpf1(db_iir1,phi_iir1,db_iir2,phi_iir2, 30 ); // disp0=30 end endfunction // -noise //---------------------------------------------------------- // // function plot_fft_y2tmx(db_fft1x,phi_fft1x,points,disp0) wb0=xget('window'); // stack old window xset('window',disp0); // create new windows [frq1,sfrq1,is1,ie1]=set_frq_fft(points); clf(); subplot(211); gainplot(frq1,db_fft1x,phi_fft1x); xtitle('frequency response of y2tmx (y2tm_noise) by fft analysis'); //subplot(212); ////plot2d(xziku1,ydat1,style=[color("blue")]); ////xtitle('waveform selected'); xset('window',wb0); // push old windows endfunction // //------------------------------------------------- function [overall_res1,db_2tube, phi_2tube ]=overall_response() [frq1,sfrq1,is1,ie1]=set_frq(1024); for v=1:sfrq1 if WT_QTY==21 then overall_res1(v)= hpf_res1(v) * ( tube2_res1(v) * yg_res1(v) + ind_noise_res1(v)); elseif WT_QTY==41 then overall_res1(v)= hpf_res1(v) * tube2_res1(v) * peq_res1(v) * yg_res1(v) ; elseif WT_QTY==51 then overall_res1(v)= hpf_res1(v) * tube2_res1(v) * stable_noise_res1(v) ; elseif WT_QTY==52 then overall_res1(v)= hpf_res1(v) * tube2_res1(v) * burst_noise_res1(v) ; else overall_res1(v)= hpf_res1(v) * tube2_res1(v) * yg_res1(v) ; end end [db_2tube,phi_2tube] =dbphi(overall_res1); endfunction //----------------------------------------------- function plot_2tube(disp0, kcode) // When kcode =0, normal, beside kcode=1, leastsq result wb0=xget('window'); // stack old window xset('window',disp0); // create new windows [frq1,sfrq1,is1,ie1]=set_frq(0); clf(); subplot(211); s23=size(db_fft1s); db_fft0=zeros(1,s23(2)); phi_fft0=zeros(1,s23(2)); v=tframe; for w=1:s23(2) db_fft0(1,w)=db_fft1s(v,w); phi_fft0(1,w)=phi_fft1s(v,w); end gainplot(frq1,db_fft0,phi_fft0); wstr1=PART_LIST0(tframe) + ' frequency response of waveform selected by fft analysis' xtitle(wstr1); [frq1,sfrq1,is1,ie1]=set_frq(1024); subplot(212); cal_overall_response(); gainplot(frq1,db_2tube',phi_2tube'); if kcode == 0 then xtitle('frequency response of tubes model by initial value for comparison to one of waveform selected by fft analysis '); elseif kcode == 1 then xtitle('frequency response of tubes model by leastsq method result from the initial value'); elseif kcode == 3 then wstr0=' frequency response of tubes model'; wstr1=sprintf('%f',l1); wstr2=sprintf('%f',r1); wstr3=sprintf('%f',l2); wstr4=sprintf('%f',r2); if WT_QTY == 3 then wstr5='l1=' + wstr1 + ' r1=' + wstr2 + ' l2=' + wstr3 + ' r2=' + wstr4 + wstr0; else // when WT_QTY == 2 then wstr5='l1=' + wstr1 + ' r1=' + wstr2 + wstr0; end xtitle(wstr5); end xset('window',wb0); // push old windows endfunction //================================================================= // // band pass filter of linear-phase FIR filter // 0 soffset then xy_ran2(loop-soffset)=yx(1,1); end end endfunction //---------------------------------------------------- function [xy_ran2]=make_noise_lpf_hpf(xlength0) //* G(z)= y[0]z0 + y[1]z-1 + y[2]z-2/ 1.0 - (x[1]z-1 + x[2]z-2) */ rand('normal'); // random generator is set to a Gaussian soffset=100; // data start offset xlength1= xlength0 + soffset; xy_ran1=rand(1, xlength1); xy_ran2=zeros(1, xlength0); yx=zeros(2,6); yx2=zeros(2,6); disp('---please wait for some time. noise calculating (hpf/lpf)'); for loop=1:xlength1 for vl=1:1 // one 1 filter // lpf yx(vl,4)=xy_ran1(loop); x0=yd_lpf(vl,1) * yx(vl,4) + yd_lpf(vl,2) * yx(vl,5) + yd_lpf(vl,3) * yx(vl,6); // x0= yi[eq0][0] * yx[eq0][3] + yi[eq0][1] * yx[eq0][4] + yi[eq0][2] * yx[eq0][5]; yx(vl,1)=x0 + yx(vl,2) * xd_lpf(vl,2) + yx(vl,3) * xd_lpf(vl,3); // yx[eq0][0] = x0 + yx[eq0][1] * xi[eq0][1] + yx[eq0][2] * xi[eq0][2]; // /* shift for next time */ yx(vl,3)=yx(vl,2); // yx[eq0][2]=yx[eq0][1]; yx(vl,2)=yx(vl,1); // yx[eq0][1]=yx[eq0][0]; yx(vl,6)=yx(vl,5); // yx[eq0][5]=yx[eq0][4]; yx(vl,5)=yx(vl,4); // yx[eq0][4]=yx[eq0][3]; // hpf yx2(vl,4)=yx(1,1); x02=yd_hpf(vl,1) * yx2(vl,4) + yd_hpf(vl,2) * yx2(vl,5) + yd_hpf(vl,3) * yx2(vl,6); // x0= yi[eq0][0] * yx[eq0][3] + yi[eq0][1] * yx[eq0][4] + yi[eq0][2] * yx[eq0][5]; yx2(vl,1)=x02 + yx2(vl,2) * xd_hpf(vl,2) + yx2(vl,3) * xd_hpf(vl,3); // yx[eq0][0] = x0 + yx[eq0][1] * xi[eq0][1] + yx[eq0][2] * xi[eq0][2]; // /* shift for next time */ yx2(vl,3)=yx2(vl,2); // yx[eq0][2]=yx[eq0][1]; yx2(vl,2)=yx2(vl,1); // yx[eq0][1]=yx[eq0][0]; yx2(vl,6)=yx2(vl,5); // yx[eq0][5]=yx[eq0][4]; yx2(vl,5)=yx2(vl,4); // yx[eq0][4]=yx[eq0][3]; end if loop > soffset then xy_ran2(loop-soffset)=yx2(1,1); end end endfunction //-------------------------------------------------------------------------- // // iir filter (bpf) design for noise proccess // function [yd_bpf, xd_bpf]=set_bpf(w0, q0, w1, q1) // g0: gain w0: frequency [rad] q0: q // gain g0 always =1. g0=1.0; // yd_bpf=zeros(2,3); xd_bpf=zeros(2,3); yd_bpf(1,1)=1.0; yd_bpf(2,1)=1.0; if (w0 <= 0.) | (w1 <= 0.) then // else //+++ 1st cal // souitizi_henkan t0=1/fs1; wout = ((w0/(2.0 * %pi)) / (1.0 / t0)) * 2.0 * %pi; //* w0-> digital w he */ wout = tan( wout /2.0); w0b = (wout * 2.0 / t0); // set_bpf_ana ya=zeros(1,3); xa=zeros(1,3); ya(1)=0.0; ya(2)= g0 * w0b / q0; ya(3)=0.0; xa(1)=1.0; xa(2)= w0b / q0; xa(3)= w0b * w0b; // set_bpf_dig yd_bpf0=zeros(1,3); xd_bpf0=zeros(1,3); //* cal X */ yd_bpf0(1)= (2.0 * ya(2) / t0) + ya(3); yd_bpf0(2)= 2.0 * ya(3); yd_bpf0(3)= (-2.0 * ya(2) / t0) + ya(3); //* cal Y */ xd_bpf0(1)= (4.0/ (t0 * t0) ) + (2.0 * xa(2) / t0) + xa(3); xd_bpf0(2)= (-8.0 / t0 /t0) + (2.0 * xa(3)); xd_bpf0(3)= (4.0/ (t0 * t0) ) + (-2.0 * xa(2) / t0) + xa(3); vl=1; // force to set xd_bpf(1) 1.0 yd_bpf(vl,1)= yd_bpf0(1)/ xd_bpf0(1); yd_bpf(vl,2)= yd_bpf0(2)/ xd_bpf0(1); yd_bpf(vl,3)= yd_bpf0(3)/ xd_bpf0(1); xd_bpf(vl,2)= -1.0 * xd_bpf0(2)/xd_bpf0(1); xd_bpf(vl,3)= -1.0 * xd_bpf0(3)/xd_bpf0(1); xd_bpf(vl,1)= xd_bpf0(1)/xd_bpf0(1); //+++ 2nd cal // souitizi_henkan t0=1/fs1; wout = ((w1/(2.0 * %pi)) / (1.0 / t0)) * 2.0 * %pi; //* w0-> digital w he */ wout = tan( wout /2.0); w0b = wout * 2.0; // set_bpf_ana ya=zeros(1,3); xa=zeros(1,3); ya(1)=0.0; ya(2)= g0 * w0b / q1; ya(3)=0.0; xa(1)=1.0; xa(2)= w0b / q1; xa(3)= w0b * w0b; // set_bpf_dig yd_bpf0=zeros(1,3); xd_bpf0=zeros(1,3); //* cal X */ yd_bpf0(1)= (2.0 * ya(2) / t0) + ya(3); yd_bpf0(2)= 2.0 * ya(3); yd_bpf0(3)= (-2.0 * ya(2) / t0) + ya(3); //* cal Y */ xd_bpf0(1)= (4.0/ (t0 * t0) ) + (2.0 * xa(2) / t0) + xa(3); xd_bpf0(2)= (-8.0 / t0 /t0) + (2.0 * xa(3)); xd_bpf0(3)= (4.0/ (t0 * t0) ) + (-2.0 * xa(2) / t0) + xa(3); vl=2; // force to set xd_bpf(1) 1.0 yd_bpf(vl,1)= yd_bpf0(1)/ xd_bpf0(1); yd_bpf(vl,2)= yd_bpf0(2)/ xd_bpf0(1); yd_bpf(vl,3)= yd_bpf0(3)/ xd_bpf0(1); xd_bpf(vl,2)= -1.0 * xd_bpf0(2)/xd_bpf0(1); xd_bpf(vl,3)= -1.0 * xd_bpf0(3)/xd_bpf0(1); xd_bpf(vl,1)= xd_bpf0(1)/xd_bpf0(1); end // if (w0 <= 0.) | w(1 <= 0.) 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 y = iir_bpf( xw, vl ) t0=1/fs1; yi = yd_bpf(vl,1) + yd_bpf(vl,2) * %e^( -1.* t0 * xw * %i) + yd_bpf(vl,3) * %e^( -2.* t0 * xw * %i); yb = 1.0 - (xd_bpf(vl,2) * %e^( -1.* t0 * xw * %i) + xd_bpf(vl,3) * %e^( -2.* t0 * xw * %i)); y = nA0 * yi / yb; endfunction function [iir_bpf1,db_iir1,phi_iir1] = iir_bpf_resp(vl,points) [frq1,sfrq1,is1,ie1]=set_frq(points); iir_bpf1=zeros(1,sfrq1); for v=1:sfrq1 xw= 2 * %pi * frq1(v); iir_bpf1(v)=iir_bpf( xw,vl ) end [db_iir1,phi_iir1] =dbphi(iir_bpf1); endfunction function plot_iir_bpf1(db_iir1,phi_iir1,db_iir2,phi_iir2,disp0) [frq1,sfrq1,is1,ie1]=set_frq(1024); w0=xget('window'); xset('window',disp0); clf(); subplot(311); gainplot(frq1,db_iir1,phi_iir1); xtitle('frequency response of iir bpf1'); subplot(312); gainplot(frq1,db_iir2,phi_iir2); xtitle('frequency response of iir bpf2'); subplot(313); gainplot(frq1,(db_iir1+db_iir2),phi_iir2); xtitle('frequency response of both and iir bpf1 and iir bpf2'); xset('window',w0); endfunction //-------------------------------------------------------------------------- // // iir filter (lpf) design for noise proccess // function [yd_lpf, xd_lpf]=set_lpf(w0, q0, xorder) // g0: gain w0: frequency [rad] q0: q // gain g0 always =1. g0=1.0; // yd_lpf=zeros(2,3); xd_lpf=zeros(2,3); yd_lpf(1,1)=1.0; yd_lpf(2,1)=1.0; if w0 <= 0. then else //+++ 1st cal // souitizi_henkan t0=1/fs1; t0x=t0; wout = ((w0/(2.0 * %pi)) / (1.0 / t0)) * 2.0 * %pi; //* w0-> digital w he */ wout = tan( wout /2.0); w0b = (wout * 2.0 / t0); y0=zeros(1,3); x0=zeros(1,3); y=zeros(1,3); x=zeros(1,3); idx=1; if xorder == 1 then // 1-st order lpf y0(0+idx)= 0.0; y0(1+idx)= g0 * w0b; y0(2+idx)= 0.0; x0(0+idx)= 1.0; x0(1+idx)= w0b; x0(2+idx)= 0.0 ; else // 2-nd order lpf y0(0+idx)= 0.0; y0(1+idx)= 0.0; y0(2+idx)= w0b * w0b * g0; x0(0+idx)=1.0; x0(1+idx)= w0b / q0; x0(2+idx)= w0b * w0b; end 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 */ vl=1; yd_lpf(vl,0+idx)=y(0+idx)/x(0+idx); yd_lpf(vl,1+idx)=y(1+idx)/x(0+idx); yd_lpf(vl,2+idx)=y(2+idx)/x(0+idx); xd_lpf(vl,1+idx)= -1.0 * x(1+idx)/x(0+idx); xd_lpf(vl,2+idx)= -1.0 * x(2+idx)/x(0+idx); xd_lpf(vl,0+idx)=x(0+idx)/x(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 = iir_lpf( xw, vl ) t0=1/fs1; yi = yd_lpf(vl,1) + yd_lpf(vl,2) * %e^( -1.* t0 * xw * %i) + yd_lpf(vl,3) * %e^( -2.* t0 * xw * %i); yb = 1.0 - (xd_lpf(vl,2) * %e^( -1.* t0 * xw * %i) + xd_lpf(vl,3) * %e^( -2.* t0 * xw * %i)); y = yi / yb; endfunction function [iir_lpf1,db_iir1,phi_iir1] = iir_lpf_resp(vl,points) [frq1,sfrq1,is1,ie1]=set_frq(points); iir_lpf1=zeros(1,sfrq1); for v=1:sfrq1 xw= 2 * %pi * frq1(v); iir_lpf1(v)=iir_lpf( xw,vl ) end [db_iir1,phi_iir1] =dbphi(iir_lpf1); endfunction function plot_iir_lpf1(db_iir1,phi_iir1,disp0) [frq1,sfrq1,is1,ie1]=set_frq(1024); w0=xget('window'); xset('window',disp0); clf(); subplot(311); gainplot(frq1,db_iir1,phi_iir1); xtitle('frequency response of iir lpf1'); xset('window',w0); endfunction // for debug lpf //w0test=2500. * 2. * %pi; //q0test=0.7; //xordertest=2; //[yd_lpf, xd_lpf]=set_lpf(w0test, q0test, xordertest); //pointstest=1024; //[iir_lpf1,db_iir1,phi_iir1] = iir_lpf_resp(1,pointstest); //plot_iir_lpf1(db_iir1,phi_iir1,21); //-------------------------------------------------------------------------- // // iir filter (hpf) design for noise proccess // function [yd_hpf, xd_hpf]=set_hpf(w0, q0, xorder) // g0: gain w0: frequency [rad] q0: q // gain g0 always =1. g0=1.0; // yd_lpf=zeros(2,3); xd_lpf=zeros(2,3); yd_lpf(1,1)=1.0; yd_lpf(2,1)=1.0; if w0 <= 0. then else //+++ 1st cal // souitizi_henkan t0=1/fs1; t0x=t0; wout = ((w0/(2.0 * %pi)) / (1.0 / t0)) * 2.0 * %pi; //* w0-> digital w he */ wout = tan( wout /2.0); w0b = (wout * 2.0 / t0); y0=zeros(1,3); x0=zeros(1,3); y=zeros(1,3); x=zeros(1,3); idx=1; if xorder == 1 then // 1-st order lpf y0(0+idx)= g0; y0(1+idx)= 0.0; y0(2+idx)= 0.0; x0(0+idx)= 1.0; x0(1+idx)= w0b; x0(2+idx)= 0.0 ; else // 2-nd order lpf y0(0+idx)= g0; y0(1+idx)= 0.0; y0(2+idx)= 0.0; x0(0+idx)=1.0; x0(1+idx)= w0b / q0; x0(2+idx)= w0b * w0b; end 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 */ vl=1; yd_hpf(vl,0+idx)=y(0+idx)/x(0+idx); yd_hpf(vl,1+idx)=y(1+idx)/x(0+idx); yd_hpf(vl,2+idx)=y(2+idx)/x(0+idx); xd_hpf(vl,1+idx)= -1.0 * x(1+idx)/x(0+idx); xd_hpf(vl,2+idx)= -1.0 * x(2+idx)/x(0+idx); xd_hpf(vl,0+idx)=x(0+idx)/x(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 = iir_hpf( xw, vl ) t0=1/fs1; yi = yd_hpf(vl,1) + yd_hpf(vl,2) * %e^( -1.* t0 * xw * %i) + yd_hpf(vl,3) * %e^( -2.* t0 * xw * %i); yb = 1.0 - (xd_hpf(vl,2) * %e^( -1.* t0 * xw * %i) + xd_hpf(vl,3) * %e^( -2.* t0 * xw * %i)); y = yi / yb; endfunction function [iir_hpf1,db_iir1,phi_iir1] = iir_hpf_resp(vl,points) [frq1,sfrq1,is1,ie1]=set_frq(points); iir_hpf1=zeros(1,sfrq1); for v=1:sfrq1 xw= 2 * %pi * frq1(v); iir_hpf1(v)=iir_hpf( xw,vl ) end [db_iir1,phi_iir1] =dbphi(iir_hpf1); endfunction function plot_iir_hpf1(db_iir1,phi_iir1,disp0) [frq1,sfrq1,is1,ie1]=set_frq(1024); w0=xget('window'); xset('window',disp0); clf(); subplot(311); gainplot(frq1,db_iir1,phi_iir1); xtitle('frequency response of iir hpf1'); xset('window',w0); endfunction // // for debug hpf //w0test=1000. * 2. * %pi; //q0test=0.7; //xordertest=2; //[yd_hpf, xd_hpf]=set_hpf(w0test, q0test, xordertest); //pointstest=1024; //[iir_hpf1,db_iir1,phi_iir1] = iir_hpf_resp(1,pointstest); //plot_iir_hpf1(db_iir1,phi_iir1,21); //w0test=1000. * 2. * %pi; //q0test=0.7; //xordertest=1; //[yd_hpf, xd_hpf]=set_hpf(w0test, q0test, xordertest); //pointstest=1024; //[iir_hpf1,db_iir1,phi_iir1] = iir_hpf_resp(1,pointstest); //plot_iir_hpf1(db_iir1,phi_iir1,22); // //-------------------------------------------------------------------------- // // iir filter (peaking equalizer ) design for "da" , for daku-on // function [yd_peq, xd_peq]=set_peq(w0, a0db, q0) // w0: frequency [rad] a0: gain[dB] q0: Q // // yd_peq=zeros(2,3); xd_peq=zeros(2,3); yd_peq(1,1)=1.0; xd_peq(2,1)=1.0; if w0 <= 0. then // else //+++ 1st cal // souitizi_henkan t0=1/fs1; wout = ((w0/(2.0 * %pi)) / (1.0 / t0)) * 2.0 * %pi; //* w0-> digital w he */ wout = tan( wout /2.0); omega = wout * 2.0; a0=10.^(a0db/20.); A=sqrt(a0); alfa=sin(omega)/(2. * q0); vl=1; //* cal X*/ xd_peq(vl,1)= 1. + alfa / A ; xd_peq(vl,2)= -2. * cos( omega); xd_peq(vl,3)= 1. - alfa / A; //* cal Y */ yd_peq(vl,1)= 1. + alfa * A; yd_peq(vl,2)= -2. * cos( omega); yd_peq(vl,3)= 1. - alfa * A; // force to set xd_bpf(1) 1.0 yd_peq(vl,1)= yd_peq(vl,1)/ xd_peq(vl,1); yd_peq(vl,2)= yd_peq(vl,2)/ xd_peq(vl,1); yd_peq(vl,3)= yd_peq(vl,3)/ xd_peq(vl,1); xd_peq(vl,2)= -1.0 * xd_peq(vl,2)/xd_peq(vl,1); xd_peq(vl,3)= -1.0 * xd_peq(vl,3)/xd_peq(vl,1); xd_peq(vl,1)= xd_peq(vl,1)/xd_peq(vl,1); end // if (w0 <= 0.) 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 y = iir_peq( xw, vl ) t0=1/fs1; yi = yd_peq(vl,1) + yd_peq(vl,2) * %e^( -1.* t0 * xw * %i) + yd_peq(vl,3) * %e^( -2.* t0 * xw * %i); yb = 1.0 - (xd_peq(vl,2) * %e^( -1.* t0 * xw * %i) + xd_peq(vl,3) * %e^( -2.* t0 * xw * %i)); y = yi / yb; endfunction function [iir_peq1,db_iir1,phi_iir1] = iir_peq_resp(vl,points) [frq1,sfrq1,is1,ie1]=set_frq(points); iir_peq1=ones(1,sfrq1); if WT_QTY == 41 then for v=1:sfrq1 xw= 2 * %pi * frq1(v); iir_peq1(v)=iir_peq( xw,vl ) end end [db_iir1,phi_iir1] =dbphi(iir_peq1); endfunction function plot_iir_peq1(db_iir1,phi_iir1,disp0) [frq1,sfrq1,is1,ie1]=set_frq(1024); w0=xget('window'); xset('window',disp0); clf(); subplot(311); gainplot(frq1,db_iir1,phi_iir1); xtitle('frequency response of iir peq1'); // xset('window',w0); endfunction //------------------ // test prg //w0=2500. * 2. * %pi; //a0db=20.; //q0=2.; //[yd_peq, xd_peq]=set_peq(w0, a0db, q0); //[iir_peq1,db_iir1,phi_iir1] = iir_peq_resp(1,1024); //plot_iir_peq1(db_iir1,phi_iir1,5); // // //--------------------------------------------------------------- //----------------------------------------------------------- // Function for disply two tubes model with "rl" noise frequency-phase response when glottis is closed // // This is maybe frequency response from "rl noise point" function y = two_tubes_rl_noise( xw ) L2= ( 1. + l1 ) * ttl_Length / 2.; L1= ttl_Length - L2; tu1=L1/c0; tu2=L2/c0; //yi = 0.5 * ( 1. + rg ) * ( 1. + r1) * ( 1. + rl ) * %e^( -1. * ( tu1 + tu2) * xw * %i); // this is for UG yi = r1 * r1 * rg * %e^( -2. * tu1 * xw * %i) + r1 * %e^( -2. * tu2 * xw * %i) + (1. - r1 * r1) * rg * %e^( -2 * (tu1 + tu2) * xw * %i); yi = (1. + rl ) * yi; yb = 1 + r1 * rg * %e^( -2 * tu1 * xw * %i) + r1 * rl * %e^( -2 * tu2 * xw * %i) + rl * rg * %e^( -2 * (tu1 + tu2) * xw * %i); // this is for UL y = yi / yb; endfunction // // Below is new revised version, changed to noise into before (1+rl) from -rl // function y = two_tubes_1rl_noise_UN( xw ) L2= ( 1. + l1 ) * ttl_Length / 2.; L1= ttl_Length - L2; tu1=L1/c0; tu2=L2/c0; //yi = 0.5 * ( 1. + rg ) * ( 1. + r1) * ( 1. + rl ) * %e^( -1. * ( tu1 + tu2) * xw * %i); // this is for UG yi = rg * r1 * (1. + rl) * %e^( -2. * tu1 * xw * %i) + (1. + rl); // this is for UN yb = 1 + r1 * rg * %e^( -2 * tu1 * xw * %i) + r1 * rl * %e^( -2 * tu2 * xw * %i) + rl * rg * %e^( -2 * (tu1 + tu2) * xw * %i); // this is for UL y = yi / yb; endfunction function y = two_tubes_1rl_noise_UG( xw ) L2= ( 1. + l1 ) * ttl_Length / 2.; L1= ttl_Length - L2; tu1=L1/c0; tu2=L2/c0; yi = 0.5 * ( 1. + rg ) * ( 1. + r1) * ( 1. + rl ) * %e^( -1. * ( tu1 + tu2) * xw * %i); // this is for UG // yi = rg * r1 * (1. + rl) * %e^( -2. * tu1 * xw * %i) + (1. + rl); // this is for UN yb = 1 + r1 * rg * %e^( -2 * tu1 * xw * %i) + r1 * rl * %e^( -2 * tu2 * xw * %i) + rl * rg * %e^( -2 * (tu1 + tu2) * xw * %i); // this is for UL y = yi / yb; endfunction //-------------------------------------------------------------- // // // // //---------------------------------------------------------------------------- function cal_overall_response() global yg_res1 global tube2_res1 L1 L2 L3 A1 A2 A3 global hpf_res1 global frq_center q0 global yd_bpf xd_bpf global yd_lpf xd_lpf global yd_hpf xd_hpf global ind_noise_res1 global stable_noise_res1 global burst_noise_res1 global yd_peq xd_peq global peq_res1 global overall_res1 db_2tube phi_2tube [yg_res1]=yg_resp(1024); [ tube2_res1, L1, L2, L3, A1, A2, A3 ]= two_tubes2_resp(); [hpf_res1]=hpf_response(1024); [frq_center,q0]= trans_waku_area( noise_waku_area); [yd_bpf, xd_bpf]=set_bpf( (frq_center * 2. * %pi), q0,(i2nd_thd_factor * frq_center * 2. * %pi), q0); [ind_noise_res1]=noise_response( 1024); [yd_lpf, xd_lpf]=set_lpf(stable_noise_lpf * 2. * %pi , q0_stable_noise,xorder_stable_noise ); [stable_noise_resl,db_iir1,phi_iir1] = iir_lpf_resp(1, 1024); [yd_hpf, xd_hpf]=set_hpf(stable_noise_hpf * 2. * %pi , q0_stable_noiseh,xorder_stable_noiseh ); [stable_noise_resh,db_iir1,phi_iir1] = iir_hpf_resp(1, 1024); stable_noise_res1=stable_noise_resl .* stable_noise_resh; [yd_bpf, xd_bpf]=set_bpf(stable_noise_lpf * 2. * %pi , q0_burst_noise, stable_noise_lpf * 2. * %pi , q0_burst_noise); [burst_noise_res1,db_iir1,phi_iir1] = iir_bpf_resp(1, 1024); [emp_frq_center,emp_q0]= trans_space_area( emp_space_area); [yd_peq, xd_peq]=set_peq( (emp_frq_center * 2. * %pi), empA0, emp_q0); [peq_res1,db_iir1,phi_iir1] = iir_peq_resp(1,1024); [overall_res1,db_2tube, phi_2tube ]=overall_response(); endfunction //----------------------------------------------------------------- limit_switch3=0.; // when limit_switch3=1., force add-offset to the calculate point value where original fft response valueis maximum will be both same. limit_switch3s=zeros(max_frames,1); //================================================================= // // +MAIN (3)program starts // about matching to 2 tubes model if T_DEMO==1 then trise=6.0; // just fit for 1st demo version tfall=0.7; // just fit for 1st demo version if SEL_CODE == 1 then // ka_sample WT_QTYs(1)=52; // 2 tube in brust noise r1s(1)=0.5; l1s(1)=0.; ttl_Lengths(1)=18.5; rls(1)=0.7; emp_space_areas(1)=5.2; empA0s(1)=20.; stable_noise_lpfs(1)= 1500./2. ; // This is center of bpf. xxx_lpf but use for bpf WT_QTYs(2)=51; // 2 tube in stable noise, unbrust r1s(2)=0.5; l1s(2)=0.; ttl_Lengths(2)=18.5; rls(2)=0.7; emp_space_areas(2)=5.2; empA0s(2)=20.; stable_noise_lpfs(2)=2000.; //--- WT_QTYs(3)=2; r1s(3)=0.5; l1s(3)=0.; ttl_Lengths(3)=18.5; rls(3)=0.7; emp_space_areas(3)=5.2; empA0s(3)=20.; stable_noise_lpfs(3)=2000.; //--- WT_QTYs(4)=2; r1s(4)=0.6; l1s(4)=0.; ttl_Lengths(4)=18.5; rls(4)=0.9; WT_QTYs(5)=2; r1s(5)=0.8; l1s(5)=0.; ttl_Lengths(5)=18.5; rls(5)=0.9; WT_QTY=WT_QTYs(tframe); r1=r1s(tframe); l1=l1s(tframe); ttl_Length=ttl_Lengths(tframe); rl=rls(tframe); emp_space_area=emp_space_areas(1); empA0=empA0s(1); stable_noise_lpf=stable_noise_lpfs(1); elseif SEL_CODE == 2 then // ?_sample WT_QTY=3; r1=-0.5; l1=0.5; r2=0.8; l2=0.; ttl_Length=23.; rl=0.9; elseif SEL_CODE == 3 then // ?_sample WT_QTY=21; // 2 tube plus independent noise r1=-0.75; l1=-0.2; ttl_Length=18.5; rl=0.9; // noise_waku_area=noise_waku_area0; i2nd_thd_factor=1.5; nA0=mix_amplitude_indep; // else // include SEL_CODE == 999 [tframe0]=edit_tframe(); [WT_QTY]= set_tube_model(); [fc,trise,tfall,tclosed]= set_2tube_para(); end //cal_overall_response() //plot_2tube(3,0); end // // -MAIN (3)program starts // //==============================================================++= // // //----------------------------------------------------------------- // global variables No. 3 // yg= zeros(1,1024); // this will be overwritten. y2tm= zeros(1,1024); // this will be overwritten. y2tm_lpf= zeros(1,1024); // this will be overwritten. y2tm_noise= zeros(1,1024); // this will be overwritten. y3out= zeros(1,1024); // this will be overwritten. y3out_least=zeros(1,1024); // this will be overwritten. // other amplitude= ones(1,2); // this will maybe be overwritten. Wy3out_filename='dummy.wav'; //--------------------------------------------------------- function [N_REPEAT]= set_section_qty() txt1=['sections quantity']; wstr1=sprintf('%f',N_REPEAT); sig1=x_mdialog('How many sections to generate',txt1, [wstr1]); N_REPEAT=int(evstr(sig1(1))); endfunction //--------------------------------------------------------- //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ //++++ // ATTENSION ! A dummy frame will be added. // fisrt( dummy frame) is same as second frame //--------------------------------------------------------- function [l0x]=make_l0x( l0s, N_REPEAT0 ) // simple linear hokan l0x=zeros(1,N_REPEAT0); count0=2; for v1=1:(QTY_FRAME-1) for v2=1:frame_lengths(v1) alfa=(l0s(v1+1) - l0s(v1)) / frame_lengths(v1); l0x(count0)=l0s(v1) + alfa * (v2 - 1); count0=count0+1; end l0x(count0)=l0s(QTY_FRAME) end l0x(1)=l0x(2); // fisrt( dummy frame) is same as second frame endfunction //--------------------------------------------------------- function [ L1z, L2z, L3z, A1z, A2z, A3z,trisez,tfallz,tclosedz,amplitudez,WT_QTY,N_REPEAT,WT_QTYz]= make_L_et_A() if WT_QTYs(1) == 3 then WT_QTY=3; else WT_QTY=2; end N_REPEAT=1+1; //A dummy frame will be added. for v=1:(QTY_FRAME-1) N_REPEAT=N_REPEAT+frame_lengths(v); end WT_QTYz=zeros(2,N_REPEAT); // WT_QTYz(1,1)=WT_QTY; // first dummy frame count0=2; for v1=1:(QTY_FRAME-1) for v2=1:frame_lengths(v1) WT_QTYz(1,count0)=WT_QTYs(v1); WT_QTYz(2,count0)=v1; count0=count0+1; end end WT_QTYz(1,count0)=WT_QTYs(QTY_FRAME); WT_QTYz(2,count0)=QTY_FRAME; // linear hokan if WT_QTY == 3 then [l1x]=make_l0x( l1s ,N_REPEAT); [r1x]=make_l0x( r1s ,N_REPEAT); [l2x]=make_l0x( l2s ,N_REPEAT); [r2x]=make_l0x( r2s ,N_REPEAT); else [l1x]=make_l0x( l1s ,N_REPEAT); [r1x]=make_l0x( r1s ,N_REPEAT); end [ttl_Lengthx]=make_l0x( ttl_Lengths ,N_REPEAT); [amplitudez]=make_l0x( amplitudes ,N_REPEAT); L1z=zeros(1,N_REPEAT); L2z=zeros(1,N_REPEAT); L3z=zeros(1,N_REPEAT); A1z=zeros(1,N_REPEAT); A2z=zeros(1,N_REPEAT); A3z=zeros(1,N_REPEAT); trisez=zeros(1,N_REPEAT); tfallz=zeros(1,N_REPEAT); tclosedz=zeros(1,N_REPEAT); for v=1:N_REPEAT if WT_QTY==3 then bunbo= l1x(v) * l2x(v) + ( l1x(v) - l2x(v) ) + 3.0; L1z(v)= ttl_Lengthx(v) * ( l2x(v) - 1.0 ) * (l1x(v) - 1.0 ) / bunbo ; L2z(v)= ttl_Lengthx(v) * ( l2x(v) - 1.0 ) * ( -1.0 * l1x(v) - 1.0 ) / bunbo ; L3z(v)= ttl_Lengthx(v) * ( l2x(v) + 1.0 ) * (l1x(v) + 1.0 ) / bunbo ; bunbo= r1x(v) * r2x(v) + ( r1x(v) - r2x(v) ) + 3.0; A1z(v)= ttl_Area * ( r2x(v) - 1.0 ) * (r1x(v) - 1.0 ) / bunbo ; A2z(v)= ttl_Area * ( r2x(v) - 1.0 ) * ( -1.0 * r1x(v) - 1.0 ) / bunbo ; A3z(v)= ttl_Area * ( r2x(v) + 1.0 ) * (r1x(v) + 1.0 ) / bunbo ; else L2z(v)= ( 1. + l1x(v) ) * ttl_Lengthx(v) / 2.; L1z(v)= ttl_Lengthx(v) - L2z(v); A2z(v)= ( 1. + r1x(v) ) * ttl_Area/ 2.; A1z(v)= ttl_Area - A2z(v); // +dummy data set L3z(v)=1; A3z(v)=1; // -dummy data set end trisez(v)=trise; tfallz(v)=tfall; tclosedz(v)=tclosed; end endfunction // function plot_L_A(L1z, L2z, L3z, A1z, A2z, A3z,WT_QTY,N_REPEAT,disp0) wb0=xget('window'); // stack old window xset('window',disp0); // create new windows clf(); dx0=zeros(1,N_REPEAT); for v=1:N_REPEAT dx0(v)=v; end if WT_QTY==3 then subplot(211); plot2d( dx0, L1z ,style=[color("red")]); plot2d( dx0, L2z ,style=[color("blue")]); plot2d( dx0, L3z ,style=[color("yellow")]); wstr0='tube length: L1(red) L2(blue) L3(yellow)'; xtitle(wstr0); subplot(212); plot2d( dx0, A1z ,style=[color("red")]); plot2d( dx0, A2z ,style=[color("blue")]); plot2d( dx0, A3z ,style=[color("yellow")]); wstr0='tube area: A1(red) A2(blue) A3(yellow) '; xtitle(wstr0); else subplot(211); plot2d( dx0, L1z ,style=[color("red")]); plot2d( dx0, L2z ,style=[color("blue")]); wstr0='tube length: L1(red) L2(blue) '; xtitle(wstr0); subplot(212); plot2d( dx0, A1z ,style=[color("red")]); plot2d( dx0, A2z ,style=[color("blue")]); wstr0='tube area: A1(red) A2(blue) '; xtitle(wstr0); end xset('window',wb0); // push old windows endfunction //---------------------------------------------------------- function plot_area_mkwave_least(disp0,code0) // disp0=4? global mess_init0 global r1s r2s l1s l2s ttl_Lengths rls noise_waku_areas i2nd_thd_factors nA0s global emp_space_areas empA0s stable_noise_lpfs global yg y2tm y3out y2tm_lpf y2tm_noise LL y3out_least mess_init0b=mess_init0; r1sb=r1s; r2sb=r2s; l1sb=l1s; l2sb=l2s; ttl_Lengthsb=ttl_Lengths; rlsb=rls; noise_waku_areasb=noise_waku_areas; i2nd_thd_factorsb=i2nd_thd_factors; nA0sb=nA0s; emp_space_areasb=emp_space_areas; empA0sb=empA0s; stable_noise_lpfsb=stable_noise_lpfs; ygb=yg; y2tmb=y2tm; y3outb=y3out; y2tm_lpfb=y2tm_lpf; y2tm_noiseb=y2tm_noise; LLb=LL; exchange8s(); mess_init0=' by leastsq method result'; if code0 == 1 then // only plot area [ L1z, L2z, L3z, A1z, A2z, A3z,trisez,tfallz,tclosedz,amplitudez,WT_QTY,N_REPEAT,WT_QTYz]= make_L_et_A(); plot_areas(L1z, L2z, L3z, A1z, A2z, A3z,WT_QTY,N_REPEAT,disp0); elseif code0 == 2 then // make waveform and plot [yg,y2tm,y3out,y2tm_lpf,y2tm_noise,LL]= make_waveform(disp0); y3out_least=y3out; plot_waveform(disp0+1); end mess_init0=mess_init0b; r1s=r1sb; r2s=r2sb; l1s=l1sb; l2s=l2sb; ttl_Lengths=ttl_Lengthsb; rls=rlsb; noise_waku_areas=noise_waku_areasb; i2nd_thd_factors=i2nd_thd_factorsb; nA0s=nA0sb; emp_space_areas=emp_space_areasb; empA0s=empA0sb; stable_noise_lpfs=stable_noise_lpfsb; yg=ygb; y2tm=y2tmb; y3out=y3outb; y2tm_lpf=y2tm_lpfb; y2tm_noise=y2tm_noiseb; LL=LLb; endfunction //++++ // ATTENSION ! A dummy frame will be added. // fisrt( dummy frame) is same as second frame // function [yg,y2tm,y3out,y2tm_lpf,y_ran3,LL]= make_waveform(disp4) global yd_lpf xd_lpf global yd_bpf xd_bpf global yd_hpf xd_hpf // y_ran3 is noise output [ L1z, L2z, L3z, A1z, A2z, A3z,trisez,tfallz,tclosedz,amplitudez,WT_QTY,N_REPEAT,WT_QTYz]= make_L_et_A(); plot_areas(L1z, L2z, L3z, A1z, A2z, A3z,WT_QTY,N_REPEAT,disp4); // disp0=4 // y_ran3=zeros(1,1); // this is only for dummy set. y2tm_lpf= zeros(1,1); // this is only for dummy set. // if WT_QTY == 3 then //Three Tubes Model of serial type tu1=L1z./c0; // delay time in 1st tube tu2=L2z./c0; // delay time in 2nd tube tu3=L3z./c0; // delay time in 3rd tube r1z=(A2z-A1z)./(A2z+A1z); // reflection coefficient between 1st tube and 2nd tube r2z=(A3z-A2z)./(A3z+A2z); // reflection coefficient between 2nd tube and 3rd tube elseif WT_QTY == 4 then //Three Tubes Model of T type tu1=L1z./c0; // delay time in 1st tube tu2=L2z./c0; // delay time in 2nd tube tu3=L3z./c0; // delay time in 3rd tube r12z=((A2z+A3z)-A1z)./(A3z+A2z+A1z); // reflection coefficient between 1st tube and others r21z=(A2z-(A1z+A3z))./(A3z+A2z+A1z); // reflection coefficient between 2nd tube and others r31z=(A3z-(A1z+A2z))./(A3z+A2z+A1z); // reflection coefficient between 3rd tube and others else //Two Tubes Model tu1=L1z./c0; // delay time in 1st tube tu2=L2z./c0; // delay time in 2nd tube r1z=(A2z-A1z)./(A2z+A1z); // reflection coefficient between 1st tube and 2nd tube end // disp('+++enter step 1 Input Glottal Volume Velocity Generation'); D_QTY=N_REPEAT; 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( tclosedz(loop) / 1000 / ts ); N2(loop)= int( trisez(loop) / 1000 / ts); N3(loop)= int( tfallz(loop) / 1000 / ts); LL(loop)= N1(loop)+N2(loop)+N3(loop); length0 = length0 + LL(loop); end // yg is Glottal Volume Velocity. rg is reflection coefficient between glottis and 1st tube // 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=rg; // 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 yg= zeros(1,length0+1); rgz= zeros(1,length0+1); tc0=1; yg(tc0)=0.; yg(tc0)=amplitudez(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; rgz(tc0)=rg_closed; elseif t0 <= (N2(loop) + N1(loop)) then yg(tc0)= 0.5 * ( 1 - cos( ( %pi / N2(loop) ) * (t0 - N1(loop)) ) ); rgz(tc0)=rg_rise; elseif t0 <= (N3(loop)+N2(loop)+N1(loop)) then yg(tc0)= cos( ( %pi / ( 2 * N3(loop) )) * ( t0 - (N2(loop)+N1(loop)) ) ); rgz(tc0)=rg_fall; else yg(tc0) =0; rgz(tc0)=rg_closed; end yg(tc0)=amplitudez(loop) * yg(tc0); end // t0 // if WT_QTYz(1,loop)==51 then [yd_lpf, xd_lpf]=set_lpf(stable_noise_lpfs(WT_QTYz(2,loop)) * 2. * %pi , q0_stable_noise,xorder_stable_noise ); //[xy_ran2s]=make_noise_lpf( int(LL(loop)) ); [yd_hpf, xd_hpf]=set_hpf(stable_noise_hpf * 2. * %pi , q0_stable_noiseh,xorder_stable_noiseh ); [xy_ran2s]=make_noise_lpf_hpf( int(LL(loop)) ); for vloop=1:int(LL(loop)) yg(tc0-int(LL(loop))+vloop)= amplitudez(loop) * xy_ran2s(vloop); end elseif WT_QTYz(1,loop)==52 then //[yd_lpf, xd_lpf]=set_lpf(stable_noise_lpfs(WT_QTYz(2,loop)) * 2. * %pi , q0_stable_noise,xorder_stable_noise ); //[xy_ran2s]=make_noise_lpf( int(LL(loop)) ); [yd_bpf, xd_bpf]=set_bpf( (stable_noise_lpfs(WT_QTYz(2,loop)) * 2. * %pi) , q0_burst_noise , (stable_noise_lpfs(WT_QTYz(2,loop)) * 2. * %pi) , q0_burst_noise ); [xy_ran2s]=make_noise_bpf( int(LL(loop)) ); for vloop=1:int(LL(loop)) yg(tc0-int(LL(loop))+vloop)= amplitudez(loop) * xy_ran2s(vloop); end end // end // loop // // if WT_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. + rgz(tc0) ) / 2.) * yg(tc0) + rgz(tc0) * ya2(M1(loop)); ya2(1)= -1. * r1z(loop) * ya1(M1(loop)) + ( 1. - r1z(loop) ) * yb2(M2(loop)); yb1(1)= ( 1 + r1z(loop) ) * ya1(M1(loop)) + r1z(loop) * yb2(M2(loop)); yb2(1)= -1. * r2z(loop) * yb1(M2(loop)) + ( 1. - r2z(loop) ) * yc2(M3(loop)); yc1(1)= ( 1 + r2z(loop) ) * yb1(M2(loop)) + r2z(loop) * yc2(M3(loop)); yc2(1)= -1. * rl * yc1(M3(loop)); y2tm(tc0)= (1 + rl) * yc1(M3(loop)); end // t0 end //loop else //Two Tubes Model disp('+++enter step 2 Two Tubes Model for Vocal Tract Calculation'); //+for noise mix if WT_QTY==21 then i_type=4; [y_ran2]=make_noise_bpf(length0+1); if i_type == 1 then disp('rl noise.'); elseif i_type == 2 then disp('1+rl noise.'); elseif ((i_type == 3) | (i_type == 4)) then disp('independent noise.'); end // 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); y_ran3=zeros(1, (length0 + 1)); end //-for noise mix //+peq y_peq0=zeros(1, length0); yx_peq=zeros(2,6); //-peq // 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); 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.; y2tm_lpf(tc0)=0.; for loop=1:D_QTY //+peq if WT_QTYz(1,loop)==41 then [emp_frq_center,emp_q0]= trans_space_area( emp_space_areas(WT_QTYz(2,loop)) ); [yd_peq, xd_peq]=set_peq( (emp_frq_center * 2. * %pi), empA0, emp_q0); else yd_peq=zeros(2,3); xd_peq=zeros(2,3); yd_peq(1,1)=1.0; xd_peq(2,1)=1.0; end //-peq 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. + rgz(tc0) ) / 2.) * yg(tc0) + rgz(tc0) * ya2(M1(loop)); ya2(1)= -1. * r1z(loop) * ya1(M1(loop)) + ( 1. - r1z(loop) ) * yb2(M2(loop)); yb1(1)= ( 1 + r1z(loop) ) * ya1(M1(loop)) + r1z(loop) * yb2(M2(loop)); yb2(1)= -1. * rl * yb1(M2(loop)); y2tm(tc0)= (1 + rl) * yb1(M2(loop)); //+peq for vl=1:1 yx_peq(vl,4)=y2tm(tc0) x0peq=yd_peq(vl,1) * yx_peq(vl,4) + yd_peq(vl,2) * yx_peq(vl,5) + yd_peq(vl,3) * yx_peq(vl,6); // x0= yi[eq0][0] * yx[eq0][3] + yi[eq0][1] * yx[eq0][4] + yi[eq0][2] * yx[eq0][5]; yx_peq(vl,1)=x0peq + yx_peq(vl,2) * xd_peq(vl,2) + yx_peq(vl,3) * xd_peq(vl,3); // yx[eq0][0] = x0 + yx[eq0][1] * xi[eq0][1] + yx[eq0][2] * xi[eq0][2]; // /* shift for next time */ yx_peq(vl,3)=yx_peq(vl,2); // yx[eq0][2]=yx[eq0][1]; yx_peq(vl,2)=yx_peq(vl,1); // yx[eq0][1]=yx[eq0][0]; yx_peq(vl,6)=yx_peq(vl,5); // yx[eq0][5]=yx[eq0][4]; yx_peq(vl,5)=yx_peq(vl,4); // yx[eq0][4]=yx[eq0][3]; end yx_peq(2,1)=0.; // uncoment when use yx_peq(2,*). y_peq0(tc0)=yx_peq(1,1)+ yx_peq(2,1); // delta5=0.05; // cross fader steps is 1.0 / delta5 if WT_QTYz(1,loop) == 41 then alfa5=1.0; if loop > 1 then if WT_QTYz(1,loop-1) <> 41 then //start... alfa5= delta5 * (t0-1); // for t0=1:int(LL(loop)) if alfa5 > 1.0 then alfa5=1.0; else // disp('(test1) cross fading'); end end end // If it's only One frame, below does not work... if loop < D_QTY then if WT_QTYz(1,loop+1) <> 41 then //...end alfa5= delta5 * (int(LL(loop)) - t0); // for t0=1:int(LL(loop)) if alfa5 > 1.0 then alfa5=1.0; else // disp('(test2) cross fading'); end end end // // cross fade mix control by alfa5. y2tm is original. y2tm(tc0)=(1.0 - alfa5) * y2tm(tc0) + alfa5 * y_peq0(tc0); end //-peq //+for noise mix 2 if WT_QTY == 21 then // 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 ((i_type == 3) | (i_type == 4)) then if y2tm_lpf(tc0) > threshold_occur_plus_id then ynoise = y_ran2(tc0); elseif y2tm_lpf(tc0) < threshold_occur_minus_id 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 if i_type == 1 then yb2(1)= yb2(1) + mix_amplitude_rl * ynoise; // mix with noise y_ran3(tc0)=mix_amplitude_rl * ynoise; elseif i_type == 2 then yb1(1)= yb1(1) + mix_amplitude_1rl * ynoise; // mix with noise y_ran3(tc0)=mix_amplitude_1rl * ynoise; elseif i_type == 3 then y2tm(tc0)= y2tm(tc0) + mix_amplitude_indep * ynoise; y_ran3(tc0)=mix_amplitude_indep * ynoise; elseif i_type == 4 then y2tm(tc0)= y2tm(tc0) + nA0 * ynoise; y_ran3(tc0)=nA0 * ynoise; end end //-for noise mix 2 end // t0 end //loop end //End of Two Tubes model // disp('+++enter step 3 Output Sound Pressure for Two Tubes Model for Vocal Tract Calculation'); // At first, low pass filter's coefficients wcz=2. * %pi * fc; al1= -1. * %e^( -1. * wcz * ts); bl0= wcz * ts; wc_org=wcz; wc_new=2. * %pi * fc; alfa=-1. * cos( (wc_org + wc_new) * ts /2. ) / cos( (wc_org - wc_new) * 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); // IIR type high pass filter y3out= zeros(1,length0+1); LI1=2; LI2=2; ya1= zeros(1,LI1); yb1= 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); end //loop disp('--- Waveform Generation finished.'); endfunction // end ofmake_waveform() //---------------------------------------------------------- function y3out_snd_play1() // 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! sy3out=size(y3out); y3outz=zeros(sy3out(2)); vm0=abs(y3out(1)); for v=1:sy3out(2) if abs(y3out(v)) > vm0 then vm0=abs(y3out(v)); end end if vm0 > 0 then for v=1:sy3out(2) y3outz(v)= y3out(v) / vm0; end end if f_win == 1 then if f412 == 1 then // for windows scilab-4.1.2 sound(y3outz' ,fs1,bit1); else // for windows scilab-3.1.1 if fs1 == 22050 then sound(y3outz,fs1,bit1); disp('This function may work, if luckily.'); elseif fs1 == 44100 then // down-sampling from 44100 to 22050 wdatw2=zeros((sy3out(2)/2)); for v=1:(sy3out(2)/2) wdatw2(v)=y3outz(2 * (v -1) +1); 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 else disp('If sound setting is good for scilab, this function maybe work.'); if f412 == 1 then // for scilab-4.1.2 sound(y3outz' ,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(y3outz,fs1,bit1); disp('This function may work, if luckily.'); elseif fs1 == 44100 then // down-sampling from 44100 to 22050 wdatw2=zeros((sy3out(2)/2)); for v=1:(sy3out(2)/2) wdatw2(v)=y3outz(2 * (v -1) +1); 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 y3out_snd_save1() // wavefilename= input(' + enter file name for saved .wav file =>',["string"]); // sy3out=size(y3out); y3outz=zeros(sy3out(2)); vm0=abs(y3out(1)); for v=1:sy3out(2) if abs(y3out(v)) > vm0 then vm0=abs(y3out(v)); end end if vm0 > 0 then for v=1:sy3out(2) y3outz(v)= y3out(v) / vm0; end end // if f412 == 1 then // for windows scilab-4.1.2 wavwrite(y3outz' ,fs1,bit1,wavefilename ); else // for windows scilab-3.1.1 wavwrite(y3outz, fs1,bit1, wavefilename); end endfunction //--------------------------------------------------------- function plot_waveform(disp0) wb0=xget('window'); // stack old window xset('window',disp0); // create new windows clf(); if WT_QTY == 21 then if (f412 == 1) | (f_win==1) then plot(yg,'b'); plot(y2tm,'y'); plot(y2tm_lpf,'g'); plot(y3out,'r'); plot(y2tm_noise,'k'); xtitle('generated wave' ,'Samples (Time)', 'Amplitude'); legend('Glottal Volume Velocity','2nd Tube Edge Volume Velocity','Smoothed Volume Velocity','Sound Pressure','Additional noise',2); // // plot additional noise fft response //[respn,db_fft1n,phi_fft1n]=y2tmx_resp(y2tm_noise,2,1024); //plot_fft_y2tmx(db_fft1n,phi_fft1n,1024,13); // plot windows #13 else ygs=size(yg); xzikur=[1:1:ygs(2)]; //plot2d(xzikur,yg,style=[color("blue")]); plot2d(xzikur,y2tm,style=[color("yellow")]); plot2d(xzikur,y2tm_lpf,style=[color("green")]); plot2d(xzikur,y3out,style=[color("red")]); plot2d(xzikur,y2tm_noise,style=[color("black")]); xtitle('generated wave by tubes model','Samples (Time)', 'Amplitude'); //legend('Glottal Volume Velocity','Sound Pressure',2); end else if f412 == 1 then plot(yg,'b'); plot(y3out,'r'); wstr1= 'generated wave by tubes model' + mess_init0; xtitle(wstr1,'Samples (Time)', 'Amplitude'); legend('Glottal Volume Velocity','Sound Pressure',2); elseif f_win == 1 then plot(yg,'b'); plot(y3out,'r'); wstr1= 'generated wave by tubes model' + mess_init0; xtitle(wstr1,'Samples (Time)', 'Amplitude'); //xtitle('generated wave by tubes model','Samples (Time)', 'Amplitude'); legend('Glottal Volume Velocity','Sound Pressure',2); else ygs=size(yg); xzikur=[1:1:ygs(2)]; //plot2d(xzikur,yg,style=[color("blue")]); plot2d(xzikur,y3out,style=[color("red")]); wstr1= 'generated wave by tubes model' + mess_init0; xtitle(wstr1,'Samples (Time)', 'Amplitude'); //xtitle('generated wave by tubes model','Samples (Time)', 'Amplitude'); //legend('Glottal Volume Velocity','Sound Pressure',2); end end xset('window',wb0); // push old windows endfunction // //================================================================= // global yg_res9=zeros(2,1); // dummy set hpf_res9=zeros(2,1); // dummy set wm8=ones(2,1); // dummy set wm8s=ones(max_frames,2); wm9=ones(2,1); // dummy set tm8=ones(2,1); // dummy set x_init8=zeros(2,1); // dummy set xopt=zeros(2,1); // dummy set xopts=zeros(max_frames,100); // first (*,1) is qty of the datas ym8=zeros(2,1); // dummy set ma0=['ttl_Length ' ; 'l1 ' ; 'r1 ' ; 'l2 ' ; 'r2 ' ; 'rl ']; ma1=['ttl_Length ' ; 'l1 ' ; 'r1 ' ; 'rl ']; ma2=['ttl_Length ' ; 'l1 ' ; 'r1 ' ; 'rl ' ;'wake_area ' ; '2nd thd '; 'amplitude ']; ma3=['ttl_Length ' ; 'l1 ' ; 'r1 ' ; 'rl ' ;'emp_space ' ; 'amplitude ']; ma4=['ttl_Length ' ; 'l1 ' ; 'r1 ' ; 'rl ' ;'noise_lpfFc ' ; ' ']; ma5=['ttl_Length ' ; 'l1 ' ; 'r1 ' ; 'rl ' ;'noise_bpfCf ' ; ' ']; //---------------------------------------------------------------- function [tm8, ym8, x_init8, wm9]=prepare8() global max_nth max_value max_wm9 [frq1,sfrq1,is1,ie1]=set_frq(fft_point1); wm9=ones(sfrq1,1); max_wm9=zeros(1,sfrq1); ym8=zeros(sfrq1,1); tm8=zeros(sfrq1,1); for v=1:sfrq1 tm8(v)= 2 * %pi * frq1(v); ym8(v)=db_fft1s(tframe,v); end [max_nth, max_value]=maxi(ym8); max_wm9(1,max_nth)=1.0; //. if WT_QTY == 3 then // when 3 tube model serial x_init8=zeros(6,1); x_init8(1,1)=ttl_Lengths(tframe); x_init8(2,1)=l1s(tframe); x_init8(3,1)=r1s(tframe); x_init8(4,1)=l2s(tframe); x_init8(5,1)=r2s(tframe); x_init8(6,1)=rls(tframe); elseif WT_QTY == 21 then // when 2 tube model plus independent noise x_init8=zeros(7,1); x_init8(1,1)=ttl_Lengths(tframe); x_init8(2,1)=l1s(tframe); x_init8(3,1)=r1s(tframe); x_init8(4,1)=rls(tframe); x_init8(5,1)=noise_waku_areas(tframe); x_init8(6,1)=i2nd_thd_factors(tframe); x_init8(7,1)=nA0s(tframe); elseif WT_QTY == 41 then // when 2 tube model plus partly emphasis x_init8=zeros(6,1); x_init8(1,1)=ttl_Lengths(tframe); x_init8(2,1)=l1s(tframe); x_init8(3,1)=r1s(tframe); x_init8(4,1)=rls(tframe); x_init8(5,1)=emp_space_areas(tframe); x_init8(6,1)=empA0s(tframe); elseif (WT_QTY == 51) | (WT_QTY == 52) then // when 2 tube model driven by noise x_init8=zeros(5,1); x_init8(1,1)=ttl_Lengths(tframe); x_init8(2,1)=l1s(tframe); x_init8(3,1)=r1s(tframe); x_init8(4,1)=rls(tframe); x_init8(5,1)=stable_noise_lpfs(tframe); else // when 2 tube model x_init8=zeros(4,1); x_init8(1,1)=ttl_Lengths(tframe); x_init8(2,1)=l1s(tframe); x_init8(3,1)=r1s(tframe); x_init8(4,1)=rls(tframe); end disp(' + prepare8 done'); endfunction //---------------------------------------------------------------- function [yg_res9, hpf_res9]=prepare9() global yg_res9 hpf_res9 yg_res9=yg_resp(fft_point1)'; hpf_res9=hpf_response(fft_point1)'; // disp(' + prepare9 done'); endfunction //----------------------------------------------------------------- function exchange8() // set leastsq result as tube paras global r1 r2 l1 l2 ttl_Length rl noise_waku_area i2nd_thd_factor nA0 global emp_space_area empA0 stable_noise_lpf wstr1=sprintf('%f',r2); wstr2=sprintf('%f',l2); wstr3=sprintf('%f',noise_waku_area); wstr4=sprintf('%f',i2nd_thd_factor); wstr5=sprintf('%f',nA0); wstr6=sprintf('%f',emp_space_area); wstr7=sprintf('%f',empA0); wstr8=sprintf('%f',stable_noise_lpf); // 2 tubes model ttl_Length=xopt(1); l1=xopt(2); r1=xopt(3); rl=xopt(4); if WT_QTY == 3 then l2=xopt(4); r2=xopt(5); rl=xopt(6); else r2=evstr(wstr1); l2=evstr(wstr2); end if WT_QTY == 21 then noise_waku_area=xopt(5); i2nd_thd_factor=xopt(6); nA0=xopt(7); else noise_waku_area=evstr(wstr3); i2nd_thd_factor=evstr(wstr4); nA0=evstr(wstr5); end if WT_QTY == 41 then emp_space_area=xopt(5); empA0=xopt(6); else emp_space_area=evstr(wstr6); empA0=evstr(wstr7); end if (WT_QTY == 51) | (WT_QTY == 52) then stable_noise_lpf=xopt(5); else stable_noise_lpf=evstr(wstr8); end endfunction //- function exchange8s() // set leastsq result as tube paras global r1s r2s l1s l2s ttl_Lengths rls noise_waku_areas i2nd_thd_factors nA0s global emp_space_areas empA0s stable_noise_lpfs ofs=1; for vloop=1:QTY_FRAME // 2 tubes model ttl_Lengths(vloop)=xopts(vloop,1+ofs); l1s(vloop)=xopts(vloop,2+ofs); r1s(vloop)=xopts(vloop,3+ofs); rls(vloop)=xopts(vloop,4+ofs); if WT_QTYs(vloop) == 3 then l2s(vloop)=xopts(vloop,4+ofs); r2s(vloop)=xopts(vloop,5+ofs); rls(vloop)=xopts(vloop,6+ofs); else r2s(vloop)=0.; l2s(vloop)=0.; end if WT_QTYs(vloop) == 21 then noise_waku_areas(vloop)=xopts(vloop,5+ofs); i2nd_thd_factors(vloop)=xopts(vloop,6+ofs); nA0s(vloop)=xopts(vloop,7+ofs); else noise_waku_areas(vloop)=0.; i2nd_thd_factors(vloop)=0.; nA0s(vloop)=0.; end if WT_QTYs(vloop) == 41 then emp_space_areas(vloop)=xopts(vloop,5+ofs); empA0s(vloop)=xopts(vloop,6+ofs); else emp_space_areas(vloop)=0.; empA0s(vloop)=0.; end if (WT_QTYs(vloop) == 51) | (WT_QTYs(vloop) == 52) then stable_noise_lpfs(vloop)=xopts(vloop,5+ofs); else stable_noise_lpfs(vloop)=0.; end end // for vloop=1:QTY_FRAME endfunction //- //function set_exchange9() // set each frame value as tube paras // global r1 r2 l1 l2 ttl_Length rl noise_waku_area i2nd_thd_factor nA0 // global WT_QTY // // WT_QTY=WT_QTYs(tframe); // ttl_Length=ttl_Lengths(tframe); // l1=l1s(tframe); // r1=r1s(tframe); // l2=l2s(tframe); // r2=r2s(tframe); // rl=rls(tframe); // noise_waku_area=noise_waku_areas(tframe); // i2nd_thd_factor=i2nd_thd_factors(tframe); // nA0=nA0s(tframe); //endfunction //- function exchange9() // set initial value as tube paras global r1 r2 l1 l2 ttl_Length rl noise_waku_area i2nd_thd_factor nA0 global stable_noise_lpf wstr1=sprintf('%f',r2); wstr2=sprintf('%f',l2); wstr3=sprintf('%f',noise_waku_area); wstr4=sprintf('%f',i2nd_thd_factor); wstr5=sprintf('%f',nA0); wstr6=sprintf('%f',emp_space_area); wstr7=sprintf('%f',empA0); wstr8=sprintf('%f',stable_noise_lpf); // 2 tubes model ttl_Length=x_init8(1); l1=x_init8(2); r1=x_init8(3); rl=x_init8(4); if WT_QTY == 3 then l2=x_init8(4); r2=x_init8(5); rl=x_init8(6); else r2=evstr(wstr1); l2=evstr(wstr2); end if WT_QTY == 21 then noise_waku_area=x_init8(5); i2nd_thd_factor=x_init8(6); nA0=x_init8(7); else noise_waku_area=evstr(wstr3); i2nd_thd_factor=evstr(wstr4); nA0=evstr(wstr5); end if WT_QTY == 41 then emp_space_area=x_init8(5); empA0=x_init8(6); else emp_space_area=evstr(wstr6); empA0=evstr(wstr7); end if (WT_QTY == 51) | (WT_QTY == 52) then stable_noise_lpf=x_init8(5); else stable_noise_lpf=evstr(wstr8); end endfunction //----------------------------------------------------------------- //(1)test1 // (1 + exp( fact0 * x0)) * (1 + exp(fact0 * x1)) // x0=x - fact1, x1= fact2 -x // // ex: (1 + exp( 50 * (x - 0.9)) * (1 + exp( 50 * (-0.9 -x)) // When fact0 = 50, fact1 = 0.9, fact2 = -0.9 // hoping that -1 < x < 1 ( fact2 ~< x ~< fact1 ) // limit_switch0=0.; // set 1 when r1r2,l1,l2 limit is enable, set 0 when limit is disable limit_switch0s=zeros(max_frames,1); fact0=20.; // // for limit of 3 tubes model fact1_3r1=-0.1; fact2_3r1=-0.9; fact1_3r2=0.9; fact2_3r2=0.1; fact1_3l1=0.9; fact2_3l1=-0.9; fact1_3l2=0.9; fact2_3l2=-0.9; fact1_3r1s=zeros(max_frames,1); fact2_3r1s=zeros(max_frames,1); fact1_3r2s=zeros(max_frames,1); fact2_3r2s=zeros(max_frames,1); fact1_3l1s=zeros(max_frames,1); fact2_3l1s=zeros(max_frames,1); fact1_3l2s=zeros(max_frames,1); fact2_3l2s=zeros(max_frames,1); // for limit of 2 tubes model fact1_2r1=0.9; fact2_2r1=-0.9; fact1_2l1=0.9; fact2_2l1=-0.9; fact1_2r1s=zeros(max_frames,1); fact2_2r1s=zeros(max_frames,1); fact1_2l1s=zeros(max_frames,1); fact2_2l1s=zeros(max_frames,1); // // // limit_switch2=1.0; // set 1 when rl limit is enable, set 0 when limit is disable limit_switch2s=zeros(max_frames,1); // for limit of rl fact1_rl=1.0; fact2_rl=0.5; fact1_rls=zeros(max_frames,1); fact2_rls=zeros(max_frames,1); // // // for noise_waku_area=noise_waku_area0; // dummy set limit_switch7=1.0; fact1_waku=5.0; fact2_waku=0.5; // for i2nd_thd_factor=0.; // dummy set limit_switch8=1.0; fact1_i2nd_thd=2.1; fact2_i2nd_thd=1.1; // for nA0=mix_amplitude_indep; // dummy set limit_switch9=1.0; fact2=200.; fact1_nA0= 0.09; fact2_nA0= 0.01; // // // for emp_space_area // dummy set limit_switch10=1.0; fact1_space=5.0; fact2_space=0.5; // for empA0=0.; // dummy set limit_switch11=1.0; fact1_empA0=30.; fact2_empA0=0.1; // for stable_noise_lpf_fc // dummy set limit_switch12=1.0; fact1_sn_lpf_fc=4000.; fact2_sn_lpf_fc=500.; // //===================================================== // // +MAIN (L)program starts // if T_DEMO==1 then if SEL_CODE == 1 then // ka_sample // part 1 brust fact1_2r1s(1)=0.9; fact2_2r1s(1)=-0.9; fact1_2l1s(1)=0.9; fact2_2l1s(1)=-0.9; fact1_rls(1)=0.85; // adjusted fact2_rls(1)=0.5; limit_switch0s(1)=0.0; // r1,r2,l1,l2 off limit_switch2s(1)=1.0; // rl on limit_switch3s(1)=0.0; // off : force add-offset to // part 2 2 tube in stable noise, unbrust fact1_2r1s(2)=0.9; fact2_2r1s(2)=-0.9; fact1_2l1s(2)=0.9; fact2_2l1s(2)=-0.9; fact1_rls(2)=0.85; // adjusted fact2_rls(2)=0.5; limit_switch0s(2)=0.0; // r1,r2,l1,l2 off limit_switch2s(2)=1.0; // rl on limit_switch3s(2)=0.0; // off : force add-offset to // part 3 fact1_2r1s(3)=0.9; fact2_2r1s(3)=-0.9; fact1_2l1s(3)=0.9; fact2_2l1s(3)=-0.9; fact1_rls(3)=0.85; // adjusted fact2_rls(3)=0.5; limit_switch0s(3)=1.0; // r1,r2,l1,l2 on limit_switch2s(3)=1.0; // rl on limit_switch3s(3)=0.0; // off : force add-offset to // part 4 fact1_2r1s(4)=0.9; fact2_2r1s(4)=-0.9; fact1_2l1s(4)=0.9; fact2_2l1s(4)=-0.9; fact1_rls(4)=0.95; fact2_rls(4)=0.5; limit_switch0s(4)=1.0; // r1,r2,l1,l2 on limit_switch2s(4)=1.0; // rl on limit_switch3s(4)=0.0; // off : force add-offset to // part 5 fact1_2r1s(5)=0.9; fact2_2r1s(5)=-0.9; fact1_2l1s(5)=0.9; fact2_2l1s(5)=-0.9; fact1_rls(5)=0.95; fact2_rls(5)=0.5; limit_switch0s(5)=1.0; // r1,r2,l1,l2 on limit_switch2s(5)=1.0; // rl on // //********************************************** // off : force add-offset to limit_switch3s(5)=0.0; // stable part of "a" must be 0 ? // if this switch set as 1, result frequency response lose peak-ness of which "a" feature has, and become flatter. // // fact1_2r1=fact1_2r1s(tframe); fact2_2r1=fact2_2r1s(tframe); fact1_2l1=fact1_2l1s(tframe); fact2_2l1=fact2_2l1s(tframe); fact1_rl=fact1_rls(tframe); fact2_rl=fact2_rls(tframe); limit_switch0=limit_switch0s(tframe); limit_switch2=limit_switch2s(tframe); limit_switch3=limit_switch3s(tframe); elseif SEL_CODE == 2 then // ?_sample fact1_3r1=-0.1; fact2_3r1=-0.9; fact1_3r2=0.9; fact2_3r2=0.1; fact1_3l1=0.9; fact2_3l1=-0.9; fact1_3l2=0.9; fact2_3l2=-0.9; fact1_rl=1.0; fact2_rl=0.5; limit_switch0=1.0; // r1,r2,l1,l2 on limit_switch2=1.0; // rl on elseif SEL_CODE == 3 then // ?_sample fact1_2r1=-0.1; fact2_2r1=-0.9; fact1_2l1=0.9; fact2_2l1=-0.9; // teiiki slop wo dasutame, teiiki ha hikaku no sample-suu ga sukunakunai-tame, kawarini weigh wo masu, Plus limit rl <=0.9 mo hituyou // this slop fitting is not well yet. fact1_rl=0.9; fact2_rl=0.5; limit_switch0=1.0; // r1,r2,l1,l2 on limit_switch2=1.0; // rl on else // include SEL_CODE == 999 // as same as initial values end end // // -MAIN (L)program starts // // //===================================================== // //--------------------------------------------- function [act0]= edit_limit0() global fact1_2r1s fact2_2r1s fact1_2l1s fact2_2l1s fact1_rls fact2_rls limit_switch0s limit_switch2s limit_switch3s global fact1_3l1s fact2_3l1s fact1_3r1s fact2_3r1s fact1_3l2s fact2_3l2s fact1_3r2s fact2_3r2s txt1=['switch on(1)/off(0)';'exp(a0*x),a0=';'l1 plus';'l1 minus';'r1 plus';'r1 minus';'l2 plus';'l2 minus';'r2 plus';'r2 minus']; txt2=['switch on(1)/off(0)';'exp(a0*x),a0=';'l1 plus';'l1 minus';'r1 plus';'r1 minus']; wstr0=sprintf('%f',limit_switch0s(tframe)); wstr1=sprintf('%f',fact0); wstr2=sprintf('%f',fact1_3l1s(tframe)); wstr3=sprintf('%f',fact2_3l1s(tframe)); wstr4=sprintf('%f',fact1_3r1s(tframe)); wstr5=sprintf('%f',fact2_3r1s(tframe)); wstr6=sprintf('%f',fact1_3l2s(tframe)); wstr7=sprintf('%f',fact2_3l2s(tframe)); wstr8=sprintf('%f',fact1_3r2s(tframe)); wstr9=sprintf('%f',fact2_3r2s(tframe)); wstr10=sprintf('%f',fact1_2l1s(tframe)); wstr11=sprintf('%f',fact2_2l1s(tframe)); wstr12=sprintf('%f',fact1_2r1s(tframe)); wstr13=sprintf('%f',fact2_2r1s(tframe)); if WT_QTY == 3 then sig1=x_mdialog('Set limit',txt1, [wstr0 ; wstr1 ; wstr2 ; wstr3 ; wstr4 ; wstr5 ; wstr6 ; wstr7 ; wstr8 ; wstr9 ]); if sig1==[] then act0=evstr(wstr1); else limit_switch0s(tframe)=evstr(sig1(1)); act0=evstr(sig1(2)); fact1_3l1s(tframe)=evstr(sig1(3)); fact2_3l1s(tframe)=evstr(sig1(4)); fact1_3r1s(tframe)=evstr(sig1(5)); fact2_3r1s(tframe)=evstr(sig1(6)); fact1_3l2s(tframe)=evstr(sig1(7)); fact2_3l2s(tframe)=evstr(sig1(8)); fact1_3r2s(tframe)=evstr(sig1(9)); fact2_3r2s(tframe)=evstr(sig1(10)); end else // 2 tube models // if WT_QTY == 2 then sig1=x_mdialog('Set limit',txt2, [wstr0 ; wstr1 ; wstr10 ; wstr11 ; wstr12 ; wstr13 ]); if sig1==[] then act0=evstr(wstr1); else limit_switch0s(tframe)=evstr(sig1(1)); act0=evstr(sig1(2)); fact1_2l1s(tframe)=evstr(sig1(3)); fact2_2l1s(tframe)=evstr(sig1(4)); fact1_2r1s(tframe)=evstr(sig1(5)); fact2_2r1s(tframe)=evstr(sig1(6)); end end endfunction //----------------------------------------------------------------- function edit_limit2() global fact1_2r1s fact2_2r1s fact1_2l1s fact2_2l1s fact1_rls fact2_rls limit_switch0s limit_switch2s limit_switch3s global fact1_3l1s fact2_3l1s fact1_3r1s fact2_3r1s fact1_3l2s fact2_3l2s fact1_3r2s fact2_3r2s txt1=['switch on(1)/off(0)';'rl plus';'rl minus']; wstr0=sprintf('%f',limit_switch2s(tframe)); wstr2=sprintf('%f',fact1_rls(tframe)); wstr3=sprintf('%f',fact2_rls(tframe)); sig1=x_mdialog('Set limit',txt1, [wstr0 ; wstr2 ; wstr3]); if sig1==[] then // else limit_switch2s(tframe)=evstr(sig1(1)); fact1_rls(tframe)=evstr(sig1(2)); fact2_rls(tframe)=evstr(sig1(3)); end endfunction //----------------------------------------------------------------- function edit_limit3() global limit_switch3s txt1=['switch on(1)/off(0)']; wstr0=sprintf('%f',limit_switch3s(tframe)); sig1=x_mdialog('Set force to add-offset control to equal maximum value',txt1, [wstr0]); if sig1==[] then // else limit_switch3s(tframe)=evstr(sig1(1)); end endfunction //----------------------------------------------------------------- function [switch7,act1_waku,act2_waku,switch8,act1_i2nd_thd,act2_i2nd_thd,switch9,act1_nA0,act2_nA0]= edit_limit7() txt1=['switch on(1)/off(0)';'waku are plus';'waku area minus';'switch on(1)/off(0)';'2nd/base freq factor plus';'2nd/base freq factor minus';'switch on(1)/off(0)';'noise amplitude plus';'noise amplitude minus']; wstr1=sprintf('%f',limit_switch7); wstr2=sprintf('%f',fact1_waku); wstr3=sprintf('%f',fact2_waku); wstr4=sprintf('%f',limit_switch8); wstr5=sprintf('%f',fact1_i2nd_thd); wstr6=sprintf('%f',fact2_i2nd_thd); wstr7=sprintf('%f',limit_switch9); wstr8=sprintf('%f',fact1_nA0); wstr9=sprintf('%f',fact2_nA0); sig1=x_mdialog('Set limit',txt1, [wstr1 ; wstr2 ; wstr3 ;str4 ; wstr5 ; wstr6; wstr7 ; wstr8 ; wstr9]); if sig1==[] then switch7=evstr(wstr1); act1_waku=evstr(wstr2); act2_waku=evstr(wstr3); switch8=evstr(wstr4); act1_i2nd_thd=evstr(wstr5); act2_i2nd_thd=evstr(wstr6); switch9=evstr(wstr7); act1_nA0=evstr(wstr8); act2_nA0=evstr(wstr9); else switch7=evstr(sig1(1)); act1_waku=evstr(sig1(2)); act2_waku=evstr(sig1(3)); switch8=evstr(sig1(4)); act1_i2nd_thd=evstr(sig1(5)); act2_i2nd_thd=evstr(sig1(6)); switch9=evstr(sig1(7)); act1_nA0=evstr(sig1(8)); act2_nA0=evstr(sig1(9)); end endfunction //----------------------------------------------------------------- function disp_limit_switch() // limit_switch0=0.0; // r1,r2,l1,l2 off if limit_switch0 >= 1. then disp(' r1,r2,l1,l2 limit control is ON'); else disp(' r1,r2,l1,l2 limit control is OFF'); end // limit_switch2=1.0; // rl on if limit_switch2 >= 1. then disp(' rl limit control is ON'); else disp(' rl limit control is OFF'); end if limit_switch3 >= 1. then disp(' force to add-offset control to equal maximum value is ON'); else disp(' force to add-offset control to equal maximum value is OFF'); end if WT_QTY==21 then if limit_switch7 >= 1. then disp(' noise_waku_area limit control is ON'); else disp(' noise_waku_area limit control is OFF'); end if limit_switch8 >= 1. then disp(' 2nd/base freq factor limit control is ON'); else disp(' 2nd/base freq factor limit control is OFF'); end if limit_switch9 >= 1. then disp(' nA0=mix_amplitude_indep limit control is ON'); else disp(' nA0=mix_amplitude_indep limit control is OFF'); end end if WT_QTY==41 then if limit_switch10 >= 1. then disp(' peq_space_area limit control is ON'); else disp(' peq_space_area limit control is OFF'); end if limit_switch11 >= 1. then disp(' peq amplitude limit control is ON'); else disp(' peq amplitude limit control is OFF'); end end if (WT_QTY==51) | (WT_QTY==52) then if limit_switch12 >= 1. then disp(' stable_noise_lpf_fc limit control is ON'); else disp(' stable_noise_lpf_fc limit control is OFF'); end end endfunction //----------------------------------------------------------------- function test_plot1exp(disp0) // disp_limit_switch(); wb0=xget('window'); // stack old window xset('window',disp0); // create new windows // dx0=[-1:0.1:1]; clf(); if WT_QTY == 3 then if limit_switch0 >= 1 then subplot(511); dx0=[(fact2_3l1s(tframe)-0.2):0.1:(fact1_3l1s(tframe)+0.2)]; plot2d( dx0,(1 + exp( fact0 * ( dx0 - fact1_3l1s(tframe)))) .* (1 + exp( fact0 * ( fact2_3l1s(tframe) - dx0))),style=[color("red")] ); wstr0='decrease value depend on l1, hoping l1 inside the limit '; xtitle(wstr0); subplot(512); dx0=[(fact2_3r1s(tframe)-0.2):0.1:(fact1_3r1s(tframe)+0.2)]; plot2d( dx0,(1 + exp( fact0 * ( dx0 - fact1_3r1s(tframe)))) .* (1 + exp( fact0 * ( fact2_3r1s(tframe) - dx0))),style=[color("red")] ); wstr0='decrease value depend on r1, hoping r1 inside the limit '; xtitle(wstr0); subplot(513); dx0=[(fact2_3l2s(tframe)-0.2):0.1:(fact1_3l2s(tframe)+0.2)]; plot2d( dx0,(1 + exp( fact0 * ( dx0 - fact1_3l2s(tframe)))) .* (1 + exp( fact0 * ( fact2_3l2s(tframe) - dx0))) ,style=[color("red")] ); wstr0='decrease value depend on l2, hoping l2 inside the limit '; xtitle(wstr0); subplot(514); dx0=[(fact2_3r2s(tframe)-0.2):0.1:(fact1_3r2s(tframe)+0.2)]; plot2d( dx0,(1 + exp( fact0 * ( dx0 - fact1_3r2s(tframe)))) .* (1 + exp( fact0 * ( fact2_3r2s(tframe) - dx0))),style=[color("red")] ); wstr0='decrease value depend on r2, hoping r2 inside the limit '; xtitle(wstr0); end // if limit_switch0 >= 1 then if limit_switch2 >= 1 then subplot(515); dx0=[(fact2_rls(tframe)-0.2):0.1:(fact1_rls(tframe)+0.2)]; plot2d( dx0,(1 + exp( fact0 * ( dx0 - fact1_rls(tframe)))) .* (1 + exp( fact0 * ( fact2_rls(tframe) - dx0))),style=[color("red")] ); wstr0='decrease value depend on rl, hoping rl inside the limit '; xtitle(wstr0); end // if limit_switch2 >= 1 then else //if WT_QTY == 2 then if limit_switch0 >= 1 then subplot(311); dx0=[(fact2_2l1s(tframe)-0.2):0.1:(fact1_2l1s(tframe)+0.2)]; plot2d( dx0,(1 + exp( fact0 * ( dx0 - fact1_2l1s(tframe)))) .* (1 + exp( fact0 * ( fact2_2l1s(tframe) - dx0))),style=[color("red")] ); wstr0='decrease value depend on l1, hoping l1 inside the limit '; xtitle(wstr0); subplot(312); dx0=[(fact2_2r1s(tframe)-0.2):0.1:(fact1_2r1s(tframe)+0.2)]; plot2d( dx0,(1 + exp( fact0 * ( dx0 - fact1_2r1s(tframe)))) .* (1 + exp( fact0 * ( fact2_2r1s(tframe) - dx0))),style=[color("red")] ); wstr0='decrease value depend on r1, hoping r1 inside the limit '; xtitle(wstr0); end //if limit_switch0 >= 1 then if limit_switch2 >= 1 then subplot(313); dx0=[(fact2_rls(tframe)-0.2):0.1:(fact1_rls(tframe)+0.2)]; plot2d( dx0,(1 + exp( fact0 * ( dx0 - fact1_rls(tframe)))) .* (1 + exp( fact0 * ( fact2_rls(tframe) - dx0))),style=[color("red")] ); wstr0='decrease value depend on rl, hoping rl inside the limit '; xtitle(wstr0); end //if limit_switch2 >= 1 then end xset('window',wb0); // push old windows endfunction //------------------------------------ function test_plot2exp(disp0) if WT_QTY==21 then if limit_switch7 >= 1. then disp(' noise_waku_area limit control is ON'); else disp(' noise_waku_area limit control is OFF'); end if limit_switch8 >= 1. then disp(' 2nd/base freq factor limit control is ON'); else disp(' 2nd/base freq factor limit control is OFF'); end if limit_switch9 >= 1. then disp(' nA0=mix_amplitude_indep limit control is ON'); else disp(' nA0=mix_amplitude_indep limit control is OFF'); end wb0=xget('window'); // stack old window xset('window',disp0); // create new windows // dx0=[-1:0.1:1]; clf(); if limit_switch7 >= 1 then subplot(311); dx0=[(fact2_waku-0.2):0.1:(fact1_waku+0.2)]; plot2d( dx0,(1 + exp( fact0 * ( dx0 - fact1_waku))) .* (1 + exp( fact0 * ( fact2_waku - dx0))),style=[color("red")] ); wstr0='decrease value depend on noise_waku_area, hoping noise_waku_area inside the limit '; xtitle(wstr0); end if limit_switch8 >= 1 then subplot(312); dx0=[(fact2_i2nd_thd-0.2):0.1:(fact1_i2nd_thd+0.2)]; plot2d( dx0,(1 + exp( fact0 * ( dx0 - fact1_i2nd_thd))) .* (1 + exp( fact0 * ( fact2_i2nd_thd - dx0))),style=[color("red")] ); wstr0='decrease value depend on 2nd/base freq, hoping 2nd/base freq inside the limit '; xtitle(wstr0); end if limit_switch9 >= 1 then subplot(313); dx0=[(fact2_nA0-0.02):0.01:(fact1_nA0+0.02)]; plot2d( dx0,(1 + exp( fact2 * ( dx0 - fact1_nA0))) .* (1 + exp( fact2 * ( fact2_nA0 - dx0))),style=[color("red")] ); wstr0='decrease value depend on nA0, hoping nA0 inside the limit '; xtitle(wstr0); end //if limit_switch9 >= 1 then xset('window',wb0); // push old windows end // if WT_QTY==21 then endfunction //---------------------------------------- // // //----------------------------------------------------------------- // y = tube2_resp9(tm8, x_init8) // // [frq1,sfrq1,is1,ie1]=set_frq(fft_point1); // clf(); // subplot(211); // gainplot(frq1,y,phi_fft1); // function y = tube2_resp9(xw, pa) // when 2 tube model L2= ( 1. + pa(2)) * pa(1) / 2.; // L2= ( 1. + l1 ) * ttl_Length / 2.; L1= pa(1) - L2; // L1= ttl_Length - L2; tu1=L1/c0; tu2=L2/c0; A2= ( 1. + pa(3)) * ttl_Area/ 2.; // A2= ( 1. + r1 ) * ttl_Area/ 2.; A1= ttl_Area - A2; yi = 0.5 * ( 1 + rg ) * ( 1 + pa(3)) * ( 1 + pa(4) ) * %e^( -1 * ( tu1 + tu2 ) .* xw * %i); yb = 1 + pa(3) * rg * %e^( -2 * tu1 .* xw * %i) + pa(3) * pa(4) * %e^( -2 * tu2 .* xw * %i) + pa(4) * rg * %e^( -2 * (tu1 + tu2) .* xw * %i); yc = yg_res9 .* hpf_res9 .* ( yi ./ yb) ; y = 20. * log10(abs(yc)); c0= max_wm9 * y; // c0 should be scalor // disp(c0); y=y - limit_switch3 * ( c0 - max_value ) * wm9; y= y - limit_switch0 * ((1 + exp( fact0 * ( pa(3) - fact1_2r1))) * (1 + exp( fact0 * ( fact2_2r1 - pa(3) ) ) )) * wm9; y= y - limit_switch0 * ((1 + exp( fact0 * ( pa(2) - fact1_2l1))) * (1 + exp( fact0 * ( fact2_2l1 - pa(2) ) ) )) * wm9; y= y - limit_switch2 * ((1 + exp( fact0 * ( pa(4) - fact1_rl))) * (1 + exp( fact0 * ( fact2_rl - pa(4) ) ) )) * wm9; endfunction //-------------------------------- // y = tube3_resp9(tm8, x_init8) // // [frq1,sfrq1,is1,ie1]=set_frq(fft_point1); // clf(); // subplot(211); // gainplot(frq1,y,phi_fft1); function y = tube3_resp9(xw, pa) // when 3 tube model serial bunbo= pa(2) * pa(4) + ( pa(2) - pa(4) ) + 3.0; // bunbo= l1 * l2 + ( l1 - l2 ) + 3.0; L1= pa(1) * ( pa(4) - 1.0 ) * (pa(2) - 1.0 ) / bunbo ; // L1= ttl_Length * ( l2 - 1.0 ) * (l1 - 1.0 ) / bunbo ; L2= pa(1) * ( pa(4) - 1.0 ) * ( -1.0 * pa(2) - 1.0 ) / bunbo ; // L2= ttl_Length * ( l2 - 1.0 ) * ( -1.0 * l1 - 1.0 ) / bunbo ; L3= pa(1) * ( pa(4) + 1.0 ) * (pa(2) + 1.0 ) / bunbo ; // L3= ttl_Length * ( l2 + 1.0 ) * (l1 + 1.0 ) / bunbo ; tu1=L1/c0; tu2=L2/c0; tu3=L3/c0; bunbo= pa(3) * pa(5) + ( pa(3) - pa(5) ) + 3.0; // bunbo= r1 * r2 + ( r1 - r2 ) + 3.0; A1= ttl_Area * ( pa(5) - 1.0 ) * (pa(3) - 1.0 ) / bunbo ; // A1= ttl_Area * ( r2 - 1.0 ) * (r1 - 1.0 ) / bunbo ; A2= ttl_Area * ( pa(5) - 1.0 ) * ( -1.0 * pa(3) - 1.0 ) / bunbo ; // A2= ttl_Area * ( r2 - 1.0 ) * ( -1.0 * r1 - 1.0 ) / bunbo ; A3= ttl_Area * ( pa(5) + 1.0 ) * (pa(3) + 1.0 ) / bunbo ; // A3= ttl_Area * ( r2 + 1.0 ) * (r1 + 1.0 ) / bunbo ; yi = 0.5 * ( 1. + rg ) * ( 1. + pa(3)) * ( 1. + pa(5)) * ( 1. + pa(6) ) * %e^( -1. * ( tu1 + tu2 + tu3) .* xw * %i); yb = 1 + rg * pa(3) * %e^( -2 * tu1 .* xw * %i) + pa(3) * pa(5) * %e^( -2 * tu2 .* xw * %i) + pa(5) * pa(6) * %e^( -2 * tu3 .* xw * %i); yb = yb + rg * pa(5) * %e^( -2 * (tu1 + tu2) .* xw * %i) + pa(3) * pa(6) * %e^( -2 * (tu2 + tu3) .* xw * %i); yb = yb + rg * pa(3) * pa(5) * pa(6) * %e^( -2 * (tu1 + tu3) .* xw * %i); yb = yb + rg * pa(6) * %e^( -2 * (tu1 + tu2 + tu3) .* xw * %i); yc = yg_res9 .* hpf_res9 .* ( yi ./ yb) ; //(1)test1 // for r1 // yc = ( 1. / ((1 + exp( fact0 * ( pa(3) - fact1_3r1))) * (1 + exp( fact0 * ( fact2_3r1 - pa(3))))) ) * yc; y = 20. * log10(abs(yc)); c0= max_wm9 * y; // c0 should be scalor // disp(c0); y=y - limit_switch3 * ( c0 - max_value ) * wm9; y= y - limit_switch0 * ((1 + exp( fact0 * ( pa(3) - fact1_3r1))) * (1 + exp( fact0 * ( fact2_3r1 - pa(3) ) ) )) * wm9; y= y - limit_switch0 * ((1 + exp( fact0 * ( pa(5) - fact1_3r2))) * (1 + exp( fact0 * ( fact2_3r2 - pa(5) ) ) )) * wm9; y= y - limit_switch0 * ((1 + exp( fact0 * ( pa(2) - fact1_3l1))) * (1 + exp( fact0 * ( fact2_3l1 - pa(2) ) ) )) * wm9; y= y - limit_switch0 * ((1 + exp( fact0 * ( pa(4) - fact1_3l2))) * (1 + exp( fact0 * ( fact2_3l2 - pa(4) ) ) )) * wm9; y= y - limit_switch2 * ((1 + exp( fact0 * ( pa(6) - fact1_rl))) * (1 + exp( fact0 * ( fact2_rl - pa(6) ) ) )) * wm9; endfunction //------------------------------------------------------------ function y = tube21_resp9(xw, pa) // when 2 tube model plus independent noise [frqx,q0x]= trans_waku_area( pa(5) ); [yd_bpf, xd_bpf]=set_bpf( (frqx * 2. * %pi), q0x,( pa(6) * frqx * 2. * %pi), q0x); t0=1/fs1; yi1 = yd_bpf(1,1) + yd_bpf(1,2) * %e^( -1.* t0 .* xw * %i) + yd_bpf(1,3) * %e^( -2.* t0 .* xw * %i); yb1 = 1.0 - (xd_bpf(1,2) * %e^( -1.* t0 .* xw * %i) + xd_bpf(1,3) * %e^( -2.* t0 .* xw * %i)); y1 = pa(7) * yi1 ./ yb1; yi2 = yd_bpf(2,1) + yd_bpf(2,2) * %e^( -1.* t0 * xw * %i) + yd_bpf(2,3) * %e^( -2.* t0 * xw * %i); yb2 = 1.0 - (xd_bpf(2,2) * %e^( -1.* t0 * xw * %i) + xd_bpf(2,3) * %e^( -2.* t0 * xw * %i)); y2 = pa(7) * yi2 ./ yb2; L2= ( 1. + pa(2)) * pa(1) / 2.; // L2= ( 1. + l1 ) * ttl_Length / 2.; L1= pa(1) - L2; // L1= ttl_Length - L2; tu1=L1/c0; tu2=L2/c0; A2= ( 1. + pa(3)) * ttl_Area/ 2.; // A2= ( 1. + r1 ) * ttl_Area/ 2.; A1= ttl_Area - A2; yi = 0.5 * ( 1 + rg ) * ( 1 + pa(3)) * ( 1 + pa(4) ) * %e^( -1 * ( tu1 + tu2 ) .* xw * %i); yb = 1 + pa(3) * rg * %e^( -2 * tu1 .* xw * %i) + pa(3) * pa(4) * %e^( -2 * tu2 .* xw * %i) + pa(4) * rg * %e^( -2 * (tu1 + tu2) .* xw * %i); yc1 = yg_res9 .* ( yi ./ yb) ; yc = hpf_res9 .* ( yc1 + y1 + y2); y = 20. * log10(abs(yc)); c0= max_wm9 * y; // c0 should be scalor // disp(c0); y=y - limit_switch3 * ( c0 - max_value ) * wm9; y= y - limit_switch0 * ((1 + exp( fact0 * ( pa(3) - fact1_2r1))) * (1 + exp( fact0 * ( fact2_2r1 - pa(3) ) ) )) * wm9; y= y - limit_switch0 * ((1 + exp( fact0 * ( pa(2) - fact1_2l1))) * (1 + exp( fact0 * ( fact2_2l1 - pa(2) ) ) )) * wm9; y= y - limit_switch2 * ((1 + exp( fact0 * ( pa(4) - fact1_rl))) * (1 + exp( fact0 * ( fact2_rl - pa(4) ) ) )) * wm9; y= y - limit_switch7 * ((1 + exp( fact0 * ( pa(5) - fact1_waku))) * (1 + exp( fact0 * ( fact2_waku - pa(5) ) ) )) * wm9; y= y - limit_switch8 * ((1 + exp( fact0 * ( pa(6) - fact1_i2nd_thd))) * (1 + exp( fact0 * ( fact2_i2nd_thd - pa(6) ) ) )) * wm9; y= y - limit_switch9 * ((1 + exp( fact2 * ( pa(7) - fact1_nA0))) * (1 + exp( fact2 * ( fact2_nA0 - pa(7) ) ) )) * wm9; endfunction //------------------------------------------------------------ function y = tube41_resp9(xw, pa) // when 2 tube model plus partly emphasis [emp_frq_center,emp_q0]= trans_space_area( pa(5)) [yd_peq, xd_peq]=set_peq( (emp_frq_center * 2. * %pi), pa(6), emp_q0); t0=1/fs1; vl=1; yi1 = yd_peq(vl,1) + yd_peq(vl,2) * %e^( -1.* t0 * xw * %i) + yd_peq(vl,3) * %e^( -2.* t0 * xw * %i); yb1 = 1.0 - (xd_peq(vl,2) * %e^( -1.* t0 * xw * %i) + xd_peq(vl,3) * %e^( -2.* t0 * xw * %i)); y1 = yi1 ./ yb1; L2= ( 1. + pa(2)) * pa(1) / 2.; // L2= ( 1. + l1 ) * ttl_Length / 2.; L1= pa(1) - L2; // L1= ttl_Length - L2; tu1=L1/c0; tu2=L2/c0; A2= ( 1. + pa(3)) * ttl_Area/ 2.; // A2= ( 1. + r1 ) * ttl_Area/ 2.; A1= ttl_Area - A2; yi = 0.5 * ( 1 + rg ) * ( 1 + pa(3)) * ( 1 + pa(4) ) * %e^( -1 * ( tu1 + tu2 ) .* xw * %i); yb = 1 + pa(3) * rg * %e^( -2 * tu1 .* xw * %i) + pa(3) * pa(4) * %e^( -2 * tu2 .* xw * %i) + pa(4) * rg * %e^( -2 * (tu1 + tu2) .* xw * %i); yc1 = yg_res9 .* ( yi ./ yb) ; yc = hpf_res9 .* yc1 .*y1 y = 20. * log10(abs(yc)); c0= max_wm9 * y; // c0 should be scalor // disp(c0); y=y - limit_switch3 * ( c0 - max_value ) * wm9; y= y - limit_switch0 * ((1 + exp( fact0 * ( pa(3) - fact1_2r1))) * (1 + exp( fact0 * ( fact2_2r1 - pa(3) ) ) )) * wm9; y= y - limit_switch0 * ((1 + exp( fact0 * ( pa(2) - fact1_2l1))) * (1 + exp( fact0 * ( fact2_2l1 - pa(2) ) ) )) * wm9; y= y - limit_switch2 * ((1 + exp( fact0 * ( pa(4) - fact1_rl))) * (1 + exp( fact0 * ( fact2_rl - pa(4) ) ) )) * wm9; y= y - limit_switch10 * ((1 + exp( fact0 * ( pa(5) - fact1_space))) * (1 + exp( fact0 * ( fact2_space - pa(5) ) ) )) * wm9; y= y - limit_switch11 * ((1 + exp( fact0 * ( pa(6) - fact1_empA0))) * (1 + exp( fact0 * ( fact2_empA0 - pa(6) ) ) )) * wm9; endfunction //------------------------------------------------------------ function y = tube51_resp9(xw, pa) // when 2 tube model driven by stable noise [yd_lpf, xd_lpf]=set_lpf( pa(5) * 2. * %pi , q0_stable_noise,xorder_stable_noise ); [yd_hpf, xd_hpf]=set_hpf( stable_noise_hpf * 2. * %pi , q0_stable_noiseh,xorder_stable_noiseh ); t0=1/fs1; vl=1; yi1 = yd_lpf(vl,1) + yd_lpf(vl,2) * %e^( -1.* t0 * xw * %i) + yd_lpf(vl,3) * %e^( -2.* t0 * xw * %i); yb1 = 1.0 - (xd_lpf(vl,2) * %e^( -1.* t0 * xw * %i) + xd_lpf(vl,3) * %e^( -2.* t0 * xw * %i)); y1 = yi1 ./ yb1; yi2 = yd_hpf(vl,1) + yd_hpf(vl,2) * %e^( -1.* t0 * xw * %i) + yd_hpf(vl,3) * %e^( -2.* t0 * xw * %i); yb2 = 1.0 - (xd_hpf(vl,2) * %e^( -1.* t0 * xw * %i) + xd_hpf(vl,3) * %e^( -2.* t0 * xw * %i)); y2 = yi2 ./ yb2; L2= ( 1. + pa(2)) * pa(1) / 2.; // L2= ( 1. + l1 ) * ttl_Length / 2.; L1= pa(1) - L2; // L1= ttl_Length - L2; tu1=L1/c0; tu2=L2/c0; A2= ( 1. + pa(3)) * ttl_Area/ 2.; // A2= ( 1. + r1 ) * ttl_Area/ 2.; A1= ttl_Area - A2; yi = 0.5 * ( 1 + rg ) * ( 1 + pa(3)) * ( 1 + pa(4) ) * %e^( -1 * ( tu1 + tu2 ) .* xw * %i); yb = 1 + pa(3) * rg * %e^( -2 * tu1 .* xw * %i) + pa(3) * pa(4) * %e^( -2 * tu2 .* xw * %i) + pa(4) * rg * %e^( -2 * (tu1 + tu2) .* xw * %i); yc1 = ( yi ./ yb) ; yc = hpf_res9 .* yc1 .* y1 .* y2 y = 20. * log10(abs(yc)); c0= max_wm9 * y; // c0 should be scalor // disp(c0); y=y - limit_switch3 * ( c0 - max_value ) * wm9; y= y - limit_switch0 * ((1 + exp( fact0 * ( pa(3) - fact1_2r1))) * (1 + exp( fact0 * ( fact2_2r1 - pa(3) ) ) )) * wm9; y= y - limit_switch0 * ((1 + exp( fact0 * ( pa(2) - fact1_2l1))) * (1 + exp( fact0 * ( fact2_2l1 - pa(2) ) ) )) * wm9; y= y - limit_switch2 * ((1 + exp( fact0 * ( pa(4) - fact1_rl))) * (1 + exp( fact0 * ( fact2_rl - pa(4) ) ) )) * wm9; y= y - limit_switch12 * ((1 + exp( fact0 * ( pa(5) - fact1_sn_lpf_fc))) * (1 + exp( fact0 * ( fact2_sn_lpf_fc - pa(5) ) ) )) * wm9; endfunction //------------------------------------------------------------ function y = tube52_resp9(xw, pa) // when 2 tube model driven by burst noise [yd_bpf, xd_bpf]=set_bpf( pa(5) * 2. * %pi , q0_burst_noise, pa(5) * 2. * %pi , q0_burst_noise); t0=1/fs1; vl=1; yi1 = yd_bpf(vl,1) + yd_bpf(vl,2) * %e^( -1.* t0 * xw * %i) + yd_bpf(vl,3) * %e^( -2.* t0 * xw * %i); yb1 = 1.0 - (xd_bpf(vl,2) * %e^( -1.* t0 * xw * %i) + xd_bpf(vl,3) * %e^( -2.* t0 * xw * %i)); y1 = yi1 ./ yb1; L2= ( 1. + pa(2)) * pa(1) / 2.; // L2= ( 1. + l1 ) * ttl_Length / 2.; L1= pa(1) - L2; // L1= ttl_Length - L2; tu1=L1/c0; tu2=L2/c0; A2= ( 1. + pa(3)) * ttl_Area/ 2.; // A2= ( 1. + r1 ) * ttl_Area/ 2.; A1= ttl_Area - A2; yi = 0.5 * ( 1 + rg ) * ( 1 + pa(3)) * ( 1 + pa(4) ) * %e^( -1 * ( tu1 + tu2 ) .* xw * %i); yb = 1 + pa(3) * rg * %e^( -2 * tu1 .* xw * %i) + pa(3) * pa(4) * %e^( -2 * tu2 .* xw * %i) + pa(4) * rg * %e^( -2 * (tu1 + tu2) .* xw * %i); yc1 = ( yi ./ yb) ; yc = hpf_res9 .* yc1 .*y1 y = 20. * log10(abs(yc)); c0= max_wm9 * y; // c0 should be scalor // disp(c0); y=y - limit_switch3 * ( c0 - max_value ) * wm9; y= y - limit_switch0 * ((1 + exp( fact0 * ( pa(3) - fact1_2r1))) * (1 + exp( fact0 * ( fact2_2r1 - pa(3) ) ) )) * wm9; y= y - limit_switch0 * ((1 + exp( fact0 * ( pa(2) - fact1_2l1))) * (1 + exp( fact0 * ( fact2_2l1 - pa(2) ) ) )) * wm9; y= y - limit_switch2 * ((1 + exp( fact0 * ( pa(4) - fact1_rl))) * (1 + exp( fact0 * ( fact2_rl - pa(4) ) ) )) * wm9; y= y - limit_switch12 * ((1 + exp( fact0 * ( pa(5) - fact1_sn_lpf_fc))) * (1 + exp( fact0 * ( fact2_sn_lpf_fc - pa(5) ) ) )) * wm9; endfunction //----------------------------------------------------------------- // e=myfun(x_init8, tm8, ym8, wm8) function e = myfun(x, tm, ym, wm) // when 2 tube model e = wm .*( tube2_resp9(tm, x) - ym ) endfunction function e = myfun2(x, tm, ym, wm) // when 3 tube model serial e = wm .*( tube3_resp9(tm, x) - ym ) endfunction function e = myfun21(x, tm, ym, wm) // when 2 tube model plus independent noise e = wm .*( tube21_resp9(tm, x) - ym ) endfunction function e = myfun41(x, tm, ym, wm) // when 2 tube model plus partly emphasis e = wm .*( tube41_resp9(tm, x) - ym ) endfunction function e = myfun51(x, tm, ym, wm) // when 2 tube model in stable noise e = wm .*( tube51_resp9(tm, x) - ym ) endfunction function e = myfun52(x, tm, ym, wm) // when 2 tube model in stable noise e = wm .*( tube52_resp9(tm, x) - ym ) endfunction //----------------------------------------------------------------- function [fopt,xopt, grdopt]=leastsq_main1s() global yg_res9 hpf_res9 global x_init8 wm9 global xopts global fact1_2r1 fact2_2r1 fact1_2l1 fact2_2l1 fact1_rl fact2_rl limit_switch0 limit_switch2 global fact1_3l1 fact2_3l1 fact1_3r1 fact2_3r1 fact1_3l2 fact2_3l2 fact1_3r2 fact2_3r2 global limit_switch3 [tm8, ym8, x_init8, wm9]=prepare8(); [yg_res9, hpf_res9]=prepare9(); // // check weight [frq1,sfrq1,is1,ie1]=set_frq(fft_point1); wm8=zeros(sfrq1,1); for v=1:sfrq1 wm8(v)=wm8s(tframe,v); end a=0.; s0=size(wm8); wm8c=ones(s0(1),1); for v=1:s0(1) a=a+abs( wm8(v)); end if a == 0. then disp('+WARNING: set weight as all 1, because wm8 are all 0'); else for v=1:s0(1) wm8c(v)=wm8(v); end end fact1_2r1=fact1_2r1s(tframe); fact2_2r1=fact2_2r1s(tframe); fact1_2l1=fact1_2l1s(tframe); fact2_2l1=fact2_2l1s(tframe); fact1_rl=fact1_rls(tframe); fact2_rl=fact2_rls(tframe); limit_switch0=limit_switch0s(tframe); limit_switch2=limit_switch2s(tframe); fact1_3r1=fact1_3r1s(tframe); fact2_3r1=fact2_3r1s(tframe); fact1_3l1=fact1_3l1s(tframe); fact2_3l1=fact2_3l1s(tframe); fact1_3r2=fact1_3r2s(tframe); fact2_3r2=fact2_3r2s(tframe); fact1_3l2=fact1_3l2s(tframe); fact2_3l2=fact2_3l2s(tframe); limit_switch3=limit_switch3s(tframe); disp_limit_switch(); // 1- the simplest call if WT_QTY == 3 then // when 3 tube model serial [fopt,xopt, grdopt] = leastsq(list(myfun2,tm8,ym8,wm8c),x_init8); elseif WT_QTY == 21 then // when 2 tube model plus independent noise [fopt,xopt, grdopt] = leastsq(list(myfun21,tm8,ym8,wm8c),x_init8); elseif WT_QTY == 41 then // when 2 tube model plus partly emphasis [fopt,xopt, grdopt] = leastsq(list(myfun41,tm8,ym8,wm8c),x_init8); elseif WT_QTY == 51 then // when 2 tube model in stable noise (WT_QTY == 51) [fopt,xopt, grdopt] = leastsq(list(myfun51,tm8,ym8,wm8c),x_init8); elseif WT_QTY == 52 then // when 2 tube model in burst noise (WT_QTY == 52) [fopt,xopt, grdopt] = leastsq(list(myfun52,tm8,ym8,wm8c),x_init8); else // when 2 tube model [fopt,xopt, grdopt] = leastsq(list(myfun,tm8,ym8,wm8c),x_init8); end //... call s0=size(xopt); s0s=s0(1); xopts(tframe,1)=s0s; for v=1:s0s xopts(tframe,v+1)=xopt(v); // first (*,1) is qty of the datas end disp( 'initial value -> result value'); for v=1:s0s wstr1=sprintf('%f',xopt(v)); wstr2=sprintf('%f',x_init8(v)); if WT_QTY == 3 then wstr3=' ' + ma0(v) + wstr2 + ' -> ' + wstr1; elseif WT_QTY == 21 then wstr3=' ' + ma2(v) + wstr2 + ' -> ' + wstr1; elseif WT_QTY == 41 then wstr3=' ' + ma3(v) + wstr2 + ' -> ' + wstr1; elseif WT_QTY == 51 then wstr3=' ' + ma4(v) + wstr2 + ' -> ' + wstr1; elseif WT_QTY == 52 then wstr3=' ' + ma5(v) + wstr2 + ' -> ' + wstr1; else wstr3=' ' + ma1(v) + wstr2 + ' -> ' + wstr1; end disp(wstr3); end disp(' '); endfunction //-----------------------+++ ---- ++ ------------------------------ function [wm8s]=reset_wm8() [frq1,sfrq1,is1,ie1]=set_frq(fft_point1); wm8s=zeros(max_frames,sfrq1); // draw all 1 disp('+set weight all 0'); endfunction //---------------- function plot_wm8(disp0) wb0=xget('window'); // stack old window xset('window',disp0); // create new windows clf(); sgamen= QTY_FRAME * 2; s23=size(db_fft1s); db_fft0=zeros(1,s23(2)); phi_fft0=zeros(1,s23(2)); for v=1:QTY_FRAME for w=1:s23(2) db_fft0(1,w)=db_fft1s(v,w); phi_fft0(1,w)=phi_fft1s(v,w); end [frq1,sfrq1,is1,ie1]=set_frq(0); subplot( sgamen,1,2*v-1); gainplot(frq1,db_fft0,phi_fft0); wstr1=PART_LIST0(v) + ' frequency response of waveform selected by fft analysis'; xtitle(wstr1); subplot( sgamen,1,2*v); [frq1,sfrq1,is1,ie1]=set_frq(fft_point1); wm8=zeros(sfrq1,1); for w=1:sfrq1 wm8(w)=wm8s(v,w); end plot2d(frq1,wm8',style=[color("red")],logflag="ln",rect=[frq1(1),-1,frq1(sfrq1),(max(wm8)+1)]); wstr0=' wmf: weight for leastsq evalution '; xtitle(wstr0); end // v=1:QTY_FRAME xset('window',wb0); // push old windows endfunction //---------------------------------------------------------------- function [wm8c]= set_wm8( arg1, arg2, arg3, xframe) [frq1,sfrq1,is1,ie1]=set_frq(fft_point1); wm8c=zeros(max_frames,sfrq1); for v2=1:max_frames for v=1:sfrq1 wm8c(v2,v)=wm8s(v2,v); end end for v=1:sfrq1 if (frq1(v) >= arg1) & (frq1(v) <= arg2) then wm8c(xframe, v)=arg3; end end endfunction //--------------------------------------------- function [wm8c]= edit_wm8() txt1=['frame no.';'from (frequency)';'to (frequency)';'weight value']; sp1=0.; ep1=0.; vl0=1.; wstr0=sprintf('%d',tframe); wstr1=sprintf('%f',sp1); wstr2=sprintf('%f',ep1); wstr3=sprintf('%f',vl0); sig1=x_mdialog('Set weight value',txt1, [wstr0 ; wstr1 ; wstr2 ; wstr3]); if sig1==[] then arg4=evstr(wstr0); arg1=evstr(wstr1); arg2=evstr(wstr2); arg3=evstr(wstr3); else arg4=evstr(sig1(1)); arg1=evstr(sig1(2)); arg2=evstr(sig1(3)); arg3=evstr(sig1(4)); end [wm8c]= set_wm8( arg1, arg2, arg3, arg4); endfunction function check_out1() ////ma0=['ttl_Length ' ; 'l1 ' ; 'r1 ' ; 'l2 ' ; 'r2 ' ; 'rl ']; ////ma1=['ttl_Length ' ; 'l1 ' ; 'r1 ' ; 'rl ']; ////ma2=['ttl_Length ' ; 'l1 ' ; 'r1 ' ; 'rl ' ;'wake_area ' ; '2nd thd '; 'amplitude ']; disp( 'result value -> initial value'); wstr1=sprintf('%f',ttl_Length); wstr3=' ' + ma0(1) + wstr1; disp(wstr3); wstr1=sprintf('%f',l1); wstr3=' ' + ma0(2) + wstr1; disp(wstr3); wstr1=sprintf('%f',r1); wstr3=' ' + ma0(3) + wstr1; disp(wstr3); if WT_QTY == 3 then wstr1=sprintf('%f',l2); wstr3=' ' + ma0(4) + wstr1; disp(wstr3); wstr1=sprintf('%f',r2); wstr3=' ' + ma0(5) + wstr1; disp(wstr3); wstr1=sprintf('%f',rl); wstr3=' ' + ma0(6) + wstr1; disp(wstr3); elseif WT_QTY == 21 then wstr1=sprintf('%f',rl); wstr3=' ' + ma2(4) + wstr1; disp(wstr3); wstr1=sprintf('%f',noise_waku_area); wstr3=' ' + ma2(5) + wstr1; disp(wstr3); wstr1=sprintf('%f',i2nd_thd_factor); wstr3=' ' + ma2(6) + wstr1; disp(wstr3); wstr1=sprintf('%f',nA0); wstr3=' ' + ma2(7) + wstr1; disp(wstr3); elseif WT_QTY == 41 then wstr1=sprintf('%f',rl); wstr3=' ' + ma3(4) + wstr1; disp(wstr3); wstr1=sprintf('%f',emp_space_area); wstr3=' ' + ma3(5) + wstr1; disp(wstr3); wstr1=sprintf('%f',empA0); wstr3=' ' + ma3(6) + wstr1; disp(wstr3); end disp(' '); endfunction //----------------------------------------------- function plot_2tubeS(disp0, kcode) // When kcode =0, normal, beside kcode=1, leastsq result wb0=xget('window'); // stack old window xset('window',disp0); // create new windows [frq1,sfrq1,is1,ie1]=set_frq(0); clf(); subplot(311); s23=size(db_fft1s); db_fft0=zeros(1,s23(2)); phi_fft0=zeros(1,s23(2)); v=tframe; for w=1:s23(2) db_fft0(1,w)=db_fft1s(v,w); phi_fft0(1,w)=phi_fft1s(v,w); end gainplot(frq1,db_fft0,phi_fft0); wstr1=PART_LIST0(tframe) + ' frequency response of waveform selected by fft analysis' xtitle(wstr1); [frq1,sfrq1,is1,ie1]=set_frq(1024); subplot(312); kcode=0; exchange9(); // set initial value as tube paras cal_overall_response(); gainplot(frq1,db_2tube',phi_2tube'); if kcode == 0 then xtitle('frequency response of tubes model by initial value for comparison to one of waveform selected by fft analysis '); elseif kcode == 1 then xtitle('frequency response of tubes model by leastsq method result from the initial value'); elseif kcode == 3 then wstr0=' frequency response of tubes model'; wstr1=sprintf('%f',l1); wstr2=sprintf('%f',r1); wstr3=sprintf('%f',l2); wstr4=sprintf('%f',r2); if WT_QTY == 3 then wstr5='l1=' + wstr1 + ' r1=' + wstr2 + ' l2=' + wstr3 + ' r2=' + wstr4 + wstr0; else // when WT_QTY == 2 then wstr5='l1=' + wstr1 + ' r1=' + wstr2 + wstr0; end xtitle(wstr5); end subplot(313); kcode=1; exchange8(); // set leastsq result as tube paras cal_overall_response(); gainplot(frq1,db_2tube',phi_2tube'); if kcode == 0 then xtitle('frequency response of tubes model by initial value for comparison to one of waveform selected by fft analysis '); elseif kcode == 1 then xtitle('frequency response of tubes model by leastsq method result from the initial value'); elseif kcode == 3 then wstr0=' frequency response of tubes model'; wstr1=sprintf('%f',l1); wstr2=sprintf('%f',r1); wstr3=sprintf('%f',l2); wstr4=sprintf('%f',r2); if WT_QTY == 3 then wstr5='l1=' + wstr1 + ' r1=' + wstr2 + ' l2=' + wstr3 + ' r2=' + wstr4 + wstr0; else // when WT_QTY == 2 then wstr5='l1=' + wstr1 + ' r1=' + wstr2 + wstr0; end xtitle(wstr5); end xset('window',wb0); // push old windows endfunction //---- function [wtframe]=edit_tframe() txt1=['frame no.']; wstr0=sprintf('%d',tframe); sig1=x_mdialog('Set present frame for analysis',txt1, [wstr0 ]); if sig1==[] then wtframe=evstr(wstr0); else wtframe=evstr(sig1(1)); end endfunction //---- function [WT_QTY]=push0(wtframe0) global r1 r2 l1 l2 ttl_Length rl noise_waku_area i2nd_thd_factor nA0 global emp_space_area empA0 global tframe global stable_noise_lpf global sp1 ep1 tframe=wtframe0; // push WT_QTY=WT_QTYs(wtframe0); ttl_Length=ttl_Lengths(wtframe0); l1=l1s(wtframe0); r1=r1s(wtframe0); l2=l2s(wtframe0); r2=r2s(wtframe0); rl=rls(wtframe0); noise_waku_area=noise_waku_areas(wtframe0); i2nd_thd_factor=i2nd_thd_factors(wtframe0); nA0=nA0s(wtframe0); emp_space_area=emp_space_areas(wtframe0); empA0=empA0s(wtframe0); stable_noise_lpf=stable_noise_lpfs(wtframe0); sp1=sp1s(wtframe0); ep1=ep1s(wtframe0); endfunction // //================================================================= // // +MAIN (4)program starts // leastsq method to estimate tube model parameters if T_DEMO==1 & (T_SKIP1 <> 1) then disp(' '); disp('a trial of leastsq method to estimate tube model parameters including rl using focuus weight.'); disp(' '); // once set wm8 as all ones [wm8s]=reset_wm8(); if SEL_CODE == 1 then // ka_sample // frame 1 ? // frame 2 ? [wm8s]= set_wm8( 430.66406 , 6000. , 1.0 ,3); [wm8s]= set_wm8( 430.66406 , 6000. , 1.0 ,4); //[wm8s]= set_wm8( 430.66406 , 6000. , 1.0 ,5); // // this method to estimate 'a' is failure one... [wm8s]= set_wm8( 430.66406 , 1550.3906 , 1.0 ,5); [wm8s]= set_wm8( 1894.921 , 3100.781 , 1.0 ,5); [wm8s]= set_wm8( 3445.3125 , 5512.5 , 1.0 ,5); plot_wm8(6); elseif SEL_CODE == 2 then // o_sample [wm8s]= set_wm8( 172.26562, 1200. , 1.0, 1 ); [wm8s]= set_wm8( 1894.9219, 3703.7109, 1.0, 1); plot_wm8(6); elseif SEL_CODE == 3 then // i_sample [wm8s]= set_wm8( 200., 7000. , 1.0, 1 ); // teiiki slop wo dasutame, teiiki ha hikaku no sample-suu ga sukunakunai-tame, kawarini weigh wo masu, Plus limit rl <=0.9 mo hituyou // this slop fitting is not well yet. disp(' Waring: This Slop fitting method of lower frequency portion is not good. Another good idea should need.'); [wm8s]= set_wm8( 200., 350. , 3.0, 1 ); plot_wm8(6); else // include SEL_CODE == 999 plot_wm8(6); [wm8s]= edit_wm8(); plot_wm8(6); [wm8s]= edit_wm8(); plot_wm8(6); end test_plot1exp(8); test_plot2exp(9); tframe0=tframe; // pop for vloop=1:QTY_FRAME tframe=vloop; WT_QTY=WT_QTYs(tframe); wstr1="--> present part is " + PART_LIST0(vloop); disp(wstr1); [fopt,xopt, grdopt]=leastsq_main1s(); // set result value plot_2tubeS(9+vloop,1); end //for vloop=1:QTY_FRAME [WT_QTY]=push0(tframe0); [ tube2_res1, L1, L2, L3, A1, A2, A3]= two_tubes2_resp(); // end // // // -MAIN (4)program starts // //==============================================================++= //================================================================= // // +MAIN (5)program starts // generation waveform if T_DEMO==1 then if SEL_CODE == 1 then // ka_sample // set each frame count length // part1 * part2 * * part3 is // frame_lengths(1)=2; // frame_lengths(2)=3; frame_lengths(1)=1; frame_lengths(2)=1; frame_lengths(3)=2; frame_lengths(4)=3; //QTY_FRAME=5 amplitudes(1)=0.2; amplitudes(2)=0.2; amplitudes(3)=1.0; amplitudes(4)=1.0; amplitudes(5)=1.0; elseif SEL_CODE == 2 then // ?_sample // elseif SEL_CODE == 3 then // ?_sample // else // include SEL_CODE == 999 // end mess_init0=' by initial value'; [yg,y2tm,y3out,y2tm_lpf,y2tm_noise,LL]= make_waveform(4); plot_waveform(5); //+for filter bank process xsp1=500; // for k_sample xep1=1500; // for k_sample if LOAD_FB_FLAG == 1 then [fb_out_init, fb_am_out_init, fb_am_out_prop_init]=load_fb2(); sp1bak=sp1; ep1bak=ep1; sp1=xsp1; ep1=xep1; fb_outbak=fb_out; fb_am_outbak=fb_am_out; fb_am_out_propbak=fb_am_out_prop; fb_out=fb_out_init; fb_am_out=fb_am_out_init; fb_am_out_prop=fb_am_out_prop_init; wdat1bak=wdat1; wdat1=y3out; xmessfb=' ka_sample original (500-1500/5)'; plot_fb(34,xmessfb); sp1=sp1bak; ep1=ep1bak; fb_out=fb_outbak; fb_am_out=fb_am_outbak; fb_am_out_prop=fb_am_out_propbak; wdat1=wdat1bak; else do_fb3_init(y3out,xsp1,xep1,34,'waveform by initial value (500-1500/5)'); if LOAD_FB_FLAG == -1 then save_fb2(); end end start_fb_no=1; end_fb_no=5; ////[y3out_replace]=am_conv0_iso(start_fb_no, end_fb_no, xsp1, xep1, 36, ' COPY org wave AM to init wave'); // from fb no.1 to fb no.5 //-for filter bank process if T_SKIP1 <> 1 then plot_area_mkwave_least(21,2); // plot area by least, mkwave=>y3out_least end end // // -MAIN (5)program starts // //==============================================================++= // // // // // // //=======$======&==========================#===================== // // Add menu buttons of which name are 'a_plotwav' and '' // ADD_MENU_PLOT=1; // set ADD_MENU 1 to add menu button of some functions // // if ADD_MENU_PLOT == 1 then delmenu('step1_plotwav'); addmenu('step1_plotwav',[ '(0)read_wav' ; '(1)set_et_plot()']); step1_plotwav=[ '[wavfile1,chs1,qty1,fs1,bit1,f412,SEL_CODE]=get_wavfile_pro(); disp_pro_wav(); [wdat1, size1]=read_one_ch_of_wav(wavfile1); sp1=1; ep1=size1; [ydat0,xziku0]=make_width_data( wdat1, size1);' ; '[sp1, ep1, tframe ]= set_sp1_ep1(size1); sp1s(tframe)=sp1; ep1s(tframe)=ep1; plot_wave1s(0);' ] ; end //if ADD_MENU_PLOT == 1 then // //----------------------------------------------------------------- // ADD_MENU_FFT=1; // set ADD_MENU 1 to add menu button of some functions // // if ADD_MENU_FFT == 1 then delmenu('step2_fft'); addmenu('step2_fft',[ '(1)set points'; '(2)fft cal'; '(3)plot smooth, peak candinate' ]); step2_fft=[ '[fft_point1]=set_fft_points()' ; '[db_fft1, phi_fft1]=do_fft_wav(); plot_fft1(1)'; ' smooths(2,1); ' ] ; end //if ADD_MENU_FFT == 1 then // //----------------------------------------------------------------- // ADD_MENU_TUBE2=1; // set ADD_MENU 1 to add menu button of some functions // // if ADD_MENU_TUBE2 == 1 then delmenu('step3_tubes'); addmenu('step3_tubes',[ '(0)set frame' ; '(1)set_tube_model()'; '(2)set_tube_initial_para()'; '(3)plot_tube_freq(3)'; 'plot_area(7)']); step3_tubes=[ ' [tframe0]=edit_tframe(); [WT_QTY]=push0(tframe0);' ; '[WT_QTY]= set_tube_model(); WT_QTYs(tframe)=WT_QTY; ' ; '[fc,trise,tfall,tclosed]= set_2tube_para();'; ' [WT_QTY]=push0(tframe); cal_overall_response(); plot_2tube(3,0);' ; 'disp(tframe,''tframe is''); [WT_QTY]=push0(tframe); [ tube2_res1, L1, L2, L3, A1, A2, A3]= two_tubes2_resp(); plot_area(7)'] ; end //if ADD_MENU_TUBE2 == 1 then // //----------------------------------------------------------------- // ADD_MENU_LEASTSQURE1=1; // set ADD_MENU 1 to add menu button of some functions // // if ADD_MENU_LEASTSQURE1 == 1 then delmenu('step4_1_limit') addmenu('step4_1_limit',[ '(0)set_weight_all_0' ; '(1)edit_plot_weight(6)' ; '(2)edit_limit' ;'(2-2)edit_limit of ind noise']); step4_1_limit=['[wm8s]=reset_wm8();' ; 'plot_wm8(6); [wm8s]= edit_wm8(); plot_wm8(6);' ; '[fact0]= edit_limit0(); edit_limit2(); edit_limit3(); test_plot1exp(8);' ; ' [limit_switch7,fact1_waku,fact2_waku,limit_switch8,fact1_i2nd_thd,fact2_i2nd_thd,limit_switch9,fact1_nA0,fact2_nA0]= edit_limit7(); test_plot2exp(9);' ]; end //----------------------------------------------------------------- // ADD_MENU_LEASTSQURE2=1; // set ADD_MENU 1 to add menu button of some functions // // if ADD_MENU_LEASTSQURE2 == 1 then delmenu('step4_2_leastsq'); addmenu('step4_2_leastsq',[ '(1)do_leastsq' ; '(2)plot_freq(9+tframe)' ; '(3)set_result_as_initial' ;'plot_area(21)' ; 'make_waveform(21)' ; 'y3out_least_snd_play1()'; 'y3out_least_snd_save1()']); step4_2_leastsq=[' [fopt,xopt, grdopt]=leastsq_main1s();' ; ' plot_2tubeS((9+tframe),1)' ; 'exchange8s(); disp(''set result value as initial done.'');' ;'plot_area_mkwave_least(21,1); ' ; 'plot_area_mkwave_least(21,2);' ; 'y3bz=y3out; y3out=y3out_least; y3out_snd_play1(); y3out=y3bz;'; 'y3bz=y3out; y3out=y3out_least; y3out_snd_save1(); y3out=y3bz;' ]; end // //----------------------------------------------------------------- // ADD_MENU_GEN=1; // set ADD_MENU 1 to add menu button of some functions // // if ADD_MENU_FFT == 1 then delmenu('step5Generat_ByInitValue'); if (f_win == 1) | (scilab_version_number >= 4) then addmenu('step5Generat_ByInitValue',[ '(1)set_section_qty()' ; '(2)make_waveform()' ; 'plot_waveform()'; 'y3out_snd_play1()'; 'y3out_snd_save1()']); step5Generat_ByInitValue=[ '[N_REPEAT]= set_section_qty()'; '[yg,y2tm,y3out,y2tm_lpf,y2tm_noise,LL]= make_waveform(4);' ; 'plot_waveform(5)'; 'y3out_snd_play1()'; 'y3out_snd_save1()';] ; else addmenu('step5Generat_ByInitValue',[ '-' ; '(2)make_waveform()' ; 'plot_waveform()'; 'y3out_snd_save1()' ]); step5Generat_ByInitValue=[ ' '; '[yg,y2tm,y3out,y2tm_lpf,LL],y2tm_noise= make_waveform(4);' ; 'plot_waveform(5)'; 'y3out_snd_save1()'] ; end end //if ADD_MENU_GEN == 1 then // // //----------------------------------------------------------------- // if f_win == 1 then ADD_MENU_SND=0; else ADD_MENU_SND=0; // set ADD_MENU 1 to add menu button of some functions end // // if ADD_MENU_SND == 1 then delmenu('sound_play'); addmenu('sound_play',[ 'snd_play1()' ; 'snd_save1()' ]); sound_play=[ 'snd_play1()' ; 'snd_save1()' ] ; end //if ADD_MENU_SND == 1 then // // // //----------------------------------------------------------------- // ADD_MENU_FB=1; // set ADD_MENU 1 to add menu button of some functions // // if ADD_MENU_FB == 1 then delmenu('step7_filter_bank'); addmenu('step7_filter_bank',[ '(1) set fb parameters(30)','(2) do filter bank','(-) plot3d_fb_am_out' ]); step7_filter_bank=[ '[st_freq, end_freq, n_stepf, fir_steps, kwindows, end_freq2]= set_fb_para(); [bpf_fb0,wft]=make_bpf_fb0(30)','[fb_out,fb_am_out,fb_am_out_prop]=do_fb(1); messfb0='' ''; plot_fb(32,messfb0)'] ; end //if ADD_MENU_FB == 1 then //----------------------------------------------------------------- // //==============================================================++= // // SCILAB's MAXIMUM STRING LENGTH is 512 ? // //////addmenu('&go'); // only windows shortcut can alt+()