diff --git a/FFT.py b/FFT.py index 59c7884..c2cef32 100644 --- a/FFT.py +++ b/FFT.py @@ -13,12 +13,12 @@ def FFT_backup(inp): def FFT_real(inp): - mode = 2 + mode = 3 if (mode == 0): inp = np.array(inp) out = FFT_np(inp) out = [np.abs(val) for val in out] - if (mode == 1): + if (mode == 1): # simple with recursion out = FFT_arr([[val, 0] for val in inp]) #out = [val for val in np.abs(out)] @@ -26,7 +26,7 @@ def FFT_real(inp): print("output:", out) out = [sqrt(val[0]**2 + val[1]**2) for val in out] # out = [2 for val in out] - if (mode == 2): + if (mode == 2): #no internal buffs tmp = [] for re in inp: @@ -43,10 +43,32 @@ def FFT_real(inp): #out = [val for val in np.abs(out)] print("FFT_calculated!") - print("output:", out) + #print("output:", out) out = [sqrt(val[0]**2 + val[1]**2) for val in out] - print(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 @@ -97,8 +119,6 @@ def FFT_arr(x): -import math - def FFT_arr_inplace(buf): """ In-place radix-2 DIT FFT для списка buf длины N=2^m, где каждый элемент — [re, im]. @@ -159,51 +179,74 @@ def FFT_arr_inplace(buf): +FP_acc = 1e3 - - - - - - - - - - -def FFT_arr_for_C(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_for_C(x[::2]) - X_odd = FFT_arr_for_C(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_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] + + # --- уровни бабочек --- + m = 2 + while m <= N: + half = m // 2 + # шаг угла Δ = 2π/m, базовый твиддл w_m = e^{-jΔ} => (cw, sw) + 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 - - - - - - \ No newline at end of file + buf = [val/ FP_acc for val in buf] + return buf + + + + + + + + + + + + + diff --git a/FP_trigonometry.py b/FP_trigonometry.py new file mode 100755 index 0000000..9476a9d --- /dev/null +++ b/FP_trigonometry.py @@ -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() + + + + + + + + + \ No newline at end of file diff --git a/__pycache__/FFT.cpython-312.pyc b/__pycache__/FFT.cpython-312.pyc index 396ab2e..4a1791b 100644 Binary files a/__pycache__/FFT.cpython-312.pyc and b/__pycache__/FFT.cpython-312.pyc differ diff --git a/__pycache__/FP_trigonometry.cpython-312.pyc b/__pycache__/FP_trigonometry.cpython-312.pyc new file mode 100644 index 0000000..7474dc9 Binary files /dev/null and b/__pycache__/FP_trigonometry.cpython-312.pyc differ diff --git a/naive_DFT.py b/naive_DFT.py index 6b3b917..a7b119f 100755 --- a/naive_DFT.py +++ b/naive_DFT.py @@ -4,8 +4,8 @@ import plotly.graph_objs as go from plotly.subplots import make_subplots from FFT import FFT_real +from FP_trigonometry import * -FP_acc = 1e3 INP_L = 1024 @@ -13,19 +13,14 @@ 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)): @@ -41,97 +36,6 @@ def DFT_naive(inp, out): 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): @@ -180,6 +84,7 @@ def FFT_tester(): 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() @@ -201,32 +106,13 @@ def DFT_tester(): 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() -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() - - +