From 0603fcb5828c336944790a9680e8934d1c8b2d6b Mon Sep 17 00:00:00 2001 From: Kizarm Date: Mon, 27 Jan 2025 16:21:54 +0100 Subject: [PATCH] add termistor --- V203F6P6/termistor/Makefile | 53 +++++++++ V203F6P6/termistor/adc.cpp | 126 ++++++++++++++++++++ V203F6P6/termistor/adc.h | 20 ++++ V203F6P6/termistor/ch32v203 | 1 + V203F6P6/termistor/common | 1 + V203F6P6/termistor/compute/Makefile | 21 ++++ V203F6P6/termistor/compute/spline.cpp | 56 +++++++++ V203F6P6/termistor/compute/termistor.py | 147 ++++++++++++++++++++++++ V203F6P6/termistor/hack.c | 21 ++++ V203F6P6/termistor/main.cpp | 77 +++++++++++++ V203F6P6/termistor/main.h | 70 +++++++++++ V203F6P6/termistor/real_fix.h | 77 +++++++++++++ V203F6P6/termistor/spline.cpp | 56 +++++++++ V203F6P6/termistor/spline.h | 120 +++++++++++++++++++ 14 files changed, 846 insertions(+) create mode 100644 V203F6P6/termistor/Makefile create mode 100644 V203F6P6/termistor/adc.cpp create mode 100644 V203F6P6/termistor/adc.h create mode 120000 V203F6P6/termistor/ch32v203 create mode 120000 V203F6P6/termistor/common create mode 100644 V203F6P6/termistor/compute/Makefile create mode 100644 V203F6P6/termistor/compute/spline.cpp create mode 100755 V203F6P6/termistor/compute/termistor.py create mode 100644 V203F6P6/termistor/hack.c create mode 100644 V203F6P6/termistor/main.cpp create mode 100644 V203F6P6/termistor/main.h create mode 100644 V203F6P6/termistor/real_fix.h create mode 100644 V203F6P6/termistor/spline.cpp create mode 100644 V203F6P6/termistor/spline.h diff --git a/V203F6P6/termistor/Makefile b/V203F6P6/termistor/Makefile new file mode 100644 index 0000000..73c56a8 --- /dev/null +++ b/V203F6P6/termistor/Makefile @@ -0,0 +1,53 @@ +TARGET?= ch32v203 + +#TOOL ?= gcc +TOOL ?= clang + +PRJ = example + +VPATH = . ./$(TARGET) ./common +BLD = ./build/ +DFLAGS = -d +LFLAGS = -g +LDLIBS = +BFLAGS = --strip-unneeded + +CFLAGS = -MMD -Wall -Wno-parentheses -ggdb -fno-exceptions -ffunction-sections -fdata-sections +CFLAGS+= -I. -I./$(TARGET) -I./common +DEL = rm -f + +# zdrojaky +OBJS = main.o usart.o adc.o hack.o +#OBJS += + +include $(TARGET)/$(TOOL).mk +BOBJS = $(addprefix $(BLD),$(OBJS)) + +all: $(BLD) $(PRJ).elf +# ... atd. +-include $(BLD)*.d +# linker +$(PRJ).elf: $(BOBJS) + -@echo [LD $(TOOL),$(TARGET)] $@ + @$(LD) $(LFLAGS) -o $(PRJ).elf $(BOBJS) $(LDLIBS) + -@echo "size:" + @$(SIZE) $(PRJ).elf + -@echo "listing:" + $(DUMP) $(DFLAGS) $(PRJ).elf > $(PRJ).lst + -@echo "OK." + $(COPY) $(BFLAGS) -O binary $(PRJ).elf $(PRJ).bin +# preloz co je potreba +$(BLD)%.o: %.c + -@echo [CC $(TOOL),$(TARGET)] $@ + @$(CC) -std=gnu99 -c $(CFLAGS) $< -o $@ +$(BLD)%.o: %.cpp + -@echo [CX $(TOOL),$(TARGET)] $@ + @$(CXX) -std=c++17 -fno-rtti -c $(CFLAGS) $< -o $@ +$(BLD): + mkdir $(BLD) +flash: $(PRJ).elf + minichlink -w $(PRJ).bin flash -b +# vycisti +clean: + $(DEL) $(BLD)* *.lst *.bin *.elf *.map *~ +.PHONY: all clean flash run diff --git a/V203F6P6/termistor/adc.cpp b/V203F6P6/termistor/adc.cpp new file mode 100644 index 0000000..3c2f29e --- /dev/null +++ b/V203F6P6/termistor/adc.cpp @@ -0,0 +1,126 @@ +#include "system.h" +#include "oneway.h" +#include "adc.h" + +static AdcDma * pInstance = nullptr; + +extern "C" { + [[gnu::interrupt]] extern void DMA1_Channel1_IRQHandler(); +} +void DMA1_Channel1_IRQHandler( void ) { + DMA1_Type::INTFR_DEF state (DMA1.INTFR); + DMA1.INTFCR.R = state.R; // clear all + if (!pInstance) return; + if (state.B.HTIF1 != RESET) pInstance->send (false); + else if (state.B.TCIF1 != RESET) pInstance->send (true); +} + +static inline void EnableClock (void) noexcept { + // Enable DMA + RCC.AHBPCENR.modify([](RCC_Type::AHBPCENR_DEF & r) -> auto { + r.B.SRAMEN = SET; + r.B.DMA1EN = SET; + return r.R; + }); + // Enable ADC + GPIOA + RCC.APB2PCENR.modify([](RCC_Type::APB2PCENR_DEF & r) -> auto { + r.B.ADC1EN = SET; + r.B.IOPAEN = SET; + return r.R; + }); + RCC.APB1PCENR.B.TIM3EN = SET; // Enable TIM3 + RCC.CFGR0.B.ADCPRE = 3u; // PCLK2 divided by 8 as ADC clock (18 MHz, ! pretaktovano 14 MHz max). + // PIN PA0 / A0, PA1 / A1 + GPIOA.CFGLR.modify([](GPIOA_Type::CFGLR_DEF & r) -> auto { + r.B.MODE0 = 0u; + r.B.CNF0 = 0u; + r.B.MODE1 = 0u; + r.B.CNF1 = 0u; + return r.R; + }); +} +static inline void Timer3Init (uint32_t us) noexcept { + TIM3.PSC.R = 143u; // 1 MHz Fs + TIM3.ATRLR.R = us - 1u; + // TRGO update for ADC + TIM3.CTLR2.B.MMS = 2u; +} +static inline void AdcCalibrate (void) noexcept { + // RESET + RCC.APB2PRSTR.B.ADC1RST = SET; + RCC.APB2PRSTR.B.ADC1RST = RESET; + // set channels + ADC1.RSQR3__CHANNEL.modify([](ADC1_Type::RSQR3__CHANNEL_DEF & r) -> auto { + r.B.SQ1__CHSEL = 1u; // CH1 + //r.B.SQ2 = 1u; // CH1 + return r.R; + }); + ADC1.RSQR1.B.L = 0u; // 1 regular conversion + ADC1.SAMPTR2_CHARGE2.modify([](ADC1_Type::SAMPTR2_CHARGE2_DEF & r) -> auto { + r.B.SMP0_TKCG0 = 7u; + r.B.SMP1_TKCG1 = 7u; + return r.R; + }); + ADC1.CTLR1.modify([](ADC1_Type::CTLR1_DEF & r) -> auto { + r.B.SCAN = SET; + //r.B.BUFEN = SET; + //r.B.PGA = 3u; // x64 + return r.R; + }); + + ADC1.CTLR2.B.ADON = SET; + ADC1.CTLR2.B.RSTCAL = SET; // Launch the calibration by setting RSTCAL + while (ADC1.CTLR2.B.RSTCAL != RESET); // Wait until RSTCAL=0 + ADC1.CTLR2.B.CAL = SET; // Launch the calibration by setting CAL + while (ADC1.CTLR2.B.CAL != RESET); // Wait until CAL=0 +} +typedef __SIZE_TYPE__ size_t; +static inline void Dma1Ch1Init (void * ptr) noexcept { + // Configure the peripheral data register address + DMA1.PADDR1.R = reinterpret_cast (& ADC1.RDATAR_DR_ACT_DCG); + // Configure the memory address + DMA1.MADDR1.R = reinterpret_cast (ptr); + // Configure the number of DMA tranfer to be performs on DMA channel 1 + DMA1.CNTR1 .R = FULL_LEN; + // Configure increment, size, interrupts and circular mode + DMA1.CFGR1.modify([] (DMA1_Type::CFGR1_DEF & r) -> auto { + r.B.PL = 3u; // highest priority + r.B.MEM2MEM = RESET; // periferal -> memory + r.B.DIR = RESET; + r.B.MINC = SET; // memory increment + r.B.MSIZE = 1u; // 16-bit + r.B.PSIZE = 1u; // 16-bit + r.B.HTIE = SET; // INT Enable HALF + r.B.TCIE = SET; // INT Enable FULL + r.B.CIRC = SET; // Circular MODE + // Enable DMA Channel 1 + r.B.EN = SET; + return r.R; + }); +} +static inline void AdcPostInit (void) noexcept { + ADC1.CTLR2.modify([](ADC1_Type::CTLR2_DEF & r) -> auto { + r.B.DMA = SET; + r.B.EXTTRIG = SET; + r.B.EXTSEL = 4u; // TRGO event of timer 3 + r.B.SWSTART = SET; + return r.R; + }); +} +//////////////////////////////////////////////////////////////////////////////////// +AdcDma::AdcDma() noexcept : pL (buffer), pH (buffer + HALF_LEN), dst (nullptr) { + pInstance = this; + EnableClock (); + Timer3Init (1000u); + NVIC.EnableIRQ (DMA1_Channel1_IRQn); + AdcCalibrate(); + Dma1Ch1Init (buffer); + AdcPostInit (); + // start timer + TIM3.CTLR1.B.CEN = SET; +} +inline void AdcDma::send(const bool b) { + if (!dst) return; + if (b) dst->Send (pH, HALF_LEN); + else dst->Send (pL, HALF_LEN); +} diff --git a/V203F6P6/termistor/adc.h b/V203F6P6/termistor/adc.h new file mode 100644 index 0000000..a8166e6 --- /dev/null +++ b/V203F6P6/termistor/adc.h @@ -0,0 +1,20 @@ +#ifndef ADCDMA_H +#define ADCDMA_H +#include +#include "oneway.h" + +static constexpr unsigned HALF_LEN = 0x80u; +static constexpr unsigned FULL_LEN = HALF_LEN * 2u; + +class AdcDma { + uint16_t * pL; + uint16_t * pH; + uint16_t buffer [FULL_LEN]; + OneWay * dst; + public: + explicit AdcDma () noexcept; + void attach (OneWay & d) { dst = & d; } + void send (const bool b); +}; + +#endif // ADCDMA_H diff --git a/V203F6P6/termistor/ch32v203 b/V203F6P6/termistor/ch32v203 new file mode 120000 index 0000000..7650c85 --- /dev/null +++ b/V203F6P6/termistor/ch32v203 @@ -0,0 +1 @@ +../ch32v203/ \ No newline at end of file diff --git a/V203F6P6/termistor/common b/V203F6P6/termistor/common new file mode 120000 index 0000000..8332399 --- /dev/null +++ b/V203F6P6/termistor/common @@ -0,0 +1 @@ +../common/ \ No newline at end of file diff --git a/V203F6P6/termistor/compute/Makefile b/V203F6P6/termistor/compute/Makefile new file mode 100644 index 0000000..fe43aa5 --- /dev/null +++ b/V203F6P6/termistor/compute/Makefile @@ -0,0 +1,21 @@ +CC = gcc +CX = g++ +LD = ld +VPATH = . .. +CFLAGS = -c -Os -fPIC -I../ + +OBJS = spline.o +SLIB = spline.so + +all: $(SLIB) + +$(SLIB): $(OBJS) + $(CX) -shared $(OBJS) -o $(SLIB) +%.o: %.c + $(CC) $(CFLAGS) $< -o $@ +%.o: %.cpp + $(CX) $(CFLAGS) $< -o $@ +spline.o: spline.cpp ../spline.h ../real_fix.h + +clean: + rm -f $(OBJS) spline.so *.png diff --git a/V203F6P6/termistor/compute/spline.cpp b/V203F6P6/termistor/compute/spline.cpp new file mode 100644 index 0000000..e31f295 --- /dev/null +++ b/V203F6P6/termistor/compute/spline.cpp @@ -0,0 +1,56 @@ +#include +#include "spline.h" +extern "C" { + double Compute (double); + void plot (); +}; +static constexpr Pair measured [] = { + { +133.9, 125.000 }, + { +152.2, 120.000 }, + { +198.3, 110.000 }, + { +260.0, 100.000 }, + { +343.6, 90.000 }, + { +456.9, 80.000 }, + { +609.9, 70.000 }, +//{ +705.1, 65.000 }, + { +814.5, 60.000 }, +//{ +940.1, 55.000 }, + { +1082.6, 50.000 }, +//{ +1242.9, 45.000 }, + { +1421.2, 40.000 }, +//{ +1616.3, 35.000 }, +//{ +1826.2, 30.000 }, + { +2047.5, 25.000 }, +//{ +2275.6, 20.000 }, +//{ +2504.7, 15.000 }, + { +2729.0, 10.000 }, +//{ +2942.4, 5.000 }, + { +3139.8, 0.000 }, +//{ +3317.3, -5.000 }, + { +3472.7, -10.000 }, +//{ +3605.2, -15.000 }, + { +3715.6, -20.000 }, +//{ +3805.4, -25.000 }, + { +3876.9, -30.000 }, + { +3960.2, -38.000 }, + { +3976.0, -40.000 }, +}; + +static const SPLINE spline (measured, false); + +double Compute (double x) { + try { + const real r = spline.interpolate (to_real(x, 16)); + return from_real(r); + } catch (const char * e) { + fprintf(stderr, "%s at x = %f\n", e, x); + }; + return 0; +} +void plot () { + printf("Cubic Spline Coefficients - počet bodů : %zd\n", array_size (measured) - 1ul); + for (auto & e: spline) { + printf ("%12g: a = %+13g, b = %+13g, c = %+13g, d = %+13g\n", +e.x, +e.a, +e.b, +e.c, +e.d); + } + printf("half ~ %f\n", Compute(2047.45)); +} diff --git a/V203F6P6/termistor/compute/termistor.py b/V203F6P6/termistor/compute/termistor.py new file mode 100755 index 0000000..3299f49 --- /dev/null +++ b/V203F6P6/termistor/compute/termistor.py @@ -0,0 +1,147 @@ +#!/usr/bin/python3 +# -*- coding: utf-8 -*- +import os +import cffi +import numpy as np +import matplotlib.pyplot as plt +from math import exp,log + +c_header = '''double Compute (double x); +void plot (); +''' +def print_table (tab): + return + K = 4095.0 + R0 = 10.0 + s = 'static constexpr Pair measured [] = {\n' + l = len(tab) + r = range (0, l) + for n in r: + m = tab [l - n - 1] + R = m[1] + u = K * R / (R + R0) + s += ' {{ {:+6.1f},{:8.3f} }},\n'.format(u, m[0]) + s += '};\n' + print (s) + +def adc (R): + K0 = 4095.0 + R0 = 10.0 + return K0 * R / (R + R0) + +def R_NTC_T (t): + T = 273.15 + t + T0 = 298.15 + R0 = 10000.0; + B = 3977.0 + R = R0 * exp (B * (1.0 / T - 1.0 / T0)) + return 0.001 * R + +def beta_diff (table): + print ('t[°C] : R comp [kΩ] | R table [kΩ] | odchylka [%]') + for m in table: + t = m[0] + R = R_NTC_T (t) + D = 100.0 * (R - m[1]) / R + print('{:+5g} : {:8.3f} | {:8.3f} | {:+.2g}%'.format(t, R, m[1], D)) + +def ComuputeStHaCoeff (table): + e = table[ 5]; a = log (1000.0 * e[1]); f = 1.0 / (e[0] + 273.15) # 0 °C + e = table[ 8]; b = log (1000.0 * e[1]); g = 1.0 / (e[0] + 273.15) # 25°C + e = table[11]; c = log (1000.0 * e[1]); h = 1.0 / (e[0] + 273.15) # 50°C + e = table[13]; d = log (1000.0 * e[1]); i = 1.0 / (e[0] + 273.15) # 70°C + M = np.matrix ([[1.0, a, a**2, a**3], + [1.0, b, b**2, b**3], + [1.0, c, c**2, c**3], + [1.0, d, d**2, d**3]]) + I = M.getI(); # inverze matice + V = np.matrix ([f, g, h, i]) + X = I * V.transpose() + letters = ['A','B','C','D']; n=0; + X = X.tolist() + R = [] + print('Steinhart-Hart koeficienty :') + for x in X: + print(' {:s} = {:+20.14e}'.format(letters[n], x[0])) + n += 1 + R.append (x[0]) + return R + +def SteinhartHart (r, coe): + R = 1000.0 * r + A = coe[0]; B = coe[1]; C = coe[2]; D = coe[3]; + L = log(R); K = L; # Ku podivu koeficienty vypočítané ze 4 bodů tabulky + W = A # Steinhart-Hart metodou dávají lepší výsledky než + W += B * K; K *= L; # výrobcem udávaný koeficient beta. + W += C * K; K *= L; + W += D * K + return 1.0 / W - 273.15 # [°C] + +def stha_diff (table, c): + print ('rozdíl teplot vypočtený podle Steinhart-Hart formule :') + for m in table: + t = SteinhartHart(m[1], c) + print ('{:+5g} => {:g}'.format(m[0], t)) +def SHCompute (u, c): + K = 4095.0; R0 = 10.0 + R = u * R0 / (K - u) + return SteinhartHart(R, c) +def graph_err (tab, C,c): + #return + xs = np.arange(134.0, 4010.0, 1.0); ys = [] + for x in xs: ys.append (SHCompute(x,c) - C.Compute(x)) + xc = []; yc = [] + for m in tab: + x = adc (m[1]) + xc.append (x) + yc.append (SHCompute(x,c) - m[0]) + fig,ap = plt.subplots (1, figsize=(6.0,4.0)) + fig.suptitle ('Chyba NTC', fontsize=15) + ap.plot (xs, ys) + ap.plot (xc, yc, '.') + ap.legend (['Odchylka od vzorce', 'Body tabulky']) + plt.grid() + plt.xlabel('ADC [0-4095]') + plt.ylabel('teplota [°C]') + plt.savefig('err.png') + #plt.show() +def graph(tab, c): + #return + xs = []; ys = []; + for m in tab: + xs.append (adc (m[1])) + ys.append (m[0]) + xc = np.arange(100.0, 4050.0, 10.0); yc = []; + for x in xc: yc.append(SHCompute(x,c)) + fig,ap = plt.subplots (1, figsize=(6.0,4.0)) + fig.suptitle ('NTC 10k', fontsize=15) + ap.plot (xs, ys, linewidth=3) + ap.plot (xc, yc) + ap.legend (['Tabulka výrobce','Steinhart-Hart vzorce']) + plt.grid() + plt.xlabel('ADC [0-4095]') + plt.ylabel('teplota [°C]') + plt.savefig('ntc.png') + #plt.show() + +def main_func(): + table = [(-40,334.202),(-38,293.759),(-30,177.797),(-25,131.380),(-20,97.923),(-15,73.612),(-10,55.805), + (-5,42.655),(0,32.869),(5,25.527),(10,19.977),(15,15.750),(20,12.507),(25,10.0),(30,8.049), + (35,6.521),(40,5.315),(45,4.358),(50,3.594),(55,2.980),(60,2.483),(65,2.080),(70,1.750),(80,1.256), + (90,0.916),(100,0.678),(110,0.509),(120,0.386),(125,0.338)] + print_table (table) + c = ComuputeStHaCoeff(table) + #stha_diff (table, c) + ffi = cffi.FFI() + ffi.cdef(c_header) + C = ffi.dlopen("./spline.so") + C.plot() + graph_err (table, C,c) + ffi.dlclose (C) + graph (table, c) + return + beta_diff (table) + +if __name__ == '__main__': + main_func() + print ("END OK") diff --git a/V203F6P6/termistor/hack.c b/V203F6P6/termistor/hack.c new file mode 100644 index 0000000..391b7bf --- /dev/null +++ b/V203F6P6/termistor/hack.c @@ -0,0 +1,21 @@ +#include +#include +typedef __SIZE_TYPE__ size_t; +size_t strlen (const char *s) { + size_t l = 0; + while (*s++) l++; + return l; +} +void *memcpy (void *dest, const void *src, size_t n) { + const char *s = (const char *) src; + char *d = (char *) dest; + int i; + for (i=0; i termring; +static const GpioClass led (GPIOB, 8); +static Usart serial (57600); +static RealOut cout; +static AdcDma adc; +static Average avg (termring); +static const SPLINE spline (measured, false); + +int main () { + led << true; + adc.attach(avg); + cout += serial; + int passcnt = 0; + for (;;) { + uint32_t value; + if (termring.Read (value)) { + // value je de facto posunuta doleva o 7, potřebujeme o 16 - rozdíl je tedy 9 + const int32_t ivalue = static_cast(value << 9); + const real temperature = spline.interpolate(real(ivalue)); + cout << "t = " << temperature << " °C\r\n"; + const bool b = passcnt & 1; + led << b; + passcnt += 1; + } + } + return 0; +} diff --git a/V203F6P6/termistor/main.h b/V203F6P6/termistor/main.h new file mode 100644 index 0000000..762f70e --- /dev/null +++ b/V203F6P6/termistor/main.h @@ -0,0 +1,70 @@ +#ifndef MAIN_H_DEF +#define MAIN_H_DEF +#include "system.h" +#include "gpio.h" +#include "usart.h" +#include "adc.h" +#include "spline.h" +static constexpr unsigned FIFOLEN = 8u; +class Average : public OneWay { + FIFO & ring; + uint32_t y, w; + public: + explicit Average (FIFO & r) : OneWay(), ring(r), y(0u), w(0u) {} + unsigned int Send(uint16_t * const ptr, const unsigned int len) override { + unsigned suma = 0u; + for (unsigned n=0u; n> 4; // průměr s postupným zapomínáním. Hodnota je násobena délkou + // bufferu, t.j. 128 (posun doleva o 7) + if (w < (134 * 128)) return len; // skip limits - nezobrazuj mimo rozsah -40 až 120 + if (w > (3976 * 128)) return len; + ring.Write (w); + return len; + } +}; +class RealOut : public BaseLayer { // Zjednodušené výpisy teploty + static constexpr unsigned BUFLEN = 64u; + static constexpr char const * decstr = "0123456789"; + unsigned dec_points; + char buf[BUFLEN]; //!< Buffer pro výstup čísla. + public: + explicit RealOut (const unsigned dp = 2) noexcept : dec_points(dp) {} + RealOut & operator << (const char * str) { + uint32_t i = 0; + while (str[i++]); // strlen + BlockDown (str, --i); + return * this; + } + RealOut & operator << (const real num) { + int mf = 1; + for (unsigned n=0u; n> 16; bool s = false; + if (p < 0) { p = -p; s = true; } + int i = BUFLEN, j = 0; buf [--i] = '\0'; + for (;;) { + const int q = p % 10; p = p / 10; j++; + buf [--i] = decstr [q]; + if (j == dec_points) { buf [--i] = '.'; j++; } + if (!p and j > (dec_points + 1u)) { break; } + } + if (s) { buf [--i] = '-'; j++; } + BlockDown(buf + i, j); + return * this; + } + protected: + uint32_t BlockDown (const char * buf, uint32_t len) { + uint32_t n, ofs = 0, req = len; + for (;;) { + // spodní vrstva může vrátit i nulu, pokud je FIFO plné + n = BaseLayer::Down (buf + ofs, req); + ofs += n; // Posuneme ukazatel + req -= n; // Zmenšíme další požadavek + if (!req) break; + } + return ofs; + } +}; +#endif // MAIN_H_DEF diff --git a/V203F6P6/termistor/real_fix.h b/V203F6P6/termistor/real_fix.h new file mode 100644 index 0000000..44d4bc5 --- /dev/null +++ b/V203F6P6/termistor/real_fix.h @@ -0,0 +1,77 @@ +#ifndef _REAL_FIX_H +#define _REAL_FIX_H +#include +namespace FIX { + static constexpr int iround (const double d) { return d > 0.0 ? int (d + 0.5) : int (d - 0.5); } + /** + * @class real + * @brief Aritmetika v pevné řádové čárce. + * + * Není to moc přesné, ale bez float koprocesoru to postačí. + * A není to kompletní, jen pro výpočet polynomu, zde to stačí. + */ + struct real { + int32_t x; // necháme veřejné + /** + * @brief Konstruktor z double + * @param _x double + */ + explicit constexpr real (const double _x) : x (iround (_x)) {} + /** + * @brief Konstruktor z celého čísla, přímo nepoužívat - explicitní převod from_int() + * Ono se to trochu pere s tím double. + * @param _x integer + */ + explicit constexpr real (const int32_t _x) : x (_x) {} + /** + * @brief Kopírovací konstruktor + */ + explicit constexpr real () : x (0) {} + // stačí přetížit jen některé operátory, výpočet je dělán poměrně úsporně + + real operator+ (const real & r) const { + return real (x + r.x); + } + real operator- (const real & r) const { + return real (x - r.x); + } + + bool operator< (const real & r) const { + return x < r.x; + } + bool operator>= (const real & r) const { + return x >= r.x; + } + real & operator+= (const real & r) { + x += r.x; + return * this; + } + // Pouze pro výpisy - zachová formátovou specifikaci pro FLT + double operator+ () const { + return double (x); + } + }; + // Násobení s posunem zde používá "plovoucí" posun, jinak to nejde. + static void mull (real & res, const real & r, const int sl = 0) { + const int64_t result = (int64_t)(res.x) * (int64_t)(r.x); +#ifdef __linux__ + // Check overflow + const uint64_t u = result < 0 ? -result : +result; + const uint32_t o = u >> (sl + 32); + if (o) {fprintf (stderr, "Overflow 0x%X\n", o); // Není to dokonalé, ale alespoň něco + // Šla by vyhodit výjimka a backtrace zjistit proč a hlavně kde. + throw ("overflow exception"); } +#endif + res.x = result >> sl; + } + //! Vytvoří double z real + static constexpr double from_real (const real & r) { + const double a = 1.0 / double (1ull << 16); + return a * double(r.x); + } + static constexpr real to_real (const double x, const int sh) { + const double a = double (1ull << sh); + return real (a * x); + } +}; // end of namespace FIX +#endif // _REAL_FIX_H diff --git a/V203F6P6/termistor/spline.cpp b/V203F6P6/termistor/spline.cpp new file mode 100644 index 0000000..e31f295 --- /dev/null +++ b/V203F6P6/termistor/spline.cpp @@ -0,0 +1,56 @@ +#include +#include "spline.h" +extern "C" { + double Compute (double); + void plot (); +}; +static constexpr Pair measured [] = { + { +133.9, 125.000 }, + { +152.2, 120.000 }, + { +198.3, 110.000 }, + { +260.0, 100.000 }, + { +343.6, 90.000 }, + { +456.9, 80.000 }, + { +609.9, 70.000 }, +//{ +705.1, 65.000 }, + { +814.5, 60.000 }, +//{ +940.1, 55.000 }, + { +1082.6, 50.000 }, +//{ +1242.9, 45.000 }, + { +1421.2, 40.000 }, +//{ +1616.3, 35.000 }, +//{ +1826.2, 30.000 }, + { +2047.5, 25.000 }, +//{ +2275.6, 20.000 }, +//{ +2504.7, 15.000 }, + { +2729.0, 10.000 }, +//{ +2942.4, 5.000 }, + { +3139.8, 0.000 }, +//{ +3317.3, -5.000 }, + { +3472.7, -10.000 }, +//{ +3605.2, -15.000 }, + { +3715.6, -20.000 }, +//{ +3805.4, -25.000 }, + { +3876.9, -30.000 }, + { +3960.2, -38.000 }, + { +3976.0, -40.000 }, +}; + +static const SPLINE spline (measured, false); + +double Compute (double x) { + try { + const real r = spline.interpolate (to_real(x, 16)); + return from_real(r); + } catch (const char * e) { + fprintf(stderr, "%s at x = %f\n", e, x); + }; + return 0; +} +void plot () { + printf("Cubic Spline Coefficients - počet bodů : %zd\n", array_size (measured) - 1ul); + for (auto & e: spline) { + printf ("%12g: a = %+13g, b = %+13g, c = %+13g, d = %+13g\n", +e.x, +e.a, +e.b, +e.c, +e.d); + } + printf("half ~ %f\n", Compute(2047.45)); +} diff --git a/V203F6P6/termistor/spline.h b/V203F6P6/termistor/spline.h new file mode 100644 index 0000000..ae90eb9 --- /dev/null +++ b/V203F6P6/termistor/spline.h @@ -0,0 +1,120 @@ +#ifndef SPLINE_H +#define SPLINE_H +/** @file + * @class SPLINE + * @brief Kubické splajny. + */ +#include "real_fix.h" +using namespace FIX; +typedef __SIZE_TYPE__ size_t; + +// Zjištění počtu prvků statického pole. +templateconstexpr size_t array_size (T (&) [N]) { return N; } +/*! + * @class Pair + * @brief Prvek vstupní tabulky pro výpočet kubických splajnů. + * */ +struct Pair { + double x, y; +}; +/*! @class SplineSet + @brief Výsledné koeficienty pro výpočet iterace pomocí polynomu 3. řádu. +*/ +struct SplineSet { + real x; //!< začátky intervalů x + real a; //!< záčátky intervalů y, koeficienty a + real b; //!< koeficienty b (x^1) + real c; //!< koeficienty c (x^2) + real d; //!< koeficienty d (x^3) +}; +//! V podstatě statický vektor pro SplineSet (bez nutnosti dymanické alokace) +template class SPLINE { + SplineSet data [N - 1]; // zde jsou schovány koeficienty (privátně) + public: + //! Konstruktor + explicit constexpr SPLINE (const Pair * const p, const bool reverse = false) noexcept { + double x [N], y [N]; // přeskupení dat - možnost inverzní funkce při reverse = true + if (reverse) {for (int i=0; i= 0; --j) { // zpětný chod + c[j] = z [j] - mu[j] * c[j+1]; + b[j] = (y[j+1] - y[j]) / h[j] - h[j] * (c[j+1] + 2*c[j]) / 3.0; + d[j] = (c[j+1] - c[j]) / (3.0 * h[j]); + } + for (int i = 0; i < n; ++i) { // závěrečný zápis koeficientů + // Pro FIX je potřeba nastavit různé posuny pro koeficienty, rozsah je příliš velký a pak by to bylo nepřesné. + data[i].x = to_real(x[i], 16); data[i].a = to_real(y[i], 16); + data[i].b = to_real(b[i], 20); data[i].c = to_real(c[i], 28); data[i].d = to_real(d[i], 36); + } + } + const SplineSet & operator[] (const int index) const { + return data [index]; + } + const size_t size () const { + return N - 1; + } + const real interpolate (const real x) const { + const unsigned n = range (x); // určíme interval ve kterém budeme počítat + const SplineSet & s = data [n]; // koeficienty polynomu v tomto intervalu + const real r = x - s.x; // souřadnice x v tomto intervalu + return cubic_poly (s, r); // vlastní výpočet + } + /***************************************************************/ + /** @class iterator + * @brief pro range-based for () */ + class iterator { + const SplineSet * ptr; + public: + iterator(const SplineSet * _ptr) : ptr (_ptr) {} + iterator operator++ () { ++ptr; return * this; } + bool operator!= (const iterator & other) const { return ptr != other.ptr; } + const SplineSet & operator* () const { return * ptr; } + }; + iterator begin () const { return iterator (data ); } + iterator end () const { return iterator (data + N - 1); } + protected: + const unsigned range (const real x) const { + int l=0, r=size() - 1u; + while (l <= r) { + const int s = (l + r) >> 1; + const real x0 = data[s + 0].x; + const real x1 = data[s + 1].x; + if (x < x0) r = s - 1; + else if (x >= x1) l = s + 1; + else return s; + } + return size() - 1u; + } + const real cubic_poly (const SplineSet & s, const real x) const { + real r(s.d); // d => 36 v konstruktoru + mull (r, x, 24); // 36+16-24 = 28 == c + r += s.c; // c => 28 + mull (r, x, 24); // 28+16-24 = 20 == b + r += s.b; // b => 20 + mull (r, x, 20); // 20+16-20 = 16 == a + r += s.a; // a=x => 16 + return r; // r=y => 16 + } +}; + +#endif // SPLINE_H