ラフな全探索と最小2乗法を使った構造フィッティングの例
─信号に含まれる２つの構造のパラメーターを推定する─
矩形波に必ず三角波を重畳して発生させる信号源Ａ（ここではこれを構造Ａと呼ぶことにする）と、必ずキャリヤ信号が伴ってＳＩＮ波形に振幅変調したものを
発生させる信号源Ｂ（こちらは構造Ｂと呼ぼう）があるとしてみよう。（具体例は下図を参照） 各構造には、周波数ｆ0、振幅height0、相対比
ratio0の３つの可変パラメーターがあり、これらのパラメーターの値を変えることで、発生する波形の形を変えることができる。
この２種類の構造から発生した信号を混合した信号を分析して、逆に、発生源の構造の各パラメーターを推定する問題を考えてみよう。
下図は、周波数f0が300Ｈｚで、振幅height0が1.2、矩形波と三角波の相対比ratio0が 0.6のときに構造Ａから生成される信号波形（サンプリング周波数は44.1KHz)
下図は、周波数f0が300Ｈｚで、振幅height0が0.3、キャリア信号と変調信号の周波数の
相対比ratio0が
11のときに構造Ｂから生成される信号波形（サンプリング周波数は44.1KHz) 高い周波数の波長が短いＳＩＮ波（キャリア）の振幅エンベロープが
ＳＩＮ波の形で変調されている。
下図は、上の２つの信号を混合合成し 更に、 ランダムなノイズを付加したもの（サンプリング周波数は44.1KHz)
まず、”構造的”に信号を発生さるということは、その発生させる物理的な仕組みに起因する制約（限度）があるはずである。この限度を、パラメーターが取ることがきる値の制約範囲(bound)として規定する。
制約範囲(bounds)
|
最小
|
中心
|
最大
|
構造Ａの周波数ｆ0
|
100
|
550
|
1000
|
構造Ａの振幅height0
|
0.1
|
1.05
|
2.0
|
構造Ａの矩形波と三角波の比率ratio0
|
0.5
|
0.7
|
0.9
|
構造Ｂの周波数ｆ0
|
100
|
550
|
1000
|
構造Ｂの振幅height0
|
0.1
|
0.55
|
1.0
|
構造Ｂのキャリアと変調の周波数の比ratio0
|
8
|
10
|
12
|
例えば、構造Ａは、矩形波を発生させると必ず三角波が同時に発生する。矩形波だけというケースは禁止で、必ず、5割りから9割の割合で三角波が付随してくることを意味している。
さて、推定値を求める方法として、最適化問題に使われる準ニュートン法などの最小2乗法のフィッティング手法を使えばよいと考えられるが、この問題の場合は上手くいかない。
その理由は、ある程度 真の解に近い、解に収束する初期値から最小2乗法を開始しないと、求めた推定値が的外れのものになって、似ても似つかぬ信号波形に
なってしまうからである。
そこで、まず、ラフな全探索を行う。
ラフな全探索では、すべてのポイントを探索していると膨大な時間がかかるため、実用的にダイナミックレンジを損なわず探索時間を少なくするため、周波数ｆ0と振幅height0は 対数logスケール上で均等に割り振ったポイントに限定して違いを計算する。そして、とりあえず、一番小さいポイントを暫定的な候補とする。
また、構造Ａと構造Ｂの2つを、同時に探索するのはなので、はじめは構造Ａだけを考えて（構造Ｂはないと仮定する） 探索し、構造Ａのめぼしがついたら、
次に、構造Ｂだけのパラメーターを可変して探索する。そして、求めた暫定的な候補を初期値として、全部のパラメーターを推定する。 つまり、個々に注目してラフな探索をして当たりをつけながら、だんだんと解に近づくように攻めて行く方法である。
これを、数値解析ソフトのGNU Octaveを使って、解いた（推定した）例を以下に示す。
octave:1> test1 ← test1.mという名前で下記のプログラムを保存しておいたものをOctave上で実行する。
Answer: A: f0=300 height0=1.2 ratio0=0.6 B: f0=300 height0=0.3 ratio0=11 ← これが正解の値
分析用のテスト信号として、構造Ａ（上段）と構造Ｂ（2段目）、更に乱数ノイズ（3段目）を加算して作成（一番下）する。
search Structure Sample A, only ← 最初は、構造Ａだけのフィッティングを行う
A: rough global search (1) ← しかも、ratio0は中心値に固定しておいて、f0とheight0の２つについて ラフな全探索を行う。
A: min f0->297.635 min height0->1.2462 value=19.253 ← ラフな検索で最小２乗誤差が一番小さかったものを仮の候補とする。
A: rough global search (2) ← f0は仮候補で固定しておいて、height0とratio0の２つについて ラフな全探索を行う。
A: min height0->1.2462 min ratio0->0.5 value=16.755 ← ラフな検索で最小２乗誤差が一番小さかったものを仮の候補とする。
A: call function leastsq (3) ← 上記で求めた仮候補を初期値としてleastsqr (Octave
のleasqrファンクション。非線形でLevenberg-Marquardt法）で探索する
A: kvg1=1 iter1=5
A: min f0->300.677 min height0->1.1865 min ratio0->0.656 ← 求めたf0, height0, ratio0の構造Ａの推定値。
下図は、構造Ａの推定値を使った波形（緑色）と分析するテスト信号（赤色）の比較である。
search Structure Sample B, while Structure Sample A is fixed ← 次に、構造Ａは上記で求めた推定値を使って、構造Ｂのフィッティングを行う。
B: rough global search (4) ← しかも、ratio0は中心値に固定しておいて、f0とheight0の２つについて ラフな全探索を行う。
B: min f0->335.982 min height0->0.1 value=8.1092 ← ラフな検索で最小２乗誤差が一番小さかったものを仮の候補とする。
B: rough global search (5) ← f0は仮候補で固定しておいて、height0とratio0の２つについて ラフな全探索を行う。
B: min height0->0.26367 min ratio0->9.7778 value=4.0194 ← ラフな検索で最小２乗誤差が一番小さかったものを仮の候補とする。
B: call function leastsq (6) ← 上記で求めた仮候補を初期値としてleasqrで探索する
B: kvg2=1 iter2=14
B: f0=305.04 height0=0.29549 ratio0=10.82 ← 求めたf0, height0, ratio0の構造Ｂの推定値。
下図は、構造Ｂの推定値を使った波形（緑色）と分析するテスト信号（赤色）の比較である。
search Structure Sample A and B, at once ← 上記で求めた仮候補を初期値としてＡとＢのずべてのパラメータを可変させて探索する。
A+B: call function leastsq (7)
A+B: kvg2=1 iter2=14
A+B: A: f0=300.679 height0=1.1866 ratio0=0.65568B: f0=305.04 height0=0.2955 ratio0=10.82← leasqrで求めた構造Ａのf0, height0,
ratio0 と構造Ｂのf0, height0, ratio0 の推定値。
下図は、構造Ａと構造Ｂの推定値を使った波形（緑色）と分析するテスト信号（赤色）の比較である。
参考までに、
ここで使ったGNU Octave用のサンプルプログラムを以下に紹介しておきます。 Octaveの外部のoptimパッケージも必要です。このサンプルプログラムの動作保証はありませんのであしからず。
>cat test1.m
%//
%// use for GNU Octave version 3.6.4 and version 3.8.1
%//
%//------------------------------------------------------
%// WARNING: There is ABSOLUTELY NO WARRANTY in this program.
%// If you use this program, everyhing will be done by your own risk.
%//
1; %// This is not a function file:
%////////////////////////////////////////////////////
%////////////////////////////////////////////////////
function y=rec_tri(f0, height0, ratio0)
%//
%// Struture Sample A: mix of rectangle wave and triangle wave
%//
%//
%// Fs is sampling frequecny [Hz]
%// f0 is the target frequecny
%// height0 is height from zero, rectangle value
%// ratio0 is ratio of triangle wave vs rectangle
%//===================================================
%// restricting condition due to this structure
%// 0.5 < ratio0 < 0.9
%//===================================================
global Fs;
global Size0
wl0= Fs/f0;
wlh0=wl0/2.0;
a=height0/wlh0;
b=height0;
y=zeros(Size0,1);
for i=1:Size0
x= (i-1);
x= x - (fix(x/wl0)) * wl0;
if (x <= wlh0) %// if x <= wlh0 then
y(i)= b + ratio0 * (a * x);
else
y(i)= -1.0 * a * (x - wlh0) + b;
endif %// end
end
endfunction
%////////////////////////////////////////////////////////
%////////////////////////////////////////////////////
%////////////////////////////////////////////////////
function y=sin_sin(f0, height0, ratio0)
%//
%// Struture Sample B: sin wave with AM modulation by sin
%//
%//
%// Fs is sampling frequecny [Hz]
%// f0 is AM modulation frequecny (per one pulse)
%// height0 is height from zero, sin
%// ratio0 is ratio of carrier frequecny vs f0
%//===================================================
%// restricting condition due to this structure
%// 8 < ratio0 < 12
%//===================================================
global Fs;
global Size0;
f0c=f0 * ratio0; %// carrier
wl0= Fs/f0c;
b=height0;
fmod0=f0 / 2; %// down to half, to be similar to triangle wave of Structure Sample A
wlm0= Fs/fmod0;
y=zeros(Size0,1);
d0=2 * pi * (f0c / Fs);
dmod0=2 * pi * (fmod0 / Fs);
for i=1:Size0
x= (i-1);
xm=x;
x= x - (fix(x/wl0)) * wl0;
xm= xm - (fix(x/wlm0)) * wlm0;
y(i)=b * sin( d0 * x) * abs( sin( dmod0 * xm));
end
endfunction
%////////////////////////////////////////////////////////
%////////////////////////////////////////////////////
%////////////////////////////////////////////////////
function y=myfunc(xin,p)
%// ignore xin, output size is always Size0
global Fs;
global Size0;
global p_global;
global para_quant;
global AB_switch0;
l0=length(p);
p_global2=p_global;
if( l0 == para_quant)
for i=1:l0
p_global2(i)=p(i);
end
else
for i=1:l0
p_global2(i + AB_switch0 * 3)=p(i); %// When
AB_switch0 is 1, B side. When AB_switch0 is 0, A side.
end
end
y=zeros(Size0,1);
%// (1) Structure sample A
f0=p_global2(1);
height0=p_global2(2);
ratio0=p_global2(3);
wl0= Fs/f0;
wlh0=wl0/2.0;
a=height0/wlh0;
b=height0;
for i=1:Size0
x= (i-1);
x= x - (fix(x/wl0)) * wl0;
if (x <= wlh0) %// if x <= wlh0 then
y(i)= b + ratio0 * (a * x);
else
y(i)= -1.0 * a * (x - wlh0) + b;
endif %// end
end
%// (2) Structure sample B
f0=p_global2(4);
height0=p_global2(5);
ratio0=p_global2(6);
f0c=f0 * ratio0; %// carrier
wl0= Fs/f0c;
b=height0;
fmod0=f0 / 2; %// down to half, to be similar to triangle wave of Structure Sample A
wlm0= Fs/fmod0;
d0=2 * pi * (f0c / Fs);
dmod0=2 * pi * (fmod0 / Fs);
for i=1:Size0
x= (i-1);
xm=x;
x= x - (fix(x/wl0)) * wl0;
xm= xm - (fix(x/wlm0)) * wlm0;
y(i)=y(i)+b * sin( d0 * x) * abs( sin( dmod0 * xm));
end
endfunction
%////////////////////////////////////////////////////////
%////////////////////////////////////////////////////////
%////////////////////////////////////////////////////////
%// COMMON DEFINITION
global Fs;
global Span0;
global Size0;
Fs=44100;
Span0=2;
global p_global;
global para_quant;
para_quant=6; %// Please set parameters quantity
p_global=zeros(para_quant,1);
%//bounds control
minstep0 = [ 100; 0.1; 0.5; 100; 0.1; 8];
maxstep0 = [1000; 2.0; 0.9; 1000; 1.0; 12];
global AB_switch0;
AB_switch0=0; %// When AB_switch0 = 0, then A side, other B side In myfunc
%////////////////////////////////////////////////////////
%////////////////////////////////////////////////////////
%////////////////////////////////////////////////////////
%// make a test signal
%//
%// (1) apply Structure sample A
f1=300;
height1=1.2;
ratio1=0.6; %// ratio1 must be in range of restricting condition due to structure sample A
Size0=fix(Span0 * Fs/f1)+1;
y1=rec_tri(f1, height1, ratio1);
%//
%// (2) aply Strutrue sample B
f2=300;
height2=0.3;
ratio2=11; %// ratio2 must be in range of restricting condition due to structure sample B
y2=sin_sin(f2, height2, ratio2);
%//
%// (3) and add random noise
seed0=0;
height3=0.1;
randn('seed', seed0);
yran1= randn(Size0,1) * height3;
%// combination: this is the test signal, analytical signal
yall=y1+y2+yran1;
%////////////////////////////////////////////////////////
%// show grpah
figure(1);
clf;
subplot(4, 1, 1);
plot( y1 );
title("mix of triangle wave and rectangle wave");
grid "on";
subplot(4, 1, 2);
plot( y2 );
title("sin wave with amplitude modulation by sin");
grid "on";
subplot(4, 1, 3);
plot( yran1 );
title("addtional random noise");
grid "on";
subplot(4, 1, 4);
plot( yall );
title("the test signal, analytical signal, above combination");
grid "on";
%////////////////////////////////////////////////////////
%////////////////////////////////////////////////////////
%////////////////////////////////////////////////////////
%// Search
%//
p_kai=[f1; height1; ratio1; f2; height2; ratio2;];
disp([ "Answer: A: f0=", num2str(p_kai(1)) ," height0=",
num2str(p_kai(2)), " ratio0=", num2str(p_kai(3)), " B: f0=",
num2str(p_kai(4)) ," height0=", num2str(p_kai(5)), " ratio0=",
num2str(p_kai(6)) ]);
xin=1:Size0;
xin=xin';
%////////////////////////////////////////////////////////
%// set tentative value
f10=(maxstep0(1)-minstep0(1))/2+minstep0(1);
height10=(maxstep0(2)-minstep0(2))/2+minstep0(2);
ratio10=(maxstep0(3)-minstep0(3))/2+minstep0(3);
f20=(maxstep0(4)-minstep0(4))/2+minstep0(4);
height20=(maxstep0(5)-minstep0(5))/2+minstep0(5);
ratio20=(maxstep0(6)-minstep0(6))/2+minstep0(6);
p_global=[f10; height10; ratio10; f20; height20; ratio20;];
%////////////////////////////////////////////////////////
%// set global search range
f1_BUNKATU=20;
lf0=log10(minstep0(1));
df0=(log10(maxstep0(1)) - log10(minstep0(1))) / (f1_BUNKATU-1);
f_list1=zeros(f1_BUNKATU,1); %// list for the A f0
for i=1:f1_BUNKATU
f_list1(i)=10^(lf0 + df0 * (i-1));
end
h1_BUNKATU=20;
lf0=log10(minstep0(2));
df0=(log10(maxstep0(2)) - log10(minstep0(2))) / (h1_BUNKATU-1);
h_list1=zeros(h1_BUNKATU,1); %// list for the A heigh0
for i=1:h1_BUNKATU
h_list1(i)=10^(lf0 + df0 * (i-1));
end
r1_BUNKATU=10;
lf0=minstep0(3);
df0=(maxstep0(3) - minstep0(3)) / (r1_BUNKATU-1);
r_list1=zeros(r1_BUNKATU,1); %// list for the A ratio0
for i=1:r1_BUNKATU
r_list1(i)=lf0 + df0 * (i-1);
end
f2_BUNKATU=20;
lf0=log10(minstep0(4));
df0=(log10(maxstep0(4)) - log10(minstep0(4))) / (f2_BUNKATU-1);
f_list2=zeros(f2_BUNKATU,1); %// list for the B f0
for i=1:f2_BUNKATU
f_list2(i)=10^(lf0 + df0 * (i-1));
end
h2_BUNKATU=20;
lf0=log10(minstep0(5));
df0=(log10(maxstep0(5)) - log10(minstep0(5))) / (h2_BUNKATU-1);
h_list2=zeros(h2_BUNKATU,1); %// list for the B heigh0
for i=1:h2_BUNKATU
h_list2(i)=10^(lf0 + df0 * (i-1));
end
r2_BUNKATU=10;
lf0=minstep0(6);
df0=(maxstep0(6) - minstep0(6)) / (r2_BUNKATU-1);
r_list2=zeros(r2_BUNKATU,1); %// list for the B ratio0
for i=1:r2_BUNKATU
r_list2(i)=lf0 + df0 * (i-1);
end
%//////////////////////////////////////////////////////////////////////
%//
%// search Structure Sample A, only
%//
%// (1)rough global search of Structure Sample A, f0 and height0, only
disp(["search Structure Sample A, only"]);
AB_switch0=0; %// set A
CAL_FLAG1=1; %//<----- SET CONTROL FLAG
if (CAL_FLAG1 == 1) %//if (CAL_FLAG = 1)
hf_value1=zeros(h1_BUNKATU , f1_BUNKATU);
p=[f10; height10; ratio10; f20; 0; ratio20;]; %// height20 is zero, due to B is zero
disp(["A: rough global search (1)"]);
for i=1:h1_BUNKATU
for j=1:f1_BUNKATU
p(1)=f_list1(j);
p(2)=h_list1(i);
y=myfunc(xin,p);
hf_value1(i,j)= sumsq( (y - yall)' );
end %// j
end %// i
save data01.mat hf_value1;
else %//if (CAL_FLAG = 1)
load data01.mat hf_value1;
end %//if (CAL_FLAG = 1)
%// find minimum
[v0, r0]=min( hf_value1);
[v1, pf0]=min(v0);
ph0=r0(pf0);
disp(["A: min f0->", num2str(f_list1(pf0)) , " min
height0->", num2str(h_list1(ph0)), " value=", num2str(v1)]);
%//////////////////////////////////
%// (2)rough global search of Structure Sample A, height0 and ratio0, when f0 is fixed
CAL_FLAG2=1; %//<----- SET CONTROL FLAG
if (CAL_FLAG2 == 1) %//if (CAL_FLAG = 1)
hr_value1=zeros(h1_BUNKATU , r1_BUNKATU);
p=[f_list1(pf0); height10; ratio10; f20; 0; ratio20;]; %// height20 is zero, due to B is zero
disp(["A: rough global search (2)"]);
for i=1:h1_BUNKATU
for j=1:r1_BUNKATU
p(2)=h_list1(i);
p(3)=r_list1(j);
y=myfunc(xin,p);
hr_value1(i,j)= sumsq( (y - yall)' );
end %// j
end %// i
save data02.mat hr_value1;
else %//if (CAL_FLAG = 1)
load data02.mat hr_value1;
end %//if (CAL_FLAG = 1)
%// find minimum
[v0, r0]=min( hr_value1);
[v1, pr0]=min(v0);
ph0=r0(pr0);
disp([ "A: min height0->", num2str(h_list1(ph0)), " min
ratio0->", num2str(r_list1(pr0)) , " value=", num2str(v1)]);
%//////////////////////////////////
%// (3)call function leastsq, start from rough global serach
CAL_FLAG3=1; %//<----- SET CONTROL FLAG
if (CAL_FLAG3 == 1) %//if (CAL_FLAG = 1)
disp(["A: call function leastsq (3)"]);
p_global=[f_list1(pf0); h_list1(ph0); r_list1(pr0); f20;
0; ratio20;]; %// height20 is zero, due to B is
zero
p_init=[f_list1(pf0); h_list1(ph0); r_list1(pr0);]; %// height20 is zero, due to B is zero
stol0=0.001; %//tolerance, scalar
niter0=20; %//max repeat kaisuu
wt1=ones(size(yall)); %// omomi
dp0 = [0.01; 0.001; 0.001;]; % tolerance of dfdp
%//bounds control
minstep1 = minstep0(1:3,1);
maxstep1 = maxstep0(1:3,1);
options.bounds= [minstep1, maxstep1];
%load filter package to call leastsqr
pkg load optim;
[fout1, pout1, kvg1, iter1] =leasqr (xin, yall, p_init, "myfunc", stol0, niter0, wt1, dp0, "dfdp", options);
disp(["A: kvg1=", num2str(kvg1), " iter1=", num2str(iter1)]); %// kvg1 = 1 , Scalar: = 1, then success
save data03.mat pout1;
else %//if (CAL_FLAG = 1)
load data03.mat pout1;
end %//if (CAL_FLAG = 1)
disp([ "A: min f0->", num2str(pout1(1)) ," min height0->",
num2str(pout1(2)), " min ratio0->", num2str(pout1(3)) ]);
%////////////////////////////////////////////////////////
%// show grpah
p_global=[pout1(1); pout1(2); pout1(3); f20; 0; ratio20;]; %// height20 is zero, due to B is zero
fout1=myfunc(xin,p_global);
figure(2);
clf;
hold on % over-write
plot( yall, "r");
plot( fout1, "g");
hold off %
title("compare (A), red is the test signal, green is estimation");
grid "on";
%//
fflush(stdout); %// flush output
%//////////////////////////////////////////////////////////////////////
%//
%// search Structure Sample B, while Structure Sample A is fixed
%//
%// (4)rough global search of Structure Sample B, f0 and height0, only
disp(["search Structure Sample B, while Structure Sample A is fixed"]);
AB_switch0=1; %// set B
CAL_FLAG4=1; %//<----- SET CONTROL FLAG
if (CAL_FLAG4 == 1) %//if (CAL_FLAG = 1)
hf_value2=zeros(h2_BUNKATU , f2_BUNKATU);
p=[pout1(1); pout1(2); pout1(3); f20; height20; ratio20;];
disp(["B: rough global search (4)"]);
for i=1:h2_BUNKATU
for j=1:f2_BUNKATU
p(4)=f_list2(j);
p(5)=h_list2(i);
y=myfunc(xin,p);
hf_value2(i,j)= sumsq( (y - yall)' );
end %// j
end %// i
save data04.mat hf_value2;
else %//if (CAL_FLAG = 1)
load data04.mat hf_value2;
end %//if (CAL_FLAG = 1)
%// find minimum
[v0, r0]=min( hf_value2);
[v1, pf1]=min(v0);
ph1=r0(pf1);
disp(["B: min f0->", num2str(f_list2(pf1)) , " min
height0->", num2str(h_list2(ph1)), " value=", num2str(v1)]);
%//////////////////////////////////
%// (5)rough global search of Structure Sample B, height0 and ratio0, when f0 is fixed
CAL_FLAG5=1; %//<----- SET CONTROL FLAG
if (CAL_FLAG5 == 1) %//if (CAL_FLAG = 1)
hr_value2=zeros(h2_BUNKATU , r2_BUNKATU);
p=[pout1(1); pout1(2); pout1(3); f_list2(pf1); height20; ratio20;];
disp(["B: rough global search (5)"]);
for i=1:h2_BUNKATU
for j=1:r2_BUNKATU
p(5)=h_list2(i);
p(6)=r_list2(j);
y=myfunc(xin,p);
hr_value2(i,j)= sumsq( (y - yall)' );
end %// j
end %// i
save data05.mat hr_value2;
else %//if (CAL_FLAG = 1)
load data05.mat hr_value2;
end %//if (CAL_FLAG = 1)
%// find minimum
[v0, r0]=min( hr_value2);
[v1, pr1]=min(v0);
ph1=r0(pr1);
disp([ "B: min height0->", num2str(h_list2(ph1)), " min
ratio0->", num2str(r_list2(pr1)) , " value=", num2str(v1)]);
%//////////////////////////////////////////////////////////////////////
%//
%// (6)call function leastsq, start from rough global serach
CAL_FLAG6=1; %//<----- SET CONTROL FLAG
if (CAL_FLAG6 == 1) %//if (CAL_FLAG = 1)
disp(["B: call function leastsq (6)"]);
p_global=[pout1(1); pout1(2); pout1(3); f_list2(pf1); h_list2(ph1); r_list2(pr1);];
p_init=[ f_list2(pf1); h_list2(ph1); r_list2(pr1);];
stol0=0.001; %//tolerance, scalar
niter0=20; %//max repeat kaisuu
wt1=ones(size(yall)); %// omomi
dp0 = [ 0.01; 0.001; 0.001;]; % tolerance of dfdp
%//bounds control
minstep1 = minstep0(4:6,1);
maxstep1 = maxstep0(4:6,1);
options.bounds= [minstep1, maxstep1];
%load filter package to call leastsqr
pkg load optim;
[fout2, pout2, kvg2, iter2] =leasqr (xin, yall, p_init, "myfunc", stol0, niter0, wt1, dp0, "dfdp", options);
disp(["B: kvg2=", num2str(kvg2), " iter2=", num2str(iter2)]); %// kvg1 = 1 , Scalar: = 1, then success
save data06.mat pout2;
else %//if (CAL_FLAG = 1)
load data06.mat pout2;
end %//if (CAL_FLAG = 1)
disp(["B: f0=",num2str(pout2(1))," height0=",num2str(pout2(2))," ratio0=",num2str(pout2(3))]);
%////////////////////////////////////////////////////////
%// show grpah
p_global=[pout1(1); pout1(2); pout1(3); pout2(1); pout2(2); pout2(3);];
fout1=myfunc(xin,p_global);
figure(3);
clf;
hold on % over-write
plot( yall, "r");
plot( fout1, "g");
hold off %
title("compare (B), red is the test signal, green is estimation");
grid "on";
%//
fflush(stdout); %// flush output
%//////////////////////////////////////////////////////////////////////
%//
%// search Structure Sample A and B, at once
%//
%//
%// (7)call function leastsq, after each leastsq
disp(["search Structure Sample A and B, at once"]);
CAL_FLAG7=1; %//<----- SET CONTROL FLAG
if (CAL_FLAG7 == 1) %//if (CAL_FLAG = 1)
disp(["A+B: call function leastsq (7)"]);
p_global=[pout1(1); pout1(2); pout1(3); pout2(1); pout2(2); pout2(3);];
p_init=[pout1(1); pout1(2); pout1(3); pout2(1); pout2(2); pout2(3);];
stol0=0.001; %//tolerance, scalar
niter0=20; %//max repeat kaisuu
wt1=ones(size(yall)); %// omomi
dp0 = [ 0.01; 0.001; 0.001; 0.01; 0.001; 0.001;]; % tolerance of dfdp
%//bounds control
minstep1 = minstep0(1:6,1);
maxstep1 = maxstep0(1:6,1);
options.bounds= [minstep1, maxstep1];
%load filter package to call leastsqr
pkg load optim;
[fout3, pout3, kvg3, iter3] =leasqr (xin, yall, p_init, "myfunc", stol0, niter0, wt1, dp0, "dfdp", options);
disp(["A+B: kvg2=", num2str(kvg2), " iter2=", num2str(iter2)]); %// kvg1 = 1 , Scalar: = 1, then success
save data07.mat pout3;
else %//if (CAL_FLAG = 1)
load data07.mat pout3;
end %//if (CAL_FLAG = 1)
disp(["A+B: A: f0=",num2str(pout3(1))," height0=",num2str(pout3(2)),"
ratio0=",num2str(pout3(3)),"B: f0=",num2str(pout3(4)),"
height0=",num2str(pout3(5))," ratio0=",num2str(pout3(6))]);
%////////////////////////////////////////////////////////
%// show grpah
p_global=[pout3(1); pout3(2); pout3(3); pout3(4); pout3(5); pout3(6);];
fout1=myfunc(xin,p_global);
figure(4);
clf;
hold on % over-write
plot( yall, "r");
plot( fout1, "g");
hold off %
title("compare (A+B), red is the test signal, green is estimation");
grid "on";
No.2b 2015年2月14日
