- Linear calibration of a pressure sensor
- Drift detection of a humidity sensor
- Temperature uncertainty of a PT100 probe
- Remove mains hum from a sensor by FFT
- Battery life of a LoRaWAN sensor node
- Health report of an IoT gateway
- Decode and acknowledge a Modbus status register
- LoRaWAN link quality report of a sensor
- A/B test of a firmware update on the radio latency
- 10-minute wind report of a weather station
Linear calibration of a pressure sensor
Raw ADC counts of a pressure sensor are recorded against a reference gauge. The script fits the gain and offset, checks the linearity and converts a new reading into bar.
#1. Raw ADC counts of a pressure sensor against a reference gauge in bar; -> "[adc]" "[412, 1235, 2050, 2871, 3690]"; -> "[ref]" "[0.00, 2.51, 5.00, 7.52, 10.01]"; #2. Linear calibration p = a + b*x; json load "rx_iot_cal" (stat fit [adc] [ref] "linear"); -> "[expr]" (json select "rx_iot_cal" /expression); #3. Convert a new reading of 3000 counts; -> "[p]" (math eval [expr] "{\"x\": 3000}"); concat "offset a = " (math round_sig (json select "rx_iot_cal" /a) 4) " bar, gain b = " (math round_sig (json select "rx_iot_cal" /b) 4) " bar/count, r2 = " (math round (json select "rx_iot_cal" /r2) 6) ", 3000 counts = " (math round [p] 3) " bar"
offset a = -1.260 bar, gain b = 0.003055 bar/count, r2 = 0.999998, 3000 counts = 7.906 bar
- Use
stat fit_modelto get the standard errors of the gain and the offset - Fit a second order polynomial with
pa polynomial_curve_fit_get_coefffor non-linear sensors - Compute the non-linearity as the largest residual in percent of the full scale
Drift detection of a humidity sensor
A humidity sensor is compared every day with a reference. The script smooths the offsets, tests for a monotonic trend and estimates the drift per month to decide a recalibration.
#1. Daily offset of a humidity sensor against a reference, in percent RH; -> "[d]" "[0.1, 0.3, 0.2, 0.4, 0.3, 0.6, 0.5, 0.7, 0.9, 0.8, 1.0, 1.1, 1.0, 1.3]"; #2. Smooth the daily noise with a 3-day moving average; -> "[m]" (dsp moving_average [d] 3); #3. Monotonic trend test and slope of the drift; json load "rx_iot_k" (stat kendall "[1,2,3,4,5,6,7,8,9,10,11,12,13,14]" [d]); json load "rx_iot_f" (stat fit "[1,2,3,4,5,6,7,8,9,10,11,12,13,14]" [d] "linear"); -> "[drift]" (math calc (concat "30*" (json select "rx_iot_f" /b))); -> "[p]" (json select "rx_iot_k" /p_value); concat "smoothed offset from " (math round (stat min [m]) 2) " to " (math round (stat max [m]) 2) " %RH, Kendall p-value = " (math round_sig [p] 3) ", drift = " (math round [drift] 2) " %RH per 30 days -> " (if (and (< [p] 0.05) (> [drift] 1)) {"recalibrate"} {"OK"})
smoothed offset from 0.2 to 1.13 %RH, Kendall p-value = 0.0000144, drift = 2.58 %RH per 30 days -> recalibrate
- Use
stat spearmanas another trend test - Remove the trend with
dsp detrendand check the remaining noise withstat stddev - Run the check every night with a job on the stored readings
Temperature uncertainty of a PT100 probe
A PT100 probe gives 109.73 ohm. The script converts it into a temperature, propagates the uncertainties of R and R0 and writes the result with an expanded uncertainty k = 2.
#1. PT100 probe: R measured with its uncertainty, R0 and alpha of the class; -> "[vars]" "{\"R\": [109.73, 0.02], \"R0\": [100, 0.01], \"alpha\": 0.00385}"; #2. Propagate the uncertainty through T = R/R0 - 1 / alpha; json load "rx_iot_u" (num uncertainty "(R/R0-1)/alpha" [vars]); #3. Expanded uncertainty k = 2, written with the usual rounding; -> "[T]" (json select "rx_iot_u" /value); -> "[U]" (math calc (concat "2*" (json select "rx_iot_u" /uncertainty))); concat "T = " (math round_uncertainty [T] [U]) " degC (k = 2), share of R0 in the variance: " (math round (math calc (concat "100*" (json select "rx_iot_u" /variance_share/R0))) 0) " %"
T = 25.27 ± 0.12 degC (k = 2), share of R0 in the variance: 23 %
- Use the Callendar-Van Dusen equation with
num solveabove 0 degC - Add the self-heating error as another variable
- Read
variance_shareto find the dominant source of uncertainty
Remove mains hum from a sensor by FFT
A vibration sensor picks up 50 Hz mains hum and its 150 Hz harmonic. With 1024 samples at 1024 Hz, the hum falls exactly into FFT bins: the script zeroes them, goes back with the inverse FFT, smooths the noise with Savitzky-Golay and draws the spectrum of the cleaned signal.
#1. Vibration sensor sampled at 1024 Hz for 1 s: slow 1.5 Hz motion, 50 Hz mains hum and its 150 Hz harmonic, a little noise; -> "[x]" (vector add (matrix col (num table "2+0.5*sin(2*pi*1.5*t)+0.3*sin(2*pi*50*t)+0.1*sin(2*pi*150*t)" "t" 0 0.9990234375 0.0009765625) 2) (stat dist_random "normal" "[0, 0.05]" 1024)); #2. FFT: 1 Hz per bin, so the hum is exactly in bins 50 and 150 and their mirrors 974 and 874; json load "rx_X" (concat "{\"X\": " (num fft [x]) "}"); json uarray "rx_X" /X 50 "0" STR; json uarray "rx_X" /X 974 "0" STR; json uarray "rx_X" /X 150 "0" STR; json uarray "rx_X" /X 874 "0" STR; #3. Back to the time domain, then Savitzky-Golay smoothing of the noise; -> "[clean]" (dsp savgol (num ifft (json select "rx_X" /X)) 21 3); #4. Spectrum of the cleaned signal: no more line at 50 and 150 Hz; plot spectrum [clean] 1024
The chart is shown in the Plot tab of the editor.
- Do the same with
dsp filterbandstop around 50 Hz when the sampling does not match the bins. - Draw raw and cleaned signals together with
plot chart.
Battery life of a LoRaWAN sensor node
A LoRaWAN node sleeps at 4 uA, measures for 50 ms every 10 s and transmits at SF12 for 400 ms every 5 minutes. Compute the average current, the battery life on two derated AA lithium cells and the date of the replacement visit, moved to a working day.
#1. LoRaWAN node: sleep 4 uA, measure 12 mA for 50 ms every 10 s, radio SF12 45 mA for 400 ms every 5 min; -> "[sleep]" (unit calc "4 uA" "mA"); -> "[meas]" (math calc "12*0.05/10"); -> "[tx]" (math calc "45*0.4/300"); #2. Average current of the duty cycle; -> "[iavg]" (math calc "[sleep]+[meas]+[tx]"); #3. Two AA lithium cells of 3000 mAh in parallel, 80 pct usable after derating; -> "[hours]" (math calc "2*3000*0.8/[iavg]"); -> "[days]" (math floor (/ [hours] 24)); #4. Replacement date from the installation day, moved to a working day; -> "[install]" "2026-11-02"; -> "[change]" (date add [install] "DAY" [days]); if (== (date is_weekend [change]) 1) {-> "[change]" (date next_day [change] "monday");}; concat "average current " (math round (* 1000 [iavg]) 1) " uA | battery life " [days] " days = " (math round (/ [days] 365.25) 2) " years | replacement visit on " [change] ", ISO week " (date week_of_year [change]) ", " (date nb_business_days [install] [change]) " business days after installation"
average current 124 uA | battery life 1612 days = 4.41 years | replacement visit on 2031-04-02, ISO week 14, 1152 business days after installation
- Compare SF7 and SF12 airtimes: the radio term drops by about 20.
- Use
date add_business_daysto plan a check visit every 60 working days.
Health report of an IoT gateway
An edge gateway running MentDB must publish a small health message every few minutes. The script reads the live CPU and available memory, extracts the number of cores and the disk usage from the system report, derives an anonymous gateway id from the hostname and builds the JSON status message.
#1. Live resources; -> "[cpu]" (metric current cpu jvm); -> "[mem]" (metric current available mem); #2. Cores and disk usage of the root file system, from the system report; -> "[sys]" (metric system); -> "[cores]" (string regexp_extract [sys] "nb_proc: ([0-9]+)" 1); -> "[disk]" (string regexp_extract [sys] "\[/\] ([0-9]+)%" 1); #3. Anonymous gateway id: first 8 hex digits of the SHA of the hostname; -> "[gid]" (string left (string sha (os hostname)) 8); #4. Status message; json load "rx_p4f_gw" "{}"; json iobject "rx_p4f_gw" / "gateway" (concat "gw-" [gid]) STR; json iobject "rx_p4f_gw" / "cores" [cores] NUM; json iobject "rx_p4f_gw" / "history_days" (string regexp_count (metric date) "[0-9]{4}-[0-9]{2}-[0-9]{2}") NUM; json iobject "rx_p4f_gw" / "status" (if (or (> [cpu] 80) (or (< [mem] 500) (> [disk] 90))) {"degraded"} {"ok"}) STR; json iobject "rx_p4f_gw" / "checks" (concat "cpu " (if (> [cpu] 80) {"high"} {"normal"}) ", memory " (if (< [mem] 500) {"low"} {"normal"}) ", disk " (if (> [disk] 90) {"full"} {"normal"})) STR; json doc "rx_p4f_gw"
{
"cores": 11.0,
"checks": "cpu normal, memory normal, disk normal",
"history_days": 31.0,
"gateway": "gw-eedfe8d8",
"status": "ok"
}- Add the free swap and the JVM heap with
metric current used mem jvm - Publish the message every 5 minutes with a job and keep the last 31 days with
metric date
Decode and acknowledge a Modbus status register
A field sensor exposes a 16-bit Modbus status register (read as hex 02C5). The script decodes the active flags with bit masks, acknowledges the two alarm bits with an XOR, sets the maintenance-mode bit with an OR and returns the register to write back in hex and binary.
#1. Register read on the bus and the meaning of its bits; -> "[reg]" (string hex_to_int "02C5"); json load "rx_p4f_bits" "[[\"running\", 0], [\"low_battery\", 2], [\"over_temp\", 3], [\"sensor_fault\", 6], [\"comm_warning\", 7], [\"door_open\", 9]]"; -> "[active]" ""; #2. Active flags; for (-> "[i]" 0) (< [i] (json count "rx_p4f_bits" /)) (++ "[i]") { if (> (math bit_and [reg] (math power 2 (json select "rx_p4f_bits" (concat "/[" [i] "]/[1]")))) 0) { -> "[active]" (concat [active] (json select "rx_p4f_bits" (concat "/[" [i] "]/[0]")) " "); }; }; json unload "rx_p4f_bits"; #3. Acknowledge the alarms of bits 6 and 7, then set maintenance mode bit 12; -> "[ack]" (math bit_xor [reg] (math bit_and [reg] 192)); -> "[new]" (math bit_or [ack] 4096); concat "0x02C5 = " (string lpad (string bin [reg]) "0" 16) " | active: " (string trim [active]) " | write back 0x" (string upper (string lpad (string int_to_hex [new]) "0" 4)) " = " (string lpad (string bin [new]) "0" 16)
0x02C5 = 0000001011000101 | active: running low_battery sensor_fault comm_warning door_open | write back 0x1205 = 0001001000000101
- Compare two successive reads with
math bit_xorto log only the flags that changed - Write the register back with a Modbus gateway REST call
LoRaWAN link quality report of a sensor
The network server gives the spreading factor of the last 20 uplinks of a sensor and the frame counters. The script builds the frequency table of the spreading factors, finds the most used one, computes the effective bit rate as the harmonic mean of the per-packet rates (equal payloads) and the packet delivery ratio.
#1. SF of the last 20 uplinks received, frame counters 1180 to 1203; -> "[sf]" "[9, 9, 10, 9, 8, 9, 9, 10, 10, 9, 7, 9, 9, 8, 9, 12, 10, 9, 9, 9]"; #2. Frequency table and share of the uplinks at SF 9 or less; json load "rx_p4f_freq" (stat frequency [sf]); -> "[sf9]" (json select "rx_p4f_freq" "/[2]/cumulative_frequency"); json unload "rx_p4f_freq"; #3. Bit rate of each uplink, BW 125 kHz, coding rate 4/5; json load "rx_p4f_rate" "[]"; json load "rx_p4f_sf" [sf]; json parse_array "rx_p4f_sf" "/" "[s]" { json iarray "rx_p4f_rate" / (* [s] (* (/ 125000 (math power 2 [s])) 0.8)) NUM; }; #4. Effective rate for equal payloads = harmonic mean; -> "[hm]" (stat harm_mean (json doc "rx_p4f_rate")); -> "[am]" (stat mean (json doc "rx_p4f_rate")); json unload "rx_p4f_rate"; json unload "rx_p4f_sf"; concat "SF mode " (stat mode [sf]) ", SF <= 9 in " (math round (* 100 [sf9]) 1) " pct | effective rate " (math round [hm] 1) " bit/s (arithmetic mean " (math round [am] 1) ") | PDR " (math round (/ (* 100 20) (- 1203 1179)) 1) " pct"
SF mode [9], SF <= 9 in 75 pct | effective rate 1319.4 bit/s (arithmetic mean 1850.6) | PDR 83.3 pct
- Enable ADR when SF 9 or less is used in more than 90 pct of the uplinks
- Compute the time on air per SF and the duty cycle used in the 1 pct sub-band
A/B test of a firmware update on the radio latency
A new firmware is deployed on part of the fleet. The script compares the end-to-end latencies of the two groups with the non-parametric Mann-Whitney test, then checks with a z-test (known jitter 4 ms from the radio datasheet) whether the new firmware meets the 125 ms target.
#1. Latencies in ms of the old and new firmware groups; -> "[old]" "[128, 135, 122, 140, 131, 126, 138, 133]"; -> "[new]" "[118, 121, 115, 129, 119, 123, 117, 120, 124]"; #2. Non-parametric comparison of the two groups; json load "rx_p4f_mw" (stat mann_whitney [old] [new]); -> "[p]" (json select "rx_p4f_mw" "/p_value"); -> "[method]" (json select "rx_p4f_mw" "/method"); #3. z-test of the new group against the 125 ms target, sigma known = 4 ms; json load "rx_p4f_z" (stat ztest [new] 125 4); -> "[z]" (json select "rx_p4f_z" "/z"); -> "[pz]" (json select "rx_p4f_z" "/p_value"); json unload "rx_p4f_mw"; json unload "rx_p4f_z"; concat "median old " (stat median [old]) " ms, new " (stat median [new]) " ms | Mann-Whitney " [method] " p = " (math round_sig [p] 3) (if (< [p] 0.05) {" -> significant change"} {" -> no proof of change"}) " | z = " (math round [z] 3) ", p = " (math round_sig [pz] 3) (if (and (< [z] 0) (< [pz] 0.05)) {" -> mean below 125 ms: roll out"} {" -> target not proven"})
median old 132 ms, new 120 ms | Mann-Whitney exact p = 0.00156 -> significant change | z = -3.25, p = 0.00115 -> mean below 125 ms: roll out
- Use
stat bootstrapstyle resampling of the medians to give a confidence interval - Check the latency distribution of each group with
stat frequency
10-minute wind report of a weather station
A weather station measures the wind direction every minute around north, where a simple average fails (350 and 10 deg average to 180). The script computes the vector mean direction with atan2, the Yamartino standard deviation of the direction and the resultant wind speed, as in WMO reports.
#1. Directions in deg and speeds in m/s, one per minute; json load "rx_p4f_dir" "[350, 355, 5, 10, 358, 2, 15, 345, 8, 0]"; json load "rx_p4f_v" "[6.1, 5.8, 6.4, 7.0, 6.2, 5.9, 7.3, 5.5, 6.6, 6.0]"; json load "rx_p4f_s" "[]"; json load "rx_p4f_c" "[]"; json load "rx_p4f_u" "[]"; json load "rx_p4f_w" "[]"; #2. Unit vector components and wind vector components; for (-> "[i]" 0) (< [i] 10) (++ "[i]") { -> "[a]" (math calc (concat (json select "rx_p4f_dir" (concat "/[" [i] "]")) "*pi/180")); -> "[v]" (json select "rx_p4f_v" (concat "/[" [i] "]")); json iarray "rx_p4f_s" / (math sin [a]) NUM; json iarray "rx_p4f_c" / (math cos [a]) NUM; json iarray "rx_p4f_u" / (* [v] (math sin [a])) NUM; json iarray "rx_p4f_w" / (* [v] (math cos [a])) NUM; }; #3. Vector mean direction and Yamartino sigma; -> "[sa]" (math avg (json doc "rx_p4f_s")); -> "[ca]" (math avg (json doc "rx_p4f_c")); -> "[dir]" (math calc (concat (math atan2 [sa] [ca]) "*180/pi")); -> "[e]" (math sqrt (- 1 (+ (* [sa] [sa]) (* [ca] [ca])))); -> "[sig]" (math calc (concat (math asin [e]) "*(1+0.1547*" [e] "^3)*180/pi")); #4. Resultant wind speed; -> "[vr]" (math hypot (math avg (json doc "rx_p4f_u")) (math avg (json doc "rx_p4f_w"))); -> "[naive]" (math avg (json doc "rx_p4f_dir")); json unload "rx_p4f_dir"; json unload "rx_p4f_v"; json unload "rx_p4f_s"; json unload "rx_p4f_c"; json unload "rx_p4f_u"; json unload "rx_p4f_w"; concat "vector mean direction " (math round (if (< [dir] 0) {+ [dir] 360} {[dir]}) 1) " deg (simple average " [naive] " deg) | Yamartino sigma " (math round [sig] 2) " deg | resultant speed " (math round [vr] 2) " m/s"
vector mean direction 0.8 deg (simple average 144.8 deg) | Yamartino sigma 8.74 deg | resultant speed 6.21 m/s
- Fit a Weibull distribution of the speeds and get the mean with
math gammaof 1 + 1/k - Flag a stuck vane when sigma stays at 0 for an hour
