ssa

1
package pmtm
2

3
import (
4
	"math/cmplx"
5

6
	"github.com/mjibson/go-dsp/fft"
7
)
8

9
func Pmtm(x []float64, n int) []float64 {
10
	/*
11
			// Размер массива
12
			N := len(x)
13
			// Расчет преобразования Фурье
14
			X := make([]float64, N)
15
			for k := 0; k < N; k++ {
16
				for i := 0; i < N; i++ {
17
					X[k] += x[i] * math.Cos(2*math.Pi*float64(k)*float64(i)/float64(N))
18
				}
19
			}
20
			// Расчет периодограммы Томсона
21
			P := make([]float64, n)
22
			for k := 0; k < n; k++ {
23
				P[k] = math.Pow(X[k], 2) / float64(N)
24
				for i := 1; i < N/2; i++ {
25
					fmt.Println("len(P)", len(P), "len(X)", len(X), k, k+i)
26
					P[k] += 2 * math.Pow(X[k+i], 2) / float64(N)
27
				}
28
			}
29
		return P
30
	*/
31
	//Fs := 1000.0 // Sampling frequency
32
	//T := 1 / Fs  // Sampling period
33
	L := 1024 // Length of signal
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, L/2+1)
45
	for i := 0; i < len(P1); i++ {
46
		P1[i] = 2 * P2[i]
47
	}
48

49
	return P2
50
}
51