238 lines
6.8 KiB
Python
Executable File
238 lines
6.8 KiB
Python
Executable File
#!/usr/bin/python3
|
|
from math import sin, cos, pi, sqrt
|
|
import plotly.graph_objs as go
|
|
from plotly.subplots import make_subplots
|
|
|
|
FP_acc = 1e3
|
|
|
|
INP_L = 1024
|
|
|
|
|
|
F_nyquist = INP_L//2
|
|
|
|
|
|
pi_FP = 1* FP_acc
|
|
|
|
|
|
def abs_f(re, im):
|
|
return sqrt(re*re + im*im)
|
|
|
|
|
|
def abs_FP(re, im):
|
|
return int(sqrt(re*re + im*im))
|
|
|
|
def sqrt_FP(val):
|
|
#print(val)
|
|
return int(sqrt(val))
|
|
|
|
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 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 = (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 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):
|
|
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 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.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.show()
|
|
|
|
|
|
|
|
|
|
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 = go.Figure()
|
|
chart.add_trace(go.Scatter(x = angs, y=res_f, name="sin_float", mode="markers+lines"))
|
|
chart.add_trace(go.Scatter(x = angs, y=[val/FP_acc for val in res_FP], name="sin_FP", mode="markers+lines"))
|
|
chart.add_trace(go.Scatter(x = angs, y=[val/FP_acc for val in res_cos_FP], name="cos_FP", mode="markers+lines"))
|
|
chart.show()
|
|
|
|
|
|
|
|
|
|
|
|
|
|
if __name__ == "__main__":
|
|
#main()
|
|
# DFT_tester()
|
|
FFT_tester()
|
|
#sin_tester()
|