ssa
1package pmtm2
3import (4"math/cmplx"5
6"github.com/mjibson/go-dsp/fft"7)
8
9func Pmtm(x []float64, n int) []float64 {10/*11// Размер массива
12N := len(x)
13// Расчет преобразования Фурье
14X := make([]float64, N)
15for k := 0; k < N; k++ {
16for i := 0; i < N; i++ {
17X[k] += x[i] * math.Cos(2*math.Pi*float64(k)*float64(i)/float64(N))
18}
19}
20// Расчет периодограммы Томсона
21P := make([]float64, n)
22for k := 0; k < n; k++ {
23P[k] = math.Pow(X[k], 2) / float64(N)
24for i := 1; i < N/2; i++ {
25fmt.Println("len(P)", len(P), "len(X)", len(X), k, k+i)
26P[k] += 2 * math.Pow(X[k+i], 2) / float64(N)
27}
28}
29return P
30*/
31//Fs := 1000.0 // Sampling frequency32//T := 1 / Fs // Sampling period33L := 1024 // Length of signal34
35// fft36fft_signal := fft.FFTReal(x)37//fmt.Println(len(fft_signal))38P2 := make([]float64, len(fft_signal))39for i := 0; i < len(P2); i++ {40// P2[i] = cmplx.Abs((fft_signal[i])*(fft_signal[i])) / 1024.041P2[i] = cmplx.Abs((fft_signal[i]) / 1024.0)42}43// Возвести в квадрат44P1 := make([]float64, L/2+1)45for i := 0; i < len(P1); i++ {46P1[i] = 2 * P2[i]47}48
49return P250}
51