ssa

1
package ssa2
2

3
import (
4
	"fmt"
5
	"log"
6
	"math"
7
	"os"
8
	"slices"
9
	"sort"
10

11
	"github.com/RB-PRO/ssa/pkg/ssa2/pmtm"
12
)
13

14
// Структура загрузчика данных
15
type Setup struct {
16
	Cad   int     // Количество кадров в секунду
17
	NPart int     // Количество долей res
18
	Win   int     // Ширина окна
19
	FMi   float64 // Частота среза для 40 уд/мин (0.6667 Гц)
20
	FMa   float64 // Частота среза для 240 уд/мин (4.0 Гц)
21
}
22

23
type SSA struct {
24
	set Setup // Структура настройки Расчётов SSA
25

26
	len int // Количество отсчётов pw
27

28
	// Временные интервалы
29
	tim []float64 // Время в секундах
30
	dt  float64   // Интервал дискретизации времени, сек
31

32
	col int // Всего колонок
33

34
	// Сигнал фотоплетизмографии
35
	pw       []float64 // Сигнал фотоплетизмографии
36
	freq     []float64
37
	Pto_fMAX []float64
38

39
	// Оценки СПМ перекрывающихся сегменов pw
40
	iGmin int
41
	iGmax int
42

43
	res int
44
}
45

46
// Создать объект SSA для работы
47
func NewSSA(pw []float64, Setting Setup) (*SSA, error) {
48
	var ssa SSA
49
	ssa.pw = pw
50
	ssa.set = Setting
51

52
	// к-во отсчётов пульсовой волны
53
	ssa.len = len(ssa.pw)
54

55
	// интервал дискретизации времени, сек
56
	ssa.dt = 1.0 / float64(ssa.set.Cad)
57

58
	// время в секундах
59
	tim := make([]float64, ssa.len)
60
	for i := 0; i < ssa.len; i++ {
61
		tim[i] = float64(i) * ssa.dt
62
	}
63
	ssa.tim = tim
64

65
	return &ssa, nil
66
}
67

68
// Расчёт параметра col
69
func (ssa *SSA) Col() (*SSA, error) {
70
	col := 1
71
	Imin := 1
72
	Imax := ssa.set.Win
73
	resFloat64 := float64(ssa.len) - float64(ssa.set.Win)*math.Floor(float64(ssa.len)/float64(ssa.set.Win))
74
	res := int(math.Floor(resFloat64 / float64(ssa.set.NPart)))
75

76
	// col - кол-во перекрывающихся сегментов в пределах len
77
	for Imax <= ssa.len {
78
		Imin = Imin + res
79
		Imax = Imax + res
80
		col = col + 1
81
	}
82
	col--
83

84
	ssa.col = col
85
	ssa.res = res
86
	return ssa, nil
87
}
88

89
// текущий сегмент pw длинною win
90
func (ssa *SSA) Spw(column int) ([]float64, error) {
91
	if ssa.col == 0 {
92
		return nil, fmt.Errorf("spw: ssa.col is 0")
93
	}
94
	if column >= ssa.col {
95
		return nil, fmt.Errorf("spw: Условие column(%d) < ssa.col(%d) не выполнено", column, ssa.col)
96
	}
97
	if len(ssa.pw) == 0 {
98
		return nil, fmt.Errorf("spw: len(pw)=0")
99
	}
100
	return ssa.pw[column*ssa.res : column*ssa.res+ssa.set.Win], nil
101
}
102

103
// Оценки СПМ перекрывающихся сегменов pw
104
func (ssa *SSA) SpwEstimation() (*SSA, error) {
105
	df := float64(ssa.set.Cad) / float64(ssa.set.Win-1)
106
	Fmin := ssa.set.FMi - 10*df
107
	Fmax := ssa.set.FMa + 10*df
108
	row := 1 + ssa.set.Win/2
109

110
	tecalF := 0.0
111
	ssa.freq = make([]float64, row)
112
	for i := 0; i < row; i++ {
113
		if math.Abs(tecalF-Fmin) <= df {
114
			ssa.iGmin = i
115
		}
116
		if math.Abs(tecalF-Fmax) <= df {
117
			ssa.iGmax = i
118
		}
119
		ssa.freq[i] = tecalF
120
		tecalF += df
121
	}
122
	SaveTXT("ssafreq.txt", ssa.freq)
123

124
	return ssa, nil
125
}
126

127
// Оценки средних частот основного тона для сегментов pw
128
func (ssa *SSA) PwEstimation() (*SSA, error) {
129
	// BlackManHar := blackmanharris.Blackmanharris(ssa.set.Win, blackmanharris.Koef4_74db)
130

131
	filePW, _ := os.Create("tests/pto.txt")
132
	defer filePW.Close()
133

134
	// Выделяем память в слайс
135
	ssa.Pto_fMAX = make([]float64, ssa.col)
136
	for i := 0; i < ssa.col; i++ {
137

138
		// Получение колонки SPW
139
		ColumnSPW, ErrSpw := ssa.Spw(i)
140
		if ErrSpw != nil {
141
			return ssa, fmt.Errorf("SSA: PwEstimation: %w", ErrSpw)
142
		}
143
		// ColumnSPW := make([]float64, len(ColumnSPW2))
144
		// copy(ColumnSPW, ColumnSPW2)
145

146
		// Получение периодограммы колонки SPW
147
		// pg_spw := periodogram.Periodogram(ColumnSPW, BlackManHar, ssa.set.Win)
148

149
		pto_spw := pmtm.Pmtm(ColumnSPW, 3, ssa.set.Win)
150

151
		for j := range pto_spw {
152
			filePW.WriteString(fmt.Sprintf("%.8f;", pto_spw[j]))
153
		}
154
		filePW.WriteString(fmt.Sprintf("\n"))
155
		pto_spw = medianFilter(pto_spw, 30)
156
		// pto_spw = sredFilter(pto_spw, 15)
157

158
		// Сотртируем, чтобы получить по возрастанию порядковвые номера в слайсе по убыванию
159
		// _, SorterIndexts_pg_spw := InsertionSort(pto_spw)
160
		SorterIndexts_pg_spw := InsertionSort2(pto_spw)
161
		// _, SorterIndexts_pg_spw := InsertionSortInt(pto_spw)
162

163
		// fmt.Println(SorterIndexts_pg_spw[0], len(pto_spw))
164
		ssa.Pto_fMAX[i] = ssa.freq[SorterIndexts_pg_spw] // float64(SorterIndexts_pg_spw)
165
		// ssa.Pto_fMAX[i] = float64(SorterIndexts_pg_spw)
166
	}
167

168
	return ssa, nil
169
}
170

171
func SaveTXT(FileName string, data []float64) {
172
	filePW, ErrOpenFile := os.Create(FileName)
173
	if ErrOpenFile != nil {
174
		panic(ErrOpenFile)
175
	}
176
	defer filePW.Close()
177
	for i := range data {
178
		if _, err := filePW.WriteString(fmt.Sprintf("%.8f\n", data[i])); err != nil {
179
			log.Println(err)
180
		}
181
	}
182
}
183

184
// Сортировка с возвратом номеров изначальных элементов
185
func InsertionSort(array []float64) ([]float64, []int) {
186
	// var indexArray int
187
	// sort.Slice(array, func(i, j int) bool {
188
	// 	if array[i] > array[j] {
189
	// 		indexArray = j
190
	// 	}
191
	// 	return array[i] > array[j]
192
	// })
193
	// fmt.Println(indexArray, array[0], array[1], array[2])
194

195
	indexArray := make([]int, len(array))
196
	for ind := range indexArray {
197
		indexArray[ind] = (ind)
198
	}
199
	for i := 1; i < len(array); i++ {
200
		j := i
201
		for j > 0 {
202
			if array[j-1] < array[j] {
203
				array[j-1], array[j] = array[j], array[j-1]
204
				indexArray[j-1], indexArray[j] = indexArray[j], indexArray[j-1]
205
			}
206
			j = j - 1
207
		}
208
	}
209
	return array, indexArray
210
}
211

212
// Сортировка с возвратом номеров изначальных элементов
213
func InsertionSortInt(array []float64) ([]float64, int) {
214
	indexArrayint := 0
215
	for i := 1; i < len(array); i++ {
216
		j := i
217
		for j > 0 {
218
			if array[j-1] < array[j] {
219
				array[j-1], array[j] = array[j], array[j-1]
220
				indexArrayint = i
221
			}
222
			j = j - 1
223
		}
224
	}
225
	return array, indexArrayint
226
}
227

228
// Сортировка с возвратом номеров изначальных элементов
229
func InsertionSort2(array []float64) int {
230
	return slices.Index(array, slices.Max(array))
231
}
232

233
func medianFilter(x []float64, n int) []float64 {
234
	// Проверка на нечетное значение n
235
	if n%2 == 0 {
236
		n++
237
	}
238

239
	// Длина входного массива
240
	length := len(x)
241

242
	// Результат фильтрации
243
	y := make([]float64, length)
244

245
	for i := 0; i < length; i++ {
246
		// Индексы для сбора значений для медианного фильтра
247
		start := i - n/2
248
		end := i + n/2
249

250
		// Гарантия, что индексы не выходят за пределы массива
251
		if start < 0 {
252
			start = 0
253
		}
254
		if end >= length {
255
			end = length - 1
256
		}
257

258
		// Извлечение значений для медианы
259
		window := x[start : end+1]
260

261
		// Сортировка окна значений и выбор медианы
262
		sortedWindow := make([]float64, len(window))
263
		copy(sortedWindow, window)
264
		sort.Float64s(sortedWindow)
265
		// Хитрый мув. При делении int(5) на int(2), получается int(2),
266
		// т.е. округление в нижнюю сторону, хотя нам нужно в старшую степень.
267
		// Поэтому из нечётного делаем чётное, а в случае получения нечётного не имеет разницы
268
		medianIndex := (len(sortedWindow) + 1) / 2
269
		// fmt.Println("medianIndex", medianIndex, "-", sortedWindow)
270
		y[i] = sortedWindow[medianIndex]
271
	}
272

273
	return y
274
}
275

276
func sredFilter(x []float64, n int) []float64 {
277
	// Проверка на нечетное значение n
278
	if n%2 == 0 {
279
		n++
280
	}
281

282
	// Длина входного массива
283
	length := len(x)
284

285
	// Результат фильтрации
286
	y := make([]float64, length)
287

288
	for i := 0; i < length; i++ {
289
		// Индексы для сбора значений для медианного фильтра
290
		start := i - n/2
291
		end := i + n/2
292

293
		// Гарантия, что индексы не выходят за пределы массива
294
		if start < 0 {
295
			start = 0
296
		}
297
		if end >= length {
298
			end = length - 1
299
		}
300

301
		// Извлечение значений для медианы
302
		window := x[start : end+1]
303

304
		// Сортировка окна значений и выбор медианы
305
		// slises
306
		sum := 0.0
307
		for _, v := range window {
308
			sum += v
309
		}
310
		y[i] = sum / float64(n)
311

312
	}
313

314
	return y
315
}
316