#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Jan 15 19:45:33 2018

@author: anne
"""
import numpy as np
from lasp import PowerSpectra, Window
import matplotlib.pyplot as plt
plt.close('all')
# def test_ps():

nfft = 8
t = np.linspace(0, 1.0, nfft, endpoint=False)

ps = PowerSpectra(nfft, Window.Rectangular)

sig = np.random.randn(nfft)
freq = 4
omg = 2*np.pi*freq
# sig = 8*np.cos(omg*t)


cps = ps.compute(sig)

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)
# test_ps()
# plt.plot(res_lasp.real-res_npy.real)
# plt.plot(res_lasp.imag-res_npy.imag)
# plt.plot(res_npy.real)
# plt.plot(res_npy.imag)
    # plt.plot(t, sig)
# print('nfft:',nfft)
# #print(nfft)
# nchannels = 2

# t = np.linspace(0,1,nfft+1)
# # print(t)
# x1 = (np.cos(4*np.pi*t[:-1])+3.2*np.sin(6*np.pi*t[:-1]))[:,np.newaxis]+10
# x = np.vstack([x1.T]*nchannels).T
# # Using transpose to get the strides right
# x = np.random.randn(nchannels,nfft).T
# print("strides: ",x.strides)
# # x.strides = (8,nfft*8)x
# # print("signal:",x)

# xms = np.sum(x**2,axis=0)/nfft
# print('Total signal power time domain: ', xms)

# X = np.fft.rfft(x,axis=0)
# # X =np.fft.fft(x)
# #X =np.fft.rfft(x)

# # print(X)
# Xs = 2*X/nfft
# Xs[np.where(np.abs(Xs) < 1e-10)] = 0
# Xs[0] /= np.sqrt(2)
# Xs[-1] /= np.sqrt(2)
# # print('single sided amplitude spectrum:\n',Xs)


# power = Xs*np.conj(Xs)/2

# # print('Frequency domain signal power\n', power)
# print('Total signal power', np.sum(power,axis=0).real)

# pstest = PowerSpectra(nfft,nchannels)
# ps = pstest.compute(x)

# fft = Fft(nfft,nchannels)
# fft.fft(x)

# ps[np.where(np.abs(ps) < 1e-10)] = 0+0j

# print('our ps: \n' , ps)

# print('Our total signal power: ',np.sum(ps,axis=0).real)


# if __name__ == '__main__':
#     nfft=2048        
#     fs = 2048
#     ps = PowerSpectra(nfft, Window.Rectangular)

#     t = np.linspace(0, 1.0, nfft, endpoint=False)
#     freq = 10
#     omg = 2*np.pi*freq
#     sig = np.sin(omg*t)
    
#     res = ps.compute(sig)
    
#     plt.plot(res[:,0,0])
#     # plt.plot(t, sig)
# print('nfft:',nfft)
# #print(nfft)
# nchannels = 2

# t = np.linspace(0,1,nfft+1)
# # print(t)
# x1 = (np.cos(4*np.pi*t[:-1])+3.2*np.sin(6*np.pi*t[:-1]))[:,np.newaxis]+10
# x = np.vstack([x1.T]*nchannels).T
# # Using transpose to get the strides right
# x = np.random.randn(nchannels,nfft).T
# print("strides: ",x.strides)
# # x.strides = (8,nfft*8)x
# # print("signal:",x)

# xms = np.sum(x**2,axis=0)/nfft
# print('Total signal power time domain: ', xms)

# X = np.fft.rfft(x,axis=0)
# # X =np.fft.fft(x)
# #X =np.fft.rfft(x)

# # print(X)
# Xs = 2*X/nfft
# Xs[np.where(np.abs(Xs) < 1e-10)] = 0
# Xs[0] /= np.sqrt(2)
# Xs[-1] /= np.sqrt(2)
# # print('single sided amplitude spectrum:\n',Xs)


# power = Xs*np.conj(Xs)/2

# # print('Frequency domain signal power\n', power)
# print('Total signal power', np.sum(power,axis=0).real)

# pstest = PowerSpectra(nfft,nchannels)
# ps = pstest.compute(x)

# fft = Fft(nfft,nchannels)
# fft.fft(x)

# ps[np.where(np.abs(ps) < 1e-10)] = 0+0j

# print('our ps: \n' , ps)

# print('Our total signal power: ',np.sum(ps,axis=0).real)