-
Notifications
You must be signed in to change notification settings - Fork 25
/
rose.go
99 lines (80 loc) · 1.87 KB
/
rose.go
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
package main
import (
"encoding/csv"
"fmt"
"os"
"strconv"
"github.com/konimarti/kalman"
"github.com/konimarti/lti"
"gonum.org/v1/gonum/mat"
)
func main() {
// load data
y := load("p133.csv")[0]
// prepare output file
file, err := os.Create("p133_out.csv")
if err != nil {
panic(err)
}
defer file.Close()
fmt.Fprintln(file, "m,f")
// define LTI system
lti := lti.Discrete{
Ad: mat.NewDense(1, 1, []float64{1}),
Bd: mat.NewDense(1, 1, nil),
C: mat.NewDense(1, 1, []float64{1}),
D: mat.NewDense(1, 1, nil),
}
// system noise / process model covariance matrix ("Systemrauschen")
Gd := mat.NewDense(1, 1, []float64{1})
ctx := kalman.Context{
// initial state
X: mat.NewVecDense(1, []float64{y[0]}),
// initial covariance matrix
P: mat.NewDense(1, 1, []float64{0}),
}
// create ROSE filter
gammaR := 9.0
alphaR := 0.5
alphaM := 0.3
filter := kalman.NewRoseFilter(lti, Gd, gammaR, alphaR, alphaM)
// no control
u := mat.NewVecDense(1, nil)
for _, row := range y {
// new measurement
y := mat.NewVecDense(1, []float64{row})
// apply filter
filter.Apply(&ctx, y, u)
// get corrected state vector
state := filter.State()
// print out input and output signals
fmt.Fprintf(file, "%3.8f,%3.8f\n", y.AtVec(0), state.AtVec(0))
}
}
func load(file string) [][]float64 {
f, err := os.Open(file)
if err != nil {
panic("could not open file")
}
defer f.Close()
r := csv.NewReader(f)
r.Comma = ','
lines, err := r.ReadAll()
if err != nil {
panic("could read data from file")
}
y := make([][]float64, len(lines))
for row, line := range lines {
tmp := make([]float64, len(line))
for i, entry := range line {
var value float64
if value, err = strconv.ParseFloat(entry, 64); err != nil {
fmt.Println("error parsing", line, i, entry)
panic("could not parse all values")
}
tmp[i] = value
}
y[row] = tmp
}
return y
}