Skip to content
Snippets Groups Projects
Commit abeb30df authored by Mads M. Pedersen's avatar Mads M. Pedersen
Browse files

fixed import errors in wetb.signal

parent a8e904c2
No related branches found
No related tags found
No related merge requests found
......@@ -4,8 +4,8 @@ Created on 13/07/2016
@author: MMPE
'''
import numpy as np
from wetb.signal_tools.filters.first_order import low_pass
from wetb.signal_tools.filters import replacer
from wetb.signal.filters.first_order import low_pass
from wetb.signal.filters import replacer
replace_by_nan = replacer.replace_by_nan
......
......@@ -4,7 +4,7 @@ Created on 10/01/2015
@author: mmpe
'''
import numpy as np
from wetb.signal_tools.filters import cy_filters
from wetb.signal.filters import cy_filters
def low_pass(input, delta_t, tau, method=1):
if isinstance(tau, (int, float)):
......
......@@ -58,12 +58,16 @@ def x2F(y, max_nfft, x=None):
F = np.r_[AB[0], (AB[1:nfft + 1] + 1j * AB[nfft + 1:]), np.zeros(nfft) ]
return F
def rx2F(y, max_fft):
def rx2F(y, max_nfft, x=None):
"""Convert non-complex signal, y, to single sided Fourier components, that satifies x(t) = sum(X(cos(iw)+sin(iw)), i=0..N)"""
d = np.arange(360)
if x is not None:
x,fit = bin_fit(x,y, d)
y = fit(d)
F = np.fft.rfft(y) / len(y)
F[1:-1] *= 2 # add negative side
F = np.conj(F)
return F[:max_fft + 1]
return F[:max_nfft + 1]
def rF2x(rF):
"""Convert single sided Fourier components, that satisfies x(t) = sum(X(cos(iw)+sin(iw)), i=0..N) to non-complex signal"""
......
......@@ -4,8 +4,13 @@ Created on 13. jan. 2017
@author: mmpe
'''
import unittest
from wetb.signal_tools.filters import first_order
from scipy import signal
import numpy as np
from wetb.signal.filters import first_order
from wetb.signal.fit._fourier_fit import F2x, x2F, rx2F
class Test_first_order_filters(unittest.TestCase):
......@@ -19,6 +24,153 @@ class Test_first_order_filters(unittest.TestCase):
plt.plot(a)
plt.plot(b)
plt.show()
def test_low_pass2(self):
t = np.linspace(0, 1.0, 2001)
xlow = np.sin(2 * np.pi * 5 * t)
xhigh = np.sin(2 * np.pi * 250 * t)
x = xlow + xhigh
b, a = signal.butter(8, 0.125)
y = signal.filtfilt(b, a, x, padlen=150)
self.assertAlmostEqual(np.abs(y - xlow).max(),0,4)
if 0:
import matplotlib.pyplot as plt
plt.plot(x)
plt.plot(y)
plt.show()
def test_low_pass3(self):
t = np.linspace(0, 1.0, 2001)
# xlow = np.sin(2 * np.pi * 5 * t)
# xhigh = np.sin(2 * np.pi * 250 * t)
# x = xlow + xhigh
x = np.sum([np.sin(t*x*2*np.pi) for x in range(1,200)],0)
cutoff = .2
b, a = signal.butter(8,cutoff)
w, h = signal.freqs(b, a)
y = signal.filtfilt(b, a, x)
F = rx2F(x, max_nfft=len(t))
tF = np.linspace(0, len(F)/t[-1],len(F))
Fy = rx2F(y, max_nfft=len(t))
if 1:
import matplotlib.pyplot as plt
# plt.plot(x)
# plt.plot(y)
# plt.show()
plt.plot(tF, np.abs(F))
plt.plot(tF, np.abs(Fy))
plt.xlim([0,260])
plt.show()
print (b,a)
#
# plt.plot(w, 20 * np.log10(abs(h)))
# plt.xscale('log')
# plt.title('Butterworth filter frequency response')
# plt.xlabel('Frequency [radians / second]')
# plt.ylabel('Amplitude [dB]')
# plt.margins(0, 0.1)
# plt.grid(which='both', axis='both')
# plt.axvline(cutoff, color='green') # cutoff frequency
# plt.show()
#
# def test_low_pass2(self):
# t = np.arange(100000)/10000
# dt = t[1]-t[0]
# hz2 = np.sin(t*2*2*np.pi) # frq = 2 hz or 2*2pi rad/s
# hz10 = np.sin(t*10*2*np.pi) # frq = 10 hz or 10*2pi rad/s
# hz20 = np.sin(t*100*2*np.pi) # frq = 100 hz or 100*2pi rad/s
# #sig = np.sum([np.sin(t*x*2*np.pi) for x in range(1,200)],0)
# sig = np.sum([np.sin(t*x*2*np.pi) for x in [1,5,10]],0)
# print (sig.shape)
# F = rx2F(sig, max_nfft=len(t))
# tF = np.linspace(0, len(F)/t[-1],len(F))
#
#
# sig_lp1 = first_order.low_pass(sig, dt, 30*dt)
# F_lp1 = rx2F(sig_lp1, max_nfft=len(t))
#
# bb, ba = signal.butter(1,5, 'low', analog=True)
# print (bb, ba)
# w1, h1 = signal.freqs(bb, ba)
# sig_lp2 = signal.filtfilt(bb,ba, sig)
# # print (sig_lp2)
# # F_lp2 = rx2F(sig_lp2, max_nfft=len(t))
# # bb, ba = signal.butter(10,10, 'low', analog=True)
# # sig_lp3 = signal.lfilter(bb,ba, sig)
# # F_lp3 = rx2F(sig_lp3, max_nfft=len(t))
# # w2, h2 = signal.freqs(bb, ba)
# # #
# # self.assertLess(b.std(), a.std())
# if 1:
# import matplotlib.pyplot as plt
# plt.plot(t, sig)
# # plt.plot(t, np.sum([1/x*np.sin(t*x*2*np.pi) for x in [1]],0))
# #
# plt.plot(t, sig_lp1)
# plt.plot(t, sig_lp2)
# #plt.plot(t, sig_lp3)
#
# plt.show()
#
# # plt.plot(tF, np.abs(F))
# # plt.plot(tF, np.abs(F_lp1))
# # # plt.plot(tF, np.abs(F_lp3))
# # # plt.plot(tF, np.abs(F_lp2))
# # plt.plot([0,1000],[0.708,.708])
# # plt.xlim([0,205])
# # plt.ylim([0,1])
# # plt.show()
# # plt.plot(w1, 20 * np.log10(abs(h1)))
# # plt.plot(w2, 20 * np.log10(abs(h2)))
# # plt.xscale('log')
# # plt.title('Butterworth filter frequency response')
# # plt.xlabel('Frequency [radians / second]')
# # plt.ylabel('Amplitude [dB]')
# # plt.margins(0, 0.1)
# # plt.grid(which='both', axis='both')
# # plt.axvline(100, color='green') # cutoff frequency
# # plt.show()
# # def test_low_pass2(self):
# # F = [0,1,0,0,0,0,0,0,0,0,0,0,0,0,.1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,.1]
# # F +=[0]*(len(F)-1)
# # x = F2x(F)
# # a = np.tile(x, 3)
# # print (x.shape)
# # # a = np.random.randint(0,100,100).astype(np.float)
# # b = first_order.low_pass(a, .1, 1)
# # bb, ba = signal.butter(10,100, 'low', analog=True)
# # #c = signal.lfilter(bb,ba, a)
# # w, h = signal.freqs(bb, ba)
# #
# # # self.assertLess(b.std(), a.std())
# # if 1:
# # import matplotlib.pyplot as plt
# # plt.plot(a)
# # #plt.ylim([-2,2])
# # plt.plot(b)
# # #plt.plot(c)
# # plt.show()
# # print (h)
# # plt.plot(w, abs(h))
# # plt.xscale('log')
# # plt.title('Butterworth filter frequency response')
# # plt.xlabel('Frequency [radians / second]')
# # plt.ylabel('Amplitude [dB]')
# # plt.margins(0, 0.1)
# # plt.grid(which='both', axis='both')
# # plt.axvline(100, color='green') # cutoff frequency
# # plt.show()
# #
# #
def test_high_pass(self):
a = np.random.randint(0,100,100).astype(np.float)
......
......@@ -5,7 +5,7 @@ Created on 12/11/2015
'''
import unittest
import numpy as np
from wetb.signal_tools.spectrum import psd
from wetb.signal.spectrum import psd
class TestSpectrum(unittest.TestCase):
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment