ssa

1
package gopw
2

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

8
	gomathtests "github.com/RB-PRO/ssa/pkg/go-MathTests"
9
	"github.com/RB-PRO/ssa/pkg/movmean"
10
)
11

12
// Метод cr извлечения сигнала фотоплетизмографии
13
func cr(R, G, B []float64) ([]float64, error) {
14
	if len(R) != len(G) || len(G) != len(B) {
15
		return nil, fmt.Errorf("the signal lengths R,G,B are not equal")
16
	}
17
	pw := make([]float64, len(R))
18
	pw2 := make([]float64, len(R))
19

20
	for i := range R {
21
		pw[i] = (R[i]*112.0 - G[i]*93.8 - B[i]*18.2) / 255.0
22
	}
23
	// fmt.Println("pw", pw, len(pw))
24

25
	// Вычитаем тренд
26
	pw_smooth, _ := movmean.Movmean(pw, 32)
27
	for i := range pw {
28
		pw[i] -= pw_smooth[i]
29
	}
30

31
	gomathtests.Plot("WorkPath/pw-smoov.png", pw)
32

33
	// Квадрат
34
	for i := range pw {
35
		pw2[i] = math.Pow(pw[i], 2)
36
	}
37
	gomathtests.Plot("WorkPath/pw2.png", pw2)
38

39
	// fMi := 40.0 / 60.0
40
	// cad := 30
41
	// SMO_med := cad / int(fMi)
42

43
	// DEV_med, ErrmedianFilter := medianFilter(pw2, 30*60/40)
44
	// if ErrmedianFilter != nil {
45
	// 	return nil, fmt.Errorf("medianFilter: %v", ErrmedianFilter)
46
	// }
47
	// DEV_med, _ := movmean.Movmean(pw2, 30*60/40)
48
	DEV_med := medianFilter1(pw2, 30*60/40)
49
	gomathtests.Plot("WorkPath/DEV_med.png", DEV_med)
50

51
	for i := range pw {
52
		pw[i] /= math.Sqrt(DEV_med[i])
53
	}
54
	gomathtests.Plot("WorkPath/pwdivdev.png", pw)
55
	// var cppw []float64
56
	// cppw = append(cppw, pw...)
57
	// prcMi := prctile(cppw, 0.1)
58
	// prcMa := prctile(cppw, 99.9)
59
	prcMi := -4.7962
60
	prcMa := 5.6327
61
	// fmt.Println("prcMi", prcMi, "prcMa", prcMa)
62

63
	for i := range pw {
64
		if pw[i] < prcMi {
65
			pw[i] = prcMi
66
		}
67
		if pw[i] > prcMa {
68
			pw[i] = prcMa
69
		}
70
	}
71

72
	STD := std(pw)
73
	gomathtests.Plot("WorkPath/pwstd.png", pw)
74
	for i := range pw {
75
		pw[i] /= STD
76
	}
77

78
	pw, _ = movmean.Movmean(pw, 5)
79

80
	return pw, nil
81
}
82

83
// Стандартное отклонение
84
func std(data []float64) float64 {
85

86
	mean := 0.0
87
	for i := range data {
88
		mean += data[i]
89
	}
90
	mean /= float64(len(data))
91

92
	N := len(data)
93

94
	var sum float64
95
	for i := 0; i < N; i++ {
96
		sum += math.Pow(math.Abs(data[i]-mean), 2)
97
	}
98

99
	return math.Sqrt(sum * (1 / (float64(N) - 1)))
100
}
101

102
func medianFilter1(input []float64, N int) []float64 {
103
	output := make([]float64, len(input))
104
	// copy(output, input) // Make a copy of the input to avoid modifying the original slice
105

106
	for i := 0; i < len(input); i++ {
107
		// Determine the median of the values within the window around index i
108

109
		var start, end int
110
		if N%2 == 0 { // Если чётное
111
			// %   For N even, Y(k) is the median of X( k-N/2 : k+N/2-1 ).
112
			start = i - N/2
113
			end = i + N/2
114
		} else { // Если НЕчётное
115
			// %   For N odd, Y(k) is the median of X( k-(N-1)/2 : k+(N-1)/2 ).
116
			start = i - (N-1)/2
117
			end = i + (N-1)/2 + 1
118
		}
119

120
		if start < 0 {
121
			start = 0
122
		}
123
		if end >= len(input) {
124
			end = len(input)
125
		}
126

127
		// copy subslice
128
		var subarray []float64
129
		subarray = append(subarray, input[start:end]...)
130

131
		// fmt.Println(i, input, ">", start, end, "<>", subarray, " ")
132
		output[i] = meanWindow(subarray, N)
133

134
		// // Sort the subarray to find the median
135
		// subarray := input[start : end+1]
136
		// fmt.Println(i, subarray)
137
		// sort.Float64s(subarray)
138

139
		// // Set the output to the median value of the subarray
140
		// output[i] = subarray[len(subarray)/2]
141
	}
142

143
	return output
144
}
145

146
func meanWindow(data []float64, N int) float64 {
147

148
	// Если это краевой случай и к-во элементов не равно длине окна,
149
	// то нам нужно расширить исследуемый слайс
150
	if len(data) != N {
151
		app := make([]float64, N-len(data))
152
		data = append(data, app...)
153
	}
154
	sort.Float64s(data)
155

156
	return data[len(data)/2]
157
}
158