ssa
1package gopw2
3import (4"fmt"5"math"6"sort"7
8gomathtests "github.com/RB-PRO/ssa/pkg/go-MathTests"9"github.com/RB-PRO/ssa/pkg/movmean"10)
11
12// Метод cr извлечения сигнала фотоплетизмографии
13func cr(R, G, B []float64) ([]float64, error) {14if len(R) != len(G) || len(G) != len(B) {15return nil, fmt.Errorf("the signal lengths R,G,B are not equal")16}17pw := make([]float64, len(R))18pw2 := make([]float64, len(R))19
20for i := range R {21pw[i] = (R[i]*112.0 - G[i]*93.8 - B[i]*18.2) / 255.022}23// fmt.Println("pw", pw, len(pw))24
25// Вычитаем тренд26pw_smooth, _ := movmean.Movmean(pw, 32)27for i := range pw {28pw[i] -= pw_smooth[i]29}30
31gomathtests.Plot("WorkPath/pw-smoov.png", pw)32
33// Квадрат34for i := range pw {35pw2[i] = math.Pow(pw[i], 2)36}37gomathtests.Plot("WorkPath/pw2.png", pw2)38
39// fMi := 40.0 / 60.040// cad := 3041// 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)48DEV_med := medianFilter1(pw2, 30*60/40)49gomathtests.Plot("WorkPath/DEV_med.png", DEV_med)50
51for i := range pw {52pw[i] /= math.Sqrt(DEV_med[i])53}54gomathtests.Plot("WorkPath/pwdivdev.png", pw)55// var cppw []float6456// cppw = append(cppw, pw...)57// prcMi := prctile(cppw, 0.1)58// prcMa := prctile(cppw, 99.9)59prcMi := -4.796260prcMa := 5.632761// fmt.Println("prcMi", prcMi, "prcMa", prcMa)62
63for i := range pw {64if pw[i] < prcMi {65pw[i] = prcMi66}67if pw[i] > prcMa {68pw[i] = prcMa69}70}71
72STD := std(pw)73gomathtests.Plot("WorkPath/pwstd.png", pw)74for i := range pw {75pw[i] /= STD76}77
78pw, _ = movmean.Movmean(pw, 5)79
80return pw, nil81}
82
83// Стандартное отклонение
84func std(data []float64) float64 {85
86mean := 0.087for i := range data {88mean += data[i]89}90mean /= float64(len(data))91
92N := len(data)93
94var sum float6495for i := 0; i < N; i++ {96sum += math.Pow(math.Abs(data[i]-mean), 2)97}98
99return math.Sqrt(sum * (1 / (float64(N) - 1)))100}
101
102func medianFilter1(input []float64, N int) []float64 {103output := make([]float64, len(input))104// copy(output, input) // Make a copy of the input to avoid modifying the original slice105
106for i := 0; i < len(input); i++ {107// Determine the median of the values within the window around index i108
109var start, end int110if N%2 == 0 { // Если чётное111// % For N even, Y(k) is the median of X( k-N/2 : k+N/2-1 ).112start = i - N/2113end = i + N/2114} else { // Если НЕчётное115// % For N odd, Y(k) is the median of X( k-(N-1)/2 : k+(N-1)/2 ).116start = i - (N-1)/2117end = i + (N-1)/2 + 1118}119
120if start < 0 {121start = 0122}123if end >= len(input) {124end = len(input)125}126
127// copy subslice128var subarray []float64129subarray = append(subarray, input[start:end]...)130
131// fmt.Println(i, input, ">", start, end, "<>", subarray, " ")132output[i] = meanWindow(subarray, N)133
134// // Sort the subarray to find the median135// 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 subarray140// output[i] = subarray[len(subarray)/2]141}142
143return output144}
145
146func meanWindow(data []float64, N int) float64 {147
148// Если это краевой случай и к-во элементов не равно длине окна,149// то нам нужно расширить исследуемый слайс150if len(data) != N {151app := make([]float64, N-len(data))152data = append(data, app...)153}154sort.Float64s(data)155
156return data[len(data)/2]157}
158