From e7ea81505316647c0ed51ea290bd474ed9bad348 Mon Sep 17 00:00:00 2001 From: Oleg Borodin Date: Sat, 10 Sep 2022 17:45:04 +0200 Subject: [PATCH] added lpf4* --- tests/{gensin1.c => lpf1.c} | 0 tests/{gensin3.c => lpf3.c} | 0 tests/{gensin3a.png => lpf3a.png} | Bin tests/{gensin3b.png => lpf3b.png} | Bin tests/lpf4.c | 156 ++++++++++++++++++++++++++++++ tests/lpf4a.png | Bin 0 -> 8401 bytes tests/lpf4b.png | Bin 0 -> 7177 bytes 7 files changed, 156 insertions(+) rename tests/{gensin1.c => lpf1.c} (100%) rename tests/{gensin3.c => lpf3.c} (100%) rename tests/{gensin3a.png => lpf3a.png} (100%) rename tests/{gensin3b.png => lpf3b.png} (100%) create mode 100644 tests/lpf4.c create mode 100644 tests/lpf4a.png create mode 100644 tests/lpf4b.png diff --git a/tests/gensin1.c b/tests/lpf1.c similarity index 100% rename from tests/gensin1.c rename to tests/lpf1.c diff --git a/tests/gensin3.c b/tests/lpf3.c similarity index 100% rename from tests/gensin3.c rename to tests/lpf3.c diff --git a/tests/gensin3a.png b/tests/lpf3a.png similarity index 100% rename from tests/gensin3a.png rename to tests/lpf3a.png diff --git a/tests/gensin3b.png b/tests/lpf3b.png similarity index 100% rename from tests/gensin3b.png rename to tests/lpf3b.png diff --git a/tests/lpf4.c b/tests/lpf4.c new file mode 100644 index 0000000..f8d289e --- /dev/null +++ b/tests/lpf4.c @@ -0,0 +1,156 @@ +/* + * Copyright 2022 Oleg Borodin + */ + + +#include +#include +#include +#include +#include +#include + +typedef struct { + double x1; + double x2; + double x3; + double rc; +} lpf_t; + + +void lpf_init(lpf_t *lpf, double freq) { + lpf->x1 = 0.0; + lpf->x2 = 0.0; + lpf->x3 = 0.0; + + double order = 3.0; + double c = 1.0 / sqrt(powf(2.0, 1.0 / order) - 1.0); + lpf->rc = 1.0 / (2.0 * c * M_PI * freq); +} + +double lpf_apply(lpf_t *lpf, double x0, double dt) { + + double k = dt / (lpf->rc + dt); + + lpf->x1 = lpf->x1 + k * (x0 - lpf->x1); + lpf->x2 = lpf->x2 + k * (lpf->x1 - lpf->x2); + lpf->x3 = lpf->x3 + k * (lpf->x2 - lpf->x3); + return lpf->x3; +} + + + +void dft(double* x1, double* y1, int m, int dir) { + long i, k; + double arg; + double cosarg, sinarg; + + double x2[m]; + double y2[m]; + + for (i = 0; i < m; i++) { + x2[i] = 0; + y2[i] = 0; + arg = -dir * 2.0 * M_PI * (double)i / (double)m; + for (k = 0; k < m; k++) { + cosarg = cos(k * arg); + sinarg = sin(k * arg); + x2[i] += (x1[k] * cosarg - y1[k] * sinarg); + y2[i] += (x1[k] * sinarg + y1[k] * cosarg); + } + } + if (dir == 1) { + for (i = 0; i < m; i++) { + x1[i] = x2[i] / (double)m; + y1[i] = y2[i] / (double)m; + } + } else { + for (i = 0; i < m; i++) { + x1[i] = x2[i]; + y1[i] = y2[i]; + } + } +} + + +int main(int argc, char** argv) { + + double t = 0; + double s = 0; + double freq0 = 20; + double dt = 0.001; + + lpf_t lpf; + + lpf_init(&lpf, 20.0); + + int c = 1024 * 2; + double xa[c]; + double xb[c]; + double y[c]; + + for (int i = 0; i < c; i++) { + t = i * dt; + for (int n = 0; n < 200; n += 20) { + double freq = freq0 + n; + s += cos(2 * M_PI * freq * t); + } + + xa[i] = s; + xb[i] = lpf_apply(&lpf, s, dt);; + } + + for (int i = 0; i < c; i++) { + y[i] = 0; + } + + double rate = 1.0 / dt; + + + FILE* fd; + + fd = popen("gnuplot -persistent", "w"); + + fprintf(fd, "set terminal png size 1600,1200\n"); + fprintf(fd, "set output \"gensin4a.png\"\n"); + fprintf(fd, "plot '-' with lines\n"); + + dft(xa, y, c, -1); + for (int i = 0; i < c; i++) { + double freq = rate * (double)i / (double)c; + double amp = sqrt(xa[i] * xa[i] + y[i] * y[i]); + + if (freq < rate / 2) { + fprintf(fd, "%e %e\n", freq, amp); + } + } + fflush(fd); + fprintf(fd, "e"); + fclose(fd); + + for (int i = 0; i < c; i++) { + y[i] = 0; + } + + + fd = popen("gnuplot -persistent", "w"); + + fprintf(fd, "set terminal png size 1600,1200\n"); + fprintf(fd, "set output \"gensin4b.png\"\n"); + fprintf(fd, "plot '-' with lines\n"); + + dft(xb, y, c, -1); + for (int i = 0; i < c; i++) { + double freq = rate * (double)i / (double)c; + double amp = sqrt(xb[i] * xb[i] + y[i] * y[i]); + + if (freq < rate / 2) { + fprintf(fd, "%e %e\n", freq, amp); + } + } + fflush(fd); + fprintf(fd, "e"); + fclose(fd); + + return 0; +} diff --git a/tests/lpf4a.png b/tests/lpf4a.png new file mode 100644 index 0000000000000000000000000000000000000000..e3e56673102b42d38213f2925fb2767103430d13 GIT binary patch literal 8401 zcmeG>2UJtpwl^_IQ4v%`K}f{00m7h9airub2nY^M5K%%Hada#YIt(SbOc6(g;Eakc zN>D_^0tlK&A|w_BA0jH95JC?vMIe-fB>x6+M*o!e-e2$k_1FJqC3l~*&)Mbdv-jEO zX8+{w>bO8dPXhpe1-qQ=Jpo`Y2mq)ewK*{6JU1#A01zA9cYEzrC=@UzJ3CtmSO5Si zSbzc$2v{IcRRwH6+Noe+00?A3iq8X}GXMZWfC5sG6!p2r#xN<^@d~EX)6?6uY15Zq zetG!tVG@aS`SRtGk`fk+_2R{g-rinVglOZ&R6**U=01_Ps6(L;aK2Me9~BA&iY5S- z1{{K)B%?D%W@mh6`e&9e4}nzyfSRP2Lx5661@Q8{z6ay3kpy`O8r4aHvA1x zWQQnNECswC1*?XZB~pm8sWp!jk7`&&@d_|ofhjL9r%)(pG@8rhGMP*~9*@Oh1xr>x zg4e7b;^Z3+0GdygKg0zSO??1B@7!g--76-8&sh-L+PL(j$KN#jehF1WjZ$_X6!nXq zej!hNvGlgsy2-nBI0MJAC)^D@Qi6b)ef)ub*5r(jH1P}AxB^MtPLdP}v@c+E9{@FR z*b}vB2~9#8Cu9f-*pA8iHek|x&`I4iT|iVJWfZ-nGHT=_o;ku`KY-MwB&nKa^lLIj3E`+wyxnCx-8#PX=4f=Q+2dFn(P zm)6@@8F1UxspxjDK9wNzTvHphYU!TPd?mZTyCxgI6Kr#6t>x8O+3T$~-+bVu)M+(F zIHSzRbJeRDf1nF;q(G_M$HLD;4{s^8QF8gWvIVqYiy9%3qHDHt)9bnN-D?nt7*A<6 zrD{xp=dHm2wS%ixAd%o|_-79Sr_|H{ZFq)14{I_0xh$Qu$dVYM>`T8|c5!OMK79lN z^BD%PN1>oo(0kq=5Mqoz4I(}n0d}E3F+$sfDK>&#jhEtRLyAxf_7%V@q_!l!*?^0D z6D=dM7jXSwM6Lriw~d+#r8Nh@C>~9Ds)QS8fD58+tIZKDc>k})IMKS{h8iFu9CfIF zJdO6m@_o05ukvi5;@NB;_URrA;bX~Qxe743bCWGILCc`SyGA>;Mi3GkZcU7{=)M9h zQOd{8=R{4)M^zJyR5YvwwGTrQijBztmj?4d>Z^@p?|k09G_Q3 znW6$6h98~Gn2l~Al@Uy{!H?ygk&jncCuyjG+9jx-BK>i@Xk#a?SQ=#6M!z)RjzFwC zcoE8Ahu*9@9RBU{u&wRyv5nBMx%p;^;8}&=EcGFV{ua9?ikyZ`INpD9=uD}mjuvFp zH-(X08EUshg~^`{cvp&Nr;jFRCD7n;-HndmIF4ERaQYjfZle#@$^`@r`%`!mUlJ0w z-hZt68sT4gyAa&Mcrwrx^Dxq)?g|p={?rsb=L>(^xnZK8egMp$!V!k)I)owzdFdzeT&?{F`gxwgU6MZC)U7`qElVo)0=bXe|aXKXg6%s~|HW z^Vr8Bh9|N{RbRW^5p9ovnk-zv6dw0FgA=rcN4|zH|3I1L^d}*54BhQBRE$CSd38V= zzgzp$&Fv4auhwfY->RTc*=m2~AfJ-lA165zI8`M>Kt1MPzx=;gn#vJOQ@N20IGHuU zkrg2B)2O0asa_LGWK{iN+*xHRYDiy?kqSaItAgWV!{J1>%)l4}+_KyQrm*>fjB}$d z$>oet3tKpN_7bw7jKNUNs#knV9_@*^V)9ib%03P_{5KcyW#OimaCqj1N*ZB29GLsc ze3fseC{#zq3++@QjXN#E?c*qj5$7Sm+OsLQR>;;AV(Y1ss;^?P5x=3ckJP0Tg-`M@Taa~ zabq;|azlDKv`&rn)D-bWgsh0~sxrVJGWjPD=;q8E3WwFJdK)S`C;1iqGl4s0JGYn< zG@h2Gmk`~3Fdj~(x7*_pUl?w&qhWFLIUxzNrvVi@;g%l*$9!x<=zClNbb>)#HIrPS z6=pwuL(<(l!9-`m7;=5&z^IhaPRc1*cy-Rhl#l;h%CW;6V_N4}=Y4Ctu;=rPDf`FF zjCrE7eQH%Ad3#Zr`PX4mKj=u$nf25f4K>ye(7pd>!1lTg@T0mD`!)pZR)fI=5itX2 zyB~n6$2mx)%lO84fzXlFn0aIkHPLOOLo9=A5dbj>l{4YnSJbHR?^0rV0rG^9GXD3$69ANyMWbr(Hg?xsN9%3k1?|)~l^Sa9Y0!Dw z@f3}=P?&fF3LZrw>sMoP$vRruW{C9(!xCgybjYJ2t>0$D?C2j%Ovr`W!c9eLNKeT= z#LuK8kiIA0BAT?;e$dM?f{+9(o2x_H-U&*`A$QnEa74$?_XQC z2Llm^RyC)AK$UA^zm5s@-&&^v3*n`HgTbVfWx@I-LR+KCJrId$MwpaL!k7o5l(!F- z2)Dh!j$E}o55bf*GglExK+!V z9?IlvYsDBK-)xm+2dW%uc{k86NlD0W+6{9P?MFO6?4y_JooqX(lv>>1zn`%h!(yse zUCl9Hmg)^_ZwD?lz#RwuzYgp$yrogasiopcX0~I^VwIi=vWlSSF=S6ZEpVT)WOCGGmyOhjUDl*%#bUZC}r3%(WyF2uUmVym~Hyr?`j?U-9Y(tuOz7iyP)cE zB@vr$qbsx>XTLuAMov7&2#MPDfy7BySpi#IXD6I09K1z%u}PYyZ_#wuH|{F0ivq83 zfr+_xDnoc4x|jdHB$1yRJvft&pI^Z$%&B@?0NB7Vz7$nVZu+1==~^5q-EB_#HA z^T2S8ux`YYpH$gS8R&7?z}F#_I?mQkHk8%P+}~X_Wu4u<&$@3+IXbN)L(}5K zS7X_eQc(k!Uan0`VjH?!GJU#B_|+4G<HoiTDniX4|;U5^{@> z*vu4+$wJtQByOgZtm#+g#mFK^F5f$)JKdX=|E?-@U*t+Cv>Eq7Or`bB5lB`FmX-6` z#!15Tg+nS)!3!l#Lc;}aa)Dc#ygK5ETQ@<@Cv|aEWdV73A7zbergBAAuPRP*m?ed?hm z5koGQKP>AN&Ao5CsV8RYJ2PHwnt4gn6}G=zQtF;q(;jUcNu%cJm>tled!Vhn=^i@5 zJyTCiChiN4gFiwU9J_jma>>-uH%v8dQp|MBr$IXnL8Yw=NvusO8jK} z;od9y;?g}=^dpDYpoSBr{FA}MQJGHrh&+uuK{;WppbyxIQ#a zt|(|ujy-VmO)`CduVpk{w)O{W*XXVO?aL1wEv$Yb3BRJ3UZ!%bjkC@0Cg(=A>;2&R z&sRj*-ItD0z4G-VW8Eo~1{f#xybxU9V9YG7<70C~p|Bw)_WQIBjTFiKZah0IECf$; z?(heW1mD$>#bnWY_RQ4&?D9E zuk*;&#+JvL4Q_JsXu0k64?K82DZR*%yY=aMmV?J&1Lm_#ijL>y8uPs;aQk$>7BqH- z5%1_nPKL!_){i8v6CCH{P2{$lIk(+_Uu@6E={E8_32xEu^FE znK!~{?IdnMzW!!O+1bs(^}Q=?;(c>jt`wKDgj-LOw>?y>jz+lnz@zEb^q^~n zJe6JT)@kQ+tZ;&C8E@)!u!NRr%()})_Qmx(BnZ5C4rQ|FPAa@_a{SGg$*q*2T15rE z-vM@qO#R3}r3bvFsE#(2Th0H)%`tAwAy4-QcIVJqo;9(AoaxND!|AQNe2jC040FBF z-5WlxnJcK{Db3eH?O3}gxSdD{g%=ul0nZuer2K^4xZ-$1%Xbrni$jQ^c-}nYc!5`c z@B&&7x1MwsOc(DJl2X QDY3hDy4n})2>9+l03{}-mH+?% literal 0 HcmV?d00001 diff --git a/tests/lpf4b.png b/tests/lpf4b.png new file mode 100644 index 0000000000000000000000000000000000000000..af9acfee4204e44721475fc8c3442465f78f8810 GIT binary patch literal 7177 zcmeHMd03N28lOZ|#2T!kC;^gMt5z#fE&+kWqgYxHpgwv}=SQd0v$L~r+_=H#^P8HQ#9}dGLT+xJE=w-X20=F z8_50oAxn=vsUU1E|9tO2KMK+Hi^s{%EJsqF(^FYjRz3V+O1c=>+m0SgL)n^j)elv6 zn{}!EZ0}S7Xso)6(ExdkIi|8jW}wf|P5-!eaE$((r1W7-+({JRDL=roU)CS*U>P#R za%x053UztJpJYA!g?(>8?wFkqL1tkJK3D5D`gJ(!>Vg`jY5Gg6cWWHaRRQ^I>8f;g zSQ|&-{5q-iPj0yZgNIi88`*xF9cmd-oMe$=9C!%|=82AvrVi_~g~yWDxGCS;)M6Pq zU-yH{t}8Ze_p_Za3TE~cqrTGC*zA&xSGQ{KJ-E|}b~$OEt5%`k?I-sO<{M2nouSfJ z2Ge7ymcew%v20bl^+3Igao~W<`wVko$|j6eMA0Ai^_<>>jVSzO<1MH-b8YR{^rYc`Sdo8a&>m1=C(CbQLZ$_pisVu%#Wf-0>Nbhw~G#3 z6+?7kmnIyM8tcSi7-_R>OO_&XtL9byOZSJiL{%#e-iK2Uu64sXDWZ>_+kXfaAzR3!1k-!UXLCzH zAMVR>rPTQhR34ZD*}vH&!VI@D4Ior+)OHaLBN&p;wd+&gdk=v_e3_Wtu=MU!?+-eC z%P0IUfQ@K0X#@Uk^410ao>d_666~rwT~qtfe7t?|YV+bg=Hc2*9EG7iZfkg|E_gae zEPsANk%10Ba?MLXklGv>?C+;yBkeg=p6aWh``P!X4sw`PwE&x}K#<~H8o6CGiqU7B zj@O7Zs*dBfcYI3d9zYhD7c?7j)Rcv;rA2Zjv78wQQRyaeNaA+s>hcXMA6C(L!ML&- z+%Dod#kDmEjzpewO7#my)HeGhj6Ax8x^$vqH`7oXajOtR35X%0H*1lt7D$Ezr7qx; zJJ}mR5(Lr_G`>yVy5L{F3QD536LT|rhALmkGJB?^KY?eA(x9I5NS9Xjgk{zrOMdnv zEQrh&9A@}audsb)ERIh}q81dFf@WXV9~s{Lhc_qJ>l$a&t!q|g;) zU0lwL;wl;~A1Pgv-W<8CF{VSdOVD5X;6{d?E}vcFL7EbbuS)M3QuW7DgT_bSh#I5| z{chBq8jQSl-MiP*p73)3{G*x^`GMiqX1!*sa z#CT3@q6c#6czt+TzYqrcOX=U5fUBdp>a)T@K&&Xj@av` zkmbL5%}U5E?9yK;mLC_ig1M+;)O4!-Sc-Jw;lQl(Vx2zWF~UNo1%cH(TGysG+pbBY z7!B8@wkHq2(N+otXG(uJAEp)$#=72@>L1@uR6mQ;(dlDme&@r3SqST**l5RZzy2wy zS~Ru;34jdY$WLHzk&sdu%LyBGX7(b0>mPYi!KpDv&Di~SANwA$mut6~e-GNpAL5@b zs@Q<~=Q9G|5De18(CB}OZY_y3`zdOKftFwRXCL2bdC(oU1GqDPn~Id&1@`abYfRta zTcH2ZOxImKwEQg89a8~;J9vDJ$wzoV&}|G^x?GpP939zp7g_KNtGB(WZa@Oma@IDy z$W%wiR~J<)^+cg!_mDz&A0`2xmFZon2*L~<$)s0RUL45YSaT$q=io&;yhBSNW{)_nW?a!)b)*$ZLpwsq%$r`H z^Ed=NovKm4L{FqPRB#86xa?f0RdJnq%R_^IicHWtx%3m)M>f7XAT zgD^Bo6K>D386kUl@Vq%H2Hl~dhUkET-BGpelVdo)fg|G!M~t%2W;S@_Gz)G-QIQ6z_ryVKNX;eZdc z-0?gi9L-`>fEkseH9m}(gaP!*-ArQcLAojd$eVJIZ4Dw)CCD{)!87xYAIE&~Ov1h_ zfkfQU6zyf*PNu(O>08Opc9K3HvcDJo{T{3$p~x9)6~S@}ZR-f7%-RGLnljHD-Ed}q2BR?*Zo5dUuH^g5sDiMHifi_ z2N4a$DN4t-t0`2gA9%R|$loUI)`H^Np(TTXEukSxeAbm)B#;kw6xYjaLRh8@s^!T- zGfFA~nS|m47S^wxZGj_;P!!+tJjNBv_HIT+*~8K2<5Ue|MuT^I3Uyso6-P?}S))?Z?mMa)H1sp^t_Kv}KZAIYvArR*`Rag}>3H!I?ylQ-ch zgY2|^hcK)&R#aMk+*t