Кусочно-линейная аппроксимация
Вчера я занимался построением зависимости отчетов в ADU от четырех каналов измеряемой температуры контроллера чиллера. Для этого замотал на термопасте четыре NTC и один платиновый терморезистор (для контроля) между двух кусочков фольгированного стеклотекстолита. Поместил это в баночку из-под детского питания, залил баночку тосолом, поставил в банку из-под кофе и залил азотом. Как только азот испарился (температура тосола опустилась до -140°C), залил в железную банку теплого тосола и стал фиксировать значения ADU и температуры по мере нагревания.
В итоге у меня получилась табличка, которую теперь нужно обработать. Каждая запись в таблице - медианное усреднение девяти последовательных измерений с интервалом в 0.1с. В среднем СКО получилось около 2.5 на все четыре канала. На низких температурах расхождение их показаний не сильно отличалось от СКО, однако, по мере уменьшения сопротивления расхождение было все больше и больше. Несмотря на то, что используемые в качестве опорных резисторов имеют паспортную точность в 0.1% (т.е. 1Ом), их сопротивления имеют значение от 994 до 998 Ом, вот, собственно, где собака и порылась!
Я не нашел, как в octave найти узлы кусочно-линейной интерполяции, поэтому пришлось велосипедить.
Итак, я занес данные в файл Temperatures.dat, где первый столбец - сопротивление терморезистора, второй - вычисленная по этому сопротивлению температура, а оставшиеся четыре - показания АЦП контроллера чиллера. Так как мне не нужна температура ниже -20°C, первые строки можно обрезать:
d = dlmread('Temperatures.dat'); d(1:17,:) = [];
T = d(:,2); % температуры
R = d(:,[3:6]); % ADU сопротивлений
Строим матрицу для кубической аппроксимации
M3=[ones(size(T)) T T.^2 T.^3];
И вычисляем коэффициенты:
km3 = M3\median(R,2)
km3 =
1002.2286743
37.1381302
0.3308417
-0.0044131
Строим новые матрицы с шагом 0.05 градусов:
Tt = [-20:0.05:20]';
Mt = [ones(size(Tt)) Tt Tt.^2 Tt.^3];
ADU = Mt * km3;
Теперь имеем обратную зависимость Tt(ADU), которой нужно подобрать кусочно-линейную аппроксимацию.
Ищем узлы:
knots = [1 getknots(ADU, Tt, 0.1) max(size(ADU))]
knots =
1 42 84 169 257 350 446 550 664 801
Строим график:
plot(ADU, Tt, ADU(knots), Tt(knots), '+')
print -dpng knots.png
Формируем таблицу:
[X0 Y0 K]=buildk(ADU, Tt, 0.1)
X0 =
427.11 467.72 514.28 622.83 753.63 909.75 1087.41 1295.45 1537.77
Y0 =
-20.0000 -17.9500 -15.8500 -11.6000 -7.2000 -2.5500 2.2500 7.4500 13.1500
K =
0.050477 0.045107 0.039150 0.033639 0.029785 0.027017 0.024996 0.023522 0.022514
Сверяем:
mR2 = median(R,2);
cT=[];for i = 1:max(size(mR2)); cT=[cT calcT(mR2(i))]; endfor
plot(R,T,'o', mR2,T,'*', mR2,cT);
print -dpng check.png
В области высоких температур все получается довольно-таки хреново:
размах показаний такой, что вытягивает почти на 1°C! Однако, найденная кусочно-линейная интерполяция вполне дает желаемую точность в ±0.5°C. Таким образом, строить индивидуальные зависимости для каждого канала не нужно! Кстати, если ограничиться размахом в 0.5°C (а не в 0.1°C), то получится всего четыре узла! Но, т.к. 9 узлов не так-то много места в памяти займут, пусть будет так!
Смеха ради посмотрел, сколько будет узлов с размахом не больше 0.01°C. Всего-то 28 штук! Т.е. данная методика вполне подходит и для аппроксимации измерений температур при помощи платиновых терморезисторов на внешнем АЦП. Можно будет ею воспользоваться, если доберусь-таки до системы точной регистрации температур (узел АЦП/мультиплексора я еще летом спаял, нужно лишь подключить терморезисторы и микроконтроллер, да написать прошивку).
Исходники утилит.
function Tout = H705(Rin)
% Converts resistance of TRD into T (degrC)
_alpha = 0.00375;
_beta = 0.16;
_delta = 1.605;
T = [-300:0.1:300];
_A = _alpha + _alpha*_delta/100.;
_B = -_alpha*_delta/1e4;
_C = zeros(size(T));
_C(find(T<0.)) = -_alpha*_beta/1e8;
rT = 1000.*(1 + _A*T + _B*T.^2 - _C.*T.^3*100. + _C.*T.^4);
%plot(T, rT);
Tout = interp1(rT, T, Rin, 'spline');
endfunction
function idx = getknots(X, Y, dYmax)
%
% idx = getknots(X, Y, dYmax) - calculate piecewise-linear approximation knots for Y(X)
% dYmax - maximal deviation
%
idx = [];
I = getnewpt(X, Y, dYmax);
if(I > 0)
L = getknots(X(1:I), Y(1:I), dYmax);
R = I-1 + getknots(X(I:end), Y(I:end), dYmax);
idx = [L I R];
endif
endfunction
function idx = getnewpt(X, Y, delt)
%
% idx = getnewpt(X, Y, delt)
% find point where Y-`linear approx` is maximal
% return index of max deviation or -1 if it is less than `delt`
%
% make approximation
[x0 y0 k] = linearapprox(X,Y);
% find new knot
adiff = abs(Y - (y0 + k*(X-x0)));
maxadiff = max(adiff);
if(maxadiff < delt)
idx = -1;
else
idx = find(adiff == max(adiff));
endif
endfunction
function [X0 Y0 K] = buildk(X, Y, dYmax)
%
% [X0 Y0 K] = buildk(X, Y, dYmax) - build arrays of knots (X0, Y0) and linear koefficients K
% for piecewise-linear approximation of Y(X) with max deviation `dYmax`
%
X0 = []; Y0 = []; K = [];
knots = [1 getknots(X, Y, dYmax) max(size(X))];
for i = 1:max(size(knots))-1
idx = knots(i):knots(i+1);
[x y k] = linearapprox(X(idx), Y(idx));
X0 = [X0 x]; Y0 = [Y0 y]; K = [K k];
endfor
endfunction
function [x0 y0 k] = linearapprox(X,Y)
%
% [a b] = linearapprox(X,Y) - find linear approximation of function Y(X) through first and last points
% y = y0 + k*(X-x0)
%
x0 = X(1); y0 = Y(1);
k = (Y(end) - y0) / (X(end) - x0);
endfunction
И функция для проверки (можно было бы построить структуру при помощи mkpp и дальше средствами octave работать, но мне интересно было проверить наиближайшее приближение к тому, что будет на МК):
function T = calcT(ADU)
%
% T = calcT(ADU) - piecewise calculation prototype
%
X0 = [427.11 467.72 514.28 622.83 753.63 909.75 1087.41 1295.45 1537.77];
Y0 = [-20.0000 -17.9500 -15.8500 -11.6000 -7.2000 -2.5500 2.2500 7.4500 13.1500];
K = [0.050477 0.045107 0.039150 0.033639 0.029785 0.027017 0.024996 0.023522 0.022514];
imax = max(size(X0)); idx = int32((imax+1)/2);
T = [];
% find index
while(idx > 0 && idx <= imax)
x = X0(idx);
half = int32(idx/2);
if(ADU < x)
%printf("%g < %g\n", ADU, x);
if(idx == 1) break; endif; % less than first index
if(ADU > X0(idx-1)) idx -= 1; break; endif; % found
idx = half; % another half
else
%printf("%g > %g\n", ADU, x);
if(idx == imax) break; endif; % more than last index
if(ADU < X0(idx+1)) break; endif; % found
idx += half;
endif
endwhile
if(idx < 1) printf("too low (<%g)!\n", X0(1)); idx = 1; endif
if(idx > imax) idx = imax; endif;
T = Y0(idx) + K(idx) * (ADU - X0(idx));
endfunction
В итоге у меня получилась табличка, которую теперь нужно обработать. Каждая запись в таблице - медианное усреднение девяти последовательных измерений с интервалом в 0.1с. В среднем СКО получилось около 2.5 на все четыре канала. На низких температурах расхождение их показаний не сильно отличалось от СКО, однако, по мере уменьшения сопротивления расхождение было все больше и больше. Несмотря на то, что используемые в качестве опорных резисторов имеют паспортную точность в 0.1% (т.е. 1Ом), их сопротивления имеют значение от 994 до 998 Ом, вот, собственно, где собака и порылась!
Я не нашел, как в octave найти узлы кусочно-линейной интерполяции, поэтому пришлось велосипедить.
Итак, я занес данные в файл Temperatures.dat, где первый столбец - сопротивление терморезистора, второй - вычисленная по этому сопротивлению температура, а оставшиеся четыре - показания АЦП контроллера чиллера. Так как мне не нужна температура ниже -20°C, первые строки можно обрезать:
d = dlmread('Temperatures.dat'); d(1:17,:) = [];
T = d(:,2); % температуры
R = d(:,[3:6]); % ADU сопротивлений
Строим матрицу для кубической аппроксимации
M3=[ones(size(T)) T T.^2 T.^3];
И вычисляем коэффициенты:
km3 = M3\median(R,2)
km3 =
1002.2286743
37.1381302
0.3308417
-0.0044131
Строим новые матрицы с шагом 0.05 градусов:
Tt = [-20:0.05:20]';
Mt = [ones(size(Tt)) Tt Tt.^2 Tt.^3];
ADU = Mt * km3;
Теперь имеем обратную зависимость Tt(ADU), которой нужно подобрать кусочно-линейную аппроксимацию.
Ищем узлы:
knots = [1 getknots(ADU, Tt, 0.1) max(size(ADU))]
knots =
1 42 84 169 257 350 446 550 664 801
Строим график:
plot(ADU, Tt, ADU(knots), Tt(knots), '+')
print -dpng knots.png
Формируем таблицу:
[X0 Y0 K]=buildk(ADU, Tt, 0.1)
X0 =
427.11 467.72 514.28 622.83 753.63 909.75 1087.41 1295.45 1537.77
Y0 =
-20.0000 -17.9500 -15.8500 -11.6000 -7.2000 -2.5500 2.2500 7.4500 13.1500
K =
0.050477 0.045107 0.039150 0.033639 0.029785 0.027017 0.024996 0.023522 0.022514
Сверяем:
mR2 = median(R,2);
cT=[];for i = 1:max(size(mR2)); cT=[cT calcT(mR2(i))]; endfor
plot(R,T,'o', mR2,T,'*', mR2,cT);
print -dpng check.png
В области высоких температур все получается довольно-таки хреново:
размах показаний такой, что вытягивает почти на 1°C! Однако, найденная кусочно-линейная интерполяция вполне дает желаемую точность в ±0.5°C. Таким образом, строить индивидуальные зависимости для каждого канала не нужно! Кстати, если ограничиться размахом в 0.5°C (а не в 0.1°C), то получится всего четыре узла! Но, т.к. 9 узлов не так-то много места в памяти займут, пусть будет так!
Смеха ради посмотрел, сколько будет узлов с размахом не больше 0.01°C. Всего-то 28 штук! Т.е. данная методика вполне подходит и для аппроксимации измерений температур при помощи платиновых терморезисторов на внешнем АЦП. Можно будет ею воспользоваться, если доберусь-таки до системы точной регистрации температур (узел АЦП/мультиплексора я еще летом спаял, нужно лишь подключить терморезисторы и микроконтроллер, да написать прошивку).
Исходники утилит.
function Tout = H705(Rin)
% Converts resistance of TRD into T (degrC)
_alpha = 0.00375;
_beta = 0.16;
_delta = 1.605;
T = [-300:0.1:300];
_A = _alpha + _alpha*_delta/100.;
_B = -_alpha*_delta/1e4;
_C = zeros(size(T));
_C(find(T<0.)) = -_alpha*_beta/1e8;
rT = 1000.*(1 + _A*T + _B*T.^2 - _C.*T.^3*100. + _C.*T.^4);
%plot(T, rT);
Tout = interp1(rT, T, Rin, 'spline');
endfunction
function idx = getknots(X, Y, dYmax)
%
% idx = getknots(X, Y, dYmax) - calculate piecewise-linear approximation knots for Y(X)
% dYmax - maximal deviation
%
idx = [];
I = getnewpt(X, Y, dYmax);
if(I > 0)
L = getknots(X(1:I), Y(1:I), dYmax);
R = I-1 + getknots(X(I:end), Y(I:end), dYmax);
idx = [L I R];
endif
endfunction
function idx = getnewpt(X, Y, delt)
%
% idx = getnewpt(X, Y, delt)
% find point where Y-`linear approx` is maximal
% return index of max deviation or -1 if it is less than `delt`
%
% make approximation
[x0 y0 k] = linearapprox(X,Y);
% find new knot
adiff = abs(Y - (y0 + k*(X-x0)));
maxadiff = max(adiff);
if(maxadiff < delt)
idx = -1;
else
idx = find(adiff == max(adiff));
endif
endfunction
function [X0 Y0 K] = buildk(X, Y, dYmax)
%
% [X0 Y0 K] = buildk(X, Y, dYmax) - build arrays of knots (X0, Y0) and linear koefficients K
% for piecewise-linear approximation of Y(X) with max deviation `dYmax`
%
X0 = []; Y0 = []; K = [];
knots = [1 getknots(X, Y, dYmax) max(size(X))];
for i = 1:max(size(knots))-1
idx = knots(i):knots(i+1);
[x y k] = linearapprox(X(idx), Y(idx));
X0 = [X0 x]; Y0 = [Y0 y]; K = [K k];
endfor
endfunction
function [x0 y0 k] = linearapprox(X,Y)
%
% [a b] = linearapprox(X,Y) - find linear approximation of function Y(X) through first and last points
% y = y0 + k*(X-x0)
%
x0 = X(1); y0 = Y(1);
k = (Y(end) - y0) / (X(end) - x0);
endfunction
И функция для проверки (можно было бы построить структуру при помощи mkpp и дальше средствами octave работать, но мне интересно было проверить наиближайшее приближение к тому, что будет на МК):
function T = calcT(ADU)
%
% T = calcT(ADU) - piecewise calculation prototype
%
X0 = [427.11 467.72 514.28 622.83 753.63 909.75 1087.41 1295.45 1537.77];
Y0 = [-20.0000 -17.9500 -15.8500 -11.6000 -7.2000 -2.5500 2.2500 7.4500 13.1500];
K = [0.050477 0.045107 0.039150 0.033639 0.029785 0.027017 0.024996 0.023522 0.022514];
imax = max(size(X0)); idx = int32((imax+1)/2);
T = [];
% find index
while(idx > 0 && idx <= imax)
x = X0(idx);
half = int32(idx/2);
if(ADU < x)
%printf("%g < %g\n", ADU, x);
if(idx == 1) break; endif; % less than first index
if(ADU > X0(idx-1)) idx -= 1; break; endif; % found
idx = half; % another half
else
%printf("%g > %g\n", ADU, x);
if(idx == imax) break; endif; % more than last index
if(ADU < X0(idx+1)) break; endif; % found
idx += half;
endif
endwhile
if(idx < 1) printf("too low (<%g)!\n", X0(1)); idx = 1; endif
if(idx > imax) idx = imax; endif;
T = Y0(idx) + K(idx) * (ADU - X0(idx));
endfunction