Кусочно-линейная аппроксимация

[eddy_em]

Вчера я занимался построением зависимости отчетов в 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