From 438c657f65a16a5c844f8c367c28b317428b00ae Mon Sep 17 00:00:00 2001 From: "J.A. de Jong - ASCEE" Date: Mon, 20 Jan 2020 12:10:24 +0100 Subject: [PATCH] First work on SLM. Seems to be working properly without pre-filtering and bandpass bank --- img/subsampling_slm.pdf | Bin 0 -> 19430 bytes lasp/c/CMakeLists.txt | 2 +- lasp/c/lasp_math_raw.h | 3 +- lasp/c/lasp_slm.c | 185 ++++++++++++++++++++++++++++++++++++++++ lasp/c/lasp_slm.h | 79 +++++++++-------- lasp/lasp_slm.py | 6 +- lasp/wrappers.pyx | 70 +++++++++++++++ 7 files changed, 306 insertions(+), 39 deletions(-) create mode 100644 img/subsampling_slm.pdf create mode 100644 lasp/c/lasp_slm.c diff --git a/img/subsampling_slm.pdf b/img/subsampling_slm.pdf new file mode 100644 index 0000000000000000000000000000000000000000..e0cdeacc13a6eaf1bbb2de6b9b4d160c2341d84c GIT binary patch literal 19430 zcmeIaWmKF?voK0)ap~mYG|i)pfD3*$Ph7QAOVFX9T1RUKJFU403Zf& z7JerL{GCWxmkC%vgMfw((Y-FP4o*4EQDpPt*vFn?43;;fU*L%^*{L@ z;LHZ*(7Lbl33CSHd#;W*bT9VJ8$%71FfcBdpzB3DO-v3%Yb4xf&M3TuEINtTF@a=* z7|&oh@t%KHRDv+6_to{nDZI68eR7e@3WiUc?&W{ob1e;k-31y)g@=b$Qo(LK^_hj( z1<26x%JYQ}boBT^fcehe(L?ghYo_1uX#rr-3~mj+D{lb2eXCf&P861S*6Y5s5|-fa zYwRO!WuVmEzjfKJ;2cz%>YBxVc5Ksg6_bGr`Mk!(tO0nw7~|+$=iG5~|HKI?v3Nd% z$bw6=?o=0vRl2O=D-WqOQIXfcZW}(BCjg0@d=Lh#;Ylku%k!^Z?K_wCzV~uRH^@AS zNLvTnEAVA*l;F6qJbnFllLAV96_R6rAr)Y?a47PS=SRY`LY-U}Cd1UXbU*0dmQ8@E|@XA-c8@s=s02sqg zoD=wwjk&sV(sjJSbc4y$v&+yS1TFf`Mozdm{sJxw=bC0NXwE}b!#{3bFX@+JECv$UUx9}j_clS z+&VZ#*)Etkqp&>zK@lmhOkZ0M&%+(6c58MxC`o~e&)vkMb)_`7bee&_dHu@q6Q-A= zdiuJYS<>0UeOWa!S3Qkt!>w|w^?+d7w|S?cE{ZXqzVpB_0&aPIIooZqJ*ck@P*`tY z5!Mnp<(_GQDaep4|85-%@L~p`d z{c(vq1{hOHgQs_fU4;;#+{JobRW59}2TJURbErb8=kfd&3&_ zKgmp6yE-z`6K8886|-zGq<$yDm4TXgwG@*!sMZt57n!!oQDC^7gc#|~j;kNo5fzVr zuN`h>-a!nJIxTxJXBChlQsV}WAWp&~E;s9TT-L{t*j-c-WFXrjQ2ZTBWSa$&3U$?< z+ASOYRv1bs@>4cdvX{YxznYf31o45>w% z>#)A2E4sBrV@KQB+!R&Em&~g|af^!aVU3v3m zb)KZ9iJdZb9V;m`VS=a76E*a?3WIhm)qYmq(|1H39<#XZUq0u%bohR!#_^&_r+2f4 zy7V{Pw5e>14BTXht#@Vh2i~wCUtt(^*M1Yut;N3iFlm$xyjBOxys>78`VcJP@bn8| z5+=?HwY47puAbV|X;U^%LjHE;26QWqh=4vh1MoA)rFSO^s>o~a^Yy{m9kZ_8&dB4w zx(IbocDXlIa5O8-`$h;moe6sH;#Q32yMvzL*hQC~wCGFi-C~Fir5wsqe-1!Uk~8X- z+2u9E-juidhLKY~;f&n7EosjP=V^rrQ#saD9?AsYQN|7k>f)WjxNeu=oQR4#69^Oy zm2UQ?^_Bma?wzD_to!{SM&q zYi;rvpkCd8hz&Z%ZqqmP3we}AudQZ1*gTJCj-wbUGa1f-;aryJC$iLm0Wd1~-tP{{ zr|yO>tI98(=-oVKmMyq=^s=thVTI-8H8tbc==|f>{B`CCDQebLxC>0TbDt6t(3B@| zYUcqvxk#;Jk>U?+1fdwz)WMev(m59|JPo%8h~D!vh~dO>7U1L=9m@u36X&jN4Ph196&FN%5g6?)1LSdKPhpy(ncv9{?saBo+MFYIw$GL27^x!(S+?D!S~@+O>5_XNrs;0 z4yH#jgYdhyk-wYt7Y@7L}p#l9opB5RGG9ai!^uEWf=yORB6o@pUa(IU$1KM@e; z?r?NfIGf~=F2WpRaOreoa6h6=$J8WnFud0EJiuKYfpIpz=p1z*S@}3dt}CFCy9$Qi z%E7FLk@-t3Z6F^9C-hhA0E5o#ZLbBq@h35^VNW$~tgs`jb4y?Zu&xUmastH~TS~mv zncw7ZeF2+LSrDh>2wQQ%N9$$FAbo*G!`SX~$5$!gfX*w+=97~`bUz0VQr#bNT5dHDLk1O?n?+$CzV(W03x<~@J% z>-waCt;xt{j&pHl9M9kaQZMbjMzbMUjn^WE4VC?zf(Jm+1<{R~hTQ=lCkC^=V~n8B zSfZ8N{du+r!HQAo-%%BNH(p~}YHacfncs*sl&fl0(Y*N3j`lfvi(U_T44LTdjU$F; zG_32@X|IYCHNr?#Q8t0EB80|7IZg0}@gg>g1*Z&&)5^tu9swu(TZN5jbq6g!PCHZj zk4!GlpMJ7OLi7BbU{DH=EB;V6uLeT+Q^@#QFHn4Vo_shqwnnTa z+t;?k_?bWFHV0>Q-d}hl8xFE#J%8RNXP;juk8;V!_hk%@FPv)fVQ1B_#{yL-3De90 zS+T+mNhBijK*G*fro~)G^C3TxxypK?&7{JUmkHknh~|Z~%n#6uTc@c+uOlV z)|Ce^?~a^cJ23Ha#AWXKmKYs0d62_~*P@kEPhG5=cPk@Ci^1j}Q7fyDqotIzq6cX1e>=}+#s!^w%Y3Y3LIUhcFPD3E_AfhbN zizD(2iS$9d#xGzqJZxm$%$HafS}ZiE5}Th>TE0Ie@>1U3c3kj_5x|Q}w5`6V3wugS ze3R@0Ii*LfPBUOXeITn#b-Ge70jt|DM(j$fRQ;UZOQ|aF`dWj5Nntqx3c>9#yjw^D zE^}*RTJK9>!Gd;Akl5r!3+v1)y-w)Y@N%_0s6t4JfQB~-u%XNpq)Wa1>OwEvsT^i@ zarcd1j!LD=kUNcT?=B{j!Xg^wIe!u38S0s~EhcFMxEr8s(&C2qT2yOF`q@KjvoxZQ zic)Zc6xxI^;k+7?si~P4_%MB^XrMFCJc6)j}#WSDplg8 zs?j1$3|2EeIlLy#-I{vq>W9J&_{WnrAxTbVz-(~!E=s}36Opez{?JSPQx&4~ls6u+ zv{?%<;Vg=6Z`GbTrG;G*%RAdn9T2V#j0pC0TP?LlTnQBzca=|vc8{vP;VCCYLpoi5 z8cB=4kUCTl@W!i6-nHZF=a)VvnKaBkxs6rd+El8J6XWMC)w(HXMozxHY3_ie-c;r4 zV%6WMLzg)aX=)}I6#l+c7yisDk%SsuRjAk-4abi}624HE$IF2dwXj~yZ%$kBgQ@=d zU^~TtYZK-fSJ27K>*AVTTz#k$troEI)5`k;<+obuxn?+pmTvyQ7_Xy`jF6~pY1stQ zU*&_{U}%SRNpX5Du=iHhr-}CyjA>Onr{wN}n~!I(94acR%I^$VurvbNNC_gPOwYUx z=yoN|r$yb~6g#Lx51UpN&IQ*|lQPp!)1tLwaFUN!dQZJC68_-FQL=nj6}jt#>gU;P zsI^?BO@m`_Vvl-MyWG(Dluo!vg;W<{Z#Sy-;_f7I}VFu0-| z+!r3$Gz^k6;1+9rv=HmuE>y;QOqRNuq>9+%T}{I z{Rz3mCTF-|XgFiHV z_Da^>jaduT&}J5B)pPzXB_jcT3w*j*PU`u3;Tose&#^48&kwsfV;7a* z=s+WSEUPi(ToG2o66e3oT~@#2^eqlJ(_3;MyJ^}ZEcWTo4}cI4^&ks2tCKR*#jewn zPMtuD6@L{4cwsCg&+4d=Ly6ktNeZ@qKU4WyJ!Qt1$#XU*VO2YebXFbsM3z5C_eL_xwA8xXTP8+oIi>1)s9y@)~PgDuBkxqjV(cuPo5z* z6B&bppLrSp5JXGPLJ~Htv6RtnJHG^1ppo3%`BF}Eq#vKgHXTsXyZMZnX3$mcR`>2y z*{XrJIlQWG(W(zWY|$#>hsh{N&R^n_iFSm5@$bx~Pzh^G5i(iUTe(Ee#y#Hyde z>kKJ<<~V%%fSy%T8f%{J5PU&aiBuRetNbcw^d*+J$)(M$dMf;-t?Ih*j#1CnN0354 zleFYo@!%KlToT4zhlcQ$r?ke^t|AdP^6;Z;=Ht>`WEy-h7&9M2p9dyO z^(wIsk*R(uVSO`F;LquyjJ>pJYL1zouHxjwV&$=W!_~~^aBjY#ZMRQuT_0(4%3grX z~BPe!cQERcf>ifyQzTYa(XJ#qU~&D5FZ_k zt?zqRKCQDPVs?n$oR#6sc1s^GAHl5>hT8PZb1~jC-PB@$u3$(;S<>f5ocb#^+ z8s^BAjx>S|U%pt(Tc<4*L|KWqH2v%t@@^WeZww!DHV$(EiU z<4+JHA*@$|w{Z2=y&(4uqx)?8(aK9@TUmC{YH3r4ps`85BQkFbyP?eF{K5echA>S~ zi$Q5BxBmt?>+` zJO z{2QDiUL*7Oc`PGs&C=Uaa@iN(J;!|;vcDFWBN|@7k1h`xHx_(RM@^<|fvvMc%r|zQ z)@dIsYJ@lyqJR@oq8ETE#5;{pMu~jt-x0~D-Ime=ZAYmMzEE9_y$mIl(BG~f&#VDk8jDGLJS{CXiQ>XJo)KX*$v(4q z-!)+X##00VS!?WM-}7!SJt@=JMN+JL)h0HD#5rjep^EntMg6mw+YaINLjzwEy{hdL z_1RbzV`2WP!&Q~g-qTouUM%42xWiPM+s6r+Y#%y!a0IvL(4zC24=qE$FFz~^Rt{fZ zsr18qmU{^k>z#c_jk=ulnu?k1a;$~+Dl(4dbe=JS3Q(zQZokZwEfQnya;!v=uVj`` z*Wdim4$qUyP|m(s>Xnf>a|``lHarEJ80~|vli1olVc1hW zLXoS1Fsr3EN^l17eDSK>T(---S9MRy2ZhBX{YH^P7jpU)Bqla0NMz4X0rh;z2mb``Io~f zo@JI=5_iSXmNi#Cu{p`xsS_I}Nt|ovG+JPDc|Bwz>EC>)zj<5;JQI^VcNl``N(@m3 z22xsr5pB5WuNf4yU45yYRr%IZSw;LpQhB=T#4SuPPo_Oy`K9GdDHiJveI)aovZ2GJ zf0k9d8#ka^ac%5;=-k+=jOX+9RC$fQuBO!mHiV^rwomwz@7*}X+Z*0-FE=EQoIXO8 zL*q5pkIa;+6tRT_Kf9Asm{U@X9}BrwGy!1b5`)RMj#!|%MzA^vm0v43&Sqa^EU zPu4BpFHzcc4tMY;vUlb&z2C;04a`rXzm#9-q9CDZ!>>?F)?dZz&>e>MPiMma65tQ^ zN@}}I!)a*AEy97u;lhE@O%VEtc)_%T^JVdYKvv|{o;le@tw56Hn6Hz0!zV!ajVNX0 zYgz$>HB-u*PiPsd6-_s^+ZY{Fp1L^oYjr(TvBM%CBe+M@Uas4SFH3kVyS!tZ#eV|D z$Ou!CJUwb(0`J?s3A-UG_{1iptPSI(AxXiWXa~U@{K24(7+d*p`Bj}0j<1U7M?~j( z_82kcwMZLsfqh|Qb<|HQ2y^6DDSWAN8n@dXr z)5}Q_VoIG3ZwV$%3Y1Z(?!^a5s~dYLXHT$7&KG+38+?F*50 znYIZdQFhhuc37r;39%+8G3MZ;>4Ig+P0TgjB5pGl*~KZwwVvf3XL8@t@Xs!E^s3mf5Jz<3p_m6Qu7$H|Ilpgawn1PVUMxCU#2(6sf}w+ zwww2@K65AS=wrdRyDvxrpnW&{`(pX+XJXPaR_NJ58AX`7E;BMUsYoqztvw&74RRM> z;J!{7iZwGJOUqB^likS}g111b*!Qpvb7-jM zd$9$dWHt_YA7bIwyLSUN2HgX<@0^@+FOp5f z)d?#fg?)Ouf(7JDPOF=l26ANa72o6=#YmJZi@1IkvC@rJvhs8sQI=|xhB zaopb)ExQxo43$`!E_5bv>QSVb%1YWV%n1b%-`fJ3nLii$?_^E_%6l-23^36dRlUJ#ZEo-*NbR;Y&WhqL+xK zEw+G$XV{L_Lj}km;*P;XJ<6V&EvaWNRJGXWZ)sAuc4NphMb{>g-xDtF5HgGL$m%G- zNO$R-EvafKk0eA_uDTrI2K#hc20Dv;-RBag!BK_Dk307j79*}n9M`bU{qjW-GutQh zBh}lhtgwCxTT)J9#^XHf#V7^KN}K8COv0W;j)AX`OFUPl^h`LU1GMT}RQPW{Gi#5q zETO!ec^z#QV07#?-hSqF-7La@y%PE6&XdNTZ(E0788R+AB`=Q(9wyBaPOb;Y;5mUd zUFL<5a?)QM*JPb^|7N0@~eCeH-=#+}-g0d2ml@G-X;!aGQE-76kt#U%!u| zgUV+Hzk8IeX7;_bVl19&chw3lnU9#mxJ*MD^$q)X&hQz z@d_WsXNReT8{pO>sMcC;qIk^xdgEt$cTYnPgfnytlAUtWIAe#VFto7QSbW&!pcMxD z_JvEFaDty@az^=aHn+jA)h3tIZb)R*<8Uw{7AV-23Z)y%e66S|;&krfUqqI(I8 zo)6|XNs9BVbDp!^LyJ$K<_>RssSQxIwh?~m$3X$%@-YgRxvkU$CvRh=fn_iO>Ql^R zX05|UpS7I}oJX*u0Nto`TV)9^hN3%|$`(S~hGfzfE*$zY81J&ZSL*N=qMEHhMWbO5uW+9_`Ng-Pt#_31u%F zd%IxKjy#Aq#*8US+2(MR_r*@;w4Zx=7)*ef_io1yqs_!ig++hC+wQ|{8C6;PND}1Z zgo0e*Rq_F{Xy;EW!d6^O)**xy!0^$a#iHU=C_E*j5d*1{I!B9JSY6Bk3aQ}qH<@E$ zA=^0XVu&XRPDa9`j8zUvMHv3>0+zL#l00>9-$53u*k+qc?#?rhO?~M0h2eLg?jeG4 zG1**1s1+B-<)?eTY}1Wl`ax+mhk3aH2egq#83nbb!tRuDORvR=1(OZXF;dEo#3B7M7jG)lAQvqkRv9-rrxl9a^>*Q`M+=3ZYrmq`!DwGnO;{MpOP zdQ5&&DKX_O_wHfi@oZFCWI=kkEJ6o!hoOte`j5yOFueoh1Y+uuPmjJVC^)f_^BM)z zaHE|0!WW&rK=U+q?2l>)B>yHd&&FQ5Gt%t0$w7ozssA-}6zc;-r-SyQ+9^XheK4xV zN!X;vdhLRnz)&J9!@euanV&d?aI}_=t&#HivspTEr`@*ISrVVS{UzC`n&Rq;iMRe~ zL7^9I!Cl@dN9XQe35}1Y(*UQfA~etQRS<;JDU7;Oq@A29&m5E4hcrpvAk4Uc75(Z^#xGpo5Zj>?W1nX*qVS7&}B3ECT7g;H?FW`%#ePF-4+jn|B zZhoA_n5vvzg*zFa(0Ji5zMwT==Nja-kv|x$vK&n7;57%N)74Q+^Wv!>dv+i~Xk3u5?%yo5LCbH_JlrFx;u67z=%jap;vg#=g) zd*&UpiK;HA#&wJeP}Btu>BQ?YcrzBXHyE8C217>eWfnF2MA!s0^(vDx6Gt5l){tSl z3CWlan?J$Hw-J2h_81^SVU^*yD7j-{o^tmei&+zbnu8{4=Kp4JzD?{!YP%g4gMN`M zG5nQ7ZlWYI0HeZO*>z@e=XJ5j+up_#((gg?rnNob0M)R7=QG7p7|ibquxQ7b3T;!^ zf#cM%$MQv`yyDs5BWg;Ru1;sZEJXD~SRC;Yqn$V83iBnKlww9=(C&%ua`xBh{+Y&>{zBoK6d?&{ zg$mm>-8^biIE#;dopuP-Q{tV`zA{?o_;X`xL^kqNAC;JsEM83@YMnOTQL#z^v8TEQyt`bp~X&H{oH zQE67>A}e>l+&hZD`n1;DRvgY>1bfqkTUteZCx60fG|*kC^(^352>(8aMsi&68D?Usf(ru&0_yqZ^!Yndv^XVVa1pl9 zdFNaIRR3AH#_i!4EVuYff4{mq5{G<^$iPEu_ztD$DQu)^wd}2ZZnST+MbBTtlB~gD zZH;Cs>kk>_>*C>yN%z1`u5p%8mjPvSbv~FcOL)aN&s#c0buw$*Il4o*PB43_6U2EU z?9=I|j5JOtPbwMd<{7?I(75+-YRJ$QF}}Elt`5j@!>LN8A(Q@i825Rq7)VIl2+FkHjoR@ewC#TYjdmNg4j#g3XA3rc4xN&I_2!ed| z>d>iUO_*mB+gv^6n;SpXloH|&PwN(S{?2?;vYZ~Z`SpV2!ZTY``X~aTC*PnQ!IXPB zW~9~{AbHJ*(~W#&q?`@fpslV?GY4_pECrd~VCL7t zKNCGwL4Ze1^sZn`S8vX;WF5;M^5kiVhME~geO4utrm1C9+x#6%!I;-ZfvYIC?nTd+ zWxN|oZ=e3k>jf6y(Ym>byh08i zMejj>2uiBn=Vel}>2s1hS{-k@->MJmfbWDCL5K$j z?7}JC2!0(h5h|J|#oS!_hVDwb2>a}L^UU;q!$&+q+<12?1T3ebnwVUaiE>F+%$t?B za<0D2AA{|kbLvOTAsG6{iON{HP*4|53y5t--`F}g^e%r; z7nQL4Fw#?%ur(;6{;s#-a9TOiye?ySW1aJztvd6mdJM%ZfgINm{+w&itD5q2Tlk#? z!Q~ov%{JHNO0-v9kiDNGkvTOZD5US5>3R*Mq+%@=Tg9zV-rnzl>n2s;(`J?5JcA+c z1+;v?4Dk=Wgx;cJaqK=n$-Py0S=yVlfmwQucnDd(R-YlBOU=8H7p8uFz*A3~Mi=q< z2Krp+*=Q3>kV|2VMr+XRt#atAS5UcWAw}2mi~QIU3F!$liLZv1!VSBsNg3jXQjsB^ z(39h1E(c4v1o5I9iLZ-SL=~}wnWhO1Enh_7AcOBE442Yfl0Fj=jlF4z@3y5oKIvjL zRn6rA^bm>#3p#jZ@O{#_-}j#Uwk=zWoLZ4`3?X<4x85AbJuQ{^W*RGL$zw5G*eEIG zX*h=f=?hP?6UQ&rT5{gViaT4Ms%RH?W(ueW#@nWOo@%RP8`r{8P@!`d4Oi#(W0c%_ zO>Hoc!5dl`{5Yoy0KI|Ez{0crJ`1bhVq*w8)M{Z)a-bVEu5&6_l3JwFKpe z%L?*~C{PJV%F%<4#DWgL8ruVyejxDkTRQ{PX_)8$Gz`pipmVaIW3sHwbeaGlpS`}J zl>>l<9fT>MYa?oCYHR{ZvM~aH54Akx(a=9+iaF?7nCkOc8Cw`WoVb;DFtk(xkpX}z zKWUg4=|M%9=-SB}IsnLlQb2v6A<)POXk_XH`n7hn1Dd$lm>60Ct$}txd!Pf*k>Uq) zpyRrRfQK``a@^d1iOv2~Y;|TfMp{MyGcz+S%g-kp;{)W!lYy0$mJPs2_n1pZ2Vnl& z7bGKqiHV7p=`n|ufq|9*!1_q@fXPnB0HUU+qi3WAHOC|PM}a?{4=6vX`Th1M79$%Q zi0;Sp0pmaN9+~`I&O`qH`j6)i$b)FCbZnre_$e04gA^=mbfAV} zVSW&dm5qs(^+C-5MizQndI0kuarXZa8Hfxh%~ z3=eJoSUo7oLjO?oZ$bb5_B;Kfu%Fs^EcnO!AKH6_{-F2?eXQ#TK zFiJU4=QXvncX;UW0D2Iugzm2dD*&i$YT#gE4|wRue{^dQ-Y-5L>j52c7k=o1Khc0P zcGmh2_M;929fk)06%3snG=EzKJ;QJEA7ZFT>zNtqgZTX=iLkYuB}kAT22T$%UH}^{ zsNVoUN~2=~fGqw2;gJ}`v!t$rohfK!cwh^p+J`3yL(b41WUl&#_75lU|Ao!}ipzgH zc>n8cJ}c{=XYd(i=)l@D9q0_-BdrY^!WwytYT7@Q=MoSn_FvPW%|)X1V0{WD&Lz%B z%;+D*^qxtA<)sjzT8A_vOK0c+?-PS4js$a9OfXDt6zqSFne-2L0)Ac~_!ss8{}ul~ zhg^E5-{k-ER%Qa-TwrAZH8VT?|7J7OGqL>H&W4fnPuqE)`Jdl}guTLw**Q&ChoW(i zw=PUmb1V;1JUnoF^B#aAir?6C4NKgUG<-?7Zx%&8rlp{UQmDhzPT6XgwbjS5N|%~> zYun+ySL7viOv_aE+~!&-h-63N`3p%J%N{ntLcW}Ks%ASIK>;lf1{ePk?mfz(6zq0{Zh-T*tjis$)&7H#P zm#Ogzw10mc%BVt2&8NBa-s#>%PoDAqs8ob{i<=7to$zj4@CwHX037}owf%PN@IYM) z8xvj7Oa>IBfaV>106I_{3e?wq`~|tdU-=JF%rAsr8H|4tf+jnUgp7bc^Zy|HBjX{8 z|5YI97UClz^KU{>{vU*YWIRO7zX)0Xq6EGCsuP6q2kDtf)7bkD~HDf0|?dB!qm$6uhfrd669zfr;eaW zr6xR31vCS+d6-GU163binc)8;aD3$NulSOVmKi`pN6W^-LI;|7{PFVO{OLfY{kZ1| z0NLfg29cm9`gMI1H1^8SF@fp%XDawJ(Sl(pLcsq^Px&2H{O`*BOXSA}_*3@(J97UH zhaMG7Psa#&>;rN?lT3gI1pwsW|L7yXo%%0DurYzA{Qq)~_Cc$^G*O(;5`jkltozy{ zWI~zxa@9E+gAl$RHPRcB%z;XU>EIz-|N9cWo=TaQ%cc$<*VnP+m?S)jo7B8z!8keXtyw>V8sW3+K)e@08? zYa*VCZ~7D|UA$dMJD&WFT6`0}K0wL8roVR4rs~UwH;efD9h)LH8J?CpZalVC9dE5= z8$_$8N%_*zRU+Zn#Xg+A4~b20T1e6WQm?pX-8fDEc)dAY7QEY7y2P$!RuX--1X;cWjJ9 zFTIlR8|mx~Jh12xO9_KIq-0r|rlw`V0;cPI!RZ>O(9yh?t4s1s$a{4)+mThYaFR_% zqdNkiskQ88KED@|c#f0YVNlYqL1Sd6GHnAHKOI1-Hn>*y;Ynx`E(_DF9%&%LkS+tJ zAuR?YyT!BA+QdFxhS-O9ljyDoVhj2hP-&7InS)&qxB?eRstU-fx3EJ2&XIBH$Q$X( z{o3r-Pu4}2IL~n1?pB(4(yR*BRXiyJD>w&^WVv5osx0=W!3iWi{88|rM&J3$3?T5+ z1Xj3{PS_p5IqvW>j&GD&<7HalokVz=wgM3@vY&X(duzY=E)D@ory!&bO4C#!F8E~; z5318&MwUP?$l`I^o37C)j^tRzp-I!RX6RU@AK4#&Y%$Rq2;Zp2?LVq>p)nr8>`X#h zkr2I1iaNo3t7e5{#vZ%Q_#wW*hL@%agO77g3w(2BG``mRi2I4;NF6<|oD?VJdS{lL z6trl+CQ?8WU1w6%Q81kbYCw|AKwnUgKV9P6mC>_r1B%qiN})N}1V*lj%i|<&$)mcY zf!Lf$UpK8k2^?Ia?$gm2f4oyZZfmRQOmz7;6UQ}M&wxK57-yK-|CQkz&~?g=LDn!S zkD<1BeK+{E3#kMfsvLJytoyCfkwo75T7%C(Xy1K0IX6k`ZJvGsZ_aiRFB`?twLA&l zv@gC-#|A9!2`p|7aT#l5s~Y2l!%%z9J0$$uHZeuv7x2QXtp zc%T4iDFU=AXAfF)d9>sI7QqWiSvEurAORweC+8Yr_Sa5qp-3j+3R4|3c#4kR zdqt@%h9T?l#rF(lLH>BFR!*p7DOCkoJL;>O1i_LPIIB>c0zJ&1xix!IKFJ5-={0nx zvW}R@;$sVmRpr)pGf=Lev#R5s55!Jp%ckc_ACqf0QB-W!Hogk&$D{1$^zb%oO;()a zO-ZD$ewihv2nJ_jMvSv(Oe$Qn&h}i7)0E}g{I1ckqyXY^h);U|1~{#B!zXu=7pj+P z$uml@hCsAUA~pSW7wjyLSe`*;+h)BJ$tUhINg`t9V#gi|zV?{40pvQ0CD~u(pJAzT zuM{Pd`J>d*+U`WlT1?p<(y=GVa#6|!OWxdAh1&kgOQQFyLZQw-uVf*&ZfRxGrQU-Ynjqm&NSTXT3GL@Hgl16%I< z`V&|%3g0AO=KC>PUL#F=rRMG;i^FTcTwM<3)a4$yLY*;G)jUrm&fv z?y(TmP}CReo7j$`Rr3c>1l0%c%!|mO zHiJ|I3L&)&|6$vHSus@pyo*}cSPC(qFgDvoGj04dAJbFi8(|H*znI0kkcyuwe;VWy z8AQ8}VcF7h+OIe$IuTrEIHT%(A{6@37LZ2V#-SdYRQst?~2A<#FS09H*0Dcc^v*z(*S?f_Aw446O7 z^RbZGy_n#DI<*4JhK_hZu^g+C7?Z1*iM5F;@kg=R7b+1FyDHULKYZWx@e5W2{Q1gd zKILr>FO#YO;Z2RuD`FeEIr|rQRF36(IW#yGj0S7zWZVc}<~9!1UNnbRRG9P}Yy$`b z-?QGgW82-=_y@cqKQk{X)!7-Uhj0=_=oaMT5Hz2&GWp(6VDWXTO$-7Njc(MpT-}c& zzS=HI7*XeA;uPfya|japo?>+C*%tDyoOQ9Zdx_4<}2 zvR?j_vkz~MYE{U9K8H9@FQ1m54cgl*K4*CTY4N4Lp!dG9^cwGK_K=egQ++yltq@R| zd<(&bK1F8bI&1EGTeN86W7;tqmrXXh_|WcStKaF2zGLwr)tjY7;bq(5A(T~Q>LOQ* zJIwM$QAYPM{SMjK13yj-dK;`o_az~;j8duo3o6mLnM zU0dqNB-k+V{y7%}YE}UBoUEs?w~(Pf3zNopnDSiv;)*2^cbtqTHtW7;OOk!S{IpP( zL4kkn-d>lew>!umtn+-b`Sp3Bl}`8haJ2EBx+*zG{H>w$)DRLNrnISFP~ zdG-N?1&&F(c4wjwmuGb0gl~xQ%X!{Su_eCV?Gn^Obt5WoX|j=Q9O$xYArCSh`8ebk zV2j{B-TLO&XoE?dqVWWNU zc5HKc2^rI9Ism5YoGu}aMS!d2*r<|X0%64k4& zekyppYuALtejxNz7uv=M?!kpi4)>Dg{QdUbV(InM7VK33t%Ft){xV>SsI$WBli()) zX+kzfts`79rm}MJeagq_Xa(GX8velh>IA5}wl9cF_{XOUh;P`CEP_+}5{@ShLOd+| zXX;G$?js~1cUoGW4DO7;+=RBz@A3@LeCPg`;05Hbf2_a%gF|CasCPN^GemZ}kjH#i%p_rACHQ=%K zKPN|^P0io32!Q^B(8H!HIfyaPIt~LZJtGCEC_8HdM}5#g9Qsf2BeJQrmEZ&Ke{xF4 zz{bb~q610skGvmJ{LviFc7{gqpdCC0c+f2W=Lhik3oruwfYGxsKg3iIAAr?gFg8Zk zzr!97K%w8m=NIim0<=W%@c9oIXpZ;~7$}DO2MnaFN7%1=LED4>gfW1EroZE{v4G;& z|D*+J@-O~=mB&EG@()@D(0=AWXqg_SxqmB} zEh95%P3G?~7RJA~Gd&~wL$CQ8Ejwrv{ckWPM%I7u$Hel_wqj!cr)6o9Xux|JW$@$)$pND0_Ci&9YFi< m50(iJe6*lPgL$wZ&^m#xox@}QrDtQMV}~ap5t0^$|9=2`eHDuU literal 0 HcmV?d00001 diff --git a/lasp/c/CMakeLists.txt b/lasp/c/CMakeLists.txt index fed41fc..fa3e715 100644 --- a/lasp/c/CMakeLists.txt +++ b/lasp/c/CMakeLists.txt @@ -19,7 +19,7 @@ add_library(lasp_lib lasp_firfilterbank.c lasp_sosfilterbank.c lasp_decimation.c - lasp_sp_lowpass.c + lasp_slm.c ) diff --git a/lasp/c/lasp_math_raw.h b/lasp/c/lasp_math_raw.h index 06cf60d..66a2901 100644 --- a/lasp/c/lasp_math_raw.h +++ b/lasp/c/lasp_math_raw.h @@ -33,6 +33,7 @@ #define d_sin sin #define d_cos cos #define d_pow pow +#define d_log10 log10 #else // LASP_DOUBLE_PRECISION not defined #define c_conj conjf @@ -47,7 +48,7 @@ #define d_sin sinf #define d_cos cosf #define d_pow powf - +#define d_log10 log10f #endif // LASP_DOUBLE_PRECISION #ifdef M_PI diff --git a/lasp/c/lasp_slm.c b/lasp/c/lasp_slm.c new file mode 100644 index 0000000..4fd3795 --- /dev/null +++ b/lasp/c/lasp_slm.c @@ -0,0 +1,185 @@ +#define TRACERPLUS (10) +#include "lasp_slm.h" +#include "lasp_assert.h" +#include "lasp_tracer.h" + +typedef struct Slm { + Sosfilterbank *prefilter; /// Pre-filter, A, or C. If NULL, not used. + Sosfilterbank *bandpass; /// Filterbank. If NULL, not used + Sosfilterbank **splowpass; /// Used for time-weighting of the squared signal + d ref_level; /// Reference value for computing decibels + us downsampling_fac; /// Every x'th sample is returned. + us cur_offset; /// Storage for offset point in input arrays + vd Leq; /// Storage for the computed equivalent levels so far. + +} Slm; + +Slm *Slm_create(Sosfilterbank *prefilter, Sosfilterbank *bandpass, const d fs, + const d tau, const d ref_level, us *downsampling_fac) { + fsTRACE(15); + assertvalidptr(downsampling_fac); + + Slm *slm = NULL; + if (tau <= 0) { + WARN("Invalid time constant"); + return NULL; + } else if (ref_level <= 0) { + WARN("Invalid reference level"); + return NULL; + } else if (fs <= 0) { + WARN("Invalid sampling frequency"); + return NULL; + } + + slm = (Slm *)a_malloc(sizeof(Slm)); + slm->ref_level = ref_level; + slm->prefilter = prefilter; + slm->bandpass = bandpass; + + /// Compute the downsampling factor. This one is chosen based on the + /// lowpass filter. Which has a -3 dB point of f = 1/(tau*2*pi). See LASP + /// documentation for the computation of its minus 20 dB point. We set the + /// reduction in its 'sampling frequency' such that its noise is at a level + /// of 20 dB less than its 'signal'. + const d fs_slm = 1 / (2 * number_pi * tau) * (1 - 0.01) / 0.01; + slm->downsampling_fac = (us)(fs / fs_slm); + slm->cur_offset = 0; + *downsampling_fac = slm->downsampling_fac; + + /// Create the single pole lowpass + us filterbank_size; + if (bandpass) { + filterbank_size = Sosfilterbank_getFilterbankSize(bandpass); + } else { + filterbank_size = 1; + } + + vd lowpass_sos = vd_alloc(6); + d b0 = 1.0 / (1 + 2 * tau * fs); + *getvdval(&lowpass_sos, 0) = b0; + *getvdval(&lowpass_sos, 1) = b0; + *getvdval(&lowpass_sos, 2) = 0; + *getvdval(&lowpass_sos, 3) = 1; + *getvdval(&lowpass_sos, 4) = (1 - 2 * tau * fs) * b0; + *getvdval(&lowpass_sos, 5) = 0; + + slm->splowpass = a_malloc(filterbank_size * sizeof(Sosfilterbank *)); + for (us ch = 0; ch < filterbank_size; ch++) { + /// Allocate a filterbank with one channel and one section. + slm->splowpass[ch] = Sosfilterbank_create(1, 1); + Sosfilterbank_setFilter(slm->splowpass[ch], 0, lowpass_sos); + } + vd_free(&lowpass_sos); + + feTRACE(15); + return slm; +} + +dmat Slm_run(Slm *slm, vd *input_data) { + fsTRACE(15); + assertvalidptr(slm); + assert_vx(input_data); + + /// First step: run the input data through the pre-filter + vd prefiltered; + if (slm->prefilter) + prefiltered = Sosfilterbank_filter(slm->prefilter, input_data); + else { + prefiltered = dmat_foreign(input_data); + } + dmat bandpassed; + if (slm->bandpass) { + bandpassed = Sosfilterbank_filter(slm->bandpass, &prefiltered); + } else { + bandpassed = dmat_foreign(&prefiltered); + } + us filterbank_size = bandpassed.n_cols; + + /// Next step: square all values. We do this in-place. Then we filter for + /// each channel. + d ref_level = slm->ref_level; + d *tmp; + + /// Pre-calculate the size of the output data + us downsampling_fac = slm->downsampling_fac; + us samples_bandpassed = bandpassed.n_rows; + iVARTRACE(15, samples_bandpassed); + us cur_offset = slm->cur_offset; + + int nsamples_output = (samples_bandpassed - cur_offset) / downsampling_fac; + while(nsamples_output*downsampling_fac + cur_offset < samples_bandpassed) + nsamples_output++; + if(nsamples_output < 0) nsamples_output = 0; + + iVARTRACE(15, nsamples_output); + iVARTRACE(15, cur_offset); + + dmat levels = dmat_alloc(nsamples_output, filterbank_size); + + for (us ch = 0; ch < bandpassed.n_cols; ch++) { + iVARTRACE(15, ch); + /// Inplace squaring of the signal + for (us sample = 0; sample < bandpassed.n_rows; sample++) { + tmp = getdmatval(&bandpassed, sample, ch); + *tmp = *tmp * *tmp; + } + + // Now that all data for the channel is squared, we can run it through + // the low-pass filter + vd chan = dmat_column(&bandpassed, ch); + + /// Apply single-pole lowpass filter for current filterbank channel + vd power_filtered = Sosfilterbank_filter(slm->splowpass[ch], &chan); + dbgassert(chan.n_rows == power_filtered.n_rows, "BUG"); + + /// Output resulting levels at a lower interval + us i = 0; + d level; + while (cur_offset < samples_bandpassed) { + iVARTRACE(10, i); + iVARTRACE(10, cur_offset); + level = 10 * d_log10(*getvdval(&power_filtered, cur_offset) / ref_level / + ref_level); + + *getdmatval(&levels, i++, ch) = level; + cur_offset = cur_offset + downsampling_fac; + } + iVARTRACE(15, cur_offset); + iVARTRACE(15, i); + dbgassert(i == nsamples_output, "BUG"); + slm->cur_offset = cur_offset - samples_bandpassed; + + vd_free(&chan); + vd_free(&power_filtered); + } + + vd_free(&prefiltered); + dmat_free(&bandpassed); + feTRACE(15); + return levels; +} + +void Slm_free(Slm *slm) { + fsTRACE(15); + assertvalidptr(slm); + if (slm->prefilter) { + Sosfilterbank_free(slm->prefilter); + } + assertvalidptr(slm->splowpass); + + us filterbank_size; + if (slm->bandpass) { + filterbank_size = Sosfilterbank_getFilterbankSize(slm->bandpass); + Sosfilterbank_free(slm->bandpass); + } else { + filterbank_size = 1; + } + for (us ch = 0; ch < filterbank_size; ch++) { + Sosfilterbank_free(slm->splowpass[ch]); + } + a_free(slm->splowpass); + + a_free(slm); + + feTRACE(15); +} diff --git a/lasp/c/lasp_slm.h b/lasp/c/lasp_slm.h index 7df37ba..67970ca 100644 --- a/lasp/c/lasp_slm.h +++ b/lasp/c/lasp_slm.h @@ -2,60 +2,71 @@ // // Author: J.A. de Jong - ASCEE // -// Description: Implementation of a real time Sound Level Meter, based on -// given slm. +// Description: Multi-purpose implementation of a real time Sound Level Meter, +// can be used for full-signal filtering, (fractional) octave band filtering, +// etc. // ////////////////////////////////////////////////////////////////////// #pragma once -#ifndef LASP_slm_H -#define LASP_slm_H +#ifndef LASP_SLM_H +#define LASP_SLM_H #include "lasp_types.h" #include "lasp_mat.h" #include "lasp_sosfilterbank.h" typedef struct Slm Slm; +#define TAU_FAST (1.0/8.0) +#define TAU_SLOW (1.0) +#define TAU_IMPULSE (35e-3) + /** - * Initializes a Sound level meter. + * Initializes a Sound level meter. NOTE: Sound level meter takes over + * ownership of pointers to the filterbanks for prefiltering and band-pass + * filtering! After passing them to the constructor of the Slm, they should not + * be touched anymore. * - * @param fb: Sosfiterbank handle, used to pre-filter data - * @param T: Single pole low pass filter time constant to use. + * @param[in] weighting: Sosfiterbank handle, used to pre-filter data. This is + * in most cases an A-weighting filter, or C-weighting. If NULL, no + * pre-filtering is done. That can be the case in situations of Z-weighting. + * @param[in] fb: Sosfiterbank handle, bandpass filters. + * @param[in] fs: Sampling frequency in [Hz], used for computing the + * downsampling factor and size of output arrays. + * @param[in] tau: Time constant of single pole low pass filter for squared data. + * Three pre-defined values can be used: TAU_FAST, for fast filtering, + * TAU_SLOW, for slow filtering, TAU_IMPULSE for impulse filtering + * @param[in] ref_level: Reference level when computing dB's. I.e. P_REF_AIR for + * sound pressure levels in air + * @param[out] downsampling_fac: Here, the used downsampling factor is stored + * which is used on returning data. If the value is for example 10, the + * 'sampling' frequency of output data from `Slm_run` is 4800 is fs is set to + * 48000. This downsampling factor is a function of the used time weighting. + * + * @return Slm: Handle of the sound level meter, NULL on error. */ -slm* slm_create(Sosfilterbank* fb); +Slm* Slm_create(Sosfilterbank* weighting,Sosfilterbank* bandpass, + const d fs, + const d tau, + const d ref_level, + us* downsampling_fac); /** - * Initialize the filter coeficients in the slm + * Run the sound level meter on a piece of time data. * - * @param fb: slm handle - * @param filter_no: Filter number in the bank - * @param coefss: Array of filter coefficients. Should have a length of - * nsections x 6, for each of the sections, it contains (b0, b1, b2, a0, - * a1, a2), where a are the numerator coefficients and b are the denominator - * coefficients. + * @param[in] slm: Slm handle + * @param[in] input_data: Vector of input data samples. * -*/ -void slm_setFilter(slm* fb,const us filter_no, - const vd coefs); - -/** - * Filters x using h, returns y - * - * @param x Input time sequence block. Should have at least one sample. - - * @return Filtered output in an allocated array. The number of - * columns in this array equals the number of filters in the - * slm. The number of output samples is equal to the number of - * input samples in x. + * @return Output result of Sound Level values in [dB], for each bank in the filter + * bank. */ -dmat slm_filter(slm* fb, - const vd* x); - +dmat Slm_run(Slm* slm, + vd* input_data); /** - * Cleans up an existing filter bank. + * Cleans up an existing Slm * * @param f slm handle */ -void slm_free(slm* f); +void Slm_free(Slm* f); -#endif // LASP_slm_H +#endif // LASP_SLM_H ////////////////////////////////////////////////////////////////////// diff --git a/lasp/lasp_slm.py b/lasp/lasp_slm.py index 930e265..8ac595c 100644 --- a/lasp/lasp_slm.py +++ b/lasp/lasp_slm.py @@ -4,7 +4,7 @@ Sound level meter implementation @author: J.A. de Jong - ASCEE """ -from .wrappers import SPLowpass +from .wrappers import Slm as pyxSlm import numpy as np from .lasp_common import (TimeWeighting, P_REF) @@ -24,11 +24,11 @@ class SLM: """ Multi-channel Sound Level Meter. Input data: time data with a certain sampling frequency. Output: time-weighted (fast/slow) sound pressure - levels in dB(A/C/Z). + levels in dB(A/C/Z). Possibly in octave bands. """ - def __init__(self, fs, tw=TimeWeighting.default): + def __init__(self, fs, tw=TimeWeighting.default, ): """ Initialize a sound level meter object. Number of channels comes from the weighcal object. diff --git a/lasp/wrappers.pyx b/lasp/wrappers.pyx index e05a2ef..52d2103 100644 --- a/lasp/wrappers.pyx +++ b/lasp/wrappers.pyx @@ -455,6 +455,76 @@ cdef extern from "lasp_decimation.h": dmat Decimator_decimate(c_Decimator* dec,const dmat* samples) nogil void Decimator_free(c_Decimator* dec) nogil + +cdef extern from "lasp_slm.h": + ctypedef struct c_Slm "Slm" + d TAU_FAST, TAU_SLOW, TAU_IMPULSE + c_Slm* Slm_create(c_Sosfilterbank* prefilter, + c_Sosfilterbank* bandpass, + d fs, d tau, d ref_level, + us* downsampling_fac) + dmat Slm_run(c_Slm* slm,vd* input_data) + void Slm_free(c_Slm* slm) + + +tau_fast = TAU_FAST +tau_slow = TAU_SLOW +tau_impulse = TAU_IMPULSE + +cdef class Slm: + cdef: + c_Slm* c_slm + us downsampling_fac + + def __cinit__(self, d[:, ::1] sos_prefilter, + d[:, ::1] sos_bandpass, + d fs, d tau, d ref_level): + + cdef: + us prefilter_nsections = sos_prefilter.shape[1] // 6 + us bandpass_nsections = sos_bandpass.shape[1] // 6 + us bandpass_nchannels = sos_bandpass.shape[0] + c_Sosfilterbank* prefilter = NULL + c_Sosfilterbank* bandpass = NULL + + if sos_prefilter is not None: + if sos_prefilter.shape[0] != 1: + raise ValueError('Prefilter should have only one channel') + prefilter = Sosfilterbank_create(1,prefilter_nsections) + if prefilter is NULL: + raise RuntimeError('Error creating pre-filter') + + if sos_bandpass is not None: + bandpass = Sosfilterbank_create(bandpass_nchannels, + bandpass_nsections) + if bandpass == NULL: + if prefilter: + Sosfilterbank_free(prefilter) + raise RuntimeError('Error creating bandpass filter') + + self.c_slm = Slm_create(prefilter, bandpass, + fs, tau, ref_level, + &self.downsampling_fac) + if self.c_slm is NULL: + Sosfilterbank_free(prefilter) + Sosfilterbank_free(bandpass) + raise RuntimeError('Error creating sound level meter') + + def run(self, d[::1] data): + cdef vd data_vd = dmat_foreign_data(data.shape[0], 1, + &data[0], False) + cdef dmat res = Slm_run(self.c_slm, &data_vd) + result = dmat_to_ndarray(&res,True) + return result + + + def __dealloc__(self): + if self.c_slm: + Slm_free(self.c_slm) + + + + cdef class Decimator: cdef: c_Decimator* dec