moved python3 prototyping to the python folder
This commit is contained in:
326
python/FFT.py
Normal file
326
python/FFT.py
Normal file
@ -0,0 +1,326 @@
|
||||
import numpy as np
|
||||
from math import sin, cos, pi, sqrt
|
||||
|
||||
|
||||
|
||||
def FFT_backup(inp):
|
||||
inp = np.array(inp)
|
||||
out = np.fft.fft(inp)
|
||||
out = [val for val in np.abs(out)]
|
||||
print(out)
|
||||
return out
|
||||
|
||||
|
||||
def FFT_real(inp):
|
||||
|
||||
mode = 3
|
||||
if (mode == 0):
|
||||
inp = np.array(inp)
|
||||
out = FFT_np(inp)
|
||||
out = [np.abs(val) for val in out]
|
||||
if (mode == 1): # simple with recursion
|
||||
|
||||
out = FFT_arr([[val, 0] for val in inp])
|
||||
#out = [val for val in np.abs(out)]
|
||||
print("FFT_calculated!")
|
||||
print("output:", out)
|
||||
out = [sqrt(val[0]**2 + val[1]**2) for val in out]
|
||||
# out = [2 for val in out]
|
||||
if (mode == 2): #no internal buffs
|
||||
|
||||
tmp = []
|
||||
for re in inp:
|
||||
tmp.append(re)
|
||||
tmp.append(0)
|
||||
|
||||
out = FFT_arr_inplace(tmp)
|
||||
|
||||
tmp = out.copy()
|
||||
out = []
|
||||
for i in range(len(tmp)//2):
|
||||
out.append([tmp[2*i], tmp[2*i+1]])
|
||||
|
||||
|
||||
#out = [val for val in np.abs(out)]
|
||||
print("FFT_calculated!")
|
||||
#print("output:", out)
|
||||
out = [sqrt(val[0]**2 + val[1]**2) for val in out]
|
||||
|
||||
if (mode == 3): #fixed-point
|
||||
|
||||
tmp = []
|
||||
for re in inp:
|
||||
tmp.append(re)
|
||||
tmp.append(0)
|
||||
|
||||
out = FFT_arr_FP(tmp)
|
||||
|
||||
tmp = out.copy()
|
||||
out = []
|
||||
for i in range(len(tmp)//2):
|
||||
out.append([tmp[2*i], tmp[2*i+1]])
|
||||
|
||||
|
||||
#out = [val for val in np.abs(out)]
|
||||
print("FFT_calculated!")
|
||||
#print("output:", out)
|
||||
out = [sqrt(val[0]**2 + val[1]**2) for val in out]
|
||||
|
||||
|
||||
|
||||
#print(out)
|
||||
return out
|
||||
|
||||
|
||||
def FFT(inp):
|
||||
return FFT_arr([[val, 0] for val in inp])
|
||||
|
||||
def FFT_np(inp):
|
||||
print("inp:", inp)
|
||||
inp = np.array(inp)
|
||||
N = inp.shape[0]
|
||||
if N == 1:
|
||||
return inp
|
||||
X_even = FFT_np(inp[::2])
|
||||
X_odd = FFT_np(inp[1::2])
|
||||
k = np.arange(N // 2)
|
||||
tw = np.exp(-2j * np.pi * k / N) * X_odd
|
||||
return np.concatenate((X_even + tw, X_even - tw))
|
||||
|
||||
def FFT_arr(x):
|
||||
print("x:", x)
|
||||
|
||||
N = len(x)
|
||||
print("len(x):", N)
|
||||
if N == 1:
|
||||
print("As len(x) == 1, return it`s value")
|
||||
return x
|
||||
X_even = FFT_arr(x[::2])
|
||||
X_odd = FFT_arr(x[1::2])
|
||||
print("X_even:",X_even)
|
||||
|
||||
tw = []
|
||||
for k in range(len(X_odd)):
|
||||
a = cos(2*pi * k/N) #real
|
||||
b = -sin(2*pi * k/N) #imag
|
||||
|
||||
c = X_odd[k][0] # real
|
||||
d = X_odd[k][1] # imag
|
||||
|
||||
print("a,b,c,d:", a,b,c,d)
|
||||
|
||||
tw.append([a*c - b*d, b*c + a*d]) #(a + ib)(c + id) = (ac - bd) + i(bc + ad)
|
||||
res = [0 for i in range(len(X_even)*2)]
|
||||
for i in range(len(X_even)):
|
||||
res[i] = [X_even[i][0] + tw[i][0], X_even[i][1] + tw[i][1]]
|
||||
res[i+len(X_even)] = [X_even[i][0] - tw[i][0], X_even[i][1] - tw[i][1]]
|
||||
return res
|
||||
|
||||
|
||||
|
||||
|
||||
def FFT_arr_inplace(buf):
|
||||
"""
|
||||
In-place radix-2 DIT FFT для списка buf длины N=2^m, где каждый элемент — [re, im].
|
||||
Без массивов twiddle: твиддл на уровне обновляется рекуррентно.
|
||||
"""
|
||||
|
||||
N = len(buf)//2
|
||||
# --- bit-reverse перестановка (чтобы бабочки шли последовательно) ---
|
||||
j = 0
|
||||
for i in range(1, N):
|
||||
bit = N >> 1
|
||||
while j & bit:
|
||||
j ^= bit
|
||||
bit >>= 1
|
||||
j |= bit
|
||||
if i < j:
|
||||
buf[i*2], buf[i*2+1], buf[j*2], buf[j*2+1] = buf[j*2], buf[j*2+1], buf[i*2], buf[i*2+1]
|
||||
|
||||
# --- уровни бабочек ---
|
||||
m = 2
|
||||
while m <= N:
|
||||
half = m // 2
|
||||
# шаг угла Δ = 2π/m, базовый твиддл w_m = e^{-jΔ} => (cw, sw)
|
||||
c_delta = cos(2.0 * math.pi / m)
|
||||
s_delta = -sin(2.0 * math.pi / m)
|
||||
|
||||
for start in range(0, N, m):
|
||||
# w = e^{-j*0*Δ} = 1 + j0
|
||||
wr, wi = 1.0, 0.0
|
||||
for k in range(half):
|
||||
u_re, u_im = buf[(start + k)*2], buf[(start + k)*2 +1]
|
||||
v_re, v_im = buf[(start + k + half)*2], buf[(start + k + half)*2 +1]
|
||||
|
||||
# t = w * v
|
||||
print("wr, v_re, wi, v_im:",wr, v_re, wi, v_im)
|
||||
t_re = wr * v_re - wi * v_im
|
||||
t_im = wr * v_im + wi * v_re
|
||||
|
||||
# верх/низ
|
||||
buf[(start + k)*2 +0] = u_re + t_re
|
||||
buf[(start + k)*2 +1] = u_im + t_im
|
||||
buf[(start + k + half)*2 + 0] = u_re - t_re
|
||||
buf[(start + k + half)*2 + 1] = u_im - t_im
|
||||
|
||||
# w *= w_m (поворот на Δ с помощью рекуррентной формулы)
|
||||
# (wr + j wi) * (cΔ + j sΔ)
|
||||
wr, wi = wr * c_delta - wi * s_delta, wr * s_delta + wi * c_delta
|
||||
|
||||
m <<= 1
|
||||
return buf
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
FP_acc = 2<<16
|
||||
|
||||
def FFT_arr_FP(buf):
|
||||
"""
|
||||
In-place radix-2 DIT FFT для списка buf длины N=2^m, где каждый элемент — [re, im].
|
||||
Без массивов twiddle: твиддл на уровне обновляется рекуррентно.
|
||||
"""
|
||||
|
||||
buf = [int(val*FP_acc) for val in buf]
|
||||
|
||||
N = len(buf)//2
|
||||
# --- bit-reverse перестановка (чтобы бабочки шли последовательно) ---
|
||||
j = 0
|
||||
for i in range(1, N):
|
||||
bit = N >> 1
|
||||
while j & bit:
|
||||
j ^= bit
|
||||
bit >>= 1
|
||||
j |= bit
|
||||
if i < j:
|
||||
buf[i*2], buf[i*2+1], buf[j*2], buf[j*2+1] = buf[j*2], buf[j*2+1], buf[i*2], buf[i*2+1]
|
||||
|
||||
# --- уровни бабочек ---
|
||||
|
||||
c_delta_FP_arr = [cos(2.0 * pi / m) * FP_acc for m in range(2,N+2)]
|
||||
s_delta_FP_arr = [-sin(2.0 * pi / m) * FP_acc for m in range(2,N+2)]
|
||||
m = 2
|
||||
while m <= N:
|
||||
half = m // 2
|
||||
# шаг угла Δ = 2π/m, базовый твиддл w_m = e^{-jΔ} => (cw, sw)
|
||||
print("N, m:", N,m)
|
||||
|
||||
c_delta_FP = c_delta_FP_arr[m -2]
|
||||
s_delta_FP = s_delta_FP_arr[m -2]
|
||||
#c_delta_FP = cos(2.0 * pi / m) * FP_acc
|
||||
#s_delta_FP = -sin(2.0 * pi / m) * FP_acc
|
||||
|
||||
for start in range(0, N, m):
|
||||
# w = e^{-j*0*Δ} = 1 + j0
|
||||
wr, wi = 1.0 * FP_acc, 0.0 * FP_acc
|
||||
for k in range(half):
|
||||
u_re, u_im = buf[(start + k)*2], buf[(start + k)*2 +1]
|
||||
v_re, v_im = buf[(start + k + half)*2], buf[(start + k + half)*2 +1]
|
||||
|
||||
# t = w * v
|
||||
print("wr, v_re, wi, v_im:", wr, v_re, wi, v_im)
|
||||
t_re = (wr * v_re - wi * v_im) / FP_acc
|
||||
t_im = (wr * v_im + wi * v_re) / FP_acc
|
||||
|
||||
# верх/низ
|
||||
buf[(start + k)*2 +0] = u_re + t_re
|
||||
buf[(start + k)*2 +1] = u_im + t_im
|
||||
buf[(start + k + half)*2 + 0] = u_re - t_re
|
||||
buf[(start + k + half)*2 + 1] = u_im - t_im
|
||||
|
||||
|
||||
# w *= w_m (поворот на Δ с помощью рекуррентной формулы)
|
||||
# (wr + j wi) * (cΔ + j sΔ)
|
||||
wr, wi = (wr * c_delta_FP - wi * s_delta_FP)/ FP_acc, (wr * s_delta_FP + wi * c_delta_FP)/ FP_acc
|
||||
|
||||
|
||||
m <<= 1
|
||||
|
||||
buf = [val/ FP_acc for val in buf]
|
||||
return buf
|
||||
|
||||
#global arrays. calculated once at compile time. They store sin and cos values as fixed-point
|
||||
c_delta_FP_arr = [int(cos(2.0 * pi / m) * FP_acc) for m in range(2,inp_L//2 + 2)]
|
||||
s_delta_FP_arr = [int(-sin(2.0 * pi / m) * FP_acc) for m in range(2, inp_L//2 + 2)]
|
||||
|
||||
def FFT_arr_FP_for_C(inp, inp_L, buf): #version for translation to C directly
|
||||
"""
|
||||
In-place radix-2 DIT FFT для списка buf длины N=2^m, где каждый элемент — [re, im].
|
||||
Без массивов twiddle: твиддл на уровне обновляется рекуррентно.
|
||||
"""
|
||||
|
||||
#buf = [int(val*FP_acc) for val in buf]
|
||||
for i in range(inp_L):
|
||||
buf[i*2] = inp[i] #take data from input as real and store it to the buf as Re, Im pairs.
|
||||
buf[i*2 + 1] = 0
|
||||
|
||||
N = inp_L//2
|
||||
# --- bit-reverse перестановка (чтобы бабочки шли последовательно) ---
|
||||
j = 0
|
||||
for i in range(1, N):
|
||||
bit = N >> 1
|
||||
while j & bit:
|
||||
j ^= bit
|
||||
bit >>= 1
|
||||
j |= bit
|
||||
if i < j:
|
||||
buf[i*2], buf[i*2+1], buf[j*2], buf[j*2+1] = buf[j*2], buf[j*2+1], buf[i*2], buf[i*2+1]
|
||||
|
||||
# --- уровни бабочек ---
|
||||
|
||||
c_delta_FP_arr = [cos(2.0 * pi / m) * FP_acc for m in range(2,N+2)]
|
||||
s_delta_FP_arr = [-sin(2.0 * pi / m) * FP_acc for m in range(2,N+2)]
|
||||
m = 2
|
||||
while m <= N:
|
||||
half = m // 2
|
||||
# шаг угла Δ = 2π/m, базовый твиддл w_m = e^{-jΔ} => (cw, sw)
|
||||
#print("N, m:", N,m)
|
||||
|
||||
c_delta_FP = c_delta_FP_arr[m -2]
|
||||
s_delta_FP = s_delta_FP_arr[m -2]
|
||||
#c_delta_FP = cos(2.0 * pi / m) * FP_acc
|
||||
#s_delta_FP = -sin(2.0 * pi / m) * FP_acc
|
||||
|
||||
for start in range(0, N, m):
|
||||
# w = e^{-j*0*Δ} = 1 + j0
|
||||
wr, wi = 1.0 * FP_acc, 0.0 * FP_acc
|
||||
for k in range(half):
|
||||
u_re, u_im = buf[(start + k)*2], buf[(start + k)*2 +1]
|
||||
v_re, v_im = buf[(start + k + half)*2], buf[(start + k + half)*2 +1]
|
||||
|
||||
# t = w * v
|
||||
#print("wr, v_re, wi, v_im:", wr, v_re, wi, v_im)
|
||||
t_re = (wr * v_re - wi * v_im) / FP_acc
|
||||
t_im = (wr * v_im + wi * v_re) / FP_acc
|
||||
|
||||
# верх/низ
|
||||
buf[(start + k)*2 +0] = u_re + t_re
|
||||
buf[(start + k)*2 +1] = u_im + t_im
|
||||
buf[(start + k + half)*2 + 0] = u_re - t_re
|
||||
buf[(start + k + half)*2 + 1] = u_im - t_im
|
||||
|
||||
|
||||
# w *= w_m (поворот на Δ с помощью рекуррентной формулы)
|
||||
# (wr + j wi) * (cΔ + j sΔ)
|
||||
wr, wi = (wr * c_delta_FP - wi * s_delta_FP)/ FP_acc, (wr * s_delta_FP + wi * c_delta_FP)/ FP_acc
|
||||
|
||||
|
||||
m <<= 1
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
149
python/FP_trigonometry.py
Executable file
149
python/FP_trigonometry.py
Executable file
@ -0,0 +1,149 @@
|
||||
#!/usr/bin/python3
|
||||
from math import sin, cos, pi, sqrt
|
||||
import plotly.graph_objs as go
|
||||
from plotly.subplots import make_subplots
|
||||
|
||||
FP_acc = 2<<16
|
||||
pi_FP = 1* FP_acc
|
||||
|
||||
|
||||
|
||||
def abs_FP(re, im):
|
||||
return int(sqrt(re*re + im*im))
|
||||
|
||||
def sqrt_FP(val):
|
||||
#print(val)
|
||||
return int(sqrt(val))
|
||||
|
||||
|
||||
|
||||
def abs_FP(re, im):
|
||||
# return sqrt(re*re + im*im)
|
||||
return int(sqrt(re*re + im*im)/FP_acc)
|
||||
|
||||
trigon_debug = 0
|
||||
|
||||
def sin_FP(phi_fp):
|
||||
if (trigon_debug):
|
||||
print("sin_FP========")
|
||||
print("phi:", phi_fp)
|
||||
if phi_fp < 0:
|
||||
if (trigon_debug):
|
||||
print("phi < 0. recursive inversion...")
|
||||
return -1 *sin_FP(-1*phi_fp)
|
||||
|
||||
while phi_fp >= 2*pi_FP:
|
||||
if (trigon_debug):
|
||||
print("phi is bigger than 2Pi. Decreasing...")
|
||||
phi_fp -= 2*pi_FP
|
||||
if (trigon_debug):
|
||||
print("phi:", phi_fp)
|
||||
|
||||
if phi_fp >= pi_FP:
|
||||
if (trigon_debug):
|
||||
print("phi > pi_FP. recursive inversion...")
|
||||
print(phi_fp, pi_FP)
|
||||
return -1*sin_FP(phi_fp - pi_FP)
|
||||
|
||||
if phi_fp == pi_FP/2:
|
||||
return 1*FP_acc
|
||||
if phi_fp == 0:
|
||||
return 0
|
||||
|
||||
if phi_fp > pi_FP/2:
|
||||
if (trigon_debug):
|
||||
print("phi > pi_FP/2. recursive inversion...")
|
||||
return sin_FP(pi_FP - phi_fp)
|
||||
|
||||
#now phi should be inside [0, Pi/2). checking...
|
||||
if phi_fp < 0:
|
||||
raise ValueError('error in sin_FP. after all checks phi < 0')
|
||||
if phi_fp >= pi_FP/2:
|
||||
raise ValueError('error in sin_FP. after all checks phi > pi_FP/2')
|
||||
|
||||
#now phi is inside [0, Pi/2). So, cos(phi) > 1 always
|
||||
if (trigon_debug):
|
||||
print("phi:", phi_fp)
|
||||
return sin_FP_constrained(phi_fp)
|
||||
|
||||
|
||||
sin_05_debug = 0
|
||||
|
||||
def sin_FP_constrained(phi_fp):
|
||||
phi_trh = pi_FP/16
|
||||
if (trigon_debug):
|
||||
print("sin_FP_constrained===========")
|
||||
print("phi:", phi_fp)
|
||||
print("check is phi inside [0, Pi/2)")
|
||||
if phi_fp < 0:
|
||||
raise ValueError('error in sin_FP. after all checks phi < 0')
|
||||
if phi_fp >= pi_FP/2:
|
||||
raise ValueError('error in sin_FP. after all checks phi > pi_FP/2')
|
||||
if (trigon_debug):
|
||||
print("Ok")
|
||||
if (phi_fp > phi_trh):
|
||||
if (sin_05_debug):
|
||||
print("phi_fp:", phi_fp," >",phi_trh,"... recusion")
|
||||
print("phi_fp/2:", phi_fp/2)
|
||||
sin_phi_05 = sin_FP_constrained(phi_fp/2)
|
||||
sin_phi = int((2*sin_phi_05 * sqrt_FP(FP_acc**2 - sin_phi_05*sin_phi_05))/FP_acc)
|
||||
if (sin_05_debug):
|
||||
print("sin_phi_05:", sin_phi_05)
|
||||
print("sin_phi:",sin_phi)
|
||||
return sin_phi
|
||||
else:
|
||||
res = int((phi_fp * 3.141592653589793238462643383279502884197169399375105820974944 * FP_acc)/FP_acc - ((phi_fp * 3.141592653589793238462643383279502884197169399375105820974944 * FP_acc)/FP_acc)**3/(6*FP_acc**2))
|
||||
# res = int((phi_fp * 1.0 * FP_acc)/FP_acc - ((phi_fp * 1.0 * FP_acc)/FP_acc)**3/(6*FP_acc**2))
|
||||
if (trigon_debug):
|
||||
print("calculating sin(x) as x:",phi_fp, res)
|
||||
|
||||
return res
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
# return sin(pi*(phi_fp/FP_acc)/(pi_FP/FP_acc))
|
||||
|
||||
|
||||
def cos_FP(phi_fp):
|
||||
return sin_FP(phi_fp - pi_FP/2)
|
||||
|
||||
|
||||
def sin_tester():
|
||||
N = 4000
|
||||
angs = [(i - N/2)/1000 for i in range(N)]
|
||||
res_f = []
|
||||
res_FP = []
|
||||
res_cos_FP = []
|
||||
for phi in angs:
|
||||
res_f.append(sin(phi*pi))
|
||||
# print(phi, phi*FP_acc*pi_FP/FP_acc)
|
||||
val_fp = sin_FP(phi*FP_acc*pi_FP/FP_acc)
|
||||
val_cos_fp = cos_FP(phi*FP_acc*pi_FP/FP_acc)
|
||||
print("angle, sin, cos:",phi, val_fp, val_cos_fp)
|
||||
res_FP.append(val_fp)
|
||||
res_cos_FP.append(val_cos_fp)
|
||||
chart = make_subplots(rows=2, cols=1)
|
||||
|
||||
chart.add_trace(go.Scatter(x = angs, y=res_f, name="sin_float", mode="markers+lines"), row=1, col=1)
|
||||
chart.add_trace(go.Scatter(x = angs, y=[val/FP_acc for val in res_FP], name="sin_FP", mode="markers+lines"), row=1, col=1)
|
||||
chart.add_trace(go.Scatter(x = angs, y=[val/FP_acc for val in res_cos_FP], name="cos_FP", mode="markers+lines"), row=1, col=1)
|
||||
chart.add_trace(go.Scatter(x = angs, y=[a - b/FP_acc for a,b in zip(res_f, res_FP)], name="error", mode="markers+lines"), row=2, col=1)
|
||||
chart.update_xaxes(matches="x1", row=2, col=1)
|
||||
|
||||
chart.show()
|
||||
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
sin_tester()
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
124
python/naive_DFT.py
Executable file
124
python/naive_DFT.py
Executable file
@ -0,0 +1,124 @@
|
||||
#!/usr/bin/python3
|
||||
from math import sin, cos, pi, sqrt
|
||||
import plotly.graph_objs as go
|
||||
from plotly.subplots import make_subplots
|
||||
|
||||
from FFT import FFT_real
|
||||
from FP_trigonometry import *
|
||||
|
||||
|
||||
INP_L = 1024
|
||||
|
||||
|
||||
F_nyquist = INP_L//2
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
def abs_f(re, im):
|
||||
return sqrt(re*re + im*im)
|
||||
|
||||
|
||||
|
||||
|
||||
def DFT_naive(inp, out):
|
||||
for f in range(len(out)):
|
||||
val_re = 0
|
||||
val_im = 0
|
||||
for n in range(len(inp)):
|
||||
phi = 2*pi*f*n/INP_L
|
||||
val_re += inp[n] * sin(phi) /INP_L
|
||||
val_im += inp[n] * cos(phi) /INP_L
|
||||
|
||||
val_abs = abs_f(val_re, val_im)
|
||||
#print("F, val_abs:",f, val_abs)
|
||||
out[f] = val_abs
|
||||
|
||||
|
||||
|
||||
|
||||
def DFT_naive_FP(inp_float, out):
|
||||
inp = [val*FP_acc for val in inp_float]
|
||||
for f in range(len(out)):
|
||||
val_re = 0
|
||||
val_im = 0
|
||||
for n in range(len(inp)):
|
||||
phi = 2*pi_FP*f*n/INP_L
|
||||
phi_sin = sin_FP(phi)
|
||||
phi_cos = cos_FP(phi)
|
||||
#print(phi, phi_sin, phi_cos)
|
||||
val_re += inp[n] * phi_sin /INP_L
|
||||
val_im += inp[n] * phi_cos /INP_L
|
||||
|
||||
val_abs = abs_FP(val_re, val_im)
|
||||
#print("F, val_abs:",f, val_abs)
|
||||
out[f] = val_abs
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
def FFT_naive(inp, out):
|
||||
fft_out = FFT_real(inp)
|
||||
for i in range(len(fft_out)):
|
||||
val = fft_out[i]/len(inp)
|
||||
if (i < len(out)):
|
||||
out[i] = val
|
||||
else:
|
||||
pass
|
||||
#out.append(val)
|
||||
|
||||
def FFT_tester():
|
||||
inp = [-1 + 0.01*i + sin(2*pi*i/10) + cos(2*pi*i/20) + sin(2*pi*i/250) + sin(2*pi*i/2.001) for i in range(INP_L)]
|
||||
# inp = [sin(2*pi*i/2.001)for i in range(INP_L)]
|
||||
out_DFT = [0 for i in range(F_nyquist + 1)]
|
||||
out_FFT = [0 for val in range(F_nyquist + 1)]
|
||||
DFT_naive(inp, out_DFT)
|
||||
FFT_naive(inp, out_FFT)
|
||||
Fourier_error = []
|
||||
for a,b in zip(out_FFT, out_DFT):
|
||||
Fourier_error.append(a - b)
|
||||
chart = make_subplots(rows=3, cols=1)
|
||||
chart.add_trace(go.Scatter(x=[i for i in range(len(inp))], y=inp, name="inp", mode="markers+lines"), row=1, col=1)
|
||||
chart.add_trace(go.Scatter(x=[i for i in range(len(out_DFT))], y=out_DFT, name="out_DFT", mode="markers+lines"), row=2, col=1)
|
||||
chart.add_trace(go.Scatter(x=[i for i in range(len(out_FFT))], y=out_FFT, name="out_FFT", mode="markers+lines"), row=2, col=1)
|
||||
chart.add_trace(go.Scatter(x=[i for i in range(len(Fourier_error))], y=Fourier_error, name="error", mode="markers+lines"), row=3, col=1)
|
||||
chart.update_xaxes(matches="x2", row=3, col=1)
|
||||
chart.show()
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
def DFT_tester():
|
||||
inp = [-1 + 0.01*i + sin(2*pi*i/10) + cos(2*pi*i/20) + sin(2*pi*i/250) + sin(2*pi*i/2.001) for i in range(INP_L)]
|
||||
# inp = [sin(2*pi*i/2.001)for i in range(INP_L)]
|
||||
out_float = [0 for i in range(F_nyquist + 1)]
|
||||
out_FP = [0 for val in out_float]
|
||||
DFT_naive(inp, out_float)
|
||||
DFT_naive_FP(inp, out_FP)
|
||||
FP_error = []
|
||||
for a,b in zip(out_float, out_FP):
|
||||
FP_error.append(a - b/FP_acc)
|
||||
chart = make_subplots(rows=3, cols=1)
|
||||
chart.add_trace(go.Scatter(x=[i for i in range(len(inp))], y=inp, name="inp", mode="markers+lines"), row=1, col=1)
|
||||
chart.add_trace(go.Scatter(x=[i for i in range(len(out_float))], y=out_float, name="out_float", mode="markers+lines"), row=2, col=1)
|
||||
chart.add_trace(go.Scatter(x=[i for i in range(len(out_FP))], y=[val/FP_acc for val in out_FP], name="out_FP", mode="markers+lines"), row=2, col=1)
|
||||
chart.add_trace(go.Scatter(x=[i for i in range(len(out_FP))], y=FP_error, name="FP_error", mode="markers+lines"), row=3, col=1)
|
||||
chart.update_xaxes(matches="x2", row=3, col=1)
|
||||
chart.show()
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
#main()
|
||||
# DFT_tester()
|
||||
FFT_tester()
|
||||
#sin_tester()
|
||||
44
python/plotter.py
Executable file
44
python/plotter.py
Executable file
@ -0,0 +1,44 @@
|
||||
#!/usr/bin/pypy3
|
||||
import plotly.graph_objs as go
|
||||
from decimal import *
|
||||
from sys import argv
|
||||
|
||||
|
||||
def main():
|
||||
chart = go.Figure()
|
||||
chart.add_trace(go.Scatter(x=[1,2,3], y=[1,2,3]))
|
||||
chart.show()
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
if len(argv) == 1:
|
||||
main()
|
||||
else:
|
||||
chart = go.Figure()
|
||||
filenames = argv[1:]
|
||||
all_data = {}
|
||||
for filename in filenames:
|
||||
f = open(filename, "rt")
|
||||
data = {"X":[], "Y":[]}
|
||||
all_data[filename] = data
|
||||
for line in f:
|
||||
try:
|
||||
line_splt = line.split(",")
|
||||
print(line_splt)
|
||||
x = line_splt[0]
|
||||
y = line_splt[1]
|
||||
#print(x,y)
|
||||
x = float(x)
|
||||
y = float(y)
|
||||
data["X"].append(x)
|
||||
data["Y"].append(y)
|
||||
except ValueError:
|
||||
pass
|
||||
except IndexError:
|
||||
pass
|
||||
f.close()
|
||||
print("data samples:",len(data["X"]))
|
||||
chart.add_trace(go.Scatter(x=data["X"], y=data["Y"], name=filename))
|
||||
# chart.update_layout(title=argv[1])
|
||||
chart.show()
|
||||
|
||||
Reference in New Issue
Block a user