ssa
1package tg2
3import (4"fmt"5"math"6"sort"7
8"github.com/RB-PRO/ssa/pkg/movmean"9"github.com/RB-PRO/ssa/pkg/oss"10)
11
12func CalcPW(RGBs []RGB_float64, Path string) (pw []float64, Err error) {13pw = make([]float64, len(RGBs))14pw2 := make([]float64, len(RGBs))15
16// CR17for i := range RGBs {18pw[i] = (RGBs[i].R*112.0 -19RGBs[i].G*93.8 -20RGBs[i].B*18.2) / 255.021
22// if _, err := filePW.WriteString(fmt.Sprintf("%.8f\n", pw[i])); err != nil {23// log.Println(err)24// }25}26
27// Вычитаем тренд28pw_smooth, _ := movmean.Movmean(pw, 32)29for i := range pw {30pw[i] -= pw_smooth[i]31}32
33// Квадрат34for i := range pw {35pw2[i] = math.Pow(pw[i], 2)36}37
38// fMi := 40.0 / 60.039// cad := 3040// SMO_med := cad / int(fMi)41
42DEV_med := medianFilter(pw2, 30*60/40)43createLineChart([]float64{}, DEV_med, Path+"DEV_med.png")44// filePWmedianFilter, _ := os.Create(Path + "medianFilter.txt")45// for i := range DEV_med {46// if _, err := filePWmedianFilter.WriteString(fmt.Sprintf("%.8f\n", DEV_med[i])); err != nil {47// log.Println(err)48// }49// }50
51for i := range pw {52pw[i] /= math.Sqrt(DEV_med[i])53}54prcMi := oss.Prctile(pw, 0.1)55prcMa := oss.Prctile(pw, 99.9)56prcMi = -0.019757prcMa = 0.020758fmt.Println("prcMi", prcMi, "prcMa", prcMa)59createLineChart([]float64{}, pw, Path+"pw2.png")60for i := range pw {61if pw[i] < prcMi {62pw[i] = prcMi63}64if pw[i] > prcMa {65pw[i] = prcMa66}67}68createLineChart([]float64{}, pw, Path+"pw.png")69
70return pw, nil71}
72
73func medianFilter(x []float64, n int) []float64 {74// Проверка на нечетное значение n75if n%2 == 0 {76n++77}78
79// Длина входного массива80length := len(x)81
82// Результат фильтрации83y := make([]float64, length)84
85for i := 0; i < length; i++ {86// Индексы для сбора значений для медианного фильтра87start := i - n/288end := i + n/289
90// Гарантия, что индексы не выходят за пределы массива91if start < 0 {92start = 093}94if end >= length {95end = length - 196}97
98// Извлечение значений для медианы99window := x[start : end+1]100
101// Сортировка окна значений и выбор медианы102sortedWindow := make([]float64, len(window))103copy(sortedWindow, window)104sort.Float64s(sortedWindow)105// Хитрый мув. При делении int(5) на int(2), получается int(2),106// т.е. округление в нижнюю сторону, хотя нам нужно в старшую степень.107// Поэтому из нечётного делаем чётное, а в случае получения нечётного не имеет разницы108medianIndex := (len(sortedWindow) + 1) / 2109// fmt.Println("medianIndex", medianIndex, "-", sortedWindow)110y[i] = sortedWindow[medianIndex]111}112
113return y114}
115