Kizarm/RISC-V
1
1
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity

add termistor

Browse source
This commit is contained in:
Kizarm 2025-01-27 16:21:54 +01:00
parent b122c8a79b
commit 0603fcb582
14 changed files with 846 additions and 0 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

53
V203F6P6/termistor/Makefile Normal file
View file

@ -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

126
V203F6P6/termistor/adc.cpp Normal file
View file

@ -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<size_t> (& ADC1.RDATAR_DR_ACT_DCG);
// Configure the memory address
DMA1.MADDR1.R = reinterpret_cast<size_t> (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);
}

20
V203F6P6/termistor/adc.h Normal file
View file

@ -0,0 +1,20 @@
#ifndef ADCDMA_H
#define ADCDMA_H
#include <stdint.h>
#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<uint16_t> * dst;
public:
explicit AdcDma () noexcept;
void attach (OneWay<uint16_t> & d) { dst = & d; }
void send (const bool b);
};
#endif // ADCDMA_H

1
V203F6P6/termistor/ch32v203 Symbolic link
View file

@ -0,0 +1 @@
../ch32v203/

1
V203F6P6/termistor/common Symbolic link
View file

@ -0,0 +1 @@
../common/

21
V203F6P6/termistor/compute/Makefile Normal file
View file

@ -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

56
V203F6P6/termistor/compute/spline.cpp Normal file
View file

@ -0,0 +1,56 @@
#include <cstdio>
#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<array_size(measured)> 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));
}

147
V203F6P6/termistor/compute/termistor.py Executable file
View file

@ -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")

21
V203F6P6/termistor/hack.c Normal file
View file

@ -0,0 +1,21 @@
#include <stdint.h>
#include <stdarg.h>
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<n; i++) d[i] = s[i];
return dest;
}
void *memset (void *s, int c, size_t n) {
char *p = (char *) s;
int i;
for (i=0; i<n; i++) p[i] = c;
return s;
}

77
V203F6P6/termistor/main.cpp Normal file
View file

@ -0,0 +1,77 @@
#include "main.h"
/*********************************************************************************
* Měření teploty s nějakou lepší přesností není zase tak jednoduché. Nakonec byl
* jako čidlo zvolen NTC termistor polské výroby, kde výrobce zaručuje přesnost
* odporu a beta parametru 1%. Ten beta parametr je celkem použitelný v malém rozsahu
* teplot (10 - 70°C), z tabulky se ale dají spočítat koeficienty Steinhart-Hartových
* vzorců, které sedí ve větším rozsahu (viz adresář compute).
* Vlastní měření je celkem jednoduché - stačí udělat odporový dělič z napájecího
* napětí na zem pomocí přesného rezistoru (zde také 1%) 10k a NTC 10k (při 25°C)
* a jeho střed připojit na vstup ADC. Pak lze spočítat závislost teploty na
* hodnotě ADC jako funkci, pro inverzní funkci udělat aproximaci pomocí kubických splajnů
* a máme fakticky hotovo. Výhoda je, že to nezávisí na napájecím (referenčním) napětí.
*
* Protože tento procesor nemá koprocesor pohyblivé řádové čárky, dodělal jsem do toho
* potřebnou aritmetiku v pevné čárce (tedy v celých číslech). Protože koeficienty
* polynomů se liší o mnoho řádů, bylo potřeba udělat bitové posuny poměrně variabilní.
* A celé to otestovat, jestli to někde nepřetéká. Pak lze zobazit teplotu na RS485
* na dvě desetinná místa a protože se teplota mění pomalu a je použit plovoucí průměr,
* zobrazení je stabilní.
*
* Aniž bych to nějak kalibroval, zdá se, že teplota celkem sedí na +/- 0.5°C.
* Ve vodě s ledem to ukazuje 0.50°C, moje tělesná teplota je pak 36.6°C.
*
* NOTE : Budič RS485 je pětivoltový, takže jeho výstup Rx data jsem připojil pro jistotu
* k procesoru přes odpor 1k. S tím, že záchytné diody přepětí nějak požerou.
* Ukazuje se, že to sice funguje, ale ADC pak dost kecá - nějak mu tam to zvýšené
* napětí proniká. Stačí pak přizemnit Rx procesoru (PA3) přes odpor 1.5 až 2k a problém
* zmizí.
**********************************************************************************/
static constexpr Pair measured [] = { // pár hodnota ADC, teplota ve °C
{ +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 },
{ +814.5, 60.000 },
{ +1082.6, 50.000 },
{ +1421.2, 40.000 },
{ +2047.5, 25.000 },
{ +2729.0, 10.000 },
{ +3139.8, 0.000 },
{ +3472.7, -10.000 },
{ +3715.6, -20.000 },
{ +3876.9, -30.000 },
{ +3960.2, -38.000 },
{ +3976.0, -40.000 },
};
/*********************************************************************************/
static FIFO<uint32_t, FIFOLEN> termring;
static const GpioClass led (GPIOB, 8);
static Usart serial (57600);
static RealOut cout;
static AdcDma adc;
static Average avg (termring);
static const SPLINE<array_size(measured)> 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<int32_t>(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;
}

70
V203F6P6/termistor/main.h Normal file
View file

@ -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<uint16_t> {
FIFO<uint32_t, FIFOLEN> & ring;
uint32_t y, w;
public:
explicit Average (FIFO<uint32_t, FIFOLEN> & 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<len; n++) {
suma += ptr [n];
}
y += suma - w; // ustálení teploty trvá dost dlouho, takže lze posílat klouzavý
w = y >> 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<dec_points; n++) mf *= 10;
int p = (mf * num.x) >> 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

77
V203F6P6/termistor/real_fix.h Normal file
View file

@ -0,0 +1,77 @@
#ifndef _REAL_FIX_H
#define _REAL_FIX_H
#include <stdint.h>
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

56
V203F6P6/termistor/spline.cpp Normal file
View file

@ -0,0 +1,56 @@
#include <cstdio>
#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<array_size(measured)> 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));
}

120
V203F6P6/termistor/spline.h Normal file
View file

@ -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.
template<class T, size_t N>constexpr 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<const int N> 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<N; i++) { x[i] = p[i].y; y[i] = p[i].x; }}
else {for (int i=0; i<N; i++) { x[i] = p[i].x; y[i] = p[i].y; }}
/* Tenhle příšerně složitý konstruktor je převzat ze stackoverflow.com (původní odkaz už zmizel)
* Není to zas taková sranda - koeficienty se počítají řešením tridiagonální matice řádu N.
* V těch indexech se snadno zabloudí, tohle kupodivu funguje (přirozené splajny) ač je to dost krátké.
* Je zajímavé, že tohle najdete spíš ve Fortranu než v C, matematici jsou zřejmě hodně konzervativní.
* */
const int n = N - 1;
double h [n];
for (int i = 0; i < n; ++i) { h[i] = x[i+1]-x[i]; }
double alpha [n]; alpha [0] = 0.0;
for (int i = 1; i < n; ++i) {
alpha[i] = 3.0 * (y[i+1] - y[i]) / h[i] - 3.0 * (y[i] - y[i-1]) / h[i-1];
}
double c [n+1], l [n+1], mu [n+1], z [n+1]; l [0] = 1.0; mu[0] = 0.0; z [0] = 0.0;
for (int i = 1; i < n; ++i) { // přímý chod
l [i] = 2.0 * (x[i+1] - x[i-1]) - h[i-1] * mu[i-1];
mu[i] = h[i] / l[i];
z [i] = (alpha[i] - h[i-1] * z[i-1]) / l[i];
}
l[n] = 1.0; z[n] = 0.0; c[n] = 0.0;
double b [n], d [n];
for (int j = n-1; j >= 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