ホームページ
<フィルターバンク(FILTER BAND)による波形の分解と合成>
- フィルター バンク(FILTER BAND)による波形の分解と合成のためのSCILAB用のサンプルプログラムa_wavefile_edit_70.sci
- 波形に2箇所の区間を設定して、
FFTによる周波数特性とフィルターバンクの出力を計算して、2箇所の区間を比較するためのSCILAB用のサンプルプログラムa_wavefile_edit_80b.sci
- SCILABのevent handlerイベントハンドラーを設定することで、キー入力（キーd(down)とキーu(up))を使って
フィ
ルターバンクの出力を３バンド単位で移動表示し、キー入力（キーf(表示のon/off))で周波数の変化を折れ線で表示できるようにしたSCILAB
用のサンプルデモプログラムa_wavefile_edit_84c.sci
フィルター
バンク(FILTER BAND)による波形の分解と合成のためのSCILAB用のサンプルプログラム
a_wavefile_edit_70.sci
参考に、ここで使ったSCILABのプログラム a_wavefile_edit_70.sci を載せておきます。
このプログラムで使う音声データの wavファイルのサンプル
を、SCILABの bin または
home ディレクトリーにコピーしておいてから、SCILAB上でこのプログラムをexecしてみてください。ＳＣＩＬＡＢについては、SCILABのホームページを見てください。
T_DEMOを１に設定すると半自動動作のデモ モードで起動します。また、
T_DEMOを０に設定するとはじめからマニュアルで操作するように起動します。はじめに、分析したい wav
ファイルのありかを指定します。 scilabの画面の上面にある選択メニュー
に新たにメニューが追加されますので、それらを使って 操作することができます。
ここでは直線位相のFIR型のバンド パス
フィルターを使って フィルターバンクを構成しています。直線位相のFIR型のフィルターの遅延時間に合わせて、バンド パス フィルターを計算した
出力の波形を遅延時間分ずらして波形を比
較するためです。メニューの中の set fb
parameters を使って、フィルターバンクの周波数範囲と、分割数と、FIR型のフィルターの次数を設定することもできます。下図は、バンド
パス
フィルターの周波数特性の例で、フィルターバンクの周波数範囲を８個に分割したバン ドパ スフィルターで構成していることを示しています。
その音素にとってファイルーバンクで分割された成分の中のどれが主要因であるかを聞いて識別することを目的に、メニューの中にある
replacementを使って、オリジナルの波形の部分を、選択したバンド パス フィルターの出力の波形の和に置き換えたものを作ることができます。
それを
wavファイルに書き出し(save)、聞くことができます。OS側のサウンドの設定が上手くいっていれば、SCILAB上から直接 合成した音を
play 再生することもできるかもしれません。
上図の例では、フィルターバンクの５番目と６番目と７番目と８番目の波形を足し合わせて、オリジナルの波形のその部分を置き換えることを示しています。
このサンプルを実際に動作させるためには、
ご使用されるプログラミング環境に合わせて修正・変更・追加などが必要かもしれません。
また、バグが含まれている可能性があります。
万一、このサンプルを動作させる場合は、あなたの責任でおこなってくださいね。
//------------------------------------------------------------------------------------------
//
// Linear-Phase FIR filter bank, AM detection, and 2D
display
//
// for window scilab-4.1.2
//
//
//
.........................................................................................
// PLEASE PAY ATTENTION that this program may have some bugs and also
you may modify program
// to fit your circumstances. If you use this program, everything done
by your own risk.
// Thanks.
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
T_DEMO=1; // set T_DEMO=1 for demonstration, beside set
T_DEMO=0 when normal
T_N_STEPF=0; // set T_N_STEPF=1 for manual set in DEMO,
beside set T_N_STEPF=0 when normal
//
// select sample speech and its portion
if T_DEMO==1 then
SEL_CODE=x_choose(['ka_sample';'ta_sample ';'sa_sample'
;'ka2_fa'; 'ta2_fa'; 'sa2_fa'],['Please select one'],'default')
if SEL_CODE <= 0 then
SEL_CODE=1;
end
end
//===================================================================
// 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.
wdat2=zeros(1,10);// data of synthesized .wav, 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.
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
// for others
f412=0; // flag to detect scilab 4.1.2
defch1=1; // default channel 1
f_win=1; // flag if windows
//---------------------------------------
//
//
//--------------------------------------
function [wavfile1,chs1,qty1,fs1,bit1,f412]=get_wavfile_pro()
// choose wave file
if T_DEMO == 1 then
if SEL_CODE==3 then
[x,ierr]=fileinfo('sa_sample.wav');
xs=size(x);
if xs(1)==0 then
disp('Choose sa_sample.wav file');
wavfile1=tk_getfile(Title="Choose
sa_sample.wav file");
else
wavfile1='sa_sample.wav';
end
elseif SEL_CODE==2 then
[x,ierr]=fileinfo('ta_sample.wav');
xs=size(x);
if xs(1)==0 then
disp('Choose ta_sample.wav file');
wavfile1=tk_getfile(Title="Choose
ta_sample.wav file");
else
wavfile1='ta_sample.wav';
end
elseif SEL_CODE==4 then
[x,ierr]=fileinfo('ka2_fa.wav');
xs=size(x);
if xs(1)==0 then
disp('Choose ka2_fa.wav file');
wavfile1=tk_getfile(Title="Choose
ka2_fa.wav file");
else
wavfile1='ka2_fa.wav';
end
elseif SEL_CODE==5 then
[x,ierr]=fileinfo('ta2_fa.wav');
xs=size(x);
if xs(1)==0 then
disp('Choose ta2_fa.wav file');
wavfile1=tk_getfile(Title="Choose
ta2_fa.wav file");
else
wavfile1='ta2_fa.wav';
end
elseif SEL_CODE==6 then
[x,ierr]=fileinfo('sa2_fa.wav');
xs=size(x);
if xs(1)==0 then
disp('Choose sa2_fa.wav file');
wavfile1=tk_getfile(Title="Choose
sa2_fa.wav file");
else
wavfile1='sa2_fa.wav';
end
else
[x,ierr]=fileinfo('ka_sample.wav');
xs=size(x);
if xs(1)==0 then
disp('Choose ka_sample.wav file');
wavfile1=tk_getfile(Title="Choose
ka_sample.wav file");
else
wavfile1='ka_sample.wav';
end
end
else
wavfile1=tk_getfile(Title="Choose xxx.wav file name");
end
// read channels and samples
cs1=wavread(wavfile1,'size');
chs1=cs1(2); // channel
qty1=cs1(1); // samples
//...
if chs1 > 2 then // invert data for scilab-4.1.2
chs1=cs1(1);
qty1=cs1(2);
f412=1;
else
f412=0;
end
//...
[y,fs1,bit1]=wavread(wavfile1,[1 1]);
endfunction
//--------------------------------------
// function 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()
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');
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();
endfunction
//----------------------------------------------------------------
function [wsp1, wep1]= reset_sp1_ep1(wsize1)
wsp1=1;
wep1=wsize1;
endfunction
//----------------------------------------------------------------
function [wsp1, wep1]= set_sp1_ep1(wsize1)
txt1=['start point';'end point'];
wstr1=sprintf('%d',sp1);
wstr2=sprintf('%d',ep1);
sig1=x_mdialog('Input start point and end point for portion
display.',txt1, [wstr1 ; wstr2 ]);
if sig1==[] then
arg1=evstr(wstr1);
arg2=evstr(wstr2);
else
arg1=evstr(sig1(1));
arg2=evstr(sig1(2));
end
//
//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;
disp(wep1,'ep1',wsp1,'sp1');
end
end
end
endfunction
//
//----------------------------------------------------------------
function snd_play2()
// this sound doesnot work well on windows scilab-3.1.1 and linux
scilab-3.1
// but, this sound works on windows scilab-4.1.2. But, linux
scilab-4.1.2 doesnot!
// for wdat2
sz=size(wdat2);
size2=sz(2);
if f_win == 1 then
if f412 == 1 then // for windows scilab-4.1.2
sound(wdat2 ,fs1,bit1);
else
// for windows scilab-3.1.1
if fs1 == 22050 then
sound(wdat2',fs1,bit1);
disp('This function may work, if luckily.');
elseif fs1 == 44100 then // down-sampling from
44100 to 22050
wdatw2=zeros(1,(int(size2/2)));
for v=1:(int(size2/2))
wdatw2(v)=wdat2(2 * v );
end
disp('down-sampling from 44100 to 22050');
sound(wdatw2',22050,bit1);
disp('This function may work, if luckily.');
else
//sound(wdatw,fs1,bit1);
disp('This function does not work.');
end
end
else
disp('If sound setting is good for scilab, this function maybe
work.');
if f412 == 1 then // for scilab-4.1.2
sound(wdat2 ,fs1,bit1);
else
// for scilab-3.1.1
if fs1 == 22050 then
sound(wdat2',fs1,bit1);
disp('This function may work, if luckily.');
elseif fs1 == 44100 then // down-sampling from
44100 to 22050
wdatw2=zeros(1,(int(size2/2)));
for v=1:(int(size2/2))
wdatw2(v)=wdat2(2 * v );
end
disp('down-sampling from 44100 to 22050');
sound(wdatw2',22050,bit1);
disp('This function may work, if luckily.');
else
//sound(wdatw,fs1,bit1);
disp('This function does not work.');
end
end
end
endfunction
//
//----------------------------------------------------------------
function snd_save2()
// for wdat2
//
wavefilename= input(' + enter file name for saved .wav file
=>',["string"]);
//
if f412 == 1 then // for windows scilab-4.1.2
wavwrite(wdat2 ,fs1,bit1,wavefilename );
else
// for windows scilab-3.1.1
wavwrite(wdat2',fs1,bit1, wavefilename);
end
endfunction
//----------------------------------------------------
function [frq1,sfrq1,is1,ie1]=set_frq( spec1024 )
// spec1024, tube model response no syuhasuu bunkainou wo ageru
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
//-----------------------------------------------
//--- 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]=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);
if T_DEMO==1 then
if SEL_CODE==3 then // sa_sample
sp1=6000;
ep1=9500;
elseif SEL_CODE==2 then // ta_sample
sp1=270;
ep1=1440;
elseif SEL_CODE==4 then // ka_2_fa
sp1=270;
ep1=2600;
elseif SEL_CODE==5 then // ta_2_fa
sp1=270;
ep1=1050;
elseif SEL_CODE==6 then // sa_2_fa
sp1=300;
ep1=3200;
elseif SEL_CODE==10 then //
ta2_from_ka2_1750hzhpf
sp1=300;
ep1=1600;
else
// ka_sample
sp1=260;
ep1=1900;
end
end
[ydat1,xziku1]=make_width_data( wdat1, size1);
//
// -MAIN (1)program starts
//
//==============================================================++=
//
//
//===================================================================
// global variables No. 2
// for filter bank
st_freq=500.0;
end_freq=st_freq * 8;
n_stepf=8;
bpf_fb0=zeros(512,3); // *,1 center frequency.
this will be overwritten.
// *,2 lower frequency of bpf band. this will be overwritten.
// *,3 higher frequency of bpf band. this will be overwritten.
fir_steps=513; // This must be odd. this will
be overwritten.
kwindows=1;
wft=zeros(n_stepf,fir_steps); // this will be overwritten.
wfm=zeros(n_stepf,fir_steps); // this will be overwritten.
fr=zeros(fir_steps); // this will be overwritten.
fb_out=zeros(n_stepf,1); // this will be overwritten.
fb_am_out=zeros(n_stepf,1); // this will be overwritten.
nmax_am_prop=5;
fb_am_out_prop=zeros(n_stepf,nmax_am_prop); // *,1
average value about fb_am_out
// *,2 peak value about fb_am_out
cfb_list=ones(1,2); // this will be overwritten.
//-------------------------------------------------
function [st_freq, end_freq, n_stepf, fir_steps, kwindows]=
set_fb_para()
// codex defines which term....
txt1=['start freq'; 'end freq'; 'freq steps'; 'steps of fir
filter (odd)'; 'priority: frq-resolution (1), drop value (2)'];
wstr1=sprintf('%f',st_freq);
wstr2=sprintf('%f',end_freq);
wstr3=sprintf('%d',n_stepf);
wstr4=sprintf('%d',fir_steps);
wstr5=sprintf('%d',kwindows);
sig1=x_mdialog('set fir filter bank parameters ',txt1, [wstr1;
wstr2; wstr3; wstr4; wstr5]);
if sig1==[] then
arg1=evstr(wstr1);
arg2=evstr(wstr2);
arg3=evstr(wstr3);
arg4=evstr(wstr4);
arg5=evstr(wstr5);
else
arg1=evstr(sig1(1));
arg2=evstr(sig1(2));
arg3=evstr(sig1(3));
arg4=evstr(sig1(4));
arg5=evstr(sig1(5));
end
st_freq=arg1;
end_freq=arg2;
n_stepf=arg3;
fir_steps=arg4;
kwindows=arg5;
endfunction
//-----------------------------------------------------------
function [bpf_fb0,wft]=make_bpf_fb0(disp0)
// log gauge divide
a1=log(st_freq);
b1=log(end_freq);
c1=(b1-a1)/n_stepf;
bpf_fb0=zeros(n_stepf,3);
for v=1:n_stepf
a2=exp(a1+ c1 * (v-1));
b2=exp(a1+ c1 * v);
c2=(b2-a2)/2. + a2;
bpf_fb0(v,1)=c2;
bpf_fb0(v,2)=a2;
bpf_fb0(v,3)=b2;
if disp0 >= 1 then
disp( [ v bpf_fb0(v,1) bpf_fb0(v,2)
bpf_fb0(v,3)]);
end
end // for V=1:n_stepf
//
// ++ mk_fir_coeff
//
w0=xget('window'); // stack old window
wft=zeros(n_stepf,fir_steps);
for vloop=1:n_stepf
cfreq(1)= bpf_fb0(vloop,2) / fs1;
cfreq(2)= bpf_fb0(vloop,3) / fs1;
fpar(1)=0.1; // dummy data
fpar(2)=-0.1; // dummy data
if kwindows == 2 then
[wft2,wfm2,fr2]=wfir('bp',fir_steps,cfreq,'hn',fpar); // gensui
riritu yuusen
else
[wft2,wfm2,fr2]=wfir('bp',fir_steps,cfreq,'hm',fpar); // bunkai
nou yuusen
end
for v2=1:fir_steps
wft(vloop,v2)=wft2(v2);
end
if disp0 >=2 then
//++ only one will be done
if vloop==1 then
xset('window',disp0); // create new windows
clf();
[frq1,sfrq1,is1,ie1]=set_frq(0);
fr0=zeros(2); // freq list selected
frq1f=zeros(2); // freq responce selected
phf0=zeros(2); // phase selected
sv0=size(fr2);
end // if vloop==1 then
//-- only one will be done
sc0=0;
for v=1:sv0(2)
fx0= fr2(v) * fs1;
if (fx0 >= frq1(1)) &(fx0 <=
frq1(sfrq1)) then
sc0=sc0+1;
fr0(sc0)=fx0;
frq1f(sc0)= 20. * log10( wfm2(v));
phf0(sc0)= -1.0 * ((fir_steps - 1)
/ 2.0 ) * ( 2. * %pi * v / fir_steps);
end
end // for v=1:sv0(2)
gainplot(fr0,frq1f,phf0);
end // vloop+1:n_stepf
end // if disp0 >=2 then
xset('window',w0); // push old windows
//
// -- mk_fir_coeff
//
endfunction
//-----------------------------------------------------------
function plot_fb(disp_code)
w0=xget('window'); // stack old window
if disp_code >= 1 then
xset('window',disp_code); // create new windows
clf();
wsize2=size1 - sp1 + 1;
wsize3=ep1-sp1+1;
zdata1=zeros(1,width0+1);
zdata1p=zeros(1,width0+1);
ydata1=zeros(1,width0+1);
for vloop=1:n_stepf
if wsize3 > (width0+1) then
wstep1= (ep1-sp1) / width0;
for v=1:(width0+1)
zdat1(v)= fb_out(vloop,1 + int(wstep1 * (v -
1)));
zdat1p(v)= fb_am_out(vloop, 1 + int(wstep1 *
(v - 1)));
xziku1(v)= sp1 + int(wstep1 * (v - 1));
if vloop==1 then
ydat1(v)= wdat1( sp1 + int(wstep1 * (v -
1)));
end
end
else
for v=1:(width0+1)
if v <= wsize3 then
zdat1(v)=fb_out(vloop,v);
zdat1p(v)=fb_am_out(vloop,v);
if vloop==1 then
ydat1(v)= wdat1((sp1-1)+v);
end
else
zdat1(v)=0.;
zdat1p(v)=0.;
if vloop==1 then
ydat1(v)= 0.;
end
end
xziku1(v)=(sp1-1)+v;
end
end
if vloop==1 then
subplot((n_stepf+1),1,1);
plot(xziku1,ydat1,'b');
wstr1=sprintf('waveform selected, original');
xtitle(wstr1);
end
subplot((n_stepf+1),1,(vloop+1));
plot(xziku1,zdat1,'b');
plot(xziku1,zdat1p,'g');
wstr1=sprintf('Band Pass filtered waveform selected (%dHz-
%dHz)',int(bpf_fb0(vloop,2)),int(bpf_fb0(vloop,3)));
mv0=fb_am_out_prop(vloop,1);
if mv0 > 0. then
mv1= int(20. * log10(mv0));
else
mv1=-120;
end
pk0=fb_am_out_prop(vloop,2);
if pk0 > 0. then
pk1= int( 20. * log10(pk0));
else
pk1=-120;
end
wstr2=sprintf('%d/%d dB',mv1,pk1); // average / peak
[dB]
xtitle( wstr1,'',wstr2 );
end // for vloop=1:n_stepf
end // if disp_code >= 1 then
xset('window',w0); // push old windows
endfunction
//-----------------------------------------------------------
function windowsize0(width0,height0)
dim0=xget("wdim");
disp(dim0);
xset("wdim",width0,height0); // real width and real height
is limited in physical display pixes
endfunction
function showwindowsize0()
dim0=xget("wdim");
disp(dim0);
endfunction
//-----------------------------------------------------------
function plot3d_fb_am_out(disp0)
//
wb0=xget('window'); // stack old window
if disp0 >= 1 then
xset('window',disp0); // create new windows
clf();
w0=ep1-sp1+1;
x1=zeros(n_stepf,1);
bai0=10;
for v=1:n_stepf
x1(v)=bpf_fb0(v,1) / bai0; //
value shoud be near w0
end
x1s=bpf_fb0(1,2) / bai0; // lower freq
x1e=bpf_fb0(n_stepf,3) / bai0; // higher freq
y1=zeros(w0,1);
for v=1:w0
y1(v)=sp1+(v-1);
end
pk0=fb_am_out_prop(1,2);
for v=1:n_stepf
if fb_am_out_prop(v,2) > pk0 then
pk0=fb_am_out_prop(v,2);
end
end
if pk0 <= 0. then
gain0=1000.;
else
gain0=1000./pk0;
end
fb_am_out2=zeros(n_stepf,w0);
for v=1:n_stepf
for v2=1:w0
fb_am_out2(v,v2)=fb_am_out(v,v2) * gain0;
end
end
clf();
plot3d(x1,y1, fb_am_out2 ,theta=130, alpha=75,
leg="x10Hz@point@1000",flag=[5,1,4],ebox=[x1s,x1e,y1(1),y1(w0),0.,1000.0]);
xtitle('filter bank am detect output');
end // if disp0 >= 1 then
xset('window',wb0); // push old windows
endfunction
//-----------------------------------------------------------
function [fb_out,fb_am_out,fb_am_out_prop]=do_fb(disp0)
//
w0=(ep1-sp1)+1;
fb_out=zeros(n_stepf,w0);
dt0= int((fir_steps-1) / 2);
ep1x=ep1+dt0;
if ep1x > size1 then
ep1x = size1;
end
sp1x=sp1+dt0;
if sp1x < fir_steps then
sp1x = fir_steps;
end
//+filtering+
for vloop=1:n_stepf
if disp0 >= 1 then
wstr1=sprintf('+entering filter bank
no.%d/%d',vloop,n_stepf);
disp(wstr1);
end
for loop=sp1x:ep1x
ydz = 0.;
for loop2=1:fir_steps
ydz = ydz + wft(vloop,loop2) *
wdat1( 1 - loop2 + loop );
end
//fir_out1(loop - dt0)=ydz;
fb_out(vloop,loop - dt0 -sp1 + 1)=ydz;
end
end // if vloop=1:n_stepf
//-filtering-
//
//+am+
nsize_max=5000;
g_work2=zeros(10,nsize_max);
fb_am_out=zeros(n_stepf,w0);
fb_am_out_prop=zeros(n_stepf,nmax_am_prop);
fir_ou1=zeros(1,size1);
fir_out1p=zeros(1,size1);
for vloop=1:n_stepf
if disp0 >= 1 then
wstr1=sprintf('+entering AM de-modulation, peak
trace no.%d/%d',vloop,n_stepf);
disp(wstr1);
end
for v2=sp1:ep1
fir_out1(v2)= fb_out(vloop,v2-sp1+1);
fir_out1p(v2)=0.;
end
//+from trace peak
//+
spx=-1;
spy=-1;
epx=-1;
epy=-1;
c0=-1;
w0=-1.;
area0=-1.;
if (ep1-sp1) > 4 then
for v=(sp1+1):(ep1-1)
if fir_out1(v) > 0. then // only
postive side
if (fir_out1(v) > fir_out1(v-1))
& (fir_out1(v) > fir_out1(v+1)) then
spx=epx;
spy=epy;
epx=v;
epy=fir_out1(v);
c0=c0+1;
//
if g_work2(1,1) < nsize_max then
if g_work2(1,1) <= 0 then
g_work2(1,1)=1; // first start
g_work2(3,1)=epx;
g_work2(4,1)=epy;
g_work2(5,1)=epx;
g_work2(6,1)=epy;
else
d0=int(g_work2(1,1));
if epy >= g_work2(6,d0) then // rising or equal
g_work2(2,d0)=0; //current
g_work2(5,d0)=epx;
g_work2(6,d0)=epy;
else
// falling
if g_work2(2,d0)==0 then // complete data
g_work2(2,d0)=1; //complete
g_work2(7,d0)= g_work2(5,d0) - g_work2(3,d0);
g_work2(8,d0)= g_work2(6,d0) - g_work2(4,d0);
g_work2(9,d0)= g_work2(7,d0) * g_work2(8,d0) / 2.0 ;
if (g_work2(7,d0) > 0.) & (g_work2(8,d0) > 0.) then
if g_work2(9,d0) > area0 then
area0=g_work2(9,d0);
w0=g_work2(8,d0) / g_work2(7,d0);
end
end
// atan
//g_work2(10,d0)= atan( (g_work2(8,d0) / g_work2(7,d0))) / 2.0 / %pi *
360.0;
//
g_work2(1,1)=g_work2(1,1)+1; // set as new start point
d1=int(g_work2(1,1));
g_work2(3,d1)=epx;
g_work2(4,d1)=epy;
g_work2(5,d1)=epx;
g_work2(6,d1)=epy;
elseif g_work2(2,d0) <= -1 then // restart
g_work2(3,d0)=epx;
g_work2(4,d0)=epy;
g_work2(5,d0)=epx;
g_work2(6,d0)=epy;
end
end
end // if g_work2(1,1) <= 0 then
else
disp(' error: g_g_work2(1,1) >= nsize_max , in do_fb() ');
end // if g_work2(1,1) < nsize_max then
//
//
if c0 >= 1 then
a0=
(epy - spy) / (epx - spx);
b0=
epy - a0 * epx;
for
v2=spx:epx
fir_out1p(v2)= a0 * v2 + b0;
end
end // if
c0 >= 1 then
//
end // if (fir_out1(v) >
fir_out1(v-1)) & (fir_out1(v) > fir_out1(v+ 1)) then
end // if fir_out1(v) > 0. then
end //v=(sp1+1):(ep1-1)
end
if c0 < 1 then
epx=sp1-1;
end
for v=(epx+1):ep1
fir_out1p(v)=0;
end
pk0=0.;
mv0=0.;
for loop=sp1:ep1
fb_am_out(vloop,loop - sp1 + 1)= fir_out1p(loop);
if fir_out1p(loop) > pk0 then
pk0=fir_out1p(loop);
end
mv0=mv0+fir_out1p(loop);
end
fb_am_out_prop(vloop,1)= mv0 / (ep1-sp1+1);
fb_am_out_prop(vloop,2)= pk0;
//+
//-from trace peak
end // for vloop=1:n_stepf
//-am-
endfunction
//-----------------------------------------------------------------
function [wdat2,rslt]=do_replace(disp0)
sz0=size(cfb_list);
cfb_list2=ones(n_stepf,2);
for v=1:n_stepf
cfb_list2(v,1)=v;
if v <= sz0(1) then
cfb_list2(v,2)=cfb_list(v,2);
end
end
[rslt]=x_matrix('Mark 1 as Combine element, in filter bank no.
list(no, 1 or 0 )' , cfb_list2);
if rslt==[] then // reset all 1
rslt=ones(n_stepf,2);
for v=1:n_stepf
rslt(v,1)=v;
end
end
wdat2=zeros(1,size1);
for v=1:size1
wdat2(v)=wdat1(v);
end
//+replacement
for v=sp1:ep1
vs0=0.;
for v2=1:n_stepf
if rslt(v2,2)==1 then // if the fb
no is Combine element
vs0= vs0 + fb_out(v2,
(v-sp1+1))
end
end
wdat2(v)=vs0;
end
//-replacement
if disp0 >= 1 then
ydat1r=zeros(1,(width0+1));
if size1 > (width0+1) then
wstep1= size1 / (width0+1);
for v=1:(width0+1)
ydat1r(v)= wdat2( 1 + int(wstep1 * (v -
1)));
// xziku1(v)= 1 + int(wstep1 * (v - 1));
end
else
for v=1:(width0+1)
if v <= size1 then
ydat1r(v)=wdat2(v);
else
ydat1r(v)=0.;
end
// xziku1(v)=(sp1-1)+v;
end
end
//
clf();
subplot(211);
plot(xziku0,ydat0,'b');
//plot2d(xziku0,ydat0,style=[color("blue")]);
xtitle('waveform all');
// When linux scilab-3.1 2nd subplot does not work well.
subplot(212);
plot(xziku0,ydat1r,'b');
//plot2d(xziku1,ydat1,style=[color("blue")]);
xtitle('replacement waveform');
end // id disp0 >= 1 thne
endfunction
//=================================================================
//
// +MAIN (2) program starts
// filter bank
//
if T_DEMO==1 then
if T_N_STEPF==1 then
n_stepf=input("-> How many filters of filter bank ?");
end
[bpf_fb0,wft]=make_bpf_fb0(3);
[fb_out,fb_am_out,fb_am_out_prop]=do_fb(1);
plot_fb(2);
end
//
// -MAIN (2)program starts
//
//==============================================================++=+
//
//
// ++ MENU MENU MENU =====================
//
// 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',[ '(1)set_et_plot()'; 'reset_plot()']);
step1_plotwav=[ '[sp1, ep1]= set_sp1_ep1(size1);
[ydat1,xziku1]=make_width_data( wdat1, size1);' ; '[sp1, ep1]=
reset_sp1_ep1(size1); [ydat1,xziku1]=make_width_data( wdat1, size1);'] ;
end //if ADD_MENU_PLOT == 1 then
//
//-----------------------------------------------------------------
//
ADD_MENU_FB=1; // set ADD_MENU 1 to
add menu button of some functions
//
//
if ADD_MENU_FB == 1 then
delmenu('step2_filter_bank');
addmenu('step2_filter_bank',[ '(1) set fb parameters','(2) do filter
bank','(-) plot3d_fb_am_out','(3) replacement', '(4) save replacement
as a file', '(-) play replacement sound' ]);
step2_filter_bank=[ '[st_freq, end_freq, n_stepf, fir_steps, kwindows]=
set_fb_para();
[bpf_fb0,wft]=make_bpf_fb0(3)','[fb_out,fb_am_out,fb_am_out_prop]=do_fb(1);
plot_fb(2)',
'plot3d_fb_am_out(1)','[wdat2,cfb_list]=do_replace(1)','snd_save2()','snd_play2()']
;
end //if ADD_MENU_FB == 1 then
//
//==============================================================++=
//
//
ADD_MENU_STATUS=1;
if ADD_MENU_STATUS == 1 then
function status1()
ver0= getversion();
wstr1=sprintf('scilab version is %s',ver0);
disp(wstr1);
endfunction
function status2()
sz=stacksize();
gsz=gstacksize();
wstr1=sprintf('current stack size is %d of %d ,(%f)',
sz(2),sz(1), sz(2)/sz(1) );
wstr2=sprintf('current global stack size is %d of %d
,(%f)', gsz(2),gsz(1), gsz(2)/gsz(1));
disp( wstr1);
disp( wstr2);
endfunction
delmenu('status_');
addmenu('status_',[ 'scilab version' ; 'check stack memory size' ]);
status_=[ 'status1()' ; 'status2()' ] ;
end //if ADD_MENU_STATUS == 1 then
//
//-------------------------------------------------------------------------------//
以上がa_wavefile_edit_70.sciです。
波形に2箇所の区間を設定して、
FFTによる周波数特性とフィルターバンクの出力を計算して、2箇所の区間を比較するためのSCILAB用のサンプルプログラム
a_wavefile_edit_80b.sci
//------------------------------------------------------------------------------------------
//
// (1)Linear-Phase FIR filter bank, AM detection, and
2D display
// (2)Compare 2 Sections on FFT frequency response
and Linear-Phase FIR filter bank out
//
// for window scilab-4.1.2
//
//
//
.........................................................................................
// PLEASE PAY ATTENTION that this program may have some bugs and also
you may modify program
// to fit your circumstances. If you use this program, everything done
by your own risk.
// Thanks.
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
T_DEMO=1; // set T_DEMO=1 for demonstration, beside set
T_DEMO=0 when normal
//
// select sample speech and its portion
if T_DEMO==1 then
SEL_CODE=x_choose(['ha1'],['Please select one'],'default')
if SEL_CODE <= 0 then
SEL_CODE=1;
end
end
//===================================================================
// 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.
wdat2=zeros(1,10);// data of synthesized .wav, 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.
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
ydat2=zeros(1,width0+1); // for portion
xziku0=[0:1:width0]; // for all data
xziku1=[0:1:width0]; // for portion
xziku2=[0:1:width0]; // for portion
sp1=1; // start point for actual
display of section 1
ep1=1; // end point for actual
display of section 1
sec2_on_off=1; // when 1, section 2 is enable. when 0, section 2 is
disenable.
sp2=1; // start point for actual
display of section 2
ep2=1; // end point for actual
display of section 2
// for fft
fft_point1=512;
db_fft1=zeros(1,fft_point1+1); // for section 1
db_fft2=zeros(1,fft_point1+1); // for section 2
phi_fft1=zeros(1,fft_point1+1); // for section 1
phi_fft2=zeros(1,fft_point1+1); // for section 2
// for others
f412=0; // flag to detect scilab 4.1.2
defch1=1; // default channel 1
f_win=1; // flag if windows
//---------------------------------------
//
//
//--------------------------------------
function [wavfile1,chs1,qty1,fs1,bit1,f412]=get_wavfile_pro()
// choose wave file
if T_DEMO == 1 then
if SEL_CODE==1 then
[x,ierr]=fileinfo('ha1.wav');
xs=size(x);
if xs(1)==0 then
disp('Choose ha1.wav file');
wavfile1=tk_getfile(Title="Choose
ha1.wav file");
else
wavfile1='ha1.wav';
end
end
else
wavfile1=tk_getfile(Title="Choose xxx.wav file name");
end
// read channels and samples
cs1=wavread(wavfile1,'size');
chs1=cs1(2); // channel
qty1=cs1(1); // samples
//...
if chs1 > 2 then // invert data for scilab-4.1.2
chs1=cs1(1);
qty1=cs1(2);
f412=1;
else
f412=0;
end
//...
[y,fs1,bit1]=wavread(wavfile1,[1 1]);
endfunction
//--------------------------------------
// function 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()
clf();
subplot(311);
plot(xziku0,ydat0,'b');
//plot2d(xziku0,ydat0,style=[color("blue")]);
wstr1='waveform all : ' + wavfile1;
xtitle( wstr1 );
// When linux scilab-3.1 2nd subplot does not work well.
subplot(312);
plot(xziku1,ydat1,'r');
//plot2d(xziku1,ydat1,style=[color("blue")]);
xtitle('waveform selected (section 1)');
az0=size(xziku0);
xziku0b=zeros(2);
ydat0b=zeros(2);
j=0;
for i=1:az0(1)
if (xziku0(i) >= sp1) & ( xziku0(i)
<= ep1) then
j=j+1;
xziku0b(j)=xziku0(i);
ydat0b(j)=ydat0(i);
end
end
if j >= 2 then
if (sp1 == 1) & ( ep1 == size1) then
// no write
else
subplot(311);
plot(xziku0b,ydat0b,'r');
end
end
//
if sec2_on_off == 1 then
subplot(313);
plot(xziku2,ydat2,'g');
//plot2d(xziku1,ydat1,style=[color("blue")]);
xtitle('waveform selected (section 2)');
az0=size(xziku0);
xziku0b=zeros(2);
ydat0b=zeros(2);
j=0;
for i=1:az0(1)
if (xziku0(i) >= sp2) & ( xziku0(i)
<= ep2) then
j=j+1;
xziku0b(j)=xziku0(i);
ydat0b(j)=ydat0(i);
end
end
if j >= 2 then
if (sp2 == 1) & ( ep2 == size1) then
// no write
else
subplot(311);
plot(xziku0b,ydat0b,'g');
end
end
//
end
endfunction
//--------------------------------------
function [ydat1,xziku1]=make_width_data( work1, wsize1,sp1,ep1)
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();
endfunction
//----------------------------------------------------------------
function [wsp1, wep1, wsp2, wep2]= reset_sp1_ep1(wsize1)
wsp1=1;
wep1=wsize1;
wsp2=1;
wep2=wsize1;
endfunction
//----------------------------------------------------------------
function [wsp1, wep1, wsec2_on_off, wsp2, wep2]= set_sp1_ep1(wsize1)
txt1=['start point(sec.1)';'end point(sec.1)';'sec.2
on(1)/off(0)' ;'start point(sec.2)';'end point(sec.2)'];
wstr1=sprintf('%d',sp1);
wstr2=sprintf('%d',ep1);
wstr3=sprintf('%d',sec2_on_off);
wstr4=sprintf('%d',sp2);
wstr5=sprintf('%d',ep2);
sig1=x_mdialog('Input start point and end point for portion
display.',txt1, [wstr1 ; wstr2; wstr3; wstr4; wstr5 ]);
if sig1==[] then
arg1=evstr(wstr1);
arg2=evstr(wstr2);
arg3=evstr(wstr3);
arg4=evstr(wstr4);
arg5=evstr(wstr5);
else
arg1=evstr(sig1(1));
arg2=evstr(sig1(2));
arg3=evstr(sig1(3));
arg4=evstr(sig1(4));
arg5=evstr(sig1(5));
end
//
//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;
disp(wep1,'ep1',wsp1,'sp1');
end
end
end
//
if arg3 == 1 then
wsec2_on_off=1;
else
wsec2_on_off=0;
end
wsp2=sp2;
wep2=ep2;
if arg4 >= 1 then
if arg5 <= wsize1 then
if arg4 < arg5 then
wsp2=arg4;
wep2=arg5;
if wsec2_on_off==1 then
disp(wep2,'ep2',wsp2,'sp2');
end
end
end
end
//
endfunction
//
//----------------------------------------------------------------
function snd_play2()
// this sound doesnot work well on windows scilab-3.1.1 and linux
scilab-3.1
// but, this sound works on windows scilab-4.1.2. But, linux
scilab-4.1.2 doesnot!
// for wdat2
sz=size(wdat2);
size2=sz(2);
if f_win == 1 then
if f412 == 1 then // for windows scilab-4.1.2
sound(wdat2 ,fs1,bit1);
else
// for windows scilab-3.1.1
if fs1 == 22050 then
sound(wdat2',fs1,bit1);
disp('This function may work, if luckily.');
elseif fs1 == 44100 then // down-sampling from
44100 to 22050
wdatw2=zeros(1,(int(size2/2)));
for v=1:(int(size2/2))
wdatw2(v)=wdat2(2 * v );
end
disp('down-sampling from 44100 to 22050');
sound(wdatw2',22050,bit1);
disp('This function may work, if luckily.');
else
//sound(wdatw,fs1,bit1);
disp('This function does not work.');
end
end
else
disp('If sound setting is good for scilab, this function maybe
work.');
if f412 == 1 then // for scilab-4.1.2
sound(wdat2 ,fs1,bit1);
else
// for scilab-3.1.1
if fs1 == 22050 then
sound(wdat2',fs1,bit1);
disp('This function may work, if luckily.');
elseif fs1 == 44100 then // down-sampling from
44100 to 22050
wdatw2=zeros(1,(int(size2/2)));
for v=1:(int(size2/2))
wdatw2(v)=wdat2(2 * v );
end
disp('down-sampling from 44100 to 22050');
sound(wdatw2',22050,bit1);
disp('This function may work, if luckily.');
else
//sound(wdatw,fs1,bit1);
disp('This function does not work.');
end
end
end
endfunction
//
//----------------------------------------------------------------
function snd_save2()
// for wdat2
//
wavefilename= input(' + enter file name for saved .wav file
=>',["string"]);
//
if f412 == 1 then // for windows scilab-4.1.2
wavwrite(wdat2 ,fs1,bit1,wavefilename );
else
// for windows scilab-3.1.1
wavwrite(wdat2',fs1,bit1, wavefilename);
end
endfunction
//
//----------------------------------------------------------------
function [fft_point1]=set_fft_points()
l1=list('points',3,['128','256','512','1024']);
wrep=x_choices('Select FFT points',list(l1));
//
fft_point1=512; // default
if wrep== 1 then
fft_point1=128;
elseif wrep== 2 then
fft_point1=256;
elseif wrep== 3 then
fft_point1=512;
elseif wrep== 4 then
fft_point1=1024;
end
endfunction
//----------------------------------------------------
function [frq1,sfrq1,is1,ie1]=set_frq( spec1024 )
// spec1024, tube model response no syuhasuu bunkainou wo ageru
tame
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 plot_fft1(disp0)
[frq1,sfrq1,is1,ie1]=set_frq(0);
//
wb0=xget('window'); // stack old window
if disp0 >= 1 then
xset('window',disp0); // create new windows
clf();
subplot(211);
gainplot(frq1,db_fft1,phi_fft1);
//plot(frq1,db_fft1);
wstr1 = 'frequency response of waveform selected (section 1): '
+ wavfile1;
wstr1 = wstr1 + ' start=' + string(sp1) + ' points=' +
string( fft_point1);
xtitle( wstr1 );
if sec2_on_off==1 then
subplot(212);
gainplot(frq1,db_fft2,phi_fft2);
//plot(frq1,db_fft1);
wstr1 = 'frequency response of waveform selected (section
2): ' + wavfile1;
wstr1 = wstr1 + ' start=' + string(sp2) + ' points='
+ string( fft_point1);
xtitle( wstr1 );
end
end // if disp0 >= 1 then
xset('window',wb0); // push old windows
endfunction
//
//----------------------------------------------------
function [db_fft1,phi_fft1]=do_fft_wav(sp1,ep1)
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
//-----------------------------------------------
//--- 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]=get_wavfile_pro();
disp_pro_wav();
[wdat1, size1]=read_one_ch_of_wav(wavfile1);
sp1=1;
ep1=size1;
sp2=1;
ep2=size1;
[ydat0,xziku0]=make_width_data( wdat1, size1,sp1,ep1);
;
if T_DEMO==1 then
if SEL_CODE==1 then // ha1
sp1=2500;
ep1=3500;
sp2=8000;
ep2=9000;
end
end
[ydat1,xziku1]=make_width_data( wdat1, size1,sp1,ep1);
[ydat2,xziku2]=make_width_data( wdat1, size1,sp2,ep2);
plot_wave1();
//
// fft
if T_DEMO==1 then
[db_fft1,phi_fft1]=do_fft_wav(sp1,ep1);
[db_fft2,phi_fft2]=do_fft_wav(sp2,ep2);
plot_fft1(5)
end
//
// -MAIN (1)program starts
//
//==============================================================++=
//
//
//===================================================================
// global variables No. 2
// for filter bank
st_freq=500.0;
end_freq=st_freq * 10;
n_stepf=10;
bpf_fb0=zeros(512,3); // *,1 center frequency.
this will be overwritten.
// *,2 lower frequency of bpf band. this will be overwritten.
// *,3 higher frequency of bpf band. this will be overwritten.
fir_steps=513; // This must be odd. this will
be overwritten.
kwindows=1;
wft=zeros(n_stepf,fir_steps); // this will be overwritten.
wfm=zeros(n_stepf,fir_steps); // this will be overwritten.
fr=zeros(fir_steps); // this will be overwritten.
fb_out=zeros(n_stepf,1); // this will be
overwritten. for section 1
fb_out2=zeros(n_stepf,1); // this will be
overwritten. for section 2
fb_am_out=zeros(n_stepf,1); // this will be overwritten. for
section 1
fb_am_out2=zeros(n_stepf,1); // this will be overwritten. for section 2
nmax_am_prop=5;
fb_am_out_prop=zeros(n_stepf,nmax_am_prop); // *,1
average value about fb_am_out
// *,2 peak value about fb_am_out
// for section 1
fb_am_out_prop2=zeros(n_stepf,nmax_am_prop); // for section 2
cfb_list=ones(1,2); // this will be overwritten.
//-------------------------------------------------
function save2()
filo=strsubst( wavfile1,'wav','dat');
save( filo
,sp1,ep1,n_stepf,bpf_fb0,fb_out,fb_am_out,fb_am_out_prop,sp2,ep2,sec2_on_off,fb_out2,fb_am_out2,fb_am_out_prop2,st_freq,end_freq,n_stepf,fir_steps,kwindows);
disp( filo, '+write' );
endfunction
function
[sp1,ep1,n_stepf,bpf_fb0,fb_out,fb_am_out,fb_am_out_prop,sp2,ep2,sec2_on_off,fb_out2,fb_am_out2,fb_am_out_prop2,st_freq,end_freq,n_stepf,fir_steps,kwindows]=load2()
filo=strsubst( wavfile1,'wav','dat');
load( filo
,'sp1','ep1','n_stepf','bpf_fb0','fb_out','fb_am_out','fb_am_out_prop','sp2','ep2','sec2_on_off','fb_out2','fb_am_out2','fb_am_out_prop2','st_freq','end_freq','n_stepf','fir_steps','kwindows');
disp( filo, '+read' );
disp('parameters are renewed. It had better to check set_et_plot
and to check set fb parameters.');
endfunction
//-------------------------------------------------
function [st_freq, end_freq, n_stepf, fir_steps, kwindows]=
set_fb_para()
// codex defines which term....
txt1=['start freq'; 'end freq'; 'freq steps'; 'steps of fir
filter (odd)'; 'priority: frq-resolution (1), drop value (2)'];
wstr1=sprintf('%f',st_freq);
wstr2=sprintf('%f',end_freq);
wstr3=sprintf('%d',n_stepf);
wstr4=sprintf('%d',fir_steps);
wstr5=sprintf('%d',kwindows);
sig1=x_mdialog('set fir filter bank parameters ',txt1, [wstr1;
wstr2; wstr3; wstr4; wstr5]);
if sig1==[] then
arg1=evstr(wstr1);
arg2=evstr(wstr2);
arg3=evstr(wstr3);
arg4=evstr(wstr4);
arg5=evstr(wstr5);
else
arg1=evstr(sig1(1));
arg2=evstr(sig1(2));
arg3=evstr(sig1(3));
arg4=evstr(sig1(4));
arg5=evstr(sig1(5));
end
st_freq=arg1;
end_freq=arg2;
n_stepf=arg3;
fir_steps=arg4;
kwindows=arg5;
endfunction
//-----------------------------------------------------------
function [bpf_fb0,wft]=make_bpf_fb0(disp0)
// log gauge divide
a1=log(st_freq);
b1=log(end_freq);
c1=(b1-a1)/n_stepf;
bpf_fb0=zeros(n_stepf,3);
for v=1:n_stepf
a2=exp(a1+ c1 * (v-1));
b2=exp(a1+ c1 * v);
c2=(b2-a2)/2. + a2;
bpf_fb0(v,1)=c2;
bpf_fb0(v,2)=a2;
bpf_fb0(v,3)=b2;
if disp0 >= 1 then
disp( [ v bpf_fb0(v,1) bpf_fb0(v,2)
bpf_fb0(v,3)]);
end
end // for V=1:n_stepf
//
// ++ mk_fir_coeff
//
w0=xget('window'); // stack old window
wft=zeros(n_stepf,fir_steps);
for vloop=1:n_stepf
cfreq(1)= bpf_fb0(vloop,2) / fs1;
cfreq(2)= bpf_fb0(vloop,3) / fs1;
fpar(1)=0.1; // dummy data
fpar(2)=-0.1; // dummy data
if kwindows == 2 then
[wft2,wfm2,fr2]=wfir('bp',fir_steps,cfreq,'hn',fpar); // gensui
riritu yuusen
else
[wft2,wfm2,fr2]=wfir('bp',fir_steps,cfreq,'hm',fpar); // bunkai
nou yuusen
end
for v2=1:fir_steps
wft(vloop,v2)=wft2(v2);
end
if disp0 >=2 then
//++ only one will be done
if vloop==1 then
xset('window',disp0); // create new windows
clf();
[frq1,sfrq1,is1,ie1]=set_frq(0);
fr0=zeros(2); // freq list selected
frq1f=zeros(2); // freq responce selected
phf0=zeros(2); // phase selected
sv0=size(fr2);
end // if vloop==1 then
//-- only one will be done
sc0=0;
for v=1:sv0(2)
fx0= fr2(v) * fs1;
if (fx0 >= frq1(1)) &(fx0 <=
frq1(sfrq1)) then
sc0=sc0+1;
fr0(sc0)=fx0;
frq1f(sc0)= 20. * log10( wfm2(v));
phf0(sc0)= -1.0 * ((fir_steps - 1)
/ 2.0 ) * ( 2. * %pi * v / fir_steps);
end
end // for v=1:sv0(2)
gainplot(fr0,frq1f,phf0);
end // vloop+1:n_stepf
end // if disp0 >=2 then
xset('window',w0); // push old windows
//
// -- mk_fir_coeff
//
endfunction
//-----------------------------------------------------------
function plot_fb(disp_code)
w0=xget('window'); // stack old window
if disp_code >= 1 then
xset('window',disp_code); // create new windows
clf();
if sec2_on_off == 1 then
bun0=2;
retu0=2;
else
bun0=1;
retu0=1;
end
wsize2=size1 - sp1 + 1;
wsize3=ep1-sp1+1;
zdata1=zeros(1,width0+1);
zdata1p=zeros(1,width0+1);
ydata1=zeros(1,width0+1);
for vloop=1:n_stepf
if wsize3 > (width0+1) then
wstep1= (ep1-sp1) / width0;
for v=1:(width0+1)
zdat1(v)= fb_out(vloop,1 + int(wstep1 * (v -
1)));
zdat1p(v)= fb_am_out(vloop, 1 + int(wstep1 *
(v - 1)));
xziku1(v)= sp1 + int(wstep1 * (v - 1));
if vloop==1 then
ydat1(v)= wdat1( sp1 + int(wstep1 * (v -
1)));
end
end
else
for v=1:(width0+1)
if v <= wsize3 then
zdat1(v)=fb_out(vloop,v);
zdat1p(v)=fb_am_out(vloop,v);
if vloop==1 then
ydat1(v)= wdat1((sp1-1)+v);
end
else
zdat1(v)=0.;
zdat1p(v)=0.;
if vloop==1 then
ydat1(v)= 0.;
end
end
xziku1(v)=(sp1-1)+v;
end
end
if vloop==1 then
subplot((n_stepf+1),bun0,1);
plot(xziku1,ydat1,'b');
// wstr1=sprintf('waveform selected,
original');
wstr1='waveform selected(section 1), original
: '+ wavfile1;
xtitle(wstr1);
end
subplot((n_stepf+1),bun0,(retu0 * vloop + 1));
plot(xziku1,zdat1,'b');
plot(xziku1,zdat1p,'g');
wstr1=sprintf('Band Pass filtered waveform
selected(section 1) (%dHz-
%dHz)',int(bpf_fb0(vloop,2)),int(bpf_fb0(vloop,3)));
mv0=fb_am_out_prop(vloop,1);
if mv0 > 0. then
mv1= int(20. * log10(mv0));
else
mv1=-120;
end
pk0=fb_am_out_prop(vloop,2);
if pk0 > 0. then
pk1= int( 20. * log10(pk0));
else
pk1=-120;
end
wstr2=sprintf('%d/%d dB',mv1,pk1); // average / peak
[dB]
xtitle( wstr1,'',wstr2 );
end // for vloop=1:n_stepf
//
// section 2
if sec2_on_off == 1 then
sp1=sp2;
ep1=ep2;
wsize2=size1 - sp1 + 1;
wsize3=ep1-sp1+1;
zdata1=zeros(1,width0+1);
zdata1p=zeros(1,width0+1);
ydata1=zeros(1,width0+1);
for vloop=1:n_stepf
if wsize3 > (width0+1) then
wstep1= (ep1-sp1) / width0;
for v=1:(width0+1)
zdat1(v)= fb_out2(vloop,1 + int(wstep1 * (v -
1)));
zdat1p(v)= fb_am_out2(vloop, 1 + int(wstep1 *
(v - 1)));
xziku1(v)= sp1 + int(wstep1 * (v - 1));
if vloop==1 then
ydat1(v)= wdat1( sp1 + int(wstep1 * (v -
1)));
end
end
else
for v=1:(width0+1)
if v <= wsize3 then
zdat1(v)=fb_out2(vloop,v);
zdat1p(v)=fb_am_out2(vloop,v);
if vloop==1 then
ydat1(v)= wdat1((sp1-1)+v);
end
else
zdat1(v)=0.;
zdat1p(v)=0.;
if vloop==1 then
ydat1(v)= 0.;
end
end
xziku1(v)=(sp1-1)+v;
end
end
sec2x=2;
if vloop==1 then
subplot((n_stepf+1),bun0, sec2x );
plot(xziku1,ydat1,'b');
// wstr1=sprintf('waveform selected,
original');
wstr1='waveform selected(section 2), original
: '+ wavfile1;
xtitle(wstr1);
end
subplot((n_stepf+1),bun0,(retu0 * vloop + sec2x));
plot(xziku1,zdat1,'b');
plot(xziku1,zdat1p,'g');
wstr1=sprintf('Band Pass filtered waveform
selected(section 2) (%dHz-
%dHz)',int(bpf_fb0(vloop,2)),int(bpf_fb0(vloop,3)));
mv0=fb_am_out_prop2(vloop,1);
if mv0 > 0. then
mv1= int(20. * log10(mv0));
else
mv1=-120;
end
pk0=fb_am_out_prop2(vloop,2);
if pk0 > 0. then
pk1= int( 20. * log10(pk0));
else
pk1=-120;
end
wstr2=sprintf('%d/%d dB',mv1,pk1); // average / peak
[dB]
xtitle( wstr1,'',wstr2 );
end // for vloop=1:n_stepf
end // if sec2_on_off == 1 then
//
end // if disp_code >= 1 then
xset('window',w0); // push old windows
endfunction
//-----------------------------------------------------------
function windowsize0(width0,height0)
dim0=xget("wdim");
disp(dim0);
xset("wdim",width0,height0); // real width and real height
is limited in physical display pixes
endfunction
function showwindowsize0()
dim0=xget("wdim");
disp(dim0);
endfunction
//-----------------------------------------------------------
function plot3d_fb_am_out(disp0)
//
wb0=xget('window'); // stack old window
if disp0 >= 1 then
xset('window',disp0); // create new windows
clf();
w0=ep1-sp1+1;
x1=zeros(n_stepf,1);
bai0=10;
for v=1:n_stepf
x1(v)=bpf_fb0(v,1) / bai0; //
value shoud be near w0
end
x1s=bpf_fb0(1,2) / bai0; // lower freq
x1e=bpf_fb0(n_stepf,3) / bai0; // higher freq
y1=zeros(w0,1);
for v=1:w0
y1(v)=sp1+(v-1);
end
pk0=fb_am_out_prop(1,2);
for v=1:n_stepf
if fb_am_out_prop(v,2) > pk0 then
pk0=fb_am_out_prop(v,2);
end
end
if pk0 <= 0. then
gain0=1000.;
else
gain0=1000./pk0;
end
fb_am_out2=zeros(n_stepf,w0);
for v=1:n_stepf
for v2=1:w0
fb_am_out2(v,v2)=fb_am_out(v,v2) * gain0;
end
end
clf();
plot3d(x1,y1, fb_am_out2 ,theta=130, alpha=75,
leg="x10Hz@point@1000",flag=[5,1,4],ebox=[x1s,x1e,y1(1),y1(w0),0.,1000.0]);
xtitle('filter bank am detect output');
end // if disp0 >= 1 then
xset('window',wb0); // push old windows
endfunction
//-----------------------------------------------------------
function [fb_out,fb_am_out,fb_am_out_prop]=do_fb(disp0,sec_no0,sp1,ep1)
//
w0=(ep1-sp1)+1;
fb_out=zeros(n_stepf,w0);
dt0= int((fir_steps-1) / 2);
ep1x=ep1+dt0;
if ep1x > size1 then
ep1x = size1;
end
sp1x=sp1+dt0;
if sp1x < fir_steps then
sp1x = fir_steps;
end
//
if disp0 >= 1 then
wstr1=sprintf('+section no.%d',sec_no0);
disp(wstr1);
end
//+filtering+
for vloop=1:n_stepf
if disp0 >= 1 then
wstr1=sprintf('+entering filter bank
no.%d/%d',vloop,n_stepf);
disp(wstr1);
end
for loop=sp1x:ep1x
ydz = 0.;
for loop2=1:fir_steps
ydz = ydz + wft(vloop,loop2) *
wdat1( 1 - loop2 + loop );
end
//fir_out1(loop - dt0)=ydz;
fb_out(vloop,loop - dt0 -sp1 + 1)=ydz;
end
end // if vloop=1:n_stepf
//-filtering-
//
//+am+
nsize_max=5000;
g_work2=zeros(10,nsize_max);
fb_am_out=zeros(n_stepf,w0);
fb_am_out_prop=zeros(n_stepf,nmax_am_prop);
fir_ou1=zeros(1,size1);
fir_out1p=zeros(1,size1);
for vloop=1:n_stepf
if disp0 >= 1 then
wstr1=sprintf('+entering AM de-modulation, peak
trace no.%d/%d',vloop,n_stepf);
disp(wstr1);
end
for v2=sp1:ep1
fir_out1(v2)= fb_out(vloop,v2-sp1+1);
fir_out1p(v2)=0.;
end
//+from trace peak
//+
spx=-1;
spy=-1;
epx=-1;
epy=-1;
c0=-1;
w0=-1.;
area0=-1.;
if (ep1-sp1) > 4 then
for v=(sp1+1):(ep1-1)
if fir_out1(v) > 0. then // only
postive side
if (fir_out1(v) > fir_out1(v-1))
& (fir_out1(v) > fir_out1(v+1)) then
spx=epx;
spy=epy;
epx=v;
epy=fir_out1(v);
c0=c0+1;
//
if g_work2(1,1) < nsize_max then
if g_work2(1,1) <= 0 then
g_work2(1,1)=1; // first start
g_work2(3,1)=epx;
g_work2(4,1)=epy;
g_work2(5,1)=epx;
g_work2(6,1)=epy;
else
d0=int(g_work2(1,1));
if epy >= g_work2(6,d0) then // rising or equal
g_work2(2,d0)=0; //current
g_work2(5,d0)=epx;
g_work2(6,d0)=epy;
else
// falling
if g_work2(2,d0)==0 then // complete data
g_work2(2,d0)=1; //complete
g_work2(7,d0)= g_work2(5,d0) - g_work2(3,d0);
g_work2(8,d0)= g_work2(6,d0) - g_work2(4,d0);
g_work2(9,d0)= g_work2(7,d0) * g_work2(8,d0) / 2.0 ;
if (g_work2(7,d0) > 0.) & (g_work2(8,d0) > 0.) then
if g_work2(9,d0) > area0 then
area0=g_work2(9,d0);
w0=g_work2(8,d0) / g_work2(7,d0);
end
end
// atan
//g_work2(10,d0)= atan( (g_work2(8,d0) / g_work2(7,d0))) / 2.0 / %pi *
360.0;
//
g_work2(1,1)=g_work2(1,1)+1; // set as new start point
d1=int(g_work2(1,1));
g_work2(3,d1)=epx;
g_work2(4,d1)=epy;
g_work2(5,d1)=epx;
g_work2(6,d1)=epy;
elseif g_work2(2,d0) <= -1 then // restart
g_work2(3,d0)=epx;
g_work2(4,d0)=epy;
g_work2(5,d0)=epx;
g_work2(6,d0)=epy;
end
end
end // if g_work2(1,1) <= 0 then
else
disp(' error: g_g_work2(1,1) >= nsize_max , in do_fb() ');
end // if g_work2(1,1) < nsize_max then
//
//
if c0 >= 1 then
a0=
(epy - spy) / (epx - spx);
b0=
epy - a0 * epx;
for
v2=spx:epx
fir_out1p(v2)= a0 * v2 + b0;
end
end // if
c0 >= 1 then
//
end // if (fir_out1(v) >
fir_out1(v-1)) & (fir_out1(v) > fir_out1(v+ 1)) then
end // if fir_out1(v) > 0. then
end //v=(sp1+1):(ep1-1)
end
if c0 < 1 then
epx=sp1-1;
end
for v=(epx+1):ep1
fir_out1p(v)=0;
end
pk0=0.;
mv0=0.;
for loop=sp1:ep1
fb_am_out(vloop,loop - sp1 + 1)= fir_out1p(loop);
if fir_out1p(loop) > pk0 then
pk0=fir_out1p(loop);
end
mv0=mv0+fir_out1p(loop);
end
fb_am_out_prop(vloop,1)= mv0 / (ep1-sp1+1);
fb_am_out_prop(vloop,2)= pk0;
//+
//-from trace peak
end // for vloop=1:n_stepf
//-am-
endfunction
//-----------------------------------------------------------------
function [wdat2,rslt]=do_replace(disp0)
sz0=size(cfb_list);
cfb_list2=ones(n_stepf,2);
for v=1:n_stepf
cfb_list2(v,1)=v;
if v <= sz0(1) then
cfb_list2(v,2)=cfb_list(v,2);
end
end
[rslt]=x_matrix('Mark 1 as Combine element, in filter bank no.
list(no, 1 or 0 )' , cfb_list2);
if rslt==[] then // reset all 1
rslt=ones(n_stepf,2);
for v=1:n_stepf
rslt(v,1)=v;
end
end
wdat2=zeros(1,size1);
for v=1:size1
wdat2(v)=wdat1(v);
end
//+replacement
for v=sp1:ep1
vs0=0.;
for v2=1:n_stepf
if rslt(v2,2)==1 then // if the fb
no is Combine element
vs0= vs0 + fb_out(v2,
(v-sp1+1))
end
end
wdat2(v)=vs0;
end
//-replacement
if disp0 >= 1 then
ydat1r=zeros(1,(width0+1));
if size1 > (width0+1) then
wstep1= size1 / (width0+1);
for v=1:(width0+1)
ydat1r(v)= wdat2( 1 + int(wstep1 * (v -
1)));
// xziku1(v)= 1 + int(wstep1 * (v - 1));
end
else
for v=1:(width0+1)
if v <= size1 then
ydat1r(v)=wdat2(v);
else
ydat1r(v)=0.;
end
// xziku1(v)=(sp1-1)+v;
end
end
//
clf();
subplot(211);
plot(xziku0,ydat0,'b');
//plot2d(xziku0,ydat0,style=[color("blue")]);
xtitle('waveform all');
// When linux scilab-3.1 2nd subplot does not work well.
subplot(212);
plot(xziku0,ydat1r,'b');
//plot2d(xziku1,ydat1,style=[color("blue")]);
xtitle('replacement waveform');
end // id disp0 >= 1 thne
endfunction
//=================================================================
//
// +MAIN (2) program starts
// filter bank
//
if T_DEMO==1 then
[bpf_fb0,wft]=make_bpf_fb0(3);
end
if T_DEMO==1 then
[fb_out,fb_am_out,fb_am_out_prop]=do_fb(1,1,sp1,ep1);
[fb_out2,fb_am_out2,fb_am_out_prop2]=do_fb(1,2,sp2,ep2);
plot_fb(2);
end
//
// -MAIN (2)program starts
//
//==============================================================++=+
//
//
// ++ MENU MENU MENU =====================
//
// 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',[ '(1) set_et_plot()'; '
reset_plot()'; '(op1) set fft points'; '(op2) do fft' ]);
step1_plotwav=[ '[sp1, ep1, sec2_on_off, sp2, ep2 ]=
set_sp1_ep1(size1); [ydat1,xziku1]=make_width_data( wdat1,
size1,sp1,ep1); [ydat2,xziku2]=make_width_data( wdat1, size1,sp2,ep2);
plot_wave1();' ; '[sp1, ep1, sp2, ep2]= reset_sp1_ep1(size1);
[ydat1,xziku1]=make_width_data( wdat1, size1);
[ydat2,xziku2]=make_width_data( wdat1, size1,sp2,ep2);
plot_wave1();'; '[fft_point1]=set_fft_points();' ;
'[db_fft1,phi_fft1]=do_fft_wav(sp1,ep1);
[db_fft2,phi_fft2]=do_fft_wav(sp2,ep2); plot_fft1(5)' ] ;
end //if ADD_MENU_PLOT == 1 then
//
//-----------------------------------------------------------------
//
ADD_MENU_FB=1; // set ADD_MENU 1 to
add menu button of some functions
//
//
if ADD_MENU_FB == 1 then
delmenu('step2_filter_bank');
addmenu('step2_filter_bank',[ '(1) set fb parameters','(2-1) do filter
bank (sec.1)','(2-2) do filter bank (sec.2)','(2-e)
plot2d_fb_am_out','(3) replacement', '(4) save replacement as a file',
'(-) play replacement sound' ]);
step2_filter_bank=[ '[st_freq, end_freq, n_stepf, fir_steps, kwindows]=
set_fb_para();
[bpf_fb0,wft]=make_bpf_fb0(3)','[fb_out,fb_am_out,fb_am_out_prop]=do_fb(1,1,sp1,ep1)',
'[fb_out2,fb_am_out2,fb_am_out_prop2]=do_fb(1,2,sp2,ep2)','plot_fb(2)','[wdat2,cfb_list]=do_replace(1)','snd_save2()','snd_play2()']
;
end //if ADD_MENU_FB == 1 then
//-----------------------------------------------------------------
//
ADD_MENU_FB2=1; // set ADD_MENU 1
to add menu button of some functions
//
//
if ADD_MENU_FB2 == 1 then
delmenu('step2_filter_bank_backup');
addmenu('step2_filter_bank_backup',[ '(-) save','(-) load' ]);
step2_filter_bank_backup=['save2()','[sp1,ep1,n_stepf,bpf_fb0,fb_out,fb_am_out,fb_am_out_prop,sp2,ep2,sec2_on_off,fb_out2,fb_am_out2,fb_am_out_prop2,st_freq,end_freq,n_stepf,fir_steps,kwindows]=load2();
plot_fb(2);'] ;
end //if ADD_MENU_FB == 1 then
//
//==============================================================++=
//
//
ADD_MENU_STATUS=1;
if ADD_MENU_STATUS == 1 then
function status1()
ver0= getversion();
wstr1=sprintf('scilab version is %s',ver0);
disp(wstr1);
endfunction
function status2()
sz=stacksize();
gsz=gstacksize();
wstr1=sprintf('current stack size is %d of %d ,(%f)',
sz(2),sz(1), sz(2)/sz(1) );
wstr2=sprintf('current global stack size is %d of %d
,(%f)', gsz(2),gsz(1), gsz(2)/gsz(1));
disp( wstr1);
disp( wstr2);
endfunction
delmenu('status_');
addmenu('status_',[ 'scilab version' ; 'check stack memory size' ]);
status_=[ 'status1()' ; 'status2()' ] ;
end //if ADD_MENU_STATUS == 1 then
//
//-------------------------------------------------------------------------------//
以上がa_wavefile_edit_80b.sciです。
SCILABのevent
handlerイベントハンドラーを設定することで、キー入力（キーd(down)とキーu(up))を使ってフィルターバンク
の出力を３バンド単位で移動表示し、キー入力（キーf(表示のon/off))で周波数の変化を折れ線で表示できるようにしたSCILAB用のサンプルデ
モプ
ログラム
a_wavefile_edit_84c.sci
以下は、メニューの中で na1と ma1と nn1 without low のいずれかを選択するときに使うサンプルのファイルです。
表示例
//------------------------------------------------------------------------------------------
//
// (1)Linear-Phase FIR filter bank, AM detection, FM
detection, and 2D display
// (2)Contrast 2 Sections
//
// for window scilab-4.1.2
// event handler
//
// From version _81_v32, save&load
data was changed.
// From version _84_v63, save&load
data was changed.
//
.........................................................................................
// PLEASE PAY ATTENTION that this program may have some bugs and also
you may modify program
// to fit your circumstances. If you use this program, everything done
by your own risk.
// Thanks.
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
T_DEMO=1; // set T_DEMO=1 for demonstration, beside set
T_DEMO=0 when normal
T_N_STEPF=0; // set T_N_STEPF=1 for manual set in DEMO,
beside set T_N_STEPF=0 when normal
T_SP_EP_SET_SKIP=1; // skip seting sp ep when T_DEMO=1
//
// select sample speech and its portion
if T_DEMO==1 then
//SEL_CODE=x_choose(['ma1';'na1';'ma1-fa';'na1-fa'],['Please
select one'],'default')
SEL_CODE=x_choose(['ma1';'na1';'nn1 without low'],['Please select
one'],'default')
if SEL_CODE <= 0 then
SEL_CODE=1;
end
end
//===================================================================
// scilab global variables
global eh_var1
global eh_var2
global eh_var3
global eh_var4
global eh_var5
//...
global eh_fm_disp // when 0, no fm detection display.
eh_fm_disp0=0;
//
//===================================================================
// 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.
wdat2=zeros(1,10);// data of synthesized .wav, 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.
Min_Freq=100; // 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
ydat2=zeros(1,width0+1); // for portion
xziku0=[0:1:width0]; // for all data
xziku1=[0:1:width0]; // for portion
xziku2=[0:1:width0]; // for portion
sp1=1; // start point for actual
display of section 1
ep1=1; // end point for actual
display of section 1
sec2_on_off=1; // when 1, section 2 is enable. when 0, section 2 is
disenable.
sp2=1; // start point for actual
display of section 2
ep2=1; // end point for actual
display of section 2
// for fft
fft_point1=512;
db_fft1=zeros(1,fft_point1+1); // for section 1
db_fft2=zeros(1,fft_point1+1); // for section 2
phi_fft1=zeros(1,fft_point1+1); // for section 1
phi_fft2=zeros(1,fft_point1+1); // for section 2
num_fft_sp1_list=10;
fft_sp1_list=zeros(num_fft_sp1_list,2); // this will be
overwritten.
// for others
f412=0; // flag to detect scilab 4.1.2
defch1=1; // default channel 1
f_win=1; // flag if windows
//---------------------------------------
//
//
//--------------------------------------
function [wavfile1,chs1,qty1,fs1,bit1,f412]=get_wavfile_pro()
// choose wave file
if T_DEMO == 1 then
if SEL_CODE==5 then
[x,ierr]=fileinfo('na1.wav');
xs=size(x);
if xs(1)==0 then
disp('Choose na1.wav file');
wavfile1=tk_getfile(Title="Choose
na1.wav file");
else
wavfile1='na1.wav';
end
elseif SEL_CODE==4 then
[x,ierr]=fileinfo('na1-fa.wav');
xs=size(x);
if xs(1)==0 then
disp('Choose na1-fa.wav file');
wavfile1=tk_getfile(Title="Choose
na1-fa.wav file");
else
wavfile1='na1-fa.wav';
end
elseif SEL_CODE==3 then
[x,ierr]=fileinfo('nn1-without-low.wav');
xs=size(x);
if xs(1)==0 then
disp('Choose nn1-without-low.wav file');
wavfile1=tk_getfile(Title="Choose
nn1-without-low.wav file");
else
wavfile1='nn1-without-low.wav';
end
elseif SEL_CODE==2 then
[x,ierr]=fileinfo('na1.wav');
xs=size(x);
if xs(1)==0 then
disp('Choose na1.wav file');
wavfile1=tk_getfile(Title="Choose
na1.wav file");
else
wavfile1='na1.wav';
end
else
[x,ierr]=fileinfo('ma1.wav');
xs=size(x);
if xs(1)==0 then
disp('Choose ma1.wav file');
wavfile1=tk_getfile(Title="Choose
ma1.wav file");
else
wavfile1='ma1.wav';
end
end
else
wavfile1=tk_getfile(Title="Choose xxx.wav file name");
end
// read channels and samples
cs1=wavread(wavfile1,'size');
chs1=cs1(2); // channel
qty1=cs1(1); // samples
//...
if chs1 > 2 then // invert data for scilab-4.1.2
chs1=cs1(1);
qty1=cs1(2);
f412=1;
else
f412=0;
end
//...
[y,fs1,bit1]=wavread(wavfile1,[1 1]);
endfunction
//--------------------------------------
// function 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()
clf();
subplot(311);
plot(xziku0,ydat0,'b');
//plot2d(xziku0,ydat0,style=[color("blue")]);
wstr1='waveform all : ' + wavfile1;
xtitle( wstr1 );
// When linux scilab-3.1 2nd subplot does not work well.
subplot(312);
plot(xziku1,ydat1,'r');
//plot2d(xziku1,ydat1,style=[color("blue")]);
xtitle('waveform selected (section 1)');
az0=size(xziku0);
xziku0b=zeros(2);
ydat0b=zeros(2);
j=0;
for i=1:az0(1)
if (xziku0(i) >= sp1) & ( xziku0(i)
<= ep1) then
j=j+1;
xziku0b(j)=xziku0(i);
ydat0b(j)=ydat0(i);
end
end
if j >= 2 then
if (sp1 == 1) & ( ep1 == size1) then
// no write
else
subplot(311);
plot(xziku0b,ydat0b,'r');
end
end
//
if sec2_on_off == 1 then
subplot(313);
plot(xziku2,ydat2,'g');
//plot2d(xziku1,ydat1,style=[color("blue")]);
xtitle('waveform selected (section 2)');
az0=size(xziku0);
xziku0b=zeros(2);
ydat0b=zeros(2);
j=0;
for i=1:az0(1)
if (xziku0(i) >= sp2) & ( xziku0(i)
<= ep2) then
j=j+1;
xziku0b(j)=xziku0(i);
ydat0b(j)=ydat0(i);
end
end
if j >= 2 then
if (sp2 == 1) & ( ep2 == size1) then
// no write
else
subplot(311);
plot(xziku0b,ydat0b,'g');
end
end
//
end
endfunction
//--------------------------------------
function [ydat1,xziku1]=make_width_data( work1, wsize1,sp1,ep1)
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();
endfunction
//----------------------------------------------------------------
function [wsp1, wep1, wsp2, wep2]= reset_sp1_ep1(wsize1)
wsp1=1;
wep1=wsize1;
wsp2=1;
wep2=wsize1;
endfunction
//----------------------------------------------------------------
function [wsp1, wep1, wsec2_on_off, wsp2, wep2]= set_sp1_ep1(wsize1)
txt1=['start point(sec.1)';'end point(sec.1)';'sec.2
on(1)/off(0)' ;'start point(sec.2)';'end point(sec.2)'];
wstr1=sprintf('%d',sp1);
wstr2=sprintf('%d',ep1);
wstr3=sprintf('%d',sec2_on_off);
wstr4=sprintf('%d',sp2);
wstr5=sprintf('%d',ep2);
sig1=x_mdialog('Input start point and end point for portion
display.',txt1, [wstr1 ; wstr2; wstr3; wstr4; wstr5 ]);
if sig1==[] then
arg1=evstr(wstr1);
arg2=evstr(wstr2);
arg3=evstr(wstr3);
arg4=evstr(wstr4);
arg5=evstr(wstr5);
else
arg1=evstr(sig1(1));
arg2=evstr(sig1(2));
arg3=evstr(sig1(3));
arg4=evstr(sig1(4));
arg5=evstr(sig1(5));
end
//
//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;
disp(wep1,'ep1',wsp1,'sp1');
end
end
end
//
if arg3 == 1 then
wsec2_on_off=1;
else
wsec2_on_off=0;
end
wsp2=sp2;
wep2=ep2;
if arg4 >= 1 then
if arg5 <= wsize1 then
if arg4 < arg5 then
wsp2=arg4;
wep2=arg5;
if wsec2_on_off==1 then
disp(wep2,'ep2',wsp2,'sp2');
end
end
end
end
//
endfunction
//
//----------------------------------------------------------------
function snd_play2()
// this sound doesnot work well on windows scilab-3.1.1 and linux
scilab-3.1
// but, this sound works on windows scilab-4.1.2. But, linux
scilab-4.1.2 doesnot!
// for wdat2
sz=size(wdat2);
size2=sz(2);
if f_win == 1 then
if f412 == 1 then // for windows scilab-4.1.2
sound(wdat2 ,fs1,bit1);
else
// for windows scilab-3.1.1
if fs1 == 22050 then
sound(wdat2',fs1,bit1);
disp('This function may work, if luckily.');
elseif fs1 == 44100 then // down-sampling from
44100 to 22050
wdatw2=zeros(1,(int(size2/2)));
for v=1:(int(size2/2))
wdatw2(v)=wdat2(2 * v );
end
disp('down-sampling from 44100 to 22050');
sound(wdatw2',22050,bit1);
disp('This function may work, if luckily.');
else
//sound(wdatw,fs1,bit1);
disp('This function does not work.');
end
end
else
disp('If sound setting is good for scilab, this function maybe
work.');
if f412 == 1 then // for scilab-4.1.2
sound(wdat2 ,fs1,bit1);
else
// for scilab-3.1.1
if fs1 == 22050 then
sound(wdat2',fs1,bit1);
disp('This function may work, if luckily.');
elseif fs1 == 44100 then // down-sampling from
44100 to 22050
wdatw2=zeros(1,(int(size2/2)));
for v=1:(int(size2/2))
wdatw2(v)=wdat2(2 * v );
end
disp('down-sampling from 44100 to 22050');
sound(wdatw2',22050,bit1);
disp('This function may work, if luckily.');
else
//sound(wdatw,fs1,bit1);
disp('This function does not work.');
end
end
end
endfunction
//
//----------------------------------------------------------------
function snd_save2()
// for wdat2
//
wavefilename= input(' + enter file name for saved .wav file
=>',["string"]);
//
if f412 == 1 then // for windows scilab-4.1.2
wavwrite(wdat2 ,fs1,bit1,wavefilename );
else
// for windows scilab-3.1.1
wavwrite(wdat2',fs1,bit1, wavefilename);
end
endfunction
//
//----------------------------------------------------------------
function [fft_point1]=set_fft_points()
l1=list('points',3,['128','256','512','1024']);
wrep=x_choices('Select FFT points',list(l1));
//
fft_point1=512; // default
if wrep== 1 then
fft_point1=128;
elseif wrep== 2 then
fft_point1=256;
elseif wrep== 3 then
fft_point1=512;
elseif wrep== 4 then
fft_point1=1024;
end
endfunction
//----------------------------------------------------
function [frq1,sfrq1,is1,ie1]=set_frq( spec1024 )
// spec1024, tube model response no syuhasuu bunkainou wo ageru
tame
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 plot_fft1(disp0)
[frq1,sfrq1,is1,ie1]=set_frq(0);
//
wb0=xget('window'); // stack old window
if disp0 >= 1 then
xset('window',disp0); // create new windows
clf();
subplot(211);
gainplot(frq1,db_fft1,phi_fft1);
//plot(frq1,db_fft1);
wstr1 = 'frequency response of waveform selected (section 1): '
+ wavfile1;
wstr1 = wstr1 + ' start=' + string(sp1) + ' points=' +
string( fft_point1);
xtitle( wstr1 );
if sec2_on_off==1 then
subplot(212);
gainplot(frq1,db_fft2,phi_fft2);
//plot(frq1,db_fft1);
wstr1 = 'frequency response of waveform selected (section
2): ' + wavfile1;
wstr1 = wstr1 + ' start=' + string(sp2) + ' points='
+ string( fft_point1);
xtitle( wstr1 );
end
end // if disp0 >= 1 then
xset('window',wb0); // push old windows
endfunction
//
//----------------------------------------------------
function [db_fft1,phi_fft1]=do_fft_wav(sp1,ep1)
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 plot_multi_fft1(disp0)
[frq1,sfrq1,is1,ie1]=set_frq(0);
//
wb0=xget('window'); // stack old window
if disp0 >= 1 then
xset('window',disp0); // create new windows
clf();
n_fft=0;
for v=1:num_fft_sp1_list
if fft_sp1_list(v,2) < 1 then
break;
elseif (fft_sp1_list(v,2)+fft_point1 - 1) > size1
then
break;
else
n_fft=n_fft+1;
end
end
//...
subplot(n_fft+1,1,1);
plot(xziku0,ydat0,'b');
//plot2d(xziku0,ydat0,style=[color("blue")]);
wstr1='waveform all : ' + wavfile1;
xtitle( wstr1 );
//....
for vloop=1:n_fft
sp1=fft_sp1_list(vloop,2);
ep1=sp1+fft_point1 - 1;
az0=size(xziku0);
xziku0b=zeros(2);
ydat0b=zeros(2);
j=0;
for i=1:az0(1)
if (xziku0(i) >= sp1) & ( xziku0(i)
<= ep1) then
j=j+1;
xziku0b(j)=xziku0(i);
ydat0b(j)=ydat0(i);
end
end
if j >= 2 then
k0=int(vloop/2);
k1=vloop- k0*2;
if k1== 1 then
plot(xziku0b,ydat0b,'r');
else
plot(xziku0b,ydat0b,'g');
end
end
end //for vloop=1:n_fft
//..
for vloop=1:n_fft
sp1=fft_sp1_list(vloop,2);
ep1=sp1+fft_point1 - 1;
[db_fft1,phi_fft1]=do_fft_wav(sp1,ep1);
subplot(n_fft+1,1,vloop+1);
gainplot(frq1,db_fft1,phi_fft1);
//plot(frq1,db_fft1);
wstr1 = 'frequency response of waveform selected: '
+ wavfile1;
wstr1 = wstr1 + ' start=' + string(sp1) + '
points=' + string( fft_point1);
xtitle( wstr1 );
end // for loop=1:n_fft
end // if disp0 >= 1 then
xset('window',wb0); // push old windows
endfunction
//--------------------------------------------------------------------------
//
function [rslt]= set_sp1s()
fft_sp1_list2=zeros(num_fft_sp1_list,2);
for v=1:num_fft_sp1_list
fft_sp1_list2(v,1)=v;
fft_sp1_list2(v,2)=fft_sp1_list(v,2);
end
[rslt]=x_matrix('Set (num,sp1), sp1 for each fft section. When
sp1 is 0, nothing done. ' , fft_sp1_list2);
if rslt==[] then // reset all 1
rslt=zeros(num_fft_sp1_list,2);
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]=get_wavfile_pro();
disp_pro_wav();
[wdat1, size1]=read_one_ch_of_wav(wavfile1);
sp1=1;
ep1=size1;
sp2=1;
ep2=size1;
[ydat0,xziku0]=make_width_data( wdat1, size1,sp1,ep1);
;
if T_DEMO==1 then
if T_SP_EP_SET_SKIP==1 then
else
if SEL_CODE==5 then
fft_sp1_list(1,2)=1500;
fft_sp1_list(2,2)=2650;
fft_sp1_list(3,2)=4050;
fft_sp1_list(4,2)=5550;
fft_sp1_list(5,2)=7800;
elseif SEL_CODE==4 then
fft_sp1_list(1,2)=500;
fft_sp1_list(2,2)=1800;
fft_sp1_list(3,2)=2550;
fft_sp1_list(4,2)=4600;
fft_sp1_list(5,2)=7800;
elseif SEL_CODE==3 then
elseif SEL_CODE==2 then
elseif SEL_CODE==1 then
sp1=2500;
ep1=3500;
sp2=8000;
ep2=9000;
end
end
end
[ydat1,xziku1]=make_width_data( wdat1, size1,sp1,ep1);
[ydat2,xziku2]=make_width_data( wdat1, size1,sp2,ep2);
plot_wave1();
//
// fft
if T_DEMO==1 then
if T_SP_EP_SET_SKIP==1 then
else
[db_fft1,phi_fft1]=do_fft_wav(sp1,ep1);
[db_fft2,phi_fft2]=do_fft_wav(sp2,ep2);
plot_fft1(5)
end
end
//
// -MAIN (1)program starts
//
//==============================================================++=
//
//
//===================================================================
// global variables No. 2
// for filter bank
st_freq=100.0;
end_freq=st_freq * 15;
n_stepf=10;
bpf_fb0=zeros(512,3); // *,1 center frequency.
this will be overwritten.
// *,2 lower frequency of bpf band. this will be overwritten.
// *,3 higher frequency of bpf band. this will be overwritten.
n_stepfn=10;
// change bpf data dimension
bpf_fbn=zeros(512,3); // change bpf data
fb_ctable=zeros(n_stepfn+1,1); // change map table: if
fb_ctable(1) is 1, then there is change, 0, no change. if -1 is not
defined yet.
// fb_ctable(2,....n_stepfn) is map index of old bpf_fbn, when no
map, it's 0
fb_ctable(1)=-1; // set that this is not defined yet.
fir_steps=1013; // This must be odd. this will
be overwritten.
kwindows=1;
wft=zeros(n_stepf,fir_steps); // this will be overwritten.
wfm=zeros(n_stepf,fir_steps); // this will be overwritten.
fr=zeros(fir_steps); // this will be overwritten.
wftn=zeros(n_stepfn,fir_steps); // for change, this will be
overwritten.
fb_out=zeros(n_stepf,1); // this will be
overwritten. for section 1
fb_out2=zeros(n_stepf,1); // this will be
overwritten. for section 2
fb_am_out=zeros(n_stepf,1); // this will be overwritten. for
section 1
fb_am_out2=zeros(n_stepf,1); // this will be overwritten. for section 2
nmax_am_prop=5;
fb_am_out_prop=zeros(n_stepf,nmax_am_prop); // *,1
average value about fb_am_out
// *,2 peak value about fb_am_out
// for section 1
fb_am_out_prop2=zeros(n_stepf,nmax_am_prop); // for section 2
fb_outn=zeros(n_stepfn,1); // for change, this will
be overwritten. for section 1
fb_out2n=zeros(n_stepfn,1); // for change, this will
be overwritten. for section 2
fb_am_outn=zeros(n_stepfn,1); // for change, this will be
overwritten. for section 1
fb_am_out2n=zeros(n_stepfn,1); // for change, this will be overwritten.
for section 2
fb_am_out_propn=zeros(n_stepfn,nmax_am_prop); // for change,
fb_am_out_prop2n=zeros(n_stepfn,nmax_am_prop); // for change,
fb_fm_out=zeros(n_stepf,1); // this will be overwritten. for
section 1
fb_fm_out2=zeros(n_stepf,1); // this will be overwritten. for section 2
fb_fm_out_prop=zeros(n_stepf,nmax_am_prop); // *,1
start point
// *,2 end point
// *,3 max value
// *,4 min value
// for section 1
fb_fm_out_prop2=zeros(n_stepf,nmax_am_prop); // for section 2
fb_fm_outn=zeros(n_stepf,1); // for change,
fb_fm_out2n=zeros(n_stepf,1); // for change,
fb_fm_out_propn=zeros(n_stepf,nmax_am_prop); // for change,
fb_fm_out_prop2n=zeros(n_stepf,nmax_am_prop); // for change,
cfb_list=ones(1,2); // this will be overwritten.
//-------------------------------------------------
function save2()
// From version _81_v32, save&load data was changed.
filo=strsubst( wavfile1,'wav','dat');
save( filo
,sp1,ep1,n_stepf,bpf_fb0,fb_out,fb_am_out,fb_am_out_prop,sp2,ep2,sec2_on_off,fb_out2,fb_am_out2,fb_am_out_prop2,st_freq,end_freq,n_stepf,fir_steps,kwindows);
disp( filo, '+write' );
endfunction
function save3()
// From version _84_v63, save&load data was changed.
filo=strsubst( wavfile1,'wav','dat');
save( filo
,sp1,ep1,n_stepf,bpf_fb0,fb_out,fb_am_out,fb_am_out_prop,sp2,ep2,sec2_on_off,fb_out2,fb_am_out2,fb_am_out_prop2,st_freq,end_freq,n_stepf,fir_steps,kwindows,fb_fm_out,fb_fm_out_prop,fb_fm_out2,fb_fm_out_prop2,wft,bpf_fb0);
disp( filo, '+write' );
endfunction
function
[sp1,ep1,n_stepf,bpf_fb0,fb_out,fb_am_out,fb_am_out_prop,sp2,ep2,sec2_on_off,fb_out2,fb_am_out2,fb_am_out_prop2,st_freq,end_freq,n_stepf,fir_steps,kwindows]=load2()
// From version _81_v32, save&load data was changed.
filo=strsubst( wavfile1,'wav','dat');
load( filo
,'sp1','ep1','n_stepf','bpf_fb0','fb_out','fb_am_out','fb_am_out_prop','sp2','ep2','sec2_on_off','fb_out2','fb_am_out2','fb_am_out_prop2','st_freq','end_freq','n_stepf','fir_steps','kwindows');
disp( filo, '+read' );
disp('parameters are renewed. It had better to check set_et_plot
and to check set fb parameters.');
endfunction
function
[sp1,ep1,n_stepf,bpf_fb0,fb_out,fb_am_out,fb_am_out_prop,sp2,ep2,sec2_on_off,fb_out2,fb_am_out2,fb_am_out_prop2,st_freq,end_freq,n_stepf,fir_steps,kwindows,fb_fm_out,fb_fm_out_prop,fb_fm_out2,fb_fm_out_prop2,wft,bpf_fb0]=load3()
// From version _84_v63, save&load data was changed.
filo=strsubst( wavfile1,'wav','dat');
//
[nams,typs,dims,vols]=listvarinfile(filo);
sa0=size(nams);
if sa0(1) == 18 then
load( filo
,'sp1','ep1','n_stepf','bpf_fb0','fb_out','fb_am_out','fb_am_out_prop','sp2','ep2','sec2_on_off','fb_out2','fb_am_out2','fb_am_out_prop2','st_freq','end_freq','n_stepf','fir_steps','kwindows');
fb_fm_out=zeros(n_stepf,1);
fb_fm_out_prop=zeros(n_stepf,nmax_am_prop);
fb_fm_out2=zeros(n_stepf,1);
fb_fm_out_prop2=zeros(n_stepf,nmax_am_prop);
wft=zeros(n_stepf,fir_steps);
bpf_fb0=zeros(512,3);
disp("Waring: The file is old version _81_v32. fb_fm_out,
fb_fm_out_prop, fb_fm_out2, fb_fm_out_prop2, wft, and bpf_fb0 are set
as empty data. ');
else
load( filo
,'sp1','ep1','n_stepf','bpf_fb0','fb_out','fb_am_out','fb_am_out_prop','sp2','ep2','sec2_on_off','fb_out2','fb_am_out2','fb_am_out_prop2','st_freq','end_freq','n_stepf','fir_steps','kwindows','fb_fm_out','fb_fm_out_prop','fb_fm_out2','fb_fm_out_prop2','wft','bpf_fb0');
end
disp( filo, '+read' );
disp('parameters are renewed. It had better to check set_et_plot
and to check set fb parameters.');
endfunction
//-------------------------------------------------
function [st_freq, end_freq, n_stepf, fir_steps, kwindows]=
set_fb_para()
// codex defines which term....
txt1=['start freq'; 'end freq'; 'freq steps'; 'steps of fir
filter (odd)'; 'priority: frq-resolution (1), drop value (2)'];
wstr1=sprintf('%f',st_freq);
wstr2=sprintf('%f',end_freq);
wstr3=sprintf('%d',n_stepf);
wstr4=sprintf('%d',fir_steps);
wstr5=sprintf('%d',kwindows);
sig1=x_mdialog('set fir filter bank parameters ',txt1, [wstr1;
wstr2; wstr3; wstr4; wstr5]);
if sig1==[] then
arg1=evstr(wstr1);
arg2=evstr(wstr2);
arg3=evstr(wstr3);
arg4=evstr(wstr4);
arg5=evstr(wstr5);
else
arg1=evstr(sig1(1));
arg2=evstr(sig1(2));
arg3=evstr(sig1(3));
arg4=evstr(sig1(4));
arg5=evstr(sig1(5));
end
st_freq=arg1;
end_freq=arg2;
n_stepf=arg3;
fir_steps=arg4;
kwindows=arg5;
endfunction
//-----------------------------------------------------------
function [bpf_fb0,wft]=make_bpf_fb0(disp0)
// log gauge divide
a1=log(st_freq);
b1=log(end_freq);
c1=(b1-a1)/n_stepf;
bpf_fb0=zeros(n_stepf,3);
for v=1:n_stepf
a2=exp(a1+ c1 * (v-1));
b2=exp(a1+ c1 * v);
c2=(b2-a2)/2. + a2;
bpf_fb0(v,1)=c2;
bpf_fb0(v,2)=a2;
bpf_fb0(v,3)=b2;
if disp0 >= 1 then
disp( [ v bpf_fb0(v,1) bpf_fb0(v,2)
bpf_fb0(v,3)]);
end
end // for V=1:n_stepf
//
// ++ mk_fir_coeff
//
w0=xget('window'); // stack old window
wft=zeros(n_stepf,fir_steps);
for vloop=1:n_stepf
cfreq(1)= bpf_fb0(vloop,2) / fs1;
cfreq(2)= bpf_fb0(vloop,3) / fs1;
fpar(1)=0.1; // dummy data
fpar(2)=-0.1; // dummy data
if kwindows == 2 then
[wft2,wfm2,fr2]=wfir('bp',fir_steps,cfreq,'hn',fpar); // gensui
riritu yuusen
else
[wft2,wfm2,fr2]=wfir('bp',fir_steps,cfreq,'hm',fpar); // bunkai
nou yuusen
end
for v2=1:fir_steps
wft(vloop,v2)=wft2(v2);
end
if disp0 >=2 then
//++ only one will be done
if vloop==1 then
xset('window',disp0); // create new windows
clf();
[frq1,sfrq1,is1,ie1]=set_frq(0);
fr0=zeros(2); // freq list selected
frq1f=zeros(2); // freq responce selected
phf0=zeros(2); // phase selected
sv0=size(fr2);
end // if vloop==1 then
//-- only one will be done
sc0=0;
for v=1:sv0(2)
fx0= fr2(v) * fs1;
if (fx0 >= frq1(1)) &(fx0 <=
frq1(sfrq1)) then
sc0=sc0+1;
fr0(sc0)=fx0;
frq1f(sc0)= 20. * log10( wfm2(v));
phf0(sc0)= -1.0 * ((fir_steps - 1)
/ 2.0 ) * ( 2. * %pi * v / fir_steps);
end
end // for v=1:sv0(2)
gainplot(fr0,frq1f,phf0);
end // vloop+1:n_stepf
end // if disp0 >=2 then
xset('window',w0); // push old windows
//
// -- mk_fir_coeff
//
endfunction
//-----------------------------------------------------------
function plot_fb(disp_code)
w0=xget('window'); // stack old window
if disp_code >= 1 then
xset('window',disp_code); // create new windows
clf();
if sec2_on_off == 1 then
bun0=2;
retu0=2;
else
bun0=1;
retu0=1;
end
wsize2=size1 - sp1 + 1;
wsize3=ep1-sp1+1;
zdata1=zeros(1,width0+1);
zdata1p=zeros(1,width0+1);
ydata1=zeros(1,width0+1);
for vloop=1:n_stepf
if wsize3 > (width0+1) then
wstep1= (ep1-sp1) / width0;
for v=1:(width0+1)
zdat1(v)= fb_out(vloop,1 + int(wstep1 * (v -
1)));
zdat1p(v)= fb_am_out(vloop, 1 + int(wstep1 *
(v - 1)));
xziku1(v)= sp1 + int(wstep1 * (v - 1));
if vloop==1 then
ydat1(v)= wdat1( sp1 + int(wstep1 * (v -
1)));
end
end
else
for v=1:(width0+1)
if v <= wsize3 then
zdat1(v)=fb_out(vloop,v);
zdat1p(v)=fb_am_out(vloop,v);
if vloop==1 then
ydat1(v)= wdat1((sp1-1)+v);
end
else
zdat1(v)=0.;
zdat1p(v)=0.;
if vloop==1 then
ydat1(v)= 0.;
end
end
xziku1(v)=(sp1-1)+v;
end
end
if vloop==1 then
subplot((n_stepf+1),bun0,1);
plot(xziku1,ydat1,'b');
// wstr1=sprintf('waveform selected,
original');
wstr1='waveform selected(section 1), original
: '+ wavfile1;
xtitle(wstr1);
end
subplot((n_stepf+1),bun0,(retu0 * vloop + 1));
plot(xziku1,zdat1,'b');
plot(xziku1,zdat1p,'g');
wstr1=sprintf('Band Pass filtered waveform
selected(section 1) (%dHz-
%dHz)',int(bpf_fb0(vloop,2)),int(bpf_fb0(vloop,3)));
mv0=fb_am_out_prop(vloop,1);
if mv0 > 0. then
mv1= int(20. * log10(mv0));
else
mv1=-120;
end
pk0=fb_am_out_prop(vloop,2);
if pk0 > 0. then
pk1= int( 20. * log10(pk0));
else
pk1=-120;
end
wstr2=sprintf('%d/%d dB',mv1,pk1); // average / peak
[dB]
xtitle( wstr1,'',wstr2 );
end // for vloop=1:n_stepf
//
// section 2
if sec2_on_off == 1 then
sp1=sp2;
ep1=ep2;
wsize2=size1 - sp1 + 1;
wsize3=ep1-sp1+1;
zdata1=zeros(1,width0+1);
zdata1p=zeros(1,width0+1);
ydata1=zeros(1,width0+1);
for vloop=1:n_stepf
if wsize3 > (width0+1) then
wstep1= (ep1-sp1) / width0;
for v=1:(width0+1)
zdat1(v)= fb_out2(vloop,1 + int(wstep1 * (v -
1)));
zdat1p(v)= fb_am_out2(vloop, 1 + int(wstep1 *
(v - 1)));
xziku1(v)= sp1 + int(wstep1 * (v - 1));
if vloop==1 then
ydat1(v)= wdat1( sp1 + int(wstep1 * (v -
1)));
end
end
else
for v=1:(width0+1)
if v <= wsize3 then
zdat1(v)=fb_out2(vloop,v);
zdat1p(v)=fb_am_out2(vloop,v);
if vloop==1 then
ydat1(v)= wdat1((sp1-1)+v);
end
else
zdat1(v)=0.;
zdat1p(v)=0.;
if vloop==1 then
ydat1(v)= 0.;
end
end
xziku1(v)=(sp1-1)+v;
end
end
sec2x=2;
if vloop==1 then
subplot((n_stepf+1),bun0, sec2x );
plot(xziku1,ydat1,'b');
// wstr1=sprintf('waveform selected,
original');
wstr1='waveform selected(section 2), original
: '+ wavfile1;
xtitle(wstr1);
end
subplot((n_stepf+1),bun0,(retu0 * vloop + sec2x));
plot(xziku1,zdat1,'b');
plot(xziku1,zdat1p,'g');
wstr1=sprintf('Band Pass filtered waveform
selected(section 2) (%dHz-
%dHz)',int(bpf_fb0(vloop,2)),int(bpf_fb0(vloop,3)));
mv0=fb_am_out_prop2(vloop,1);
if mv0 > 0. then
mv1= int(20. * log10(mv0));
else
mv1=-120;
end
pk0=fb_am_out_prop2(vloop,2);
if pk0 > 0. then
pk1= int( 20. * log10(pk0));
else
pk1=-120;
end
wstr2=sprintf('%d/%d dB',mv1,pk1); // average / peak
[dB]
xtitle( wstr1,'',wstr2 );
end // for vloop=1:n_stepf
end // if sec2_on_off == 1 then
//
end // if disp_code >= 1 then
xset('window',w0); // push old windows
endfunction
//-----------------------------------------------------------
function my_eventhandler(win,x,y,ibut)
global eh_var1
global eh_var2
global eh_var3
global eh_var4
global eh_var5
global eh_fm_disp0
if ibut==-1000 then return,end
[x,y]=xchange(x,y,'i2f')
xinfo(msprintf('Event code %d at
mouse position is (%f,%f)',ibut,x,y))
if ibut== 100 // d
if (eh_var2 + eh_var3 - 1)
< n_stepf then
eh_var2=eh_var2+1
end
plot_fb2(eh_var1, eh_var2,
eh_var3,eh_fm_disp0)
elseif ibut== 117 // u
if eh_var2 > 0 then
eh_var2=eh_var2-1
end
plot_fb2(eh_var1, eh_var2,
eh_var3,eh_fm_disp0)
elseif ibut== 102 // f
// set if fm put display or
not
if eh_fm_disp0 <= 0 then
eh_fm_disp0 = 1
disp('+enable FM out display ')
else
eh_fm_disp0 = 0
disp('-disable FM out display ');
end
plot_fb2(eh_var1, eh_var2,
eh_var3,eh_fm_disp0)
elseif ibut==3 // Left mouse button has
been clicked
eh_var4=x
elseif ibut==5 // Right mouse button has
been clicked
eh_var5=x
if eh_var5 > eh_var4 then
fs2= fs1 /
(eh_var5 - eh_var4)
if
sec2_on_off == 1 then
// x-y position is based by latest subplot area
// so, when mouse isnot in latest subplot area, x-y value is different
real one.
else
xinfo(msprintf('Event code %d at mouse position is (%f,%f) x-freq %f
Hz',ibut,x,y,fs2))
end
end
end
// seteventhandler('my_eventhandler')
// seteventhandler('') //suppress the event handler
endfunction
//-----------------------------------------------------------
function plot_fb2(disp_code, nstartline, nlines, fm_disp0)
w0=xget('window'); // stack old window
if disp_code >= 1 then
xset('window',disp_code); // create new windows
clf();
if sec2_on_off == 1 then
bun0=2;
retu0=2;
else
bun0=1;
retu0=1;
end
wsize2=size1 - sp1 + 1;
wsize3=ep1-sp1+1;
zdata1=zeros(1,width0+1);
zdata1p=zeros(1,width0+1);
ydata1=zeros(1,width0+1);
xziku1=zeros(1,width0+1);
if fm_disp0 > 0 then
zdata1fm=zeros(1,width0+1);
nlinesx=nlines * 2; // extend nlines
retu1=2;
else // no display fm out
nlinesx=nlines;
retu1=1;
end
for vloop=0:n_stepf
if ( vloop >= nstartline) & ( vloop <= (nstartline +
nlines -1)) then
if wsize3 > (width0+1) then
wstep1= (ep1-sp1) / width0;
for v=1:(width0+1)
if vloop==0 then
ydat1(v)= wdat1( sp1 + int(wstep1
* (v - 1)));
xziku1(v)= sp1 + int(wstep1 * (v -
1));
else
zdat1(v)= fb_out(vloop,1 +
int(wstep1 * (v - 1)));
zdat1p(v)= fb_am_out(vloop, 1 +
int(wstep1 * (v - 1)));
xziku1(v)= sp1 + int(wstep1 * (v -
1));
if fm_disp0 > 0 then
zdat1fm(v)=
fb_fm_out(vloop, 1 + int(wstep1 * (v - 1)));
end
end
end
else
for v=1:(width0+1)
if v <= wsize3 then
if vloop==0 then
ydat1(v)=
wdat1((sp1-1)+v);
xziku1(v)=(sp1-1)+v;
else
zdat1(v)=fb_out(vloop,v);
zdat1p(v)=fb_am_out(vloop,v);
xziku1(v)=(sp1-1)+v;
if fm_disp0 > 0 then
zdat1fm(v)= fb_fm_out(vloop, v);
end
end
else
if vloop==0 then
ydat1(v)= 0.;
xziku1(v)=(sp1-1)+v;
else
zdat1(v)=0.;
zdat1p(v)=0.;
xziku1(v)=(sp1-1)+v;
if fm_disp0 > 0 then
zdat1fm(v)= 0.;
end
end
end
end //
end //if wsize3 > (width0+1) then
//
if vloop==0 then // *if_123
else
if ((fm_disp0 > 0) & ( fb_fm_out_prop(vloop,2) >
0 )) then // if there is one more
spx1=0;
spx2=0;
for vv=1:(width0+1)
if xziku1(vv) <
fb_fm_out_prop(vloop,1) then
spx1=vv;
end
if xziku1(vv) <
(fb_fm_out_prop(vloop,2) + 1) then
spx2=vv;
end
end
lng0=spx2 - spx1;
if lng0 > 0 then
zdat1fm2=zeros(1,lng0);
ziku2=zeros(1,lng0);
for vv=spx1:spx2
zdat1fm2(vv-spx1+1)=zdat1fm(vv);
ziku2(vv-spx1+1)=xziku1(vv);
end
end
end
end // *if_123
//
if vloop==0 then
subplot(nlinesx ,bun0,1);
plot(xziku1,ydat1,'b');
// wstr1=sprintf('waveform selected,
original');
wstr1='waveform selected(section 1), original
: '+ wavfile1;
xtitle(wstr1);
else
subplot( nlinesx, bun0,(retu1 * retu0 * (vloop -
nstartline) + 1));
plot(xziku1,zdat1,'b');
plot(xziku1,zdat1p,'g');
wstr1=sprintf('Band Pass filtered waveform
selected(section 1) (%dHz-
%dHz)',int(bpf_fb0(vloop,2)),int(bpf_fb0(vloop,3)));
mv0=fb_am_out_prop(vloop,1);
if mv0 > 0. then
mv1= int(20. * log10(mv0));
else
mv1=-120;
end
pk0=fb_am_out_prop(vloop,2);
if pk0 > 0. then
pk1= int( 20. * log10(pk0));
else
pk1=-120;
end
wstr2=sprintf('%d/%d dB',mv1,pk1); // average
/ peak [dB]
xtitle( wstr1,'',wstr2 );
end
//+++ display fm out
if vloop == 0 then // *if_124
else
if ((fm_disp0 > 0) & ( fb_fm_out_prop(vloop,2) >
0 )) then
if sec2_on_off == 1 then
subplot( nlinesx, bun0,(retu1 *
retu0 * (vloop - nstartline) + 1+ retu1));
else
subplot( nlinesx, bun0,(retu1 *
retu0 * (vloop - nstartline) + 2));
end
plot(xziku1,zdat1fm,'w');
if lng0 > 0 then // lng is defined
above
plot(ziku2,zdat1fm2,'r');
end
wstr2=sprintf('cycle [Hz]'); // write
unit
xtitle( wstr1,'',wstr2 );
end
end // *if_124
//--- display fm out
end //if ( vloop >= nstartline) & ( vloop <=
(nstartline + nlines -1)) then
end // for vloop=0:n_stepf
//
// section 2
if sec2_on_off == 1 then
sp1=sp2;
ep1=ep2;
wsize2=size1 - sp1 + 1;
wsize3=ep1-sp1+1;
zdata1=zeros(1,width0+1);
zdata1p=zeros(1,width0+1);
ydata1=zeros(1,width0+1);
for vloop=0:n_stepf
if ( vloop >= nstartline) & ( vloop <= (nstartline +
nlines -1)) then
if wsize3 > (width0+1) then
wstep1= (ep1-sp1) / width0;
for v=1:(width0+1)
if vloop==0 then
ydat1(v)= wdat1( sp1 +
int(wstep1 * (v - 1)));
xziku1(v)= sp1 + int(wstep1
* (v - 1));
else
zdat1(v)= fb_out2(vloop,1 +
int(wstep1 * (v - 1)));
zdat1p(v)= fb_am_out2(vloop, 1 +
int(wstep1 * (v - 1)));
xziku1(v)= sp1 + int(wstep1 * (v -
1));
if fm_disp0 > 0 then
zdat1fm(v)=
fb_fm_out2(vloop, 1 + int(wstep1 * (v - 1)));
end
end
end
else
for v=1:(width0+1)
if v <= wsize3 then
if vloop==0 then
ydat1(v)=
wdat1((sp1-1)+v);
xziku1(v)=(sp1-1)+v;
else
zdat1(v)=fb_out2(vloop,v);
zdat1p(v)=fb_am_out2(vloop,v);
xziku1(v)=(sp1-1)+v;
if
fm_disp0 > 0 then
zdat1fm(v)= fb_fm_out2(vloop, v);
end
end
else
if vloop==0 then
ydat1(v)=
0.;
xziku1(v)=(sp1-1)+v;
else
zdat1(v)=0.;
zdat1p(v)=0.;
xziku1(v)=(sp1-1)+v;
if
fm_disp0 > 0 then
zdat1fm(v)= 0.;
end
end
end
end
end
//
if vloop==0 then // *if_125
else
if ((fm_disp0 > 0) & ( fb_fm_out_prop2(vloop,2)
> 0 )) then // if there is one more
spx1=0;
spx2=0;
for vv=1:(width0+1)
if xziku1(vv) <
fb_fm_out_prop2(vloop,1) then
spx1=vv;
end
if xziku1(vv) <
(fb_fm_out_prop2(vloop,2) + 1) then
spx2=vv;
end
end
lng0=spx2 - spx1;
if lng0 > 0 then
zdat1fm2=zeros(1,lng0);
ziku2=zeros(1,lng0);
for vv=spx1:spx2
zdat1fm2(vv-spx1+1)=zdat1fm(vv);
ziku2(vv-spx1+1)=xziku1(vv);
end
end
end
end // *if_125
//
sec2x=2;
if vloop==0 then
subplot(nlinesx ,bun0, sec2x );
plot(xziku1,ydat1,'b');
// wstr1=sprintf('waveform selected,
original');
wstr1='waveform selected(section 2), original
: '+ wavfile1;
xtitle(wstr1);
else
subplot( nlinesx ,bun0,(retu1 * retu0 * (vloop -
nstartline ) + sec2x));
plot(xziku1,zdat1,'b');
plot(xziku1,zdat1p,'g');
wstr1=sprintf('Band Pass filtered waveform
selected(section 2) (%dHz-
%dHz)',int(bpf_fb0(vloop,2)),int(bpf_fb0(vloop,3)));
mv0=fb_am_out_prop2(vloop,1);
if mv0 > 0. then
mv1= int(20. * log10(mv0));
else
mv1=-120;
end
pk0=fb_am_out_prop2(vloop,2);
if pk0 > 0. then
pk1= int( 20. * log10(pk0));
else
pk1=-120;
end
wstr2=sprintf('%d/%d dB',mv1,pk1); // average
/ peak [dB]
xtitle( wstr1,'',wstr2 );
end
//+++ display fm out
if vloop == 0 then // *if_126
else
if ((fm_disp0 > 0) & ( fb_fm_out_prop2(vloop,2)
> 0 )) then
subplot( nlinesx, bun0,(retu1 * retu0 * (vloop
- nstartline) + sec2x + 2));
plot(xziku1,zdat1fm,'w');
if lng0 > 0 then // lng is defined
above
plot(ziku2,zdat1fm2,'r');
end
wstr2=sprintf('cycle [Hz]'); // write
unit
xtitle( wstr1,'',wstr2 );
end
end // *if_126
//--- display fm out
end // if ( vloop >= nstartline) & ( vloop <=
(nstartline + nlines -1)) then
end // for vloop=1:n_stepf
end // if sec2_on_off == 1 then
//
end // if disp_code >= 1 then
xset('window',w0); // push old windows
endfunction
//-----------------------------------------------------------
function windowsize0(width0,height0)
dim0=xget("wdim");
disp(dim0);
xset("wdim",width0,height0); // real width and real height
is limited in physical display pixes
endfunction
function showwindowsize0()
dim0=xget("wdim");
disp(dim0);
endfunction
//-----------------------------------------------------------
function plot3d_fb_am_out(disp0)
//
wb0=xget('window'); // stack old window
if disp0 >= 1 then
xset('window',disp0); // create new windows
clf();
w0=ep1-sp1+1;
x1=zeros(n_stepf,1);
bai0=10;
for v=1:n_stepf
x1(v)=bpf_fb0(v,1) / bai0; //
value shoud be near w0
end
x1s=bpf_fb0(1,2) / bai0; // lower freq
x1e=bpf_fb0(n_stepf,3) / bai0; // higher freq
y1=zeros(w0,1);
for v=1:w0
y1(v)=sp1+(v-1);
end
pk0=fb_am_out_prop(1,2);
for v=1:n_stepf
if fb_am_out_prop(v,2) > pk0 then
pk0=fb_am_out_prop(v,2);
end
end
if pk0 <= 0. then
gain0=1000.;
else
gain0=1000./pk0;
end
fb_am_out2=zeros(n_stepf,w0);
for v=1:n_stepf
for v2=1:w0
fb_am_out2(v,v2)=fb_am_out(v,v2) * gain0;
end
end
clf();
plot3d(x1,y1, fb_am_out2 ,theta=130, alpha=75,
leg="x10Hz@point@1000",flag=[5,1,4],ebox=[x1s,x1e,y1(1),y1(w0),0.,1000.0]);
xtitle('filter bank am detect output');
end // if disp0 >= 1 then
xset('window',wb0); // push old windows
endfunction
//-----------------------------------------------------------
function [fb_out,fb_am_out,fb_am_out_prop]=do_fb(disp0,sec_no0,sp1,ep1)
//
w0=(ep1-sp1)+1;
fb_out=zeros(n_stepf,w0);
dt0= int((fir_steps-1) / 2);
ep1x=ep1+dt0;
if ep1x > size1 then
ep1x = size1;
end
sp1x=sp1+dt0;
if sp1x < fir_steps then
sp1x = fir_steps;
end
//
if disp0 >= 1 then
wstr1=sprintf('+section no.%d',sec_no0);
disp(wstr1);
end
//+filtering+
for vloop=1:n_stepf
if disp0 >= 1 then
wstr1=sprintf('+entering filter bank
no.%d/%d',vloop,n_stepf);
disp(wstr1);
end
for loop=sp1x:ep1x
ydz = 0.;
for loop2=1:fir_steps
ydz = ydz + wft(vloop,loop2) *
wdat1( 1 - loop2 + loop );
end
//fir_out1(loop - dt0)=ydz;
fb_out(vloop,loop - dt0 -sp1 + 1)=ydz;
end
end // if vloop=1:n_stepf
//-filtering-
//
//+am+
nsize_max=5000;
g_work2=zeros(10,nsize_max);
fb_am_out=zeros(n_stepf,w0);
fb_am_out_prop=zeros(n_stepf,nmax_am_prop);
fir_out1=zeros(1,size1);
fir_out1p=zeros(1,size1);
for vloop=1:n_stepf
if disp0 >= 1 then
wstr1=sprintf('+entering AM de-modulation, peak
trace no.%d/%d',vloop,n_stepf);
disp(wstr1);
end
for v2=sp1:ep1
fir_out1(v2)= fb_out(vloop,v2-sp1+1);
fir_out1p(v2)=0.;
end
//+from trace peak
//+
spx=-1;
spy=-1;
epx=-1;
epy=-1;
c0=-1;
w0=-1.;
area0=-1.;
if (ep1-sp1) > 4 then
for v=(sp1+1):(ep1-1)
if fir_out1(v) > 0. then // only
postive side
if (fir_out1(v) > fir_out1(v-1))
& (fir_out1(v) > fir_out1(v+1)) then
spx=epx;
spy=epy;
epx=v;
epy=fir_out1(v);
c0=c0+1;
//
if g_work2(1,1) < nsize_max then
if g_work2(1,1) <= 0 then
g_work2(1,1)=1; // first start
g_work2(3,1)=epx;
g_work2(4,1)=epy;
g_work2(5,1)=epx;
g_work2(6,1)=epy;
else
d0=int(g_work2(1,1));
if epy >= g_work2(6,d0) then // rising or equal
g_work2(2,d0)=0; //current
g_work2(5,d0)=epx;
g_work2(6,d0)=epy;
else
// falling
if g_work2(2,d0)==0 then // complete data
g_work2(2,d0)=1; //complete
g_work2(7,d0)= g_work2(5,d0) - g_work2(3,d0);
g_work2(8,d0)= g_work2(6,d0) - g_work2(4,d0);
g_work2(9,d0)= g_work2(7,d0) * g_work2(8,d0) / 2.0 ;
if (g_work2(7,d0) > 0.) & (g_work2(8,d0) > 0.) then
if g_work2(9,d0) > area0 then
area0=g_work2(9,d0);
w0=g_work2(8,d0) / g_work2(7,d0);
end
end
// atan
//g_work2(10,d0)= atan( (g_work2(8,d0) / g_work2(7,d0))) / 2.0 / %pi *
360.0;
//
g_work2(1,1)=g_work2(1,1)+1; // set as new start point
d1=int(g_work2(1,1));
g_work2(3,d1)=epx;
g_work2(4,d1)=epy;
g_work2(5,d1)=epx;
g_work2(6,d1)=epy;
elseif g_work2(2,d0) <= -1 then // restart
g_work2(3,d0)=epx;
g_work2(4,d0)=epy;
g_work2(5,d0)=epx;
g_work2(6,d0)=epy;
end
end
end // if g_work2(1,1) <= 0 then
else
disp(' error: g_g_work2(1,1) >= nsize_max , in do_fb() ');
end // if g_work2(1,1) < nsize_max then
//
//
if c0 >= 1 then
a0=
(epy - spy) / (epx - spx);
b0=
epy - a0 * epx;
for
v2=spx:epx
fir_out1p(v2)= a0 * v2 + b0;
end
end // if
c0 >= 1 then
//
end // if (fir_out1(v) >
fir_out1(v-1)) & (fir_out1(v) > fir_out1(v+ 1)) then
end // if fir_out1(v) > 0. then
end //v=(sp1+1):(ep1-1)
end
if c0 < 1 then
epx=sp1-1;
end
for v=(epx+1):ep1
fir_out1p(v)=0;
end
pk0=0.;
mv0=0.;
for loop=sp1:ep1
fb_am_out(vloop,loop - sp1 + 1)= fir_out1p(loop);
if fir_out1p(loop) > pk0 then
pk0=fir_out1p(loop);
end
mv0=mv0+fir_out1p(loop);
end
fb_am_out_prop(vloop,1)= mv0 / (ep1-sp1+1);
fb_am_out_prop(vloop,2)= pk0;
//+
//-from trace peak
end // for vloop=1:n_stepf
//-am-
endfunction
//-----------------------------------------------------------
function
[fb_outn,fb_am_outn,fb_am_out_propn]=do_fbn(disp0,sec_no0,sp1,ep1,n_stepfn,wftn)
//
n_stepf=n_stepfn; // set n_stepfn as n_stepf
//
w0=(ep1-sp1)+1;
fb_outn=zeros(n_stepf,w0);
dt0= int((fir_steps-1) / 2);
ep1x=ep1+dt0;
if ep1x > size1 then
ep1x = size1;
end
sp1x=sp1+dt0;
if sp1x < fir_steps then
sp1x = fir_steps;
end
//
if disp0 >= 1 then
wstr1=sprintf('+section no.%d',sec_no0);
disp(wstr1);
end
//+filtering+
for vloop=1:n_stepf
cmap=fb_ctable(vloop+1); // get index
if cmap > 0 then
if disp0 >= 1 then
wstr1=sprintf('+copying old filter bank data
no.%d/%d',vloop,n_stepf);
disp(wstr1);
end
for loop=1:w0
if sec_no0 == 2 then
fb_outn(vloop,loop)=fb_out2(cmap,loop);
else
fb_outn(vloop,loop)=fb_out(cmap,loop);
end
end
else
if disp0 >= 1 then
wstr1=sprintf('+entering new filter bank
no.%d/%d',vloop,n_stepf);
disp(wstr1);
end
for loop=sp1x:ep1x
ydz = 0.;
for loop2=1:fir_steps
ydz = ydz + wftn(vloop,loop2) *
wdat1( 1 - loop2 + loop );
end
//fir_out1(loop - dt0)=ydz;
fb_outn(vloop,loop - dt0 -sp1 + 1)=ydz;
end
end // if cmap > 0 then else
end // if vloop=1:n_stepf
//-filtering-
//
//+am+
nsize_max=5000;
g_work2=zeros(10,nsize_max);
fb_am_outn=zeros(n_stepf,w0);
fb_am_out_propn=zeros(n_stepf,nmax_am_prop);
fir_ou1=zeros(1,size1);
fir_out1p=zeros(1,size1);
w0s=w0; // w0 will be use as other
for vloop=1:n_stepf
cmap=fb_ctable(vloop+1); // get index
if cmap > 0 then
if disp0 >= 1 then
wstr1=sprintf('+copying old AM de-modulation, peak
trace data no.%d/%d',vloop,n_stepf);
disp(wstr1);
end
// wstr1=sprintf('s1 %d s2 %d w0 %d',sp1,ep1,w0);
// disp(wstr1);
for loop=1:w0s
if sec_no0 == 2 then
fb_am_outn(vloop,loop)= fb_am_out2(cmap,loop);
else
fb_am_outn(vloop,loop)= fb_am_out(cmap,loop);
end
end
if sec_no0 == 2 then
fb_am_out_propn(vloop,1)=
fb_am_out_prop2(cmap,1);
fb_am_out_propn(vloop,2)=
fb_am_out_prop2(cmap,2);
else
fb_am_out_propn(vloop,1)=
fb_am_out_prop(cmap,1);
fb_am_out_propn(vloop,2)=
fb_am_out_prop(cmap,2);
end
else
if disp0 >= 1 then
wstr1=sprintf('+entering new AM de-modulation, peak
trace no.%d/%d',vloop,n_stepf);
disp(wstr1);
end
for v2=sp1:ep1
fir_out1(v2)= fb_outn(vloop,v2-sp1+1);
fir_out1p(v2)=0.;
end
//+from trace peak
//+
spx=-1;
spy=-1;
epx=-1;
epy=-1;
c0=-1;
w0=-1.;
area0=-1.;
if (ep1-sp1) > 4 then
for v=(sp1+1):(ep1-1)
if fir_out1(v) > 0. then // only
postive side
if (fir_out1(v) > fir_out1(v-1))
& (fir_out1(v) > fir_out1(v+1)) then
spx=epx;
spy=epy;
epx=v;
epy=fir_out1(v);
c0=c0+1;
//
if g_work2(1,1) < nsize_max then
if g_work2(1,1) <= 0 then
g_work2(1,1)=1; // first start
g_work2(3,1)=epx;
g_work2(4,1)=epy;
g_work2(5,1)=epx;
g_work2(6,1)=epy;
else
d0=int(g_work2(1,1));
if epy >= g_work2(6,d0) then // rising or equal
g_work2(2,d0)=0; //current
g_work2(5,d0)=epx;
g_work2(6,d0)=epy;
else
// falling
if g_work2(2,d0)==0 then // complete data
g_work2(2,d0)=1; //complete
g_work2(7,d0)= g_work2(5,d0) - g_work2(3,d0);
g_work2(8,d0)= g_work2(6,d0) - g_work2(4,d0);
g_work2(9,d0)= g_work2(7,d0) * g_work2(8,d0) / 2.0 ;
if (g_work2(7,d0) > 0.) & (g_work2(8,d0) > 0.) then
if g_work2(9,d0) > area0 then
area0=g_work2(9,d0);
w0=g_work2(8,d0) / g_work2(7,d0);
end
end
// atan
//g_work2(10,d0)= atan( (g_work2(8,d0) / g_work2(7,d0))) / 2.0 / %pi *
360.0;
//
g_work2(1,1)=g_work2(1,1)+1; // set as new start point
d1=int(g_work2(1,1));
g_work2(3,d1)=epx;
g_work2(4,d1)=epy;
g_work2(5,d1)=epx;
g_work2(6,d1)=epy;
elseif g_work2(2,d0) <= -1 then // restart
g_work2(3,d0)=epx;
g_work2(4,d0)=epy;
g_work2(5,d0)=epx;
g_work2(6,d0)=epy;
end
end
end // if g_work2(1,1) <= 0 then
else
disp(' error: g_g_work2(1,1) >= nsize_max , in do_fbn() ');
end // if g_work2(1,1) < nsize_max then
//
//
if c0 >= 1 then
a0=
(epy - spy) / (epx - spx);
b0=
epy - a0 * epx;
for
v2=spx:epx
fir_out1p(v2)= a0 * v2 + b0;
end
end // if
c0 >= 1 then
//
end // if (fir_out1(v) >
fir_out1(v-1)) & (fir_out1(v) > fir_out1(v+ 1)) then
end // if fir_out1(v) > 0. then
end //v=(sp1+1):(ep1-1)
end
if c0 < 1 then
epx=sp1-1;
end
for v=(epx+1):ep1
fir_out1p(v)=0;
end
pk0=0.;
mv0=0.;
for loop=sp1:ep1
fb_am_outn(vloop,loop - sp1 + 1)= fir_out1p(loop);
if fir_out1p(loop) > pk0 then
pk0=fir_out1p(loop);
end
mv0=mv0+fir_out1p(loop);
end
fb_am_out_propn(vloop,1)= mv0 / (ep1-sp1+1);
fb_am_out_propn(vloop,2)= pk0;
//+
//-from trace peak
end // if cmap > 0 then else
end // for vloop=1:n_stepf
//-am-
endfunction
//-----------------------------------------------------------------
function [n_stepfn, rslt, wftn , fb_ctable]=replace_set_fb(disp0)
// copy old bpf setting
w_bpf_fb0=zeros(n_stepf,3);
for v=1:n_stepf
w_bpf_fb0(v,1)=v; // replace center freq
to no.
w_bpf_fb0(v,2)=bpf_fb0(v,2);
w_bpf_fb0(v,3)=bpf_fb0(v,3);
end
//
[rslt]=x_matrix('set bpf no, low freq, high freq (if value is
minus,it will be deleted)' , w_bpf_fb0);
if rslt==[] then // set old bpf setting
for v=1:n_stepf
for w=1:3
rslt(v,w)=bpf_fb0(v,w);
end
end
n_stepfn=n_stepf;
else // back to center freq
n_stepfn=0;
for v=1:n_stepf
if ((rslt(v,1) > 0) & (rslt(v,2) > 0)
& (rslt(v,3) > 0) ) then
n_stepfn=n_stepfn+1;
w_bpf_fb0(n_stepfn,1)= rslt(v,2) + (
rslt(v,3) - rslt(v,2) ) / 2.0;
w_bpf_fb0(n_stepfn,2)= rslt(v,2);
w_bpf_fb0(n_stepfn,3)= rslt(v,3);
end
end
rslt=zeros(n_stepfn,3);
// for sort by center freq. eq. value (*,1)
w_rslt=zeros(n_stepfn,3);
w_r=zeros(n_stepfn,1);
for v=1:n_stepfn
w_r(v)=-1 * w_bpf_fb0(v,1);
for w=1:3
w_rslt(v,w)=w_bpf_fb0(v,w);
end
end
[w_s,w_k]=sort(w_r);
for v=1:n_stepfn
v2=w_k(v);
for w=1:3
rslt(v,w)=w_rslt(v2,w);
end
end
end
// make map
fb_ctable=zeros(n_stepfn+1,1);
for v=1:n_stepfn
for v2=1:n_stepf
if ( ( abs(rslt(v,2) - bpf_fb0(v2,2))
< 0.01 ) & ( abs(rslt(v,3) - bpf_fb0(v2,3)) <
0.01 ) ) then
fb_ctable(v+1)=v2;
break;
end
end
if fb_ctable(v+1) == 0 then
fb_ctable(1)=1; // there is new
end
end
//
for v=1:n_stepfn
if disp0 >= 1 then
disp( [ v rslt(v,1) rslt(v,2) rslt(v,3)
fb_ctable(v+1)]);
end
end
//
//
// ++ mk_fir_coeff
//
w0=xget('window'); // stack old window
wftn=zeros(n_stepf,fir_steps);
for vloop=1:n_stepfn
cfreq(1)= rslt(vloop,2) / fs1;
cfreq(2)= rslt(vloop,3) / fs1;
fpar(1)=0.1; // dummy data
fpar(2)=-0.1; // dummy data
if kwindows == 2 then
[wft2,wfm2,fr2]=wfir('bp',fir_steps,cfreq,'hn',fpar); // gensui
riritu yuusen
else
[wft2,wfm2,fr2]=wfir('bp',fir_steps,cfreq,'hm',fpar); // bunkai
nou yuusen
end
for v2=1:fir_steps
wftn(vloop,v2)=wft2(v2);
end
if disp0 >=2 then
//++ only one will be done
if vloop==1 then
xset('window',disp0); // create new windows
clf();
[frq1,sfrq1,is1,ie1]=set_frq(0);
fr0=zeros(2); // freq list selected
frq1f=zeros(2); // freq responce selected
phf0=zeros(2); // phase selected
sv0=size(fr2);
end // if vloop==1 then
//-- only one will be done
sc0=0;
for v=1:sv0(2)
fx0= fr2(v) * fs1;
if (fx0 >= frq1(1)) &(fx0 <=
frq1(sfrq1)) then
sc0=sc0+1;
fr0(sc0)=fx0;
frq1f(sc0)= 20. * log10( wfm2(v));
phf0(sc0)= -1.0 * ((fir_steps - 1)
/ 2.0 ) * ( 2. * %pi * v / fir_steps);
end
end // for v=1:sv0(2)
gainplot(fr0,frq1f,phf0);
end // vloop+1:n_stepf
end // if disp0 >=2 then
xset('window',w0); // push old windows
//
// -- mk_fir_coeff
//
endfunction
//-----------------------------------------------------------------
function
[n_stepf,bpf_fb0,wft,fb_out,fb_am_out,fb_am_out_prop,fb_out2,fb_am_out2,fb_am_out_prop2,fb_fm_out,fb_fm_out_prop,fb_fm_out2,fb_fm_out_prop2]=ch_set()
n_stepf=n_stepfn;
bpf_fb0=bpf_fbn;
wft=wftn;
fb_out=fb_outn;
fb_am_out=fb_am_outn;
fb_am_out_prop=fb_am_out_propn;
fb_fm_out=fb_fm_outn;
fb_fm_out_prop=fb_fm_out_propn;
//
fb_out2=fb_out2n;
fb_am_out2=fb_am_out2n;
fb_am_out_prop2=fb_am_out_prop2n;
fb_fm_out2=fb_fm_out2n;
fb_fm_out_prop2=fb_fm_out_prop2n;
endfunction
//-----------------------------------------------------------------
function [wdat2,rslt]=do_replace(disp0)
sz0=size(cfb_list);
cfb_list2=ones(n_stepf,2);
for v=1:n_stepf
cfb_list2(v,1)=v;
if v <= sz0(1) then
cfb_list2(v,2)=cfb_list(v,2);
end
end
[rslt]=x_matrix('Mark 1 as Combine element, in filter bank no.
list(no, 1 or 0 )' , cfb_list2);
if rslt==[] then // reset all 1
rslt=ones(n_stepf,2);
for v=1:n_stepf
rslt(v,1)=v;
end
end
wdat2=zeros(1,size1);
for v=1:size1
wdat2(v)=wdat1(v);
end
//+replacement
for v=sp1:ep1
vs0=0.;
for v2=1:n_stepf
if rslt(v2,2)==1 then // if the fb
no is Combine element
vs0= vs0 + fb_out(v2,
(v-sp1+1))
end
end
wdat2(v)=vs0;
end
//-replacement
if disp0 >= 1 then
ydat1r=zeros(1,(width0+1));
if size1 > (width0+1) then
wstep1= size1 / (width0+1);
for v=1:(width0+1)
ydat1r(v)= wdat2( 1 + int(wstep1 * (v -
1)));
// xziku1(v)= 1 + int(wstep1 * (v - 1));
end
else
for v=1:(width0+1)
if v <= size1 then
ydat1r(v)=wdat2(v);
else
ydat1r(v)=0.;
end
// xziku1(v)=(sp1-1)+v;
end
end
//
clf();
subplot(211);
plot(xziku0,ydat0,'b');
//plot2d(xziku0,ydat0,style=[color("blue")]);
xtitle('waveform all');
// When linux scilab-3.1 2nd subplot does not work well.
subplot(212);
plot(xziku0,ydat1r,'b');
//plot2d(xziku1,ydat1,style=[color("blue")]);
xtitle('replacement waveform');
end // id disp0 >= 1 thne
endfunction
//=================================================================
//
// +MAIN (2) program starts
// filter bank
//
if T_DEMO==1 then
if T_N_STEPF==1 then
n_stepf=input("-> How many filters of filter bank ?");
end
// [bpf_fb0,wft]=make_bpf_fb0(3);
end
if T_DEMO==1 then
// [fb_out,fb_am_out,fb_am_out_prop]=do_fb(1,1,sp1,ep1);
// [fb_out2,fb_am_out2,fb_am_out_prop2]=do_fb(1,2,sp2,ep2);
// plot_fb(2);
[sp1,ep1,n_stepf,bpf_fb0,fb_out,fb_am_out,fb_am_out_prop,sp2,ep2,sec2_on_off,fb_out2,fb_am_out2,fb_am_out_prop2,st_freq,end_freq,n_stepf,fir_steps,kwindows,fb_fm_out,fb_fm_out_prop,fb_fm_out2,fb_fm_out_prop2,wft,bpf_fb0]=load3();
plot_fb(2);
[ydat1,xziku1]=make_width_data( wdat1, size1,sp1,ep1);
[ydat2,xziku2]=make_width_data( wdat1, size1,sp2,ep2);
plot_wave1();
end
//
// -MAIN (2)program starts
//
//==============================================================++=+
//
function
[fb_fm_outx,fb_fm_out_propx]=cycle_detect(disp0,sec_no0,sp1,ep1,n_stepfn,force_do,xswitch0)
// if force_do is 1, calculate all
// if xswitch0 is 1, use fb_out2n,fb_out2n. Other fb_out, fb_out2
n_stepf=n_stepfn; // set n_stepfn as n_stepf
nsize_max=5000;
g_work2=zeros(10,nsize_max);
fir_out1=zeros(1,size1);
w0=(ep1-sp1)+1;
w0s=w0; // w0 will be use as other
fb_fm_outx=zeros(n_stepf,w0);
fb_fm_out_propx=zeros(n_stepf,nmax_am_prop);
//
if disp0 >= 1 then
wstr1=sprintf('+section no.%d',sec_no0);
disp(wstr1);
end
//
for vloop=1:n_stepf
if ((fb_ctable(1) < 0) | (force_do == 1)) then // this is not
defined yet
cmap=0; // do cycle_detect
else
cmap=fb_ctable(vloop+1); // get index
end
if cmap > 0 then
if disp0 >= 1 then
wstr1=sprintf('+copying old FM de-modulation, peak
trace data no.%d/%d',vloop,n_stepf);
disp(wstr1);
end
// wstr1=sprintf('s1 %d s2 %d w0 %d',sp1,ep1,w0);
// disp(wstr1);
for loop=1:w0s
if sec_no0 == 2 then
fb_fm_outx(vloop,loop)= fb_fm_out2(cmap,loop);
else
fb_fm_outx(vloop,loop)= fb_fm_out(cmap,loop);
end
end
if sec_no0 == 2 then
fb_fm_out_propx(vloop,1)=
fb_fm_out_prop2(cmap,1);
fb_fm_out_propx(vloop,2)=
fb_fm_out_prop2(cmap,2);
fb_fm_out_propx(vloop,3)=
fb_fm_out_prop2(cmap,3);
fb_fm_out_propx(vloop,4)=
fb_fm_out_prop2(cmap,4);
else
fb_fm_out_propx(vloop,1)=
fb_fm_out_prop(cmap,1);
fb_fm_out_propx(vloop,2)=
fb_fm_out_prop(cmap,2);
fb_fm_out_propx(vloop,3)=
fb_fm_out_prop(cmap,3);
fb_fm_out_propx(vloop,4)=
fb_fm_out_prop(cmap,4);
end
else
if disp0 >= 1 then
wstr1=sprintf('+entering FM de-modulation, peak
trace no.%d/%d',vloop,n_stepf);
disp(wstr1);
end
for v2=sp1:ep1
if sec_no0 == 2 then
if xswitch0 == 1 then
fir_out1(v2)=
fb_out2n(vloop,v2-sp1+1);
else
fir_out1(v2)=
fb_out2(vloop,v2-sp1+1);
end
else
if xswitch0 == 1 then
fir_out1(v2)=
fb_outn(vloop,v2-sp1+1);
else
fir_out1(v2)=
fb_out(vloop,v2-sp1+1);
end
end
end
//+from trace peak
//+
spx=-1;
spy=-1;
epx=-1;
epy=-1;
c0=-1;
w0=-1.;
area0=-1.;
//
c1=0; // counter
cf=0; // 1st data flag
fb_fm_out_propx(vloop,1)=-1.; // start point
fb_fm_out_propx(vloop,2)=-1.; // end point
fb_fm_out_propx(vloop,3)=0.0; // max
fb_fm_out_propx(vloop,4)=fs1; // min
spa=-1.0;
spb=-1.0;
vala=-1.0;
valb=-1.0;
//
//
if (ep1-sp1) > 4 then
for v=(sp1+1):(ep1-1)
if fir_out1(v) > 0. then // only
postive side
if (fir_out1(v) > fir_out1(v-1))
& (fir_out1(v) > fir_out1(v+1)) then
spx=epx;
spy=epy;
epx=v;
epy=fir_out1(v);
c0=c0+1;
//
if g_work2(1,1) < nsize_max then
if g_work2(1,1) <= 0 then
g_work2(1,1)=1; // first start
g_work2(3,1)=epx;
g_work2(4,1)=epy;
g_work2(5,1)=epx;
g_work2(6,1)=epy;
else
d0=int(g_work2(1,1));
if epy >= g_work2(6,d0) then // rising or equal
g_work2(2,d0)=0; //current
g_work2(5,d0)=epx;
g_work2(6,d0)=epy;
else
// falling
if g_work2(2,d0)==0 then // complete data
g_work2(2,d0)=1; //complete
g_work2(7,d0)= g_work2(5,d0) - g_work2(3,d0);
g_work2(8,d0)= g_work2(6,d0) - g_work2(4,d0);
g_work2(9,d0)= g_work2(7,d0) * g_work2(8,d0) / 2.0 ;
if (g_work2(7,d0) > 0.) & (g_work2(8,d0) > 0.) then
if g_work2(9,d0) > area0 then
area0=g_work2(9,d0);
w0=g_work2(8,d0) / g_work2(7,d0);
end
end
// atan
//g_work2(10,d0)= atan( (g_work2(8,d0) / g_work2(7,d0))) / 2.0 / %pi *
360.0;
//
g_work2(1,1)=g_work2(1,1)+1; // set as new start point
d1=int(g_work2(1,1));
g_work2(3,d1)=epx;
g_work2(4,d1)=epy;
g_work2(5,d1)=epx;
g_work2(6,d1)=epy;
elseif g_work2(2,d0) <= -1 then // restart
g_work2(3,d0)=epx;
g_work2(4,d0)=epy;
g_work2(5,d0)=epx;
g_work2(6,d0)=epy;
end
end
end // if g_work2(1,1) <= 0 then
else
disp(' error: g_g_work2(1,1) >= nsize_max , in cycle_detect ');
end // if g_work2(1,1) < nsize_max then
//
//
if c0 >= 1 then
if disp0 >= 2 then
disp([ (epx - spx) spx epx ]);
end
// check 1st data is miss-peak-detection ?
if ((cf == 0) & (fir_out1(spx - 1) == 0.)) then
cf=-1; // flag de-set
if disp0 >= 1 then
disp('waring: maybe 1st data is miss peak detection, ignore one
');
end
else //+ 1st data OK or others
c1=c1+1; // find one
spa=spb; // push
vala=valb; // push
spb= int((epx - spx) /2) + spx; // start/end point is defined as
middle point
if ( epx - spx ) <= 0.01 then
disp('error (epx - spx) <= -1 in cycle_detect ');
valb= fs1; // this is error !
else
valb= fs1 / ( epx - spx );
end
end // if ((c1 == 0) & (fir_out1(spx - 1) == 0.)) then
//
//
if
c1 >= 1 then
if valb > fb_fm_out_propx(vloop,3) then
fb_fm_out_propx(vloop,3) = valb; // set max
end
if valb < fb_fm_out_propx(vloop,4) then
fb_fm_out_propx(vloop,4) = valb; // set min
end
if c1 == 1 then
fb_fm_out_propx(vloop,1) = spb; // set start point
end
fb_fm_out_propx(vloop,2) = spb; // set end point
if c1 >= 2 then // draw line between middle points
a0= (valb - vala) / (spb - spa);
b0= vala - a0 * spa;
for v2=spa:spb
fb_fm_outx(vloop,v2-sp1+1)= a0 * v2 + b0;
end
end //c1 >= 2 then
end
// if c1 >= 1 then
end // if c0
>= 1 then
//
end // if (fir_out1(v) >
fir_out1(v-1)) & (fir_out1(v) > fir_out1(v+ 1)) then
end // if fir_out1(v) > 0. then
end //v=(sp1+1):(ep1-1)
end //if (ep1-sp1) > 4 then
//
if c1 >= 2 then
for vv=sp1:(fb_fm_out_propx(vloop,1)-1) // draw same
line rest before
fb_fm_outx(vloop, vv-sp1+1) = fb_fm_outx(
vloop, fb_fm_out_propx(vloop,1) - sp1+1 );
end
for vv=(fb_fm_out_propx(vloop,2)+1):ep1 // draw line
rest after
fb_fm_outx(vloop, vv-sp1+1) =
fb_fm_outx(vloop, fb_fm_out_propx(vloop,2) - sp1+1 );
end
end
end // cmap > 0 then
end // for vloop=1:n_stepf
//-am-
endfunction
//
//==================================================================
//
// +MAIN (3)program starts
//
if T_DEMO==1 then
eh_var1=2;
eh_var2=0;
eh_var3=3;
xset('window', eh_var1);
seteventhandler('my_eventhandler');
end
//
// -MAIN (3)program starts
//
//===================================================================
//
// ++ MENU MENU MENU =====================
//
// 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',[ '(1) set_et_plot()'; '
reset_plot()'; '(op1-1) set fft points'; '(op1-2) do fft'; '(op2-1)set
multi fft sp1'; '(op2-2) do ffts' ]);
step1_plotwav=[ '[sp1, ep1, sec2_on_off, sp2, ep2 ]=
set_sp1_ep1(size1); [ydat1,xziku1]=make_width_data( wdat1,
size1,sp1,ep1); [ydat2,xziku2]=make_width_data( wdat1, size1,sp2,ep2);
xset(''window'',0); plot_wave1();' ; '[sp1, ep1, sp2, ep2]=
reset_sp1_ep1(size1); [ydat1,xziku1]=make_width_data( wdat1, size1);
[ydat2,xziku2]=make_width_data( wdat1, size1,sp2,ep2);
plot_wave1();'; '[fft_point1]=set_fft_points();' ;
'[db_fft1,phi_fft1]=do_fft_wav(sp1,ep1);
[db_fft2,phi_fft2]=do_fft_wav(sp2,ep2); plot_fft1(5)' ;'[fft_sp1_list]=
set_sp1s()';'plot_multi_fft1(5)'] ;
end //if ADD_MENU_PLOT == 1 then
//
//-----------------------------------------------------------------
//
ADD_MENU_FB=1; // set ADD_MENU 1 to
add menu button of some functions
//
//
if ADD_MENU_FB == 1 then
delmenu('step2_filter_bank');
addmenu('step2_filter_bank',[ '(1) set fb parameters','(2-1) do filter
bank (sec.1)','(2-2) do filter bank (sec.2)','(2-e)
plot2d_fb_am_out','(3) replacement', '(4) save replacement as a file',
'(-) play replacement sound' ,'(1b)change fb parameters','(2b)do change
filter bank(sec.1)','(2c)do change filter bank(sec.2)','(3e)change
fix']);
step2_filter_bank=[ '[st_freq, end_freq, n_stepf, fir_steps, kwindows]=
set_fb_para();
[bpf_fb0,wft]=make_bpf_fb0(3)','[fb_out,fb_am_out,fb_am_out_prop]=do_fb(1,1,sp1,ep1);
[fb_fm_out,fb_fm_out_prop]=cycle_detect(1,1,sp1,ep1,n_stepf,1,0);',
'[fb_out2,fb_am_out2,fb_am_out_prop2]=do_fb(1,2,sp2,ep2);
[fb_fm_out2,fb_fm_out_prop2]=cycle_detect(1,2,sp2,ep2,n_stepf,1,0);','plot_fb(2)','[wdat2,cfb_list]=do_replace(1)','snd_save2()','snd_play2()','[n_stepfn,bpf_fbn,wftn,fb_ctable]=replace_set_fb(3);','[fb_outn,fb_am_outn,fb_am_out_propn]=do_fbn(1,1,sp1,ep1,n_stepfn,wftn);
[fb_fm_outn,fb_fm_out_propn]=cycle_detect(1,1,sp1,ep1,n_stepfn,0,1);','[fb_out2n,fb_am_out2n,fb_am_out_prop2n]=do_fbn(1,2,sp2,ep2,n_stepfn,wftn);
[fb_fm_out2n,fb_fm_out_prop2n]=cycle_detect(1,2,sp2,ep2,n_stepfn,0,1);','eh_var2=0;
[n_stepf,bpf_fb0,wft,fb_out,fb_am_out,fb_am_out_prop,fb_out2,fb_am_out2,fb_am_out_prop2,fb_fm_out,fb_fm_out_prop,fb_fm_out2,fb_fm_out_prop2]=ch_set();
fb_ctable(1)=-1;'] ;
end //if ADD_MENU_FB == 1 then
//
//-----------------------------------------------------------------
//
ADD_MENU_EH=1; // set ADD_MENU 1 to add menu button of some
functions
//
//
if ADD_MENU_EH == 1 then
delmenu('step2_event_handler');
addmenu('step2_event_handler',['(1) set event handler window(2), key in
d(=down) or u(=up), f(=FM out disp on/off)','(e) suppress handler']);
step2_event_handler=[ 'eh_var1=2; eh_var2=0; eh_var3=3;
xset(''window'', eh_var1); seteventhandler(''my_eventhandler'');'
, 'seteventhandler('''') ' ];
end // if ADD_MENU_EH == 1 then
//
//-----------------------------------------------------------------
//
ADD_MENU_FB2=1; // set ADD_MENU 1
to add menu button of some functions
//
//
if ADD_MENU_FB2 == 1 then
delmenu('step2_filter_bank_backup');
addmenu('step2_filter_bank_backup',[ '(-) save','(-) load' ]);
step2_filter_bank_backup=['save3()','[sp1,ep1,n_stepf,bpf_fb0,fb_out,fb_am_out,fb_am_out_prop,sp2,ep2,sec2_on_off,fb_out2,fb_am_out2,fb_am_out_prop2,st_freq,end_freq,n_stepf,fir_steps,kwindows,fb_fm_out,fb_fm_out_prop,fb_fm_out2,fb_fm_out_prop2,wft,bpf_fb0]=load3();
plot_fb(2);'] ;
end //if ADD_MENU_FB == 1 then
//
//-----------------------------------------------------------------
//
ADD_MENU_OSC=1; // set ADD_MENU 1
to add menu button of some functions
//
//
if ADD_MENU_OSC == 1 then
delmenu('step3_cycle_detect');
addmenu('step3_cycle_detect',[ '(1) cycle detect (sec.1)' , '(2)
cycle detect (sec.2)' ]);
step3_cycle_detect=['[fb_fm_out,fb_fm_out_prop]=cycle_detect(1,1,sp1,ep1,n_stepf,1,0);',
'[fb_fm_out2,fb_fm_out_prop2]=cycle_detect(1,2,sp2,ep2,n_stepf,1,0);'];
end //if ADD_MENU_OSC == 1 then
//
//
//==============================================================++=
//
//
ADD_MENU_STATUS=1;
if ADD_MENU_STATUS == 1 then
function status1()
ver0= getversion();
wstr1=sprintf('scilab version is %s',ver0);
disp(wstr1);
endfunction
function status2()
sz=stacksize();
gsz=gstacksize();
wstr1=sprintf('current stack size is %d of %d ,(%f)',
sz(2),sz(1), sz(2)/sz(1) );
wstr2=sprintf('current global stack size is %d of %d
,(%f)', gsz(2),gsz(1), gsz(2)/gsz(1));
disp( wstr1);
disp( wstr2);
endfunction
delmenu('status_');
addmenu('status_',[ 'scilab version' ; 'check stack memory size' ]);
status_=[ 'status1()' ; 'status2()' ] ;
end //if ADD_MENU_STATUS == 1 then
//
//-------------------------------------------------------------------------------//
以上がa_wavefile_edit_84c.sciです。
No.16 2008年5月11日
