ssa

1
package pchip
2

3
import (
4
	"fmt"
5
	"math"
6
)
7

8
func exteriorSlope(d1, d2, h1, h2 float64) float64 {
9
	var s float64
10
	var signd1 float64
11
	var signs float64
12
	s = ((2.0*h1+h2)*d1 - h1*d2) / (h1 + h2)
13
	signd1 = d1
14
	if d1 < 0.0 {
15
		signd1 = -1.0
16
	} else if d1 > 0.0 {
17
		signd1 = 1.0
18
	} else {
19
		if d1 == 0.0 {
20
			signd1 = 0.0
21
		}
22
	}
23
	signs = s
24
	if s < 0.0 {
25
		signs = -1.0
26
	} else if s > 0.0 {
27
		signs = 1.0
28
	} else {
29
		if s == 0.0 {
30
			signs = 0.0
31
		}
32
	}
33
	if signs != signd1 {
34
		s = 0.0
35
	} else {
36
		signs = d2
37
		if d2 < 0.0 {
38
			signs = -1.0
39
		} else if d2 > 0.0 {
40
			signs = 1.0
41
		} else {
42
			if d2 == 0.0 {
43
				signs = 0.0
44
			}
45
		}
46
		if (signd1 != signs) && (math.Abs(s) > math.Abs(3.0*d1)) {
47
			s = 3.0 * d1
48
		}
49
	}
50
	return s
51
}
52

53
type PchipCoefs struct {
54
	A, B, C, D []float64
55
}
56

57
func Pchip(x, y, new_x []float64, x_len, new_x_len int) ([]float64, []float64, PchipCoefs) {
58
	// x_len = 12 // кал
59
	new_y := make([]float64, new_x_len)
60
	var low_ip1 int
61
	var hs float64
62
	del := make([]float64, x_len-1)
63
	slopes := make([]float64, x_len)
64
	h := make([]float64, x_len-1)
65

66
	var hs3 float64
67
	var w1 float64
68
	//var ix int
69
	pp_coefs := make([]float64, (x_len-1)+(3*(x_len-1)))
70

71
	var low_i int
72
	var high_i int
73
	var mid_i int
74
	for low_ip1 := 0; low_ip1 < x_len-1; low_ip1++ {
75
		hs = x[low_ip1+1] - x[low_ip1]
76
		del[low_ip1] = (y[low_ip1+1] - y[low_ip1]) / hs
77
		h[low_ip1] = hs
78
	}
79

80
	for low_ip1 := 0; low_ip1 < x_len-2; low_ip1++ {
81
		hs = h[low_ip1] + h[low_ip1+1]
82
		hs3 = 3.0 * hs
83
		w1 = (h[low_ip1] + hs) / hs3
84
		hs = (h[low_ip1+1] + hs) / hs3
85
		hs3 = 0.0
86
		if del[low_ip1] < 0.0 {
87
			if del[low_ip1+1] <= del[low_ip1] {
88
				hs3 = del[low_ip1] / (w1*(del[low_ip1]/del[low_ip1+1]) + hs)
89
			} else {
90
				if del[low_ip1+1] < 0.0 {
91
					hs3 = del[low_ip1+1] / (w1 + hs*(del[low_ip1+1]/del[low_ip1]))
92
				}
93
			}
94
		} else {
95
			if del[low_ip1] > 0.0 {
96
				if del[low_ip1+1] >= del[low_ip1] {
97
					hs3 = del[low_ip1] / (w1*(del[low_ip1]/del[low_ip1+1]) + hs)
98
				} else {
99
					if del[low_ip1+1] > 0.0 {
100
						hs3 = del[low_ip1+1] / (w1 + hs*(del[low_ip1+1]/del[low_ip1]))
101
					}
102
				}
103
			}
104
		}
105

106
		slopes[low_ip1+1] = hs3
107
	}
108
	if x_len <= 3 {
109
		return new_y, pp_coefs, PchipCoefs{}
110
	}
111
	fmt.Println("x_len", x_len)
112

113
	slopes[0] = exteriorSlope(del[0], del[1], h[0], h[1])
114
	slopes[x_len-1] = exteriorSlope(del[x_len-2], del[x_len-3], h[x_len-2], h[x_len-3])
115
	for low_ip1 := 0; low_ip1 < x_len-1; low_ip1++ {
116
		hs = (del[low_ip1] - slopes[low_ip1]) / h[low_ip1]
117
		hs3 = (slopes[low_ip1+1] - del[low_ip1]) / h[low_ip1]
118
		pp_coefs[low_ip1] = (hs3 - hs) / h[low_ip1]
119
		pp_coefs[low_ip1+x_len-1] = 2.0*hs - hs3
120
		pp_coefs[low_ip1+(2*(x_len-1))] = slopes[low_ip1]
121
		pp_coefs[low_ip1+(3*(x_len-1))] = y[low_ip1]
122
	}
123

124
	for ix := 0; ix < new_x_len; ix++ {
125
		low_i = 0
126
		low_ip1 = 2
127
		high_i = x_len
128
		for high_i > low_ip1 {
129
			mid_i = ((low_i + high_i) + 1) >> 1
130
			if new_x[ix] >= x[mid_i-1] {
131
				low_i = mid_i - 1
132
				low_ip1 = mid_i + 1
133
			} else {
134
				high_i = mid_i
135
			}
136
		}
137

138
		hs = new_x[ix] - x[low_i]
139
		hs3 = pp_coefs[low_i]
140
		for low_ip1 := 0; low_ip1 < 3; low_ip1++ {
141
			hs3 = hs*hs3 + pp_coefs[low_i+(low_ip1+1)*(x_len-1)]
142
		}
143

144
		new_y[ix] = hs3
145
		/*
146
			// my
147
			if new_y[ix] < 0.15 {
148
				new_y[ix] = new_y[ix-1]
149
				fmt.Println("YESS")
150
			}
151
		*/
152
	}
153

154
	var coefs PchipCoefs
155
	lenCoefs := len(pp_coefs) / 4
156
	coefs.A = make([]float64, lenCoefs)
157
	coefs.B = make([]float64, lenCoefs)
158
	coefs.C = make([]float64, lenCoefs)
159
	coefs.D = make([]float64, lenCoefs)
160
	for i := 0; i < lenCoefs; i++ {
161
		coefs.A[i] = pp_coefs[i]
162
		coefs.B[i] = pp_coefs[1*(len(pp_coefs)/4)+i]
163
		coefs.C[i] = pp_coefs[2*(len(pp_coefs)/4)+i]
164
		coefs.D[i] = pp_coefs[3*(len(pp_coefs)/4)+i]
165
	}
166

167
	return new_y, pp_coefs, coefs
168
}
169

170
// Функция, реализующая сглаживание с помощью фильтра Савицкого-Голея
171