From 59f99cc466b41bb3fc037832239b1a3248f70cfd Mon Sep 17 00:00:00 2001 From: Theodor Chikin Date: Wed, 8 Oct 2025 18:01:08 +0300 Subject: [PATCH] implemented in-place FFT c-style. (by ChatGPT) --- FFT.py | 171 +++++++++++++++++++++++++++++++- __pycache__/FFT.cpython-312.pyc | Bin 1712 -> 7502 bytes naive_DFT.py | 4 +- 3 files changed, 168 insertions(+), 7 deletions(-) diff --git a/FFT.py b/FFT.py index 25e6b14..4620b9a 100644 --- a/FFT.py +++ b/FFT.py @@ -1,5 +1,5 @@ import numpy as np - +from math import sin, cos, pi, sqrt @@ -12,17 +12,38 @@ def FFT_backup(inp): def FFT_real(inp): - inp = np.array(inp) - out = FFT(inp) - out = [val for val in np.abs(out)] + + mode = 2 + if (mode == 0): + inp = np.array(inp) + out = FFT_np(inp) + out = [np.abs(val) for val in out] + if (mode == 1): + + 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): + + + out = FFT_arr_inplace([[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] + print(out) return out def FFT(inp): - return FFT_np(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: @@ -33,3 +54,143 @@ def FFT_np(inp): 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 + + + + +import math + +def FFT_arr_inplace(buf): + """ + In-place radix-2 DIT FFT для списка buf длины N=2^m, где каждый элемент — [re, im]. + Без массивов twiddle: твиддл на уровне обновляется рекуррентно. + """ + N = len(buf) + # --- 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], buf[j] = buf[j], buf[i] + + # --- уровни бабочек --- + m = 2 + while m <= N: + half = m // 2 + # шаг угла Δ = 2π/m, базовый твиддл w_m = e^{-jΔ} => (cw, sw) + cΔ = math.cos(2.0 * math.pi / m) + sΔ = -math.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] + v_re, v_im = buf[start + k + half] + + # t = w * v + t_re = wr * v_re - wi * v_im + t_im = wr * v_im + wi * v_re + + # верх/низ + buf[start + k][0] = u_re + t_re + buf[start + k][1] = u_im + t_im + buf[start + k + half][0] = u_re - t_re + buf[start + k + half][1] = u_im - t_im + + # w *= w_m (поворот на Δ с помощью рекуррентной формулы) + # (wr + j wi) * (cΔ + j sΔ) + wr, wi = wr * cΔ - wi * sΔ, wr * sΔ + wi * cΔ + + m <<= 1 + return buf + + + + + + + + + + + + + + + + + + + + + +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 + + + + + + + + + \ No newline at end of file diff --git a/__pycache__/FFT.cpython-312.pyc b/__pycache__/FFT.cpython-312.pyc index 7ed90aa1e765975b6c3d8ef9833776a5aab6a668..c7e41bbd0320e7986aa2b19f56946fb11398e29b 100644 GIT binary patch literal 7502 zcmeHMYiv}<6`uRt`|xXR1GXV#*Y7n3WK0}Fk|v8C3M4^sOd{%(xYv8vc)`!RYXk0G zbq$Hw#g2S8LCQKo%JR>Y{(v-TY0{*S_D7{k(Opc?>efcy9K(#BDfdpm$W5(kYcb-iVO4$95mzcJKL{ zwcRJ&3)~cwy{QdwQ;dz{S~tw(ZW!*f67G`@GawB^r1_EUW|56A@Bu1DrH^7Q5LkT2 zc|M>IGVG8PmCbRdkgb=~>C4&Ofu7hU-b5rW%+1GZRKkl&O~9$a8?*d20sT4dAjRgq z>{M3v;yjoC6WNzf=x20pXZU>bcEp^>S=6_-Jp94W@eymd`W%DBL@`!+pEBSZ8&E=i z)ji<=Qv)m6?GZGrlNvT-i^o($BO1O#%Yk3OUbO-|j@8E^T%fB5jCBBv&ckZ>x-52Z zSoLc~m<}I&5MvQouS)~iMV5dWX?|ru!xxCDg7qL&1lbcm7;6Me z+Zt(KkV?P?w?x|0)pfU}`b6vWj+wStn^E)Hqo!0JX@@F1?XKjUV{Vr=CmEwX;e379 z&C=$12{_zkMN`^%i`}%qicyczb&IV`?|A68)Hc^Tz0+U}Z=6q1aWT=Jq-OUfy|epM zC3T<4^`_K@xwqZP2=xVzM;49+hpq zza${};rmMx+W1hW9ee=J>ld3k%Getfaz};m3Fq#3fc>OR?x+@iU*qm*ScxwVpeFdT z2O4%yqGLOtve9Qv*A2?=QWN3!=Ia@eVvXQrTAt;^eGxg|CB%U=0e! zCK-u6RxS)D&Q8EC@*G(kJ+`t{Mdd=-JtK4s+s&#n0bW37Nli_F&F1*(4DJhaHtyIE zR~4PsIgF#^3DWlEainfIvg=^Y9t;pblF6;U53Ef5_qv_YAA5MIE>3&e-FSzB%pUakFB_yx3~et@l7$ zv_jeOkr8Fz1od}OcAH|`LN}F$z$*I^<5PRu(j+D=N)dEKE4x_CKz0^ow#WBFhryDfMCoj{~>=puejfaG{4>Ep9~ zKzB56R1}&%#iBLOCQqZ19WS5IjARSDA=$`=-Evs^WI03>&M*KtSr2>mOsk!l3HLhxKRfnru zs4Pmcs9uZmSrmX*z+yCiFo@X`@XA&s7$QrN>JbJY)(d)9<1!0USZN0^AOUnh3eG77_(Gu@6Mvbm`+td?)VUp8u8(yOWx;o+w^HLRH z7~xTWi>|!m5CNRzlJjN83vW~bz%3z>sgEcNf*kcu_r^Mmme@X{Gj;&79Pffy5bd1S zVuuWG?69H4x~CpAg^Gyxj!+m?AUATx9Hmll!Bi zKpQkVjnQ~l>=YiGtw~Zz_eEybV;0qZ#?~#Cff!$uqjgmmY7%{KZ%Hb%Tj$FgZwsF2 z=sBk;cyMy%$q4@98yn~42bK`cJVQm4bWu6T8%R*GfyDC=!o84enLV5;s!NW~DHl&& zabDD~^j>pb^0+C(mC&?b#tSaDpFM~SK6+* zuRe0E^Xg;QtFL$7sQV*zqcPQbz^r^4&zr7^_DyYx72agmFLGFbMIK3oIMox8*mMs^ z7=I`f_lhe+v3xLN0Xp09-Hzjs0Mg6WV8?R6$6a8%I3c4a03*XJ*cLkh3Zs>f78o?7 z;w+3pXl5IwrZeNEp0|keY7P!8$Ptz>`>~55gXuAyy3|#P)Wx5ZHt{oztn~ zRR|fBR-|((gM$`8-n>^>GpJ%i*mOuahJD+rLXeTOhbdu)jbk49m#3-!Hb+2v!FHMu z;_BRGY{pRFJl}#{h${DEyc;PjFH`Eg4KZJZ(E>ObQ7^!lginmjC=3@Z_v-wLvEp!v zF6qT&kFs8p^9W%l^pZEgPcn*Ug-i89-HBVv{O~F&Med?VFI{7d+z*+DQ8q8#vK+SA zUq_MQA?c2+oWxrZ&n0C&mz)_RQi7M6KFgx%3xLNQXc>f%%bKCl1~m=mfkVZyMG#)X zD7SxTE6EF8@G&NEa_bIP$NobuNDW-(d*+9!8CPo3{GE9X>JQ8-uHLae(r;cfKT1ux z4nDr)<-ul``AeXE--Rj6UzzWvrp({CQm>gGn(v#xh2KZ1)2_cp&bwaH{LQYw;E`=) zT{8|tzjgt|mDD7#zh=H`er&$$3Y`e3>VSW@D|H%rt^sF|!)1O1tEbMSBA5nVV6~6U ztI)u&Ag1|#>NKzcaU;O?158O;ur?;jEV(K$J_<*G%ab1XPDaY&6S+P(BFzJRS`$=| zD=VS?RbIwZW-)La#J^b7F>MDXWO@T3iyE}per2HV3#`SyPtTKQ@YC}a7Yr#{$f8eZ z@C#(T5IY8Gwgq(nW{s2fIBADSI|S_@j^eWp$^OVVC-`QBG+ey|zu-TKbNWI;4k=o5 zra$7o!#lt=$t8x<@J>CtgdEt3q$N2z8iS-9H+<(B(!3+uI@u2c;LMcCK*YDmA-N#x zn`!{qgh;$>avwMlhY)Q!GnU>^m1v2NCxU?AV82+790COJ zf*9?b;=%VBRIDx1aaStJd@Q(7{?3Nk4e`g5Evd5Fd8rOQ9mvHqyJEYH_6xP|G|V={ zpG|sG<#qEowFepRI@cSc;EFom*t~?wxZP3iPVt7xzDRGhdZD=7I5hRsB~;?r9d)PI zRV3DdtEh~>9FxE~tND)Pp;VEP<&66KWCaw06iD^}7I67KSy?Mt(Hh}svdINOeXB709TonE(0ba56IccR+7_6j zHAD3^Ny9!oj;kc5x$&PS4HXIJY-^%3(K@?+UUuQEVFL-;YoV`p)-f-+zRw!|S6PF$ z6Ry%zWN|}dgCnOh-{9Mzi)5WP{lHlS`J0M-JSX2+GvC4R_gLbqiDFzr1ha`7|5A7i zzQ>OZ_k=m)JhZhxrtppl~NirC8T#i9jKdb3jO#e&+f3D$MH4U;YG$`*1^ z=#qm8gLvz$hguI_6#oZ7JjeA~SnyEVni!4a@ zbow^^FeAUSaoICWd4)b&5i00Qn$TboIe<)=pI;Ijsg8g>gbIM;IKqZC;~;n}e77z} zfnqAW0}zoOH(v2BPsU;^n(A?L>vqiLJvMy-%|Ge4@@kdBZWwqP8XtKy{djeyLu%N| zf8o(d$ZXOM-1k`#`e7Y%pEOSlWgVOc@syVe>h}KlZ?bEUOd6w0qwGQ(DM|quA28F)>0D*