2018-01-29 15:14:50 +00:00
|
|
|
#!/usr/bin/env python3
|
|
|
|
# -*- coding: utf-8 -*-
|
|
|
|
"""
|
2022-10-05 12:57:39 +00:00
|
|
|
Testing code for power spectra
|
2018-01-29 15:14:50 +00:00
|
|
|
"""
|
|
|
|
import numpy as np
|
2022-09-03 18:59:14 +00:00
|
|
|
from lasp import PowerSpectra, Window
|
2022-10-05 12:57:39 +00:00
|
|
|
|
|
|
|
def test_ps():
|
|
|
|
"""
|
|
|
|
Check Parsevall for single-sided power spectra
|
|
|
|
"""
|
|
|
|
nfft = 2048
|
|
|
|
t = np.linspace(0, 1.0, nfft, endpoint=False)
|
2022-09-03 18:59:14 +00:00
|
|
|
|
2022-10-05 12:57:39 +00:00
|
|
|
ps = PowerSpectra(nfft, Window.WindowType.Rectangular)
|
2022-09-03 18:59:14 +00:00
|
|
|
|
2022-10-05 12:57:39 +00:00
|
|
|
sig = np.random.randn(nfft)
|
|
|
|
freq = 4
|
|
|
|
omg = 2*np.pi*freq
|
|
|
|
# sig = 8*np.cos(omg*t)
|
|
|
|
|
|
|
|
|
2022-10-11 12:50:44 +00:00
|
|
|
cps = ps.compute(sig[:,None])
|
2022-10-05 12:57:39 +00:00
|
|
|
|
|
|
|
pow1 = np.sum(sig**2)/sig.size
|
|
|
|
pow2 = np.sum((cps[:,0,0]).real)
|
|
|
|
|
|
|
|
# print(pow1)
|
|
|
|
# print(pow2)
|
|
|
|
# plt.plot(cps[:,0,0])
|
|
|
|
assert np.isclose(pow2 - pow1,0, atol=1e-1)
|
|
|
|
|
|
|
|
if __name__ == '__main__':
|
|
|
|
test_ps()
|
2022-09-03 18:59:14 +00:00
|
|
|
|