From 9c30fc44059d8d9112a71cc6dc6474c99096eab4 Mon Sep 17 00:00:00 2001 From: Theodor Chikin Date: Wed, 8 Oct 2025 19:07:22 +0300 Subject: [PATCH] moved fixed-point trigonometry to a separate file --- FFT.py | 145 ++++++++++++------- FP_trigonometry.py | 149 ++++++++++++++++++++ __pycache__/FFT.cpython-312.pyc | Bin 8236 -> 9603 bytes __pycache__/FP_trigonometry.cpython-312.pyc | Bin 0 -> 6093 bytes naive_DFT.py | 126 +---------------- 5 files changed, 249 insertions(+), 171 deletions(-) create mode 100755 FP_trigonometry.py create mode 100644 __pycache__/FP_trigonometry.cpython-312.pyc 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 396ab2e2758da4a32390e6b1ddea8fcbfef67499..4a1791bcf02f97624ef7231546eaa435afad0965 100644 GIT binary patch delta 2481 zcmb7GeQZ-z6u&Mz{z{GBytc)=f0*Hf+1;zkHpnRAKGUx*~7>KVWKzU8) zXG*A{%_VUpXlBA65HLc<9~cQy6C(-6tbrgY82N*!i!rH1{(xxooVV_y6HL6x`}N-Q zb z$3=yKel9>&DJm0U@ETtUx|+mxiIWtAYSQjHLi9`QPySOXv2#L$Dj8B#Bm~6_0R(>$ z=CZkB?UM}9Q3D%kC0i$K{QqY7I2U9WL=U?xDuvG`SXK@3?3y^uGQ@!=0XzFcTwjt^ z+S=A_*{bVX+E%mWQnRfZc%|@#S8@YPnn&qfCRsW~Pbzq{k@Z+!GuQRDoI(_P@AQZH zZpu{IGwor%j@agft3dPIly<-HWreLhAb%zL+Gh%1bH4W3!Z*#n_W7(??xN-FwCof8 zAi8qkJ2^|Khn-OxsF!`O1lg~um983?WBZM#0amTf2j;74&8jINt(PNAMi3Anr$pxw zDgXi{S$a^Pj>~$483;2G8UV6%YnE=y(omLax{rOUPE|~62pt0t)T@++2d+4;@bm$T zc^=I0pd#i;F9Y})J{&Cpt}x9tkI={?_K~+-G98@EKKIsI(T6PE1=|dHXVV5&RlZcg z!`Bg_13l%RaPdZ%?t~9)9%~6vb2bcC_`#s=Z-O8{AO=}Nj?oaelZ3d<;7g1~Mx_lJ zK^M(|6?V43NpnfBKw>D1IWmVa%GGF;i-RFf<6?ZoG9sYy26wh#ZVB1E7bFCO)A%Sj z6cY{J5TjsTcbGGXBnA}VdW*O;Im|~#c|$Px-C!|HlXV7m3LR|&_Myw|0-7+61}oO@ zs;nx|n!vF!3D^nNO@v)^`Qyob0?}v;`pV#TPa-5m^^Jyvaly8h2xib6nh-G?8yC1W z^r3AmOU1d1$r1ao7BusJ>e=FCf~s{~7CD16q_cQ{Fd&bAflG=e7Heoy@j6x#vyIs> zEF%XA`nN&f4tY=vJ7xg$1&3$eB}GJFnJGd}OpTR}sTh_K#ff_3IpegQF-Js2_i<)62p16P)jsUkYpwl$oLzDG>?l2k{ok^jL8F4Qx=_=AmcVl4HweKoE@$QjYFzg z0jlAiC=tgkM+yRz(Fiy&W~3;tVpC6+D-;;mu^Dn9IO}7BGvhZ|5SuCbnAqS%lRpR~ z`L_UN2hKVmC!+Gg{QpK!HZK4>z(wf^Y#)_+m7EV03zb1Y~_a&hp)RBU!0!^_H!nkUo54A$P)UT7tr5`XXdu9SmQp%{cQR#9I&- zBeVi!mCzP#drxF5JRw;nJUZ8xqC^|QG6cMa^>&2i2rCd)BD{$35aA>LJ8$2S}?`hutCp4gPA?%j4rElaICdN@y9mR;;}Mc3T5r^`<3=jYsz=ihR; zQtQ&zzSjP-6U#G>+H>BF+L#djMSdou&SI-8V{wJJJ*jBAvA^~7q>o?DxMp7jKbkJ- zU#$2nKDZ`0xS=z%rZdyrc_Y-FS=60qJ;ME73)|A11)a;yBBGS=y}7J+ac&{8Dv32o zDc3@Y&wo?)_qh}9iAv>PLVZzwi{}U8Ao?z@s@sbbQs!vw9>q?Rj0>-?M#@ zz8+}!E`_lWb3CQ=HgWw(YOa%-yev|C-mw563iO7&q|iC}ipj8Uq0dr#?jXf~8Pa`W Ko`u>?>;C{3H~|y@ delta 820 zcmZvaOGs2v7{||d?>w&eF*Y;lJXE}+V|qsi4a-NQ4y6+z43rh*qQSbFie{Y2d@z`h z3wI%BGeWCYNk+uo2Sz~weEeLVZx+@q^4vP{YdC?Rj zyxs(sao58$_CtIG!0x+}<&gI*x(UH(-o@jWJRNLR?DAUlBrgoH z&z^+KlldsH8gHXusBF-Co^5*j6rCse_{!Y@>$|s;x0gEl-V5GUaV7DlOkGYCogShz z799@>CX0?ll+L2l#F-26s?F+wUe{JwNsH1`cPuTrp8H~HjhXT|RI*oc%w5Nmj#|62 z4G?6n)g)A~U+R!t!DWi*vMy?$MA%@b{9Qa_*2{+a_`2icw^eek;6p&AEO5_G`Nq4 zn6t%d16_6Df`7feaHi0{Ub(I4%ZJ{|p#l_Qg*eyZ>)vfyx8?(rP(eF+jp4K1*Jh`$ zFHrvV2<>J&@tR&8U3{<7!|c64=^i+4gU}>~aQ;u}`1Fi9J2l!%GklBREZ5Ca4Phe= fnv(&5q~np{h!d^?Z8rr3RL6CgyBUECr9OWEg15Ex diff --git a/__pycache__/FP_trigonometry.cpython-312.pyc b/__pycache__/FP_trigonometry.cpython-312.pyc new file mode 100644 index 0000000000000000000000000000000000000000..7474dc9c6a07fc6532adf7b864fd53084b6615c7 GIT binary patch literal 6093 zcmcIoU2GKB6~6Pcv;VtZ|6w2@>mLkbj2DwYD~CoQ0rO*U0E3g5kkxo+yk5NPH8W!y zw!6&+C5kQFX01@gNIbYI4@Tx84|z!6cu1>AZKYjyo48XWAvOKs%{A#m=tFzX%#L?3 zAf-`zrMYu|&OPVcd(Qptx%anTuZut$xboNGpZW>;Ggiz*mkSSm28EA_NJJ`564U^d zpafJ?MVu9OY9D|$z^6CyikBGoLGm7&Nb`V!6i?qi*2HAabzYD^@l zX!KZ2W0ecC+RIa#HzABjVP$-9ES^-A7$G#dPr-NySz@ zrAeCL1kzNR!hZ@{7t3_!zmMFcdT-Ji6B|`!4hW!3%aTUN5(h<29%u(ohnl4YbjLO0=J*7agb-MD5 z?14VrQW^)VDMTW&7h1~3(lytYQqO9!rx`L+>-Fn-qrTwXR&;L9)7y=SJaJJ!1y6v9 z03RdC8!ApFQl^R&v9!`@vRHOk3fxd+MxZ7)S$McD)NF-^Xj9RcdJ=8MGMT;T#C^z( zQ_)rRgl!{l0q=cJ;cYiZ>w2sknc6@|#yT^PMXU95ST|OeO3my`lZrNa+g^(%n=Skt z6-_Fw*3p#}v#x@bA_Wd>axqt0W#8`QjM;(AQY_@1_ zwU5h`2{vA{PDL(Q$)Z}#m+h0AB5j@?pj|8(?@t+Rh|*P%Rn=${<1*Hji~*|~YHyXb7r)6I`VH2RVr1x$gP zS|M7b5DGGsitei?0N=H(E7ImBv2JJ)V4(s*7HSJT1Nho>D0izEf#^x(Om6ztHk*lP z*Wt+6+7;?G7`$xN&TWzBH2EQ9r>PNyrfdby8Q=(_FJiA9`X=@0wbb5O(Vdn5Lg;X#nWF;z!BCxO&b(e=+j56(OJ*$3}ey1!J@w}O&!B&#gn`i{cFHN>|#0Q z1beXnyS2UdQJK`Uk3)rt!|<0k>NGlxA-3ueJF3J)DR_2Qd+j^YU2DW}mH+TE3MqY7BvwWllPC|A0o4Lc8<{#T>( z<@>1keGDd24G|$88IKDXwIbSXzsLs#C3vxmcgl?@8{<;>In0_M(>%y**B%2rWel>6 z&sbPHn!Ll+m$CCmZemWH|3QAoizRL!;$zQT z{XBK8p~N*U2b$-*7dn4&xER=(J-X}*M+ZpmGp?fu&2otnru^cVd7r8mO`f4C$> z3%+PRF0T+W0Vweg0(Uo+{)UiAI%Cu_SbK*)$UtA5I>-R!VdetT*DKysYv z!`d09AN!6QzaHZ^LcNZS{?{3}QMlh=vakN7d<3=5Pb?@U?)iIeU-szZu#5IW&w0=e znKi_R$Yu?24G6YwBG?*qukOe^Z5)rfW3`}_5rfJX4-RgGc=UEcQv&iqfgty3BPnF>40Qr!L&{4R#oP*FQ12;rK4BFiF zB7>z_ddAO^44t8(B*VbR!e>j<2gtdPU<{XL&ZDGtse%pZbpf-9>gFGi+G!33TD*(t zNyj#vN2ONj1HH`g8G8ktzt?z&8Zf1ZY?@2k(~d=0HB5+PPnuoYdY!O6kz(z9qsYw9 zS)$e$906%*(RPb7Q9a{KJI^CntV=3JTzWd+V3sO>ICU7&E;?>GO-)4`aCKP_`1^5X z!6RH{kAOAKxY90*MGSB7vfwWnm+1aK$hZ9$@;(21`6?bSdP*~Jo|QK-=;Z#;69M0c zeFDDn7PM#F8Bf}s_B^^exXod|)W8Zc8BfM7!b6e8Mrip*ZG30*J z)>t_?sc38@8CNtXo>W8eq@bqYVM-o>>(P!lMEwdpU{Mz}>XOEe3JFPL6G>5mu?&pG zN3FPla0p_WO~0FH^u8n?2JPg14_*)h$6fw!Q+!gN0XYpFCwULKhl^~Y)tcvkBNdRg)a&hC76veoSfvn z8r7>&@Yn|86p$(@;USH~U-t>kg?|HlH3otPje{p4xYcMb{9^}{ZINL?R^?h0PYD_| zsNmgGKfTS&ct6sQdsR|EAGr@mF2JY!UO&v2I9IlDxoO)1^~>#FI(N*Soa&sqkPA#J z(?7@^xYCh#?#Q-4w&e0;x5BL?+q&crWZRa!fmLE>nzE54*7?D)caL4I$@N{Yp6&l= zbBTTSYi~{N)$82s(d*-d+Ljyrh1wlOZ+o_TIj~`lx!!%vbz@JVVf(H6LPMw+*ac@d zMhwr)_06{~w3XOBaEbch_`AoA`xxu{imij&RPOvA;PT=Jncq|qcf;Jt`I9Brw(Qn> zwEct7yP>O1AGgi6&DSo}eRB5pJEiK^^XJ~qi&8dJqKB3qCArOg$shWxW}*4k&|=#i z{?Hf0`E&77bt12*`P7ebM0))Q*Oz=xWkgT@)IeSsD^*{}PhQOb8U! z&9|j5`irle{)#{ItLH!K{OsgHw9la6uV@Y=xmH>TjgENSg1GEq4ij&*{r@OAhDM-f3sn zw(6j$rd5un8dg0t#jiL>RUqqGWqGP)` Op7+~V7{YGS_4p64h6wop literal 0 HcmV?d00001 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() - - +