濁音の「だ」の波形をよく観察してみると、波形の初期の部分に比較的高
い周波数の波形がより大きく含まれていることがわかる。そこで、２管模型の出力にピーク イコライザー（PEAK
EQUALIZER)をつかって一部の周波数を持ち上げ強調することによって、「だ」が実現できないかどうかを実験してみることにした。
スペクトル表示でみると分かりにくいが、下図のピンク色の部分がもりあがっている。
「だ」の波形のすべてをフィッティングするのは大変なので、代表的と思われる３つの部分を抜き取って最小２乗法でフィッティングしてみることにした。
最小２乗法をおこなうとき、評価の重み付けを使用することによって 期待した形でななくオリジナルとの差が小さくなる予想外の解が出てくることを避けるた
め、オリジナルの波形の周波数特性が最大値をとる周波数においての２管模型の理論周波数特性の値がオリジナルと同じになるようなオフセットを加える方法を
導入した。
「だ」の３つの部分のフィッティングの結果を下図に示す。一番右の３番目は「あ」そのものである。
最小２乗法をつかって推定？（最適化）したパラメーターをつかって作った音を参考にリ
ンクしておこう。濁音感が強めの音になっている。（ちなみに、手動で与えた初期値を使った作った音は
こんな感じ。）
代表と思われる３つの部分のフレームだけしかパラメータは計算していないので、代表フレーム以外の間のフレームのパラメータは代表フレームからの値から直
線補間で計算した。全部で７フレーム、音を２管模型で求めている。はじまりのNo.2とNo.3のフレームはピーク イコライザーをつかって強調してあ
る。
参考に、上記を計算した SCILABのデモ プログラム
a_wavfile_edit_647.sci を載
せておきます。 入力には 音声サンプルのda_sample.wav（「だ」）を使います。
このプログラムを実
際に動作させるためには、
ご使用されるプログラミング環境に合わせて修正・変更・追加などが必要かもしれません。
また、バグが含まれている可能性があります。
万一、このデモンストレーション用プログラムを動作させる場合は、あなたの責任でおこなってくださいね。
//-----------------------------------------------------------------------------------------
// 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
//
//
// 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
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 ind_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
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 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 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_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()
// choose wave file
if T_DEMO == 1 then
//+ select
//SEL_CODE=x_choose0(['da_sample (into 3 parts)';'-'; '-';
'manual select'],['Please select one'],'default')
SEL_CODE=x_choose0(['da_sample (into 3 parts)';'-'; '-';
'-'],['Please select one'],'default')
if SEL_CODE == 4 then // when manual select
SEL_CODE = 999;
elseif SEL_CODE <= 0 then
SEL_CODE=1;
end
//- select
if SEL_CODE == 1 then // da_sample part 1
disp(' This is demonstration. Copy this file and
da_sample.wav to your scilab''s bin or home
directory.');
[x,ierr]=fileinfo('da_sample.wav');
xs=size(x);
if xs(1)==0 then
disp('Choose da_sample.wav file');
if scilab_version_number >= 5 then
wavfile1=uigetfile(["*.wav"],"","Choose da_sample.wav file");
else
wavfile1=tk_getfile(Title="Choose
da_sample.wav file");
end
else
wavfile1='da_sample.wav';
end
elseif SEL_CODE == 2 then // da_sample part 2
disp(' This is demonstration. Copy this file and
da_sample.wav to your scilab''s bin or home
directory.');
[x,ierr]=fileinfo('da_sample.wav');
xs=size(x);
if xs(1)==0 then
disp('Choose da_sample.wav file');
if scilab_version_number >= 5 then
wavfile1=uigetfile(["*.wav"],"","Choose da_sample.wav file");
else
wavfile1=tk_getfile(Title="Choose
da_sample.wav file");
end
else
wavfile1='da_sample.wav';
end
elseif SEL_CODE == 3 then // da_sample part3
disp(' This is demonstration. Copy this file and
da_sample.wav to your scilab''s bin or home
directory.');
[x,ierr]=fileinfo('da_sample.wav');
xs=size(x);
if xs(1)==0 then
disp('Choose da_sample.wav file');
if scilab_version_number >= 5 then
wavfile1=uigetfile(["*.wav"],"","Choose da_sample.wav file");
else
wavfile1=tk_getfile(Title="Choose
da_sample.wav file");
end
else
wavfile1='da_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 // da_sample
QTY_FRAME=3;
sp1s(1)=1435;
ep1s(1)=2000;
sp1s(2)=3129;
ep1s(2)=3700;
sp1s(3)=7810;
ep1s(3)=8400;
//
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)
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(0); // 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);
//////////////////////////////////////////////////////
smooths(2);
////[sn,swnd,sm1out,sm2out,npk,pklist]=smooth1(5);
////plot_fft_sm1(5);
/////////////////////////////////////////////////////
end
//
// -MAIN (2)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.
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
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
//
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.
iir_bpf1_res1=zeros(1,1); // this will be overwritten.
iir_bpf2_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']);
//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;
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
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));
//
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
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) ;
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<cfreq(1),cfreq(2)<.5 sampling frequency = 1
//----------------------------------------------------------
function [xy_ran2]=make_noise(xlength0)
cfreq(1) = (frq_center - ( freq_band / 2)) / fs1;
cfreq(2) = (frq_center + ( freq_band / 2)) / fs1;
fpar(1)=1; // dummy data
fpar(2)=0; // dummy data
[wft,wfm,fr]=wfir( 'bp' , nstep_fir ,cfreq, 'hm' ,fpar);
rand('normal'); // random generator is set to a Gaussian
//xlength0=512; // dummy set
xlength1= xlength0 + nstep_fir;
xy_ran1=rand(1, xlength1);
xy_ran2=zeros(1, xlength0);
disp('---please wait for some time. noise calculating');
for loop=1:xlength0
ydz = 0.;
for loop2=1:nstep_fir
ydz = ydz + wft(loop2) * xy_ran1( nstep_fir - loop2 + loop);
end
xy_ran2(loop)=ydz;
end
endfunction
//----------------------------------------------------
function [xy_ran2]=make_noise_bpf(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
xlength1= xlength0 + 10;
xy_ran1=rand(1, xlength1);
xy_ran2=zeros(1, xlength0);
yx=zeros(2,6);
disp('---please wait for some time. noise calculating');
for loop=1:xlength0
for vl=1:2
yx(vl,4)=xy_ran1(loop)
x0=yd_bpf(vl,1) * yx(vl,4) + yd_bpf(vl,2) * yx(vl,5) +
yd_bpf(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_bpf(vl,2) + yx(vl,3) *
xd_bpf(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];
end
xy_ran2(loop)=yx(1,1) + i2nd_thd_factor * yx(2,1);
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 (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); // tthis 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); // tthis 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); // tthis 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 ind_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);
[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 // da_sample
WT_QTYs(1)=41;
r1s(1)=0.1;
l1s(1)=-0.3;
ttl_Lengths(1)=18.5;
rls(1)=0.7;
limit_switch3s(1)=1.0;
emp_space_areas(1)=5.2;
empA0s(1)=20.;
//---
WT_QTYs(2)=2;
r1s(2)=0.4;
l1s(2)=0.;
ttl_Lengths(2)=18.5;
rls(2)=0.9;
limit_switch3s(2)=1.0;
WT_QTYs(3)=2;
r1s(3)=0.8;
l1s(3)=0.;
ttl_Lengths(3)=18.5;
rls(3)=0.9;
limit_switch3s(3)=0.0;
//
WT_QTY=WT_QTYs(tframe);
r1=r1s(tframe);
l1=l1s(tframe);
ttl_Length=ttl_Lengths(tframe);
rl=rls(tframe);
limit_switch3=limit_switch3s(1);
emp_space_area=emp_space_areas(1);
empA0=empA0s(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.
// 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
//----------------------------------------------------------
//++++
// 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()
// 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,4); // 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
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 '];
//----------------------------------------------------------------
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);
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
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);
// 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
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
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
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
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);
// 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
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;
//
//=====================================================
//
// +MAIN (L)program starts
//
if T_DEMO==1 then
if SEL_CODE == 1 then // a_sample
// part 1
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
// part 2
fact1_2r1s(2)=0.9;
fact2_2r1s(2)=-0.9;
fact1_2l1s(2)=0.9;
fact2_2l1s(2)=-0.9;
fact1_rls(2)=0.95;
fact2_rls(2)=0.5;
limit_switch0s(2)=0.0; // r1,r2,l1,l2 off
limit_switch2s(2)=1.0; // rl on
// 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.95;
fact2_rls(3)=0.5;
limit_switch0s(3)=0.0; // r1,r2,l1,l2 off
limit_switch2s(3)=1.0; // rl on
//
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);
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
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
//-----------------------------------------------------------------
// 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 [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);
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;
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
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);
endfunction
//
//=================================================================
//
// +MAIN (4)program starts
// leastsq method to estimate tube model parameters
if T_DEMO==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 // da_sample
[wm8s]= set_wm8( 344.53125 , 861.32812 , 1.0 ,1);
[wm8s]= set_wm8( 1205.8594 , 2842.3828 , 1.0 ,1);
[wm8s]= set_wm8( 430.66406 , 1000. , 1.0 ,2);
[wm8s]= set_wm8( 1200. , 1800. , 1.0 ,2);
[wm8s]= set_wm8( 2100. , 2670.1172 , 1.0 ,2);
[wm8s]= set_wm8( 430.66406 , 1722.6562 , 1.0 ,3);
[wm8s]= set_wm8( 2000. , 3500. , 1.0 ,3);
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 // da_sample
// set each frame count length
// part1 * part2 * * part3 is
// frame_lengths(1)=2;
// frame_lengths(2)=3;
frame_lengths(1)=2;
frame_lengths(2)=3;
//QTY_FRAME=3
amplitudes(1)=1.0;
amplitudes(2)=1.0;
amplitudes(3)=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();
plot_waveform(5);
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' ; 'load_wfft()'; 'save_wfft()']);
step2_fft=[ '[fft_point1]=set_fft_points()' ; '[db_fft1,
phi_fft1]=do_fft_wav(); plot_fft1(1)';
'[sn,swnd,sm1out,sm2out,npk,pklist]=smooth1(1); plot_fft_sm1(2)'
;'[fft_point1,db_fft1,phi_fft1]=load_wfft(); plot_fft1(1)' ;
'save_wfft()'] ;
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' ]);
step4_2_leastsq=[' [fopt,xopt, grdopt]=leastsq_main1s();' ;
' plot_2tubeS((9+tframe),1)' ; 'exchange8s(); mess_init0='' by leastsq
method result'';' ];
end
//
//-----------------------------------------------------------------
//
ADD_MENU_GEN=1; // set ADD_MENU 1
to add menu button of some functions
//
//
if ADD_MENU_FFT == 1 then
delmenu('step5_generation');
if (f_win == 1) | (scilab_version_number >= 4) then
addmenu('step5_generation',[ '(1)set_section_qty()' ;
'(2)make_waveform()' ; 'plot_waveform()'; 'y3out_snd_play1()';
'y3out_snd_save1()' ]);
step5_generation=[ '[N_REPEAT]= set_section_qty()';
'[yg,y2tm,y3out,y2tm_lpf,y2tm_noise,LL]= make_waveform();' ;
'plot_waveform(5)'; 'y3out_snd_play1()'; 'y3out_snd_save1()'] ;
else
addmenu('step5_generation',[ '-' ; '(2)make_waveform()' ;
'plot_waveform()'; 'y3out_snd_save1()' ]);
step5_generation=[ ' '; '[yg,y2tm,y3out,y2tm_lpf,LL],y2tm_noise=
make_waveform();' ; '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
//
//==============================================================++=
//
// SCILAB's MAXIMUM STRING LENGTH is 512 ?
//
//////addmenu('&go'); // only windows shortcut
can alt+()
No.3c 2009年8月15日