ssa

1
package pmtm
2

3
import "math"
4

5
// Pxx = PMTM(X,NW,NFFT) указывает длину БПФ, используемую для вычисления оценок PSD.
6
// Для реального X Pxx имеет (NFFT/2+1) строк, если NFFT четное, и (NFFT + 1)/2 строк,
7
// если NFFT нечетное. Для сложного X значение Pxx всегда имеет длину NFFT.
8
// Если значение NFFT указано как пустое, значение NFFT устанавливается равным либо 256,
9
// либо следующей степени 2, превышающей длину X, в зависимости от того, что на % больше.
10
func Pmtm(x []float64, NW int, NFFT int) []float64 {
11

12
	// // Размер массива
13
	// N := len(x)
14
	// // Расчет преобразования Фурье
15
	// X := make([]float64, N)
16
	// for k := 0; k < N; k++ {
17
	// 	for i := 0; i < N; i++ {
18
	// 		X[k] += x[i] * math.Cos(2*math.Pi*float64(k)*float64(i)/float64(N))
19
	// 	}
20
	// }
21
	// // Расчет периодограммы Томсона
22
	// P := make([]float64, (NFFT/2 + 1))
23
	// for k := 0; k < (NFFT/2 + 1); k++ {
24
	// 	P[k] = math.Pow(X[k], 2) / float64(N)
25
	// 	for i := 1; i < N/2; i++ {
26
	// 		// fmt.Println("len(P)", len(P), "len(X)", len(X), k, k+i)
27
	// 		P[k] += 2 * math.Pow(X[k+i], 2) / float64(N)
28
	// 	}
29
	// }
30
	// return P
31

32
	//Fs := 1000.0 // Sampling frequency
33
	//T := 1 / Fs  // Sampling period
34

35
	// // fft
36
	// fft_signal := fft.FFTReal(x)
37
	// //fmt.Println(len(fft_signal))
38
	// P2 := make([]float64, len(fft_signal))
39
	// for i := 0; i < len(P2); i++ {
40
	// 	// P2[i] = cmplx.Abs((fft_signal[i])*(fft_signal[i])) / 1024.0
41
	// 	P2[i] = cmplx.Abs((fft_signal[i]) / 1024.0)
42
	// }
43
	// // Возвести в квадрат
44
	// P1 := make([]float64, (NFFT/2 + 1))
45
	// for i := 0; i < len(P1); i++ {
46
	// 	P1[i] = 2 * P2[i]
47
	// }
48
	// return P2
49

50
	// Вычислите размер Pxx в зависимости от типа nfft.
51
	var PxxLen int
52
	if NFFT%2 == 0 {
53
		PxxLen = NFFT/2 + 1
54
	} else {
55
		PxxLen = (NFFT + 1) / 2
56
	}
57

58
	// Создайте слайс Pxx с правильной длиной.
59
	Pxx := make([]float64, PxxLen)
60
	alpha := 2.0
61

62
	for k := 0; k < PxxLen; k++ {
63
		freq := float64(k+1) / float64(NFFT)
64
		sinSum := 0.0
65
		cosSum := 0.0
66

67
		for n := 0; n < len(x); n++ {
68
			sinSum += x[n] * math.Sin(2*math.Pi*freq*float64(n))
69
			cosSum += x[n] * math.Cos(2*math.Pi*freq*float64(n))
70
		}
71

72
		Sk := (sinSum*sinSum + cosSum*cosSum) / float64(len(x))
73

74
		Pxx[k] = Sk / (alpha * float64(NW))
75
	}
76

77
	return Pxx
78
}
79