common.R (4999B)
1 library(DBI) 2 3 # https://stackoverflow.com/questions/13023274/how-to-do-printf-in-r 4 printf <- function(...) cat(sprintf(...)) 5 6 reference_value <- 2.045 # measured with commercial anemometer 7 8 calc_origin <- function (filename="anemometer-no-wind") { 9 data <- read.table(filename, sep="", header=FALSE) 10 names(data) <- c("t", "x", "y", "z", "sx", "sy", "sz") 11 origin = list(x=mean(data$x), y=mean(data$y), z=mean(data$z)) 12 origin 13 } 14 15 abs_sqrt <- function (x) { sign(x) * sqrt(abs(x)) } 16 17 calc_coefficient <- function (filename, axis) { 18 origin <- calc_origin(sprintf("%s-origin",filename)) 19 data <- read.table(filename, sep="", header=FALSE) 20 names(data) <- c("t", "x", "y", "z", "sx", "sy", "sz") 21 sign <- data[1,sprintf("s%s",axis)] 22 sign*abs_sqrt(mean(data[[axis]])-origin[[axis]])/reference_value 23 } 24 25 calibration_coefficients <- function(debug=FALSE) { 26 # *1 --- negative, *2 --- positive, *reverse --- positive 27 cx1 <- calc_coefficient("build/x", "x") 28 cx2 <- calc_coefficient("build/x-reverse", "x") 29 cy1 <- calc_coefficient("build/y", "y") 30 cy2 <- calc_coefficient("build/y-reverse", "y") 31 cz1 <- calc_coefficient("build/z", "z") 32 cz2 <- calc_coefficient("build/z-reverse", "z") 33 if (debug) { 34 printf("-1 m/s is x %f y %f z %f\n", cx1, cy1, cz1) 35 printf("1 m/s is x %f y %f z %f\n", cx2, cy2, cz2) 36 } 37 list(x1=cx1, x2=cx2, y1=cy1, y2=cy2, z1=cz1, z2=cz2) 38 #weight_0 = (calc_origin("build/z-fan")$z - calc_origin("build/z-fan-origin")$z) 39 #weight_1 = (calc_origin("build/z-fan-and-wind")$z - calc_origin("build/z-fan-and-wind-origin")$z) 40 #print(weight_0) 41 #print(weight_1) 42 #print(abs_sqrt(weight_1-weight_0)/reference_value) 43 } 44 45 # convert load cell anemometer sensor value to wind speed in m/s 46 sample_to_speed <- function(x, c1, c2) { 47 # remove linear trend 48 t <- c(1:length(x)) 49 reg <- lm(x~t) 50 x <- x - reg$fitted.values 51 #x <- x - mean(x) # convert from unsigned to signed sensor values 52 x <- sign(x)*sqrt(abs(x)) # convert from force to velocity 53 # scale sensor values to wind speed using calibration coefficients 54 x[x<0] = x[x<0] / -c1 55 x[x>0] = x[x>0] / c2 56 x 57 } 58 59 speed <- function(x,y,z=0) { sqrt(x**2 + y**2 + z**2) } 60 direction <- function(x,y) { atan2(y,x) } 61 62 select_samples <- function (timestamp, time_delta, url="samples/load-cell.sqlite3", 63 table="samples") { 64 db <- dbConnect(RSQLite::SQLite(), url) 65 velocity <- NULL 66 if (is.null(timestamp)) { 67 velocity = dbGetQuery(db, 68 sprintf("SELECT x,y,z FROM %s ORDER BY timestamp, t", table)) 69 velocity = dbGetQuery(db, 70 sprintf("SELECT timestamp,x,y,z FROM %s 71 WHERE ABS(x-:xmean)<10000 AND ABS(y-:ymean)<10000 AND ABS(z-:zmean)<10000 72 ORDER BY timestamp, t", table), 73 params = list(xmean=mean(velocity$x), ymean=mean(velocity$y), zmean=mean(velocity$z))) 74 } else { 75 velocity = dbGetQuery(db, 76 sprintf("SELECT x,y,z FROM %s WHERE timestamp BETWEEN :t0 AND (:t0 + :dt) ORDER BY timestamp, t", table), 77 params = list(t0 = timestamp, dt=time_delta)) 78 velocity = dbGetQuery(db, 79 sprintf("SELECT timestamp,x,y,z FROM %s 80 WHERE timestamp BETWEEN :t0 AND (:t0 + :dt) 81 AND ABS(x-:xmean)<10000 AND ABS(y-:ymean)<10000 AND ABS(z-:zmean)<10000 82 ORDER BY timestamp, t", table), 83 params = list(t0 = timestamp, dt=time_delta, 84 xmean=mean(velocity$x), ymean=mean(velocity$y), zmean=mean(velocity$z))) 85 } 86 if (nrow(velocity) == 0) { return(data.frame()) } 87 # convert sensor values to wind speed in m/s 88 coefficients <- calibration_coefficients() 89 velocity$x = sample_to_speed(velocity$x, coefficients$x1, coefficients$x2) 90 velocity$y = sample_to_speed(velocity$y, coefficients$y1, coefficients$y2) 91 velocity$z = sample_to_speed(velocity$z, coefficients$z1, coefficients$z2) 92 # compute wind speed and direction 93 velocity$speed = speed(velocity$x,velocity$y,velocity$z) 94 velocity$direction = direction(velocity$x,velocity$y) 95 dbDisconnect(db) 96 velocity 97 } 98 99 datetime_format <- "%FT%H:%M" 100 101 time_string_to_timestamp <- function(str) { 102 as.integer(as.POSIXct(strptime(str, datetime_format))) 103 } 104 105 make_directory <- function (path) { 106 dir.create(path, recursive=TRUE, showWarnings=FALSE) 107 } 108 109 # https://stackoverflow.com/questions/32370485/convert-radians-to-degree-degree-to-radians 110 rad2deg <- function(rad) {(rad * 180) / (pi)} 111 deg2rad <- function(deg) {(deg * pi) / (180)} 112 113 # https://stackoverflow.com/questions/27455948/in-r-whats-the-simplest-way-to-scale-a-vector-to-a-unit-vector 114 unit <- function(x) { 115 length = sqrt(sum(x**2)) 116 ifelse(length == 0, x, x/length) 117 } 118 119 # normalized rmse 120 rmse <- function(estimated, actual) { 121 s_max <- max(actual) 122 s_min <- min(actual) 123 sqrt(mean((estimated - actual)**2)) / (s_max - s_min) 124 }