FreeCADのバージョンを確認していると0.13になっていたので，新しいバージョンに変えることにした。
Windows 8.1にインストールしてみるが，どうも使用できない．
0.14も確認したが，使用できない．
参った
ある実験に関して次のような関数があるとする（未知）
a_ = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0];
hFunc = @(p1, p2, p3)(a_(1)*p1.^2+a_(2)*p2.^2+a_(3)*p3.^2+a_(4)*p1.*p2+a_(5)*p2.*p3+a_(6)*p3.*p1+a_(7)*p1+a_(8)*p2+a_(9)*p3+a_(10));
%パラメータを決めて実験をし、ノイズを含む実験結果
%hFunc(p(:, 1), p(:, 2), p(:, 3)) + 0.005 * rand(size(a))
%を得る．
p = rand(10, 3);
p = [p;p;p;p];
% 例えばのパラメータの組み合わせ
a = p(:, 1);
b = p(:, 2);
c = p(:, 3);
% 関数系を推定しフィッティングを行う．
large_X = [a.^2,b.^2,c.^2,a.*b,b.*c,c.*a,a,b,c,ones(size(a))];
beta = regress(hFunc(p(:, 1), p(:, 2), p(:, 3)) + 0.005 * rand(size(a)), large_X)
betaはフィッティングされた関数の係数列
regressという関数の利用方法を検討中．
例えば２次の項までを持った関数
a_ = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0];
hFunc = @(p1, p2, p3)(a_(1)*p1.^2+a_(2)*p2.^2+a_(3)*p3.^2+a_(4)*p1.*p2+a_(5)*p2.*p3+a_(6)*p3.*p1+a_(7)*p1+a_(8)*p2+a_(9)*p3+a_(10));
を定義する．a_は係数列．この関数を近似する．
試しの引数を
p = rand(10, 3);
a = p(:, 1);
b = p(:, 2);
c = p(:, 3);
とする．（a, b, cが引き数列）
係数の推定は以下のようにする．
large_X = [a.^2,b.^2,c.^2,a.*b,b.*c,c.*a,a,b,c,ones(size(a))];
beta = regress(hFunc(p(:, 1), p(:, 2), p(:, 3)), large_X);
betaが推定された係数．
最初は
pkg install -forge statistics
pkg load statistics
普通のC/C++で記述する場合には全く問題にならないが、
先日のプログラムでOctaveで実行すると少し長いデータでは時間が掛かり気味．
以下はまだ実験中のコード。（これでも全然速くはならない．）
x_buf = zeros(1, 3);
y_buf = zeros(1, 2);
y = zeros(size(x));
for k = 1:length(x);
x_buf(2:3) = x_buf(1:2);
x_buf(1) = x(k);
y(k) = b(1)*x_buf(1) + b(2)*x_buf(2) + b(3)*x_buf(3) - a(2) * y_buf(1) - a(3) * y_buf(2);
y_buf(2) = y_buf(1);
y_buf(1) = y(k);
end
y2 = zeros(size(x));
for k = 1:length(x);
x_buf(2:3) = x_buf(1:2);
x_buf(1) = x(k);
y2(k) = sum(b .* x_buf) - sum(a(2:3) .* y_buf);
y_buf(2) = y_buf(1);
y_buf(1) = y2(k);
end
filter関数の中身が少し分かったので、開始直後の特性についてもいじることができるようになった。
時定数の大きなフィルタを作る場合には、フィルタの動作開始から安定するまでの間 時間が掛かる。
そこを少し手を加えることができるようになる。（下はa(1)=1を前提としている実装．違う場合a, b をa(1)で割る）
t = 0:0.1:10; % Sampling freq. 10 Hz
% Sampling freq. 10 Hz
t = 0:0.1:10; % Sampling freq. 10 Hz
x = sin(t * 3.0 * 2 * pi) + sin(t * 1.0 * 2 * pi) + 5.;
[b, a] = butter(2, 0.4);
filt_x = filter(b, a, x);
% org. filter
x_buf = zeros(1, 3);
y_buf = zeros(1, 2);
y = zeros(size(x));
for k = 1:length(x);
x_buf(2:3) = x_buf(1:2);
x_buf(1) = x(k);
y(k) = b(1)*x_buf(1) + b(2)*x_buf(2) + b(3)*x_buf(3) - a(2) * y_buf(1) - a(3) * y_buf(2);
y_buf(2) = y_buf(1);
y_buf(1) = y(k);
end
% org. filter2 最初だけいじってみる
y2 = zeros(size(x));
x_buf = zeros(1, 3);
y_buf = zeros(1, 2);
x_buf = x(1:3);
for k = 1:length(x);
x_buf(2:3) = x_buf(1:2);
x_buf(1) = x(k);
y2(k) = b(1)*x_buf(1) + b(2)*x_buf(2) + b(3)*x_buf(3) - a(2) * y_buf(1) - a(3) * y_buf(2);
y_buf(2) = y_buf(1);
y_buf(1) = y2(k);
end