Create a gist now

Instantly share code, notes, and snippets.

What would you like to do?
# An R script to estimate MQ gas sensors correlation curve and compute Ro, min and max Rs/Ro
#
# Copyright (c) Davide Gironi, 2016
#
# Released under GPLv3.
# Please refer to LICENSE file for licensing information.
# How to use this script:
# 1) set limits as datasheet curve ("xlim" and "ylim")
# ex.
# xlim = c(10, 1000)
# ylim = c(0.1, 10)
# 2) find out datasheet curve points, and write it out (to "pointsdata")
# each line it's a point on cartesian coordinate system
# the useful WebPlotDigitizer app can help you extract points from the graph
# ex.
# pointsdata = "
# 10.052112405371744, 2.283698378106183
# 20.171602728600178, 1.8052797165878915
# 30.099224396434586, 1.5715748803154423
# 50.09267987761949, 1.3195287228519417
# 80.38812026903305, 1.1281218760133969
# 90.12665922665023, 1.0815121769656304
# 100.52112405371739, 1.0430967861855598
# 199.62996638292853, 0.8000946404902397
# "
# 3) optional for Ro estimation: measure the sensor resistance (set it to "mres" ohm value) at a know amount of gas
# set it to 0 if you do not need the Ro estimation
# ex.
# mres = 26954
# 4) optional for Ro estimation: set the know amount of gas for the resistance measure of the previous step (to "mppm")
# set it to 0 if you do not need the Ro estimation
# ex.
# mppm = 392
# 5) optional for min-max Rs/Ro estimation: set the minand max amount of gas the sensor will react to (as "minppm" and "maxppm")
# set it to 0 if you do not need the min-max Rs/Ro estimation
# ex.
# minppm = 10
# maxppm = 200
library(data.table)
#remove old variables
rm(list=ls())
#set input values
xlim = c(0.1, 10)
ylim = c(0.1, 10)
minppm = 0
maxppm = 0
mres = 0
mppm = 0
pointsdata = "
YOUR_POINTS_HERE
"
#load points using fread
setnames(points <- fread(pointsdata, sep=",", sep2="\n"), c("x","y"))
#set named list of points, and swapped list of points
#points will be used to plot and compute values as datasheet figure
#pointsrev will be used to plot and compute values for the correlation function, it's the datasheet figure with swapped axis
x <- as.vector(points[,x])
y <- as.vector(points[,y])
points = list(x=x, y=y)
pointsrev = list(x=y, y=x)
#the nls (Nonlinear Least Squares) it's used to perform the power regression on points
#in order to work, nls needs an estimation of staring values
#we use log-log slope estimation to find intitial values
#estimate fit curve initial values
xfirst = head(points$x, n=1)
xlast = tail(points$x, n=1)
yfirst = head(points$y, n=1)
ylast = tail(points$y, n=1)
bstart= log(ylast/yfirst)/log(xlast/xfirst)
astart = yfirst/(xfirst^bstart)
#perform the fit
fit <- nls("y~a*x^b", start=list(a=astart,b=bstart), data=points)
#estimate fitref curve initial values
xfirstrev = head(pointsrev$x, n=1)
xlastrev = tail(pointsrev$x, n=1)
yfirstrev = head(pointsrev$y, n=1)
ylastrev = tail(pointsrev$y, n=1)
bstartrev = log(ylastrev/yfirstrev)/log(xlastrev/xfirstrev)
astartrev = yfirstrev/(xfirstrev^bstartrev)
fitrev <- nls("y~a*x^b", start=list(a=astartrev,b=bstartrev), data=pointsrev)
#plot fit curve (log-log scale)
fiteq = function(x){coef(fit)["a"]*x^(coef(fit)["b"])}
plot(points, log="xy", col="blue", xlab="ppm", ylab="Rs/Ro", xlim=xlim, ylim=ylim, panel.first=grid(equilogs=FALSE))
curve(fiteq, col="red", add=TRUE)
#plot fitrev curve (log-log scale)
fiteqrev = function(x){coef(fitrev)["a"]*x^(coef(fitrev)["b"])}
plot(pointsrev, log="xy", col="blue", xlab="Rs/Ro", ylab="ppm", xlim=ylim, ylim=xlim, panel.first=grid(equilogs=FALSE))
curve(fiteqrev, col="red", add=TRUE)
#plot fit curve (linear scale)
fiteq = function(x){coef(fit)["a"]*x^(coef(fit)["b"])}
plot(points, col="blue", xlab="ppm", ylab="Rs/Ro", panel.first=grid(equilogs=FALSE))
curve(fiteq, col="red", add=TRUE)
#plot fitrev curve (linear scale)
fiteqrev = function(x){coef(fitrev)["a"]*x^(coef(fitrev)["b"])}
plot(pointsrev, col="blue", xlab="Rs/Ro", ylab="ppm", panel.first=grid(equilogs=FALSE))
curve(fiteqrev, col="red", add=TRUE)
#estimate min Rs/Ro
cat("\nCorrelation function coefficients")
cat("\nEstimated a\n")
cat(" ")
cat(coef(fitrev)["a"])
cat("\nEstimated b\n")
cat(" ")
cat(coef(fitrev)["b"])
cat("\n")
#estimate min Rs/Ro
if (minppm != 0) {
minRsRo = (maxppm/coef(fitrev)["a"])^(1/coef(fitrev)["b"])
cat("\nEstimated min Rs/Ro\n")
cat(" ")
cat(minRsRo)
cat("\n")
}
#estimate max Rs/Ro
if (maxppm != 0) {
maxRsRo = (minppm/coef(fitrev)["a"])^(1/coef(fitrev)["b"])
cat("\nEstimated max Rs/Ro\n")
cat(" ")
cat(maxRsRo)
cat("\n")
}
#estimate Ro
if (mppm != 0 && mres != 0) {
Ro = mres*(coef(fitrev)["a"]/mppm)^(1/coef(fitrev)["b"])
cat("\nEstimated Ro\n")
cat(" ")
cat(Ro)
cat("\n")
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment