ssa

1
package ssaApp
2

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

8
	"github.com/RB-PRO/ssa/pkg/blackmanharris"
9
	"github.com/RB-PRO/ssa/pkg/oss"
10
	"github.com/RB-PRO/ssa/pkg/periodogram"
11
	"github.com/RB-PRO/ssa/pkg/ssa"
12
	"gonum.org/v1/gonum/mat"
13
)
14

15
func SSS_spw2(pw, fmp []float64) {
16
	s := ssa.New("File")
17
	s.Graph = true // Создавать графики
18
	s.Xlsx = true  // Сохранять в Xlsx
19
	s.Var(pw, fmp)
20
	s.Spw_Form(pw) // Создать spw
21

22
	// # 1, 2, 3, 4
23
	s.SET_Form() // SSA - анализ сегментов pw
24

25
	// # 5
26
	// Оценка АКФ сингулярных троек для сегментов pw
27
	// Визуализация АКФ сингулярных троек для сегментов pw
28
	s.AKF_Form() // Оценка АКФ сингулярных троек для сегментов pw
29

30
	// # 6, 7
31
	// Огибающие АКФ сингулярных троек sET12 сегментов pw
32
	// Нормированные АКФ сингулярных троек sET12 сегментов pw
33
	s.Envelope()
34

35
	// # 8
36
	// Мгновенная частота нормированной АКФ сингулярных троек sET12 для сегментов pw
37
	s.MomentFrequency()
38

39
	// # 9
40
	// Визуализация СПМ сингулярных троек сегменов pw
41
	s.VisibleSPM()
42

43
	// # 10
44
	// Агрегирование сегментов очищенной пульсовой волны cpw
45
	s.AggregationPW()
46

47
}
48

49
/*
50
	func SSA_spw(pw, fmp []float64) {
51
		defer timer("SSA_spw")()
52
		// Сегменты отсчётов pw
53
		N := len(pw) // Количество отсчетов pw
54
		win := 1024
55
		res := N - win*int(math.Floor(float64(N)/float64(win)))
56
		nPart := 20 // Количество долей res
57
		res = int(math.Floor(float64(res) / float64(nPart)))
58
		//overlap := (float64(win) - float64(res)) / float64(win)
59
		S := 1
60
		Imin := 1
61
		Imax := win
62

63
		var ns []float64
64
		for Imax <= N {
65
			ns = append(ns, float64(S)) // номер текущего сегмента pw
66
			Imin = Imin + res
67
			Imax = Imax + res
68
			S++
69
		}
70
		S-- // кол-во перекрывающихся сегментов pw в пределах N
71

72
		spw := mat.NewDense(win, S, nil)
73
		//fmt.Println("Размеры spw:", win, S)
74
		for j := 0; j < S; j++ {
75
			for i := 0; i < win; i++ {
76
				k := (j) * res
77
				spw.Set(i, j, pw[k+i])
78
			}
79
		}
80

81
		// Set general parameters
82
		cad := 30              // 30 кадров/сек
83
		dt := 1 / float64(cad) // интервал дискретизации времени, сек
84
		tim := make([]float64, N)
85
		for index := 1; index < N; index++ {
86
			tim[index] = tim[index-1] + dt
87
		}
88

89
		ns2 := make([]int, S)
90
		for index := range ns2 {
91
			ns2[index] = (index + 1)
92
		}
93

94
		L := make([]float64, S)
95
		for index := range L { // цикл по сегментам pw
96
			L[index] = math.Floor(float64(cad) / fmp[index]) // кол-во отсчетов основного тона pw
97
		}
98

99
		K := 5
100
		M := int(float64(K) * oss.Max(L)) // параметр вложения в траекторное пространство
101
		// SSA - анализ сегментов pw
102
		seg := 100 // номер сегмента pw для визуализации
103
		nET := 4   // кол-во сингулярных троек для сегментов pw
104

105
		//var sET12 mat.Dense
106
		sET12_sum2 := mat.NewDense(win, 2, nil) // НЕ ФАКТ, ЧТО К-во строк win
107
		sET34_sum2 := mat.NewDense(win, 2, nil) // НЕ ФАКТ, ЧТО К-во строк win
108
		sET12 := mat.NewDense(win, S, nil)      // НЕ ФАКТ, ЧТО К-во строк win
109
		sET34 := mat.NewDense(win, S, nil)      // НЕ ФАКТ, ЧТО К-во строк win
110

111
		for j := 0; j < S; j++ {
112
			C, LBD, RC := ssa.SSA(win, M, spw.ColView(j), nET)
113
			//fmt.Println(j, S)
114
			RC_T := mat.DenseCopyOf(RC.T())
115

116
			sET12_sum2.SetCol(0, RC_T.RawRowView(0))
117
			sET12_sum2.SetCol(1, RC_T.RawRowView(1))
118
			sET12.SetCol(j, oss.Sum2(*sET12_sum2))
119
			sET12_sum2.Zero()
120

121
			sET34_sum2.SetCol(0, RC_T.RawRowView(2))
122
			sET34_sum2.SetCol(1, RC_T.RawRowView(3))
123
			sET34.SetCol(j, oss.Sum2(*sET34_sum2))
124
			sET34_sum2.Zero()
125

126
			if j == seg {
127
				graph.Imagesc(C, "C")
128
				oss.Matlab_mat_Dense(&C, 1, "C", "")
129
				log.Println("Covariance matrix")
130
				graph.MakeGraphOfArray(LBD, "LBD")
131

132
				oss.Matlab_arr_float(LBD, 2, "LBD", "")
133
				log.Println("Eigenvalues")
134

135
				err_makeGraphYX_sET12 := graph.MakeGraphYX_VecDense(
136
					*mat.NewVecDense(win, tim[0:win]),
137
					*(mat.VecDenseCopyOf(spw.ColView(j))),
138
					*(mat.NewVecDense(len(oss.Vec_in_ArrFloat(sET12.ColView(j))), oss.Vec_in_ArrFloat(sET12.ColView(j)))),
139
					"origin", "sET12")
140
				oss.Matlab_arr_float(tim, 3, "tim", "")
141
				oss.Matlab_mat_Dense(spw, 3, "spw", "")
142
				oss.Matlab_mat_Dense(sET12, 3, "sET12", "")
143
				log.Println("Original time series and reconstruction sET12")
144

145
				err_makeGraphYX_sET34 := graph.MakeGraphYX_VecDense(
146
					*mat.NewVecDense(win, tim[0:win]),
147
					*(mat.VecDenseCopyOf(spw.ColView(j))),
148
					*(mat.NewVecDense(len(oss.Vec_in_ArrFloat(sET34.ColView(j))), oss.Vec_in_ArrFloat(sET34.ColView(j)))),
149
					"origin", "sET34")
150
				oss.Matlab_arr_float(tim, 4, "tim", "")
151
				oss.Matlab_mat_Dense(spw, 4, "spw", "")
152
				oss.Matlab_mat_Dense(sET34, 4, "sET34", "")
153
				log.Println("Original time series and reconstruction sET34")
154

155
				if err_makeGraphYX_sET12 != nil {
156
					fmt.Println(err_makeGraphYX_sET12)
157
				}
158
				if err_makeGraphYX_sET34 != nil {
159
					fmt.Println(err_makeGraphYX_sET34)
160
				}
161
			}
162
		}
163

164
		oss.SafeToXlsxMatrix(sET12, "sET12")
165
		oss.SafeToXlsxMatrix(sET34, "sET34")
166

167
		// *****************
168
		// Оценка АКФ сингулярных троек для сегментов pw
169
		lag := int(math.Floor(float64(win) / 10.0)) // % наибольший лаг АКФ <= win/10
170
		lagS := 2 * lag
171
		Acf_sET12 := ssa.ACF_estimation_of_singular_triples(lagS, win, S, *sET12)
172
		oss.SafeToXlsxM(*Acf_sET12, "Acf_sET12")
173
		// *****************
174
		// Визуализация АКФ сингулярных троек для сегментов pw
175
		lgl := make([]float64, lag)
176
		for m := 0; m < len(lgl); m++ {
177
			lgl[m] = float64(m + 1)
178
		}
179
		time := make([]float64, lag)
180
		for m := 1; m < len(time); m++ {
181
			time[m] = time[m-1] + dt
182
		}
183
		oss.Matlab_arr_float(ns, 5, "ns", "")
184
		oss.Matlab_arr_float(time, 5, "time", "")
185
		oss.Matlab_mat_Dense(Acf_sET12, 5, "Acf_sET12", "")
186
		log.Println("Визуализация АКФ сингулярных троек для сегментов pw")
187

188
		// *****************
189
		// Огибающая по критерию локальных максимумов abs(acf_sET12)
190
		//power := 0.75 // параметр спрямляющего преобразования
191
		EnvAcf_sET12 := *mat.NewDense(lag, S, nil)
192
		AcfNrm_sET12 := *mat.NewDense(lag, S, nil)
193
		//for j := 16; j <= 16; j++ { // цикл по сегментам АКФ
194
		for j := 0; j < S; j++ { // цикл по сегментам АКФ
195
			Acf_sET12_col := *mat.VecDenseCopyOf(Acf_sET12.ColView(j))
196
			absTS := oss.AbsVector(Acf_sET12_col)
197
			at1 := absTS.AtVec(0)
198
			at2 := absTS.AtVec(1)
199

200
			maxTS := *mat.NewVecDense(lag, nil)
201
			maxTS.SetVec(0, at1)
202

203
			maxN := *mat.NewVecDense(lag, nil)
204
			maxN.SetVec(0, 1)
205

206
			var Nmax int
207

208
			for m := 2; m < lag; m++ {
209
				at3 := absTS.AtVec(m)
210
				if (at1 <= at2) && (at2 >= at3) {
211
					Nmax++                        // номер очередного узла интерполяции (счетчик максимумов)
212
					maxN.SetVec(Nmax, float64(m)) // номер очередного максимума для ряда absTS
213
					maxTS.SetVec(Nmax, at2)       // отсчет очередного узла интерполяции
214
				}
215
				at1 = at2
216
				at2 = at3
217
			}
218
			Nmax++                               // количество узлов интерполяции
219
			maxN.SetVec(Nmax, float64(lag))      // номер отсчета absTS финального узла интерполяции
220
			maxTS.SetVec(Nmax, absTS.AtVec(lag)) // отсчет absTS финального узла интерполяции
221
			NumMax := maxN.SliceVec(0, Nmax+1)
222

223
			// Интерполяция огибающей АКФ
224

225
			// acfEnvelope := pchip.Pchip(Vec_in_ArrFloat(NumMax),
226
			// 	oss.Vec_in_ArrFloat(maxTS.SliceVec(0, Nmax+1)),
227
			// 	(lgl),
228
			// 	NumMax.Len(), len(lgl))
229

230
			acfEnvelope, _, _ := pchip.Pchip(oss.Vec_in_ArrFloat(NumMax), oss.Vec_in_ArrFloat(maxTS.SliceVec(0, Nmax+1)), lgl, NumMax.Len(), len(lgl))
231

232
			//fmt.Println(len(lgl), NumMax.Len())
233
			//acfEnvelope := pchip.Pchip2(lgl, oss.Vec_in_ArrFloat(NumMax))
234

235
			//fmt.Println(Nmax, len(lgl))
236

237
			EnvAcf_sET12.SetCol(j, acfEnvelope)
238

239
			// нормированные АКФ
240
			AcfNrm_sET12.SetCol(j, oss.VecDense_in_float64(oss.Vector_DivElemVec((Acf_sET12.Slice(0, lag, j, j+1)), EnvAcf_sET12.ColView(j))))
241

242
			//fmt.Println(AcfNrm_sET12.At(lag-1, j)) // Тут возникает 850+ для 16-ти
243
		}
244

245
		// Обход ошибки вывода с 856, заменив последнюю строку
246
		EnvAcf_sET12 = oss.EditLastRow(EnvAcf_sET12)
247
		AcfNrm_sET12 = oss.EditLastRow(AcfNrm_sET12)
248

249
		// *****************
250
		oss.SafeToXlsxM(EnvAcf_sET12, "EnvAcf_sET12")
251
		oss.SafeToXlsxM(AcfNrm_sET12, "AcfNrm_sET12")
252

253
		// 6 - Огибающие АКФ сингулярных троек sET12 сегментов pw
254
		oss.Matlab_arr_float(ns, 6, "ns", "")
255
		oss.Matlab_arr_float(time, 6, "time", "")
256
		oss.Matlab_mat_Dense(&EnvAcf_sET12, 6, "EnvAcf_sET12", "")
257
		log.Println("Огибающие АКФ сингулярных троек sET12 сегментов pw")
258
		graph.SaveDat_2(EnvAcf_sET12, "File_For_MatLab"+oss.OpSystemFilder+strconv.Itoa(6)+oss.OpSystemFilder+"EnvAcf_sET12"+".dat")
259
		graph.SaveDat(ns, "File_For_MatLab"+oss.OpSystemFilder+strconv.Itoa(6)+oss.OpSystemFilder+"ns"+".dat")
260
		graph.SaveDat(time, "File_For_MatLab"+oss.OpSystemFilder+strconv.Itoa(6)+oss.OpSystemFilder+"time"+".dat")
261

262
		// 7 - Нормированные АКФ сингулярных троек sET12 сегментов pw
263
		oss.Matlab_arr_float(ns, 7, "ns", "")
264
		oss.Matlab_arr_float(time, 7, "time", "")
265
		oss.Matlab_mat_Dense(&AcfNrm_sET12, 7, "AcfNrm_sET12", "")
266
		Folder7 := "File_For_MatLab" + oss.OpSystemFilder + strconv.Itoa(7) + oss.OpSystemFilder
267
		graph.SaveDat_2(AcfNrm_sET12, Folder7+"AcfNrm_sET12"+".dat")
268
		graph.SaveDat(ns, Folder7+"ns"+".dat")
269
		graph.SaveDat(time, Folder7+"time"+".dat")
270
		graph.SplotMatrixFromFile(graph.Option3D{ // Задаём настройки 3D графика
271
			FileNameDat: Folder7 + "AcfNrm_sET12.dat",
272
			FileNameOut: Folder7 + "AcfNrm_sET12.png",
273
			Titile:      "Нормированные АКФ сингулярных троек sET12 сегментов pw",
274
			Xlabel:      "ns",
275
			Ylabel:      "lag,s",
276
			Zlabel:      "Acf_Nrm",
277
		}) // Делаем график
278
		log.Println("Нормированные АКФ сингулярных троек sET12 сегментов pw")
279

280
		// ********************************************************************
281
		// Мгновенная частота нормированной АКФ сингулярных троек sET12 для сегментов pw
282
		insFrc_AcfNrm, insFrc_AcfNrmErr := Instantaneous_frequency_of_normalized_ACF_sET12(AcfNrm_sET12, S, lag, dt, lgl)
283
		if insFrc_AcfNrmErr != nil {
284
			log.Println(insFrc_AcfNrmErr)
285
		}
286

287
		// filter savitzky-goley
288
		filter, savitzky_goley_Error := savitzkygolay.NewFilterWindow(53)
289
		if savitzky_goley_Error != nil {
290
			log.Fatalln(savitzky_goley_Error)
291
		}
292
		smo_insFrc_AcfNrm, _ := filter.Process(insFrc_AcfNrm, lgl)
293
		//smo_insFrc_AcfNrm := insFrc_AcfNrm
294

295
		//smo_insFrc_AcfNrm := savitzky_goley(insFrc_AcfNrm, 33, 2)
296

297
		oss.Matlab_arr_float(ns, 8, "ns", "")
298
		oss.Matlab_arr_float(insFrc_AcfNrm, 8, "insFrc_AcfNrm", "")
299
		oss.Matlab_arr_float(smo_insFrc_AcfNrm, 8, "smo_insFrc_AcfNrm", "")
300
		err_insFrc_AcfNrm := graph.MakeGraphYX_float64(
301
			insFrc_AcfNrm,
302
			ns,
303
			"insFrc_AcfNrm")
304
		if err_insFrc_AcfNrm != nil {
305
			fmt.Println(err_insFrc_AcfNrm)
306
		}
307
		err_insFrc_AcfNrm = graph.MakeGraphYX_float64(
308
			smo_insFrc_AcfNrm,
309
			ns,
310
			"smo_insFrc_AcfNrm")
311
		if err_insFrc_AcfNrm != nil {
312
			fmt.Println(err_insFrc_AcfNrm)
313
		}
314
		err_insFrc_AcfNrm = graph.MakeGraphYX_VecDense(
315
			*mat.NewVecDense(len(ns), ns),
316
			*mat.NewVecDense(len(insFrc_AcfNrm), insFrc_AcfNrm),
317
			*mat.NewVecDense(len(smo_insFrc_AcfNrm), smo_insFrc_AcfNrm),
318
			"origin", "insFrc_AcfNrm")
319
		if err_insFrc_AcfNrm != nil {
320
			fmt.Println(err_insFrc_AcfNrm)
321
		}
322

323
		// Оценки СПМ сингулярных троек для сегменов pw
324
		var iGmin, iGmax int
325
		smopto := 3 // параметр сглаживания периодограммы Томсона
326
		// Визуализация СПМ сингулярных троек сегменов pw
327
		fmi := 40.0 / 60.0                                            // частота среза для 40 уд/мин (0.6667 Гц)
328
		fma := 240.0 / 60.0                                           // частота среза для 240 уд/мин (4.0 Гц)
329
		Nf := 1 + win/2                                               // кол-во отсчетов частоты
330
		df := float64(cad) / float64(win-1)                           // интервал дискретизации частоты, Гц
331
		Fmin := fmi - float64(10*df)                                  // частота в Гц, min
332
		Fmax := fma + float64(10*df)                                  // частота в Гц, max
333
		pto_sET12 := pto_sET12_init(*sET12, *spw, smopto, win, Nf, S) // Расчёт оценки СПМ сингулярных троек для сегменов pw
334
		oss.SafeToXlsxMatrix(pto_sET12, "pto_sET12")                  // Сохранить в Xlsx матрицу оценки СПМ
335

336
		f := make([]float64, Nf)
337
		for i := 2; i < Nf; i++ {
338
			f[i] = f[i-1] + df // частота в герцах
339
			if math.Abs(f[i]-Fmin) <= df {
340
				iGmin = i
341
			}
342
			if math.Abs(f[i]-Fmax) <= df {
343
				iGmax = i
344
			}
345
		}
346
		fG := make([]float64, iGmax)
347
		for i := 0; i < iGmax; i++ {
348
			fG[i] = f[i] // сетка частот 3D-графика
349
		}
350
		oss.Matlab_arr_float(ns, 9, "ns", "")
351
		oss.Matlab_arr_float(fG, 9, "fG", "")
352
		oss.Matlab_mat_Dense(pto_sET12, 9, "pto_sET12", "")
353
		oss.Matlab_variable(iGmin, 9, "iGmin", "")
354
		oss.Matlab_variable(iGmax, 9, "iGmax", "")
355

356
		// Оценки средних частот основного тона сингулярных троек сегментов pw
357
		pto_fMAX12 := make([]float64, S)
358
		for index := range pto_fMAX12 {
359
			_, I := oss.MaxArrFloat64(oss.Vec_in_ArrFloat(pto_sET12.ColView(index))) // Поиск индекса максимального значения массива
360
			pto_fMAX12[index] = f[I]
361
		}
362
		oss.Matlab_arr_float(ns, 10, "ns", "")
363
		oss.Matlab_arr_float(pto_fMAX12, 10, "pto_fMAX12", "")
364
		err_pto_fMAX12 := graph.MakeGraphYX_float64(
365
			pto_fMAX12,
366
			ns,
367
			"pto_fMAX12")
368
		if err_pto_fMAX12 != nil {
369
			fmt.Println(err_pto_fMAX12)
370
		}
371
		Folder9 := "File_For_MatLab" + oss.OpSystemFilder + strconv.Itoa(9) + oss.OpSystemFilder
372
		graph.SaveDat_2(*pto_sET12, Folder9+"pto_sET12"+".dat")
373
		graph.SaveDat(ns, Folder9+"ns"+".dat")
374
		graph.SaveDat(time, Folder9+"time"+".dat")
375

376
		graph.SplotMatrixFromFile(graph.Option3D{ // Задаём настройки 3D графика
377
			FileNameDat: Folder9 + "pto_sET12.dat",
378
			FileNameOut: Folder9 + "pto_sET12.png",
379
			Titile:      "Периодограмма Блэкмана–Харриса sET12 сегментов pw",
380
			Xlabel:      "ns",
381
			Ylabel:      "f,Hz",
382
			Zlabel:      "P(f)",
383
		}) // Делаем график
384

385
		log.Println("Нормированные АКФ сингулярных троек sET12 сегментов pw")
386

387
		// oss.SafeToXlsx(f, "f")
388

389
		// ***
390
		// Агрегирование сегментов очищенной пульсовой волны cpw
391
		NSF := win + res*(S-1) // номер финального отсчета финального сегмента <= N
392
		NumS, cpw_avr, cpw_med, cpw_iqr := wav(NSF, S, win, res, *sET12)
393
		oss.SafeToXlsx(NumS, "NumS")
394

395
		oss.Matlab_variable(NSF, 10, "NSF", "")
396
		oss.Matlab_arr_float(tim, 10, "tim", "")
397
		oss.Matlab_arr_float(cpw_avr, 10, "cpw_avr", "")
398
		oss.Matlab_arr_float(cpw_med, 10, "cpw_med", "")
399
		oss.Matlab_arr_float(cpw_iqr, 10, "cpw_iqr", "")
400
	}
401
*/
402
func timer(name string) func() {
403
	start := time.Now()
404
	return func() {
405
		fmt.Printf("%s took %v\n", name, time.Since(start))
406
	}
407
}
408
func wav(N, S, W, res int, sET mat.Dense) ([]float64, []float64, []float64, []float64) {
409
	NS := make([]float64, N)
410
	w_avr := make([]float64, N)
411
	w_med := make([]float64, N)
412
	w_iqr := make([]float64, N)
413

414
	ET := mat.NewDense(N, S, nil)
415

416
	for j := 0; j < S; j++ { // цикл по сегментам
417
		for i := 0; i < W; i++ {
418
			k := (j) * res
419
			//fmt.Print(i, k, j, "=")
420
			//fmt.Println(sET.At(i, j))
421
			ET.Set(i+k, j, sET.At(i, j)) // сдвинутый сегмент ET(:,j)
422
		}
423
	}
424

425
	Smp := make([]float64, N*S)
426
	for i := 0; i < N; i++ {
427
		var nSi int
428
		for j := 0; j < S; j++ {
429
			if ET.At(i, j) != 0.0 {
430
				nSi++
431
				Smp[nSi] = ET.At(i, j)
432
			}
433
		}
434
		NS[i] = float64(nSi)                      // кол-во сегментов для текущего i
435
		w_avr[i] = oss.Mean(Smp[:nSi])            // выборочная средняя
436
		w_med[i] = oss.Median_floatArr(Smp[:nSi]) // медиана
437
		w_iqr[i] = (oss.Prctile(Smp[:nSi], 75) - oss.Prctile(Smp[:nSi], 25)) / 2.0
438
	}
439

440
	return NS, w_avr, w_med, w_iqr
441
}
442

443
// Формирование оценки СПМ сингулярных троек для сегменов pw
444
func pto_sET12_init(sET12 mat.Dense, spw mat.Dense, smopto, win, Nf, S int) *mat.Dense {
445
	pto_sET12 := mat.NewDense(Nf, S, nil)
446

447
	// Расчёт окна Блэкмана_Харриса шириной win
448
	// и с заданными коэффициентами
449
	BlackManHar := blackmanharris.Blackmanharris(win, blackmanharris.Koef4_74db)
450

451
	for j := 0; j < S; j++ {
452
		// Периодограмма Блэкмана_Харриса
453
		// pto_sET12(:,j) = periodogram(spw(:,j),blackmanharris(win),win); % Блэкмана-Харриса
454
		pto_sET12.SetCol(j, periodogram.Periodogram(oss.Vec_in_ArrFloat(spw.ColView(j)), BlackManHar, win))
455

456
		// Периодограмма Томсона
457
		// pto_sET12.SetCol(j, pmtmMy(sET12.ColView(j), smopto, win))
458
	}
459
	return pto_sET12
460
}
461

462
// func pmtmMy(sET12 mat.Vector, smopto, win int) []float64 {
463
// 	return pmtm.Pmtm(oss.Vec_in_float64(sET12), 1024)
464
// }
465

466
// Расчёты вектора PhaAcfNrm, модуль от Акосинуса.
467
func MakePhaAcfNrm(vect mat.Vector) *mat.VecDense {
468
	output := mat.VecDenseCopyOf(vect)
469
	for i := 0; i < output.Len(); i++ {
470
		output.SetVec(i, math.Abs(math.Acos(output.AtVec(i))))
471
	}
472
	return output
473
}
474