*frequency response of band pass filters of the filter bank
*a filter bank program, This program is just for a reference
//------------------------------------------------------------------------------------------
//
// 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
//
//-------------------------------------------------------------------------------//
*wav file sample data
*another scilab program defines 2 sections on a waveform, calculates
fft frequency response and filter bank output, and compares each on
them. This program is just for a reference.
//------------------------------------------------------------------------------------------
//
// (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
//
//-------------------------------------------------------------------------------//
*wav file sample data