ssa

1
// # periodogram
2
// Оценка спектральной плотности мощности периодограммой
3
//
4
// В этом пакете мы объявляем функцию periodogram,
5
// которая принимает на вход сигнал x в виде среза []float64,
6
// окно window в виде среза []float64 и количество точек nfft для быстрого преобразования Фурье.
7
// Функция возвращает периодограмму в виде среза []float64.
8
//
9
// Мы начинаем с создания среза pxx для хранения периодограммы итерационно перебираем сигнал x в блоках размера nfft.
10
// Для каждого из этих блоков мы применяем окно window покомпонентно, умножая его на элементы блока x.
11
// Затем мы вычисляем быстрое преобразование Фурье (FFT) для окончательного преобразования блока x в частотный домен.
12
// Мы вычисляем квадрат модуля каждого элемента FFT и добавляем его к соответствующему элементу pxx.
13
//
14
// Также мы объявляем функцию fft, которая рекурсивно вычисляет быстрое преобразование Фурье для сигнала x.
15
// Эта функция работает за O(n log n) и использует формулу Баттерворта-Тьюки.
16
// Мы начинаем с проверки базового случая, когда n равно 1, и возвращаем значение x в качестве комплексного числа.
17
// В противном случае мы разделяем входной сигнал на две половины even и odd и рекурсивно вычисляем FFT для каждой половины.
18
// Затем мы объединяем два FFT с помощью преобразования Баттерворта-Тьюки и возвращаем результат.
19
package periodogram
20

21
import (
22
	"math"
23
	"math/cmplx"
24

25
	"github.com/mjibson/go-dsp/fft"
26
)
27

28
// Расчёт модифицированной периодограммы с окном
29
func Periodogram(x []float64, window []float64, nfft int) []float64 {
30
	pxx := make([]float64, nfft/2+1)
31
	for i := 0; i <= len(x)-nfft; i += nfft {
32
		block := x[i : i+nfft-1]
33

34
		//fmt.Println(i, len(x)-nfft, len(block), len(window))
35
		for j := range block {
36
			block[j] *= window[j]
37
		}
38
		FFT := fft.FFTReal(block)
39

40
		// fmt.Println(len(FFT), FFT)
41

42
		for j := range pxx {
43
			abs := cmplx.Abs(FFT[j])
44
			pxx[j] += math.Pow(abs, 2) / float64(nfft)
45
		}
46
	}
47
	return pxx
48
}
49