From d652669bcd85f93fe8490975469938968d44b595 Mon Sep 17 00:00:00 2001 From: Moliner Eloi Date: Mon, 30 Aug 2021 18:42:12 +0300 Subject: [PATCH] uploading more shit --- PEMOQ/PEMOQ.m | 5 + PEMOQ/audioqual.p | Bin 0 -> 7564 bytes PEMOQ/audioqual_hi.p | Bin 0 -> 8257 bytes PEMOQ/pemo1_SL.mexa64 | Bin 0 -> 13080 bytes PEMOQ/pemo2.mexa64 | Bin 0 -> 13592 bytes PEMOQ/pemo3.mexa64 | Bin 0 -> 12896 bytes PEMOQ/pemo4.mexa64 | Bin 0 -> 13200 bytes PEMOQ/toeplitzC.c | 128 +++++++++++++++ PQevalAudio/CB/PQCB.m | 99 +++++++++++ PQevalAudio/CB/PQDFTFrame.m | 60 +++++++ PQevalAudio/CB/PQeval.m | 112 +++++++++++++ PQevalAudio/CB/PQgroupCB.m | 63 +++++++ PQevalAudio/CB/PQspreadCB.m | 66 ++++++++ PQevalAudio/MOV/PQavgMOVB.m | 263 ++++++++++++++++++++++++++++++ PQevalAudio/MOV/PQframeMOV.m | 65 ++++++++ PQevalAudio/MOV/PQloudTest.m | 29 ++++ PQevalAudio/MOV/PQmovBW.m | 45 +++++ PQevalAudio/MOV/PQmovEHS.m | 135 +++++++++++++++ PQevalAudio/MOV/PQmovModDiffB.m | 43 +++++ PQevalAudio/MOV/PQmovNLoudB.m | 33 ++++ PQevalAudio/MOV/PQmovNMRB.m | 36 ++++ PQevalAudio/MOV/PQmovPD.m | 42 +++++ PQevalAudio/MOV/PQprtMOV.m | 35 ++++ PQevalAudio/MOV/PQprtMOVCi.m | 40 +++++ PQevalAudio/Misc/PQHannWin.m | 10 ++ PQevalAudio/Misc/PQIntNoise.m | 10 ++ PQevalAudio/Misc/PQRFFT.m | 33 ++++ PQevalAudio/Misc/PQRFFTMSq.m | 13 ++ PQevalAudio/Misc/PQWOME.m | 13 ++ PQevalAudio/Misc/PQdataBoundary.m | 106 ++++++++++++ PQevalAudio/Misc/PQgetData.m | 59 +++++++ PQevalAudio/Misc/PQinitFMem.m | 13 ++ PQevalAudio/Misc/PQtConst.m | 13 ++ PQevalAudio/Misc/PQwavFilePar.m | 32 ++++ PQevalAudio/PEAQTest.m | 9 + PQevalAudio/PQevalAudio.m | 177 ++++++++++++++++++++ PQevalAudio/PQnNet.m | 100 ++++++++++++ PQevalAudio/Patt/PQadapt.m | 107 ++++++++++++ PQevalAudio/Patt/PQloud.m | 53 ++++++ PQevalAudio/Patt/PQmodPatt.m | 34 ++++ conf/conf.yaml | 112 +++++++++++++ conf/dset/dataset.yaml | 13 ++ experiments/trained_model/noti | 0 43 files changed, 2206 insertions(+) create mode 100644 PEMOQ/PEMOQ.m create mode 100644 PEMOQ/audioqual.p create mode 100644 PEMOQ/audioqual_hi.p create mode 100644 PEMOQ/pemo1_SL.mexa64 create mode 100644 PEMOQ/pemo2.mexa64 create mode 100644 PEMOQ/pemo3.mexa64 create mode 100644 PEMOQ/pemo4.mexa64 create mode 100644 PEMOQ/toeplitzC.c create mode 100644 PQevalAudio/CB/PQCB.m create mode 100644 PQevalAudio/CB/PQDFTFrame.m create mode 100644 PQevalAudio/CB/PQeval.m create mode 100644 PQevalAudio/CB/PQgroupCB.m create mode 100644 PQevalAudio/CB/PQspreadCB.m create mode 100644 PQevalAudio/MOV/PQavgMOVB.m create mode 100644 PQevalAudio/MOV/PQframeMOV.m create mode 100644 PQevalAudio/MOV/PQloudTest.m create mode 100644 PQevalAudio/MOV/PQmovBW.m create mode 100644 PQevalAudio/MOV/PQmovEHS.m create mode 100644 PQevalAudio/MOV/PQmovModDiffB.m create mode 100644 PQevalAudio/MOV/PQmovNLoudB.m create mode 100644 PQevalAudio/MOV/PQmovNMRB.m create mode 100644 PQevalAudio/MOV/PQmovPD.m create mode 100644 PQevalAudio/MOV/PQprtMOV.m create mode 100644 PQevalAudio/MOV/PQprtMOVCi.m create mode 100644 PQevalAudio/Misc/PQHannWin.m create mode 100644 PQevalAudio/Misc/PQIntNoise.m create mode 100644 PQevalAudio/Misc/PQRFFT.m create mode 100644 PQevalAudio/Misc/PQRFFTMSq.m create mode 100644 PQevalAudio/Misc/PQWOME.m create mode 100644 PQevalAudio/Misc/PQdataBoundary.m create mode 100644 PQevalAudio/Misc/PQgetData.m create mode 100644 PQevalAudio/Misc/PQinitFMem.m create mode 100644 PQevalAudio/Misc/PQtConst.m create mode 100644 PQevalAudio/Misc/PQwavFilePar.m create mode 100644 PQevalAudio/PEAQTest.m create mode 100644 PQevalAudio/PQevalAudio.m create mode 100644 PQevalAudio/PQnNet.m create mode 100644 PQevalAudio/Patt/PQadapt.m create mode 100644 PQevalAudio/Patt/PQloud.m create mode 100644 PQevalAudio/Patt/PQmodPatt.m create mode 100644 conf/conf.yaml create mode 100644 conf/dset/dataset.yaml create mode 100644 experiments/trained_model/noti diff --git a/PEMOQ/PEMOQ.m b/PEMOQ/PEMOQ.m new file mode 100644 index 0000000..f9a45ca --- /dev/null +++ b/PEMOQ/PEMOQ.m @@ -0,0 +1,5 @@ +function [ODG]=PEMOQ(ref,test) +[refa, fs]=audioread(ref); +[testa, fs]=audioread(test); +[PSM, PSMt, ODG, PSM_inst] = audioqual(refa, testa, fs); + diff --git a/PEMOQ/audioqual.p b/PEMOQ/audioqual.p new file mode 100644 index 0000000000000000000000000000000000000000..d0709ec7bfa67986cfc64049bea1f584f973f455 GIT binary patch literal 7564 zcmV;79dqJ#FfcAKFm^C7E-)|v2H+fn6JNCe04X^D03B=q07inPocNA} zN0*|$`b&vwLKgCy)LfY{+H`d^lZ-zONb z1)?;_>V?_?6+FC{85wvUHLL+i1?Yof4?`lIf&z^-rqs4oM069@^oa%Y~Mch=~n9Z*iTc6ZKjH4NoCQV`wNHi6or z9`bj#D9UyyLM~O}{gd*=Bu2kC5Co{W#CA0@&~R&5oX1B!jCe?m5%Q5N6`biGx{p-V zugENe5D3{O2y8)|eCT)|J8CLy#K}6wiJg12<9iX*7i%U>*VcWtG?Q)Gj^J=HD@U0L zN>y@&E5C!kN@Fgwaob8%i|P&=?)e@$1@CWPh`g|eOG26qVX$cGkv_&r^@?xkRd~oO znCC0F+;p%|*a}S}BaqkzWt) zr+*uS-OPb9@)>+r19=;h&Pnqni|Tmw6}o1;Nl(GvQ7fxPc6ZomqVabma9OXPV>v#w zrlaOFOaKGUqInzd>37uYlu`zWxdM^4Dpu_G??Vk7z4%KJID%DA6crAxFucdf=e&K{ zHDOWp5?Nt8w!D60GP5zANv~FEfs1~>R5;M=RW_{=7Mxi=$Un1~wRG3v= z58icPVfXh05xRWnZo9o?H%!^fKA(fbAHS%^EX$MdbU?&UEfd9#KKq38P<)ee3xcL$iS9 z6u|veRlr4SMXlu`z1Tv6WID8)? zJllO?hagK24$UJm|J22uyCrHH$adsr7Df|A;9PvsHeOn+aP?;V;69OLBkUpm&V2lH zuLO|^FZ1Ozh^Os=RRjgXYjeqQseyN-S=7e+G=vaU`AF~p1golJLOBt_80m5wXCfVc zQ;>y{BI~U_=ABo`=LubrdCict_E>op7(~w{xI-;MPG&cpEA21lm~ z9r&x?Y391%Id&0p`s%K^o)x(3jh|QzF;9#(62}vVrFBN)kpYWtMpAI_F@Nr z7yaNSPI+p|xaH3wK8h}iu)4C#0n#I`k`{)o>dc3CNR9U!u~|iCe>7w)7s8ImJH4v&;&e>fuUus5`Z|iYju69R+%xGjZ_rHO#UZ&4m zOGuX=L3o`>)?wLFlDWwBfp$%rv_u+_F%V~z`Os)?xtMS#Ud0^um<#fy_^1%TYF(f$ zS+2)Dmz5~p2i%B8PqDuDCga%$s@wyCV%QOp{)k}LDm-?8UP2CMzvIvder50R#1)}c z#bSEc^uc&Tz}2yE(e2pc+U|&(ry$3Hl=bO0-z@xw*WU3kgmSLLiMaO7m6@asAe9#4 zcfd#<%8@o!)zO-M;xwGh!$D)jwV^3zpq=47v1a{zuL{m5NPYdum2`bN_!sikjcchm`k25@nZ1ZpJE?ddnmX;K7 zVAS<0p&5Env@7j%(W@v9pb2rsZl1G^HuJ|IDK?m)o)k1IjqnV1QYdvF0Bu;hij!dl zQG@esFA~ZC^Jxtv?)u=jzUS_x#jw-7Y0#7Iw#xId=Y(@+aAbF!1(BL$7J4?bY$p!U z#6o08c^oKqEf}R)G)8Z?bc6)K)_yY?{C}zaU2!tLVI;}uroN_Te+81B@WN;#O-n-$ z+k#?IN{q10$_B#U=df}sX!`J{Chafz$mI4303N)v#Fy2eTjS3>uV*XURIwM2pZhv) zxUB7R*`;(b7F068?MKZhNYT=I-TE8qZ_#X~w<`ac8^w#>eUCMifm)afv-kWiR%YW; zv^-30JVnjKwfJRnFdaXiZ|Qql*|+L@&mREU8zWb{YB3$ z?*(2UL9UB(%9&2$`4jzoo`l;^zJ`xewka4CSM&aUW<*Yr$2o832%ny983hrJz#QFv zUGPD`O4CCKi~xETMpAoiF2dfTE91yijKsk$?Q>~0J~G7MacjomK1s5gO9#(j5tHnK z;>nsOK+(*5PYGCaxwhS}OvqP{?hG~|Z{uLzPy;r)@TJ0tDnZT%$ReqfC0?2xbg@@Q zZ?3|IHlT636v6Dvwc)XQdK&yE-Q%je9Y+ic)UbG6_802H=lcp=S^nshp`TUTEuDb= zKRCwl^73y?fMlpQcfg{MuA``NG(9u?4&4rwTY0=x*%ycDhY|IYAsxNI_OY{P#vBZf zS}jvz{>GbW`FpX^JGzn*e)IR}Pyf;i2OhK;LFV_9tO1KD=t*I3{w4MoqTRkK zVPNk}SOI+|ac3?6akyHnWN0~oJg4KH+a!_e zMzMYAoevVlOCAF6_W1knt)4H?^!deXZ|ik3h0!NtU(Qyu3Lz43F2-7d@`M5VH@M?l z0d(Q%4vye}L7k%6D$+`b+R}_4J>Bh!@M-rRGDDYK>!qrP5hKYOaK~cavizVpB+g2* zDa4kS{|uMEmo2WFPn>8jGf34Y)T7bBifCog@msCeTq%g)vWXq5W{d=uI-3M7c`cz9 z9kQRwYn)()CYWYCf*U$8J~Vl4qc=dceZ7|_CGv3f&_ooZs1JNh+?-Ca*9BNnC63n> zXD?NHotPifdN#0OG|0#tVmYQ8zq%t2XxavWd4?mm*cSf70!nU$IQ|8IpU4%rkJpe+ zJ*DoW!_*0&pWX%yX$)4o243~pH~Y624E_i8{H<_9GT-;deu?}^A=hL>{m0pdF8~`T z%%8Ny9ntYL6W>Tv@s-OqlWTIIF?qm@4!Ui&5&MxE6h=vr0DJENZs2ahU;|R8A7Ad` zi+c(E+B1aG=ora3HKkXg;C(M?8`1a%i*R-A|!Ow1D&d5CDRy>Yx2->667m z*8u_17GpNP+pY8s5@GYX1+w_tP8lyDmrM&jDK5r22aDInj!lB1h`0qk7EmWynSPVN$712?jxUg6}8(Ttxlv5voa0Iz(^Cg<9}a6ZjB zY{7s1=S_a2HEvL~vC3wwY`QI7&qOe?oH)A&*y*m6$V9ECwJ19{to$j|X0T9q-2;sae$JQ63TbKFf|&=@cf zd)FbwRW}+X{79c92VwHGPY~@M+PFX#iqy>Z`h&YZ_@Rs9B<&lCi+pXYTjO?v;@*dN z4((c;KTMY$GGoV5`+3yMhva~9Y8}6+)qGYxlekuYhJNFkIo+d12zlsqUB&bH z+>dfS6zb6FI5u zzyE8q5!g-!>sn!ht6xdlp3qi@zo)NhYQ2*QVvAvdn@k?57&f^?>jKV`JS&7bt-sDB2K1vvH??ir!FH`2w6ir&7#Gr zu(=(@Nsi-{(W2%CCHA0O5ey|WJuU))_gPT<43&?Dtpimcc9R0#fw&fL`${@Z%$2Ki z?6ks5kZw~=!jas4NQhO9^8pg`?{Nh>{QL5f+rpG<=i z3|D=Mps4d=TGGELGeF)Y67C5H5*<1f;)|V9|D#DQ>w%1>JuN`tkVuk9Lx|9+x}P!6 zx+t20*jEEc>38lwR$%_7LVw?25=bIKTR78GMCnxKsesc>AwMIFcBwg)Ov?|n@La3P z6D?4giIyC;R&n-s16|FhbD5tY6m+dqpbZ_0VpZoG862X%_>clRszC_{~fM`-jw{H41*o$EQEefX104)z2<#(c23y#YXt_Fyj0Xx(5OvA2L`3a~;Ca*Bxye_E zUGU0{(yrA;Xw=+1uLR@DOaCP zer$v^fs!SsI*MshU~*yQ2Hr$e$q8abS!)V3IF`<%a8Q~D*a#w5hDtj2Q3Q~Q|0lBg zmS^XiK~ihcT{Bv7KD9lCg7n!(lL39z`3+)}X1-w=`@{(84uWy1Kd=2Z?yJ<_#&J!u zc*D#htygfhN|Iuqn?FrLafHVR?Yj%0lN-{_B_gI@<6S++vGQM?A-F}{Pa`s|A(d_+dy)+E3qL7fX8TA= zwRBi>rhEbW4pfh_Tv0K7UkYPS+P8K(_$ufK%!*yEk|(-sq6J&$N0*2u^+X-INFL zw%#gQ#MZl;2IoB0;J;>)6Y|{&q;XBRAD&Cmr>#IaFj_cuGkP+u*#+9$mQuv`qUVZX zL+fFYyDPUOoi>g!e0jaa!S*ix@$jYG+_1SrpSYlr-I?GQ@524q1i@9~H4lyV+U2qa z8kJ}yLbfHnEYHg5@4@t>w%X$72>C!uhqw8y-h>mX((*O;Rp%GE3ZY~8dE^;MY}6P( zTlbYSWVLcrO-9y6m0)69Gz)3%?!&Tq-@2tj=QPj51GV9>7A3pqrCg7xU;I+~w>q9@ zl4vCUSwVHfq0RhltE1EH#}^w@w=-y$1BD0(E_zl)B>|g>p0_;_3hmIrgoS$uwggwo!9|k`6Mu0WMHhMuJeDjZ%NAK4i*R9TGq~fE*A>9 zoZnc>Dn2hSmmf0y)PRB})G$Z^AwM|<%Q3crjJTIpp4~zbUqw|5jr-04@b|BAEIR3M zYtuMOj*?KeZo<+)!k~Z#-(UZqaMipyrzxyM!b(i}RAm-9@`{5$E&uLJ;CL66y1G7h z(Zq99ybE}+Vt7_j+H}xUP@xRxD6t;yTv)IGxF9Ly3<}R*!aAcfE3AIW`R#z%>OpAw zQxfFZWbDxybn@%;k0}8Xu1aH|7H8@chCd|qw?=2-y$xbpqzV+9BjC^agI9XF0leIc zDd69nCihsukB~rwJ$S1RbVW!U`whWFaL~{sxaMQ-K~lC>9Sx8R@gqq7gDu6`Lw|kw zAH7H`e4j7Uks>L1C=5V>;j2hQwAS!#7ERpip52#}!$NX)Z9!QOV|ST~Tj4b>px+_d zS7#{1fKPAjv5aOu3CL_UgBqz~>F~I=n5Ob*2|t`;04QL8C2J3JK5vv8&<dU3&^6p;J8o-isFrV)!S~4`(^pz19x_KAaX+>ggrmAzC@ANaip>6T7(H+Jb1Ca; z>?Ue4Nvm&^D%ExZbs())gnym8xzWUfh}=fx|4QAm69k2*n@8*6;}PB%*Z}(n}^p z{>k1y59HUvPpjct_wes3AloDC-YXQl+)iY9ov;=fZ-E~hG;szRw6Gb8$V;6hdcKiC zWhKPJ9sDQyy$F}z!KwX#eSU!vzZ<8zVHsriatH8E>=NCzUSz?vk(M_4I?IT^$7SF} zD`PR*Y{K6tMl%VwL{J)Zh~$L$U%%Dj$d9fQym zvD-?G?v}`+CAbsx4=5OX|BRmWRW-^J_yQ|B`&_S2Z$5Fn9~ySu0yjhbMi*p*^L2a+ z5i{S&a|t2Xnft>2$M=l8dwM>dJ&Twe{{LJr=eQy8Q=Lr~l6Q8ISkZ0W9!dd(tlCL! zlW$43$;n&YZ4%_za7gZKf_P95M4%kUhE(OFPvSgtzR1k1$?z3r(V|6itTls(2-FR* z>Wow2JTcn*6AfyswPtxg!x7%OD$CWieSy~G)k<2GwG3qV*$x{4z{&k1s(^0}-W7I#s_8HJvY||{eOAc=wcHJ{RNV(2(UlYu;xjf zT678amKFP*!7NG74@EPzUXkhH`2d4KOk#eBe@s|02_a~O-ngMjK87yQe66u_6OR>B zw%eR7ss8DQK>LEzftISrO+fTW)~hR{kmL}m3Kx@rQZn_Ol)hd0ChAim0_x#e7Q7YO z?jNZ-ZO7M{7}#+0Uo!Bh2!|CeIK9sNG2#Wg#FG99WQ#wGYQp7J5BV!ysloMri_+!S zsvQtPCvToT@5xOo^`1ZE^T0bwbMLq3?CFbD%b)H6f&y6RmWeDFrjBH587kPBCk#tk z@BtaR<1Q5<9)NiNb>1$4;booS94&v=I7eO|q9wFP<2=_XYn4?h4#KhNHh{28zuQ!|x%r~7!}j+&jsjeEcuep`9*7rffcd&hF+}1%^)Imk%G$V;TvvWi=@bT{w0gf}z zHTJjpNN+=9Fe|>t*s?GrIYYc^AVe4X(Ca?DMW0^LYTn-Gx(e~l?4Vrb> z9C{BQf8iO>=am$%u*wcwCbM&KboO&WPde zK~w{>uiudiVQd`}^3Jv(V2gAU$6~XC!=1aG6RYI_y%^yerjmsuGUj~tq<5aONlWwT zTCJzdTyJ<-6!~s~V58G4gtCkt!r)0tP|i;-#pWk9wxf({8iEaX)We&5UkCi1Lu#1% zpAw9<86`vJH6WaxAVHhvtGmGFDDbuMeb-ujUIgLc1o7FEXV4^yTaFs(7;*Sg7~dvE z`2R2WSun&@^g>pz`uScgK%r_1)6$5p0y`(*J>E3ygU3Uo^RfV4)c49_jeq=z-UFeL zl41)bgm6<;r$B*I^W`{)CagWR_R1gSQ#J&`NIJ_S_|yjjr(bkkafz$@|~NO->zpABgqK#9B=vEXcQ>WdL_5caj z3pCA>b`>VR3vwXCmdR>s^sX!`cp8yW!Y?VA)>mXSjPPpkjPEZbug77il|Os`Ra+}i zqX(wSuR9~}Nc@R`49h>ewN*dlM8Wc9yG*Zck! zO6wk=VWGSELC0a@m8SbBpX@s7jz4+eHZUDOr&Lsx`FwfxQZt{pW3flQ9}(pu=+=n< z?Ja&PF^QrE&K}L}j*ms9qpz6x?ML@_b8x7|h<;aokLu_NrM39(dsZ!uLqnzQal82%8{VDwfKF*Y`fUk=ufUak`ALg6rqo7 zP9I{VI7t*OZtoko0UCga!O03@uCYN52&iEVD?0hHGx+@{o469DtUFB93ao>82);*g zGT9XKd4@JbUuxc$Ck+$Ky67JW7K1W8D6|lMBRQ!rvHe(XylJD-=5qlUnsBY9oN_r<4g+1i|Hn2 zT$A701t8KvT6Wy;oK+4Rni2-Dr!-_7S%Z#Qq9&vSzb4;z7HVjt1_0j}&d(ve%f?X} zhl@65$uS5J^2pH%M9TIibXccmw*%L7X~J;~q*Z?Puyoi4%4{&?x8J!gsLnn>_NM-L zwkf#aBVT!vqo8>bVr7ol*ygoa`KRF#P^#$O0PC8~Pz%qV&M?WwlnWqvN1106q~4gK z#x3~qcN3$4A><@Au;R5YF>5OO>xS<`(qNns3xFC`R5f7do4gnH5hKGsQAkCdiW17k zzum9x=vj78_ZmoW#@~vC9+qu#R_jS^+?`%)DuWrQ0!tt|GhMA6%vK zc9p*Cm75tfP$(IZtX2{%MUhRAhYn@>m9hz_r(L_FE_F^QIgUKxDGA$iqCZHWKolFO zO&e7^n~FjIlFyQ+o)@(G%nNGEIGI}zJES=DW33||F{rm; z?&u+>yKp^T$)-Q6nc~znDvGc6mNZAa=lw6&(xpp>H1x91__~Y`Fp+Ir9opq&W6sv> z`N$K#I_qBacUvJw1&Zt_m|trsX5!VhnMeKASDPUw31ytu7B|=rP{*JazK7&s-;+h= zQj-AH;#4fB_Qo6y3Th^`8axlJ>%nJ29@{?sQSaB`d5Q~nWzUxR7kQ9k~G)javx>N8{ZBM*J4rzbk=&ZcQHu}2#fn&bh z{ojO1uNTE|3ux&2gfE>#4#6UfW(`k0-}%VnNs2&l!%!6t1sn4)?qQ+90;YVdrB3tT zJ>ZmPp23y6OPsrBwnAu5WSCFd)x%#9ADBGwML=7j?GY|%hVc!$7_gkXQ&dg18_$OY zrQ|irtjad1#!A)=_vnDX?niPo6Z;nq$I;#Vv#r2PJ?NsIv73xibQ=cXoq+)X|3Y_f zrX5j`+Rf;K#N?S?3LJrZL-FkKe>V!n>K|cAgZ;!}(c~ zwC-SL+-M9BYxIH2xhzD+pUj!B3M!A3cEjOQY1;cD&qakHzy&gam#cc#Tv(hvRGk6B zHSbJmAM5&L)7k_xgZI!gL`b&rs`7J@w_k98PM?0<5VVv=KA-MiEYGAP2&i`_5Ol9& zY06Q;yr;SZ^Ayl>g(SvaT_vdWS3DGvG$tnC-IU5|>vVpKz_U0A_E*0uyHb{|=OeLR zDDNRj8@9l&6DwABa98Qh*s)3ayqHwsj3P+MtvQeu;;j+=-7!e-4CeAAy)Zp}DvOXR zmT7#7XjaV-9z@g}c~-Tr!h_ASAqk*NC45CNN$0HUAWjjgtGZA9Px|?>fMg}A`BkPm zqS(WacLI0DVi1yMqZ60)^(OM}lYsx}H@UucUH5+ouxZ(a;W?;-VH-i8T6OW6!Z)Z&>wcG&P4+DmiBF7O znc6|Ia7p6-e~qmrCHTmwJ{zijR{J>HK$LFP`rZo<9nuAocX@lX^F5*5aA8ANC1_vg zRKLRrAzr)z_W{hbe8{Qkmz6x#Llr|S)DmowQh{VBT5pyAg-KR1@^7N!2twKfz1xlg zJ>~HbAZlmx8iGS;?<{;|r8$*U$pHF?&sE9awklNj3p!rueH9q*+=^&gE!fl}zNHIW ziHA%`_$a|t`sJ?3>1e|dD2ETqOUt__a(X<%ybX0FEwKV}QV@xbvTfoRl8IFMa*sYqpYVc?ph{JZVs`BYDBiQjc_edeijz&cH5>%NdeF^ zi3T+QU2Lh>=U?`my^E_MxC|l2mz89V#HIVAhjvd@lImHo=V)J@&`V;{b>RV9ZX{uH zr&`6T27)TFmGqB$a;#$_C2l`Q9Pu?qopf~+zCd2(s)W}1C^kcl)r)zMk>*oov2TKe zzBV(r&%;Kjy^%$T%k`coWTFY|tCL6=n*5DSQ@VNu22*WBIJh6)HyH@}m4JHFYw!=j zDuqvRnz|VTg1&@!P9}(34toC+Bdu(BWK^VVw4CxSiA7ib+kBo2I4r?=shW;MqOg~@ zft!9=D&<;1y+c=Otm@NNxLMpMCTwm$X@b*bQI1E7eR3#Ot?G9)lI5rvt%0tx#g&a# zmTY!h!&%Ww&2uWfOLZvf%YUl%e$s`p0~2QuS$3s&i_ZXyb@(i+vMF(_jUz^UfJPm8 z_4A`K!8!_k%?`cio(_&yFRDn(E(jFB%uWUkUwZuQ#UrKgCPf-T1x%Ec>Rq|3M9(oEMAlY;w)E7j|Hp87i$9a+upE-@9JD5vhG!MdDqf%* z5Jq&eA5PS6?eN}sux1bQ>69g$v`(KaglUU@?RP?L4YC4fGp*5`smhd(n_;f{cVkUb zEGW-^Ne=Ueba(Nk>tZcpx7=+S>auybXsZq2zv;hWjvNfZU51C3$-5XGmuLNVPr%o^ z(CR3}H|4c0(6QOBIG6P5@f*B<6aB;hG_mF&f&Spbpl7H;UN9ExPHAG@mQ~bBWyerH z(lWMZ|7A?EO=n7JM|h=T9vZsmhAg@NCLO-0@PwLn@av`HF= zSHyU7;mF;dL{m`g!sIt8mwW4o$)osRz#k;vo`p6POY_Qd|A@EsZ}evS9}+wMvRNxd zQ>h!XM3>)eEuQ&fJrI@o<%%x7-y|k;T9-Z(iTsbs^ld!sMCx zExnl>dp5oqis+T$GqBp4jv)0yq}Sbc=s>Q&C80`LtwpCL894G*e$)nP2m!zQx|>T7 z%pR!dA&B^$?RBg!ylW?Nfao?er$R(c0>(RUq{x>>aLOI~Hkmi30#OSmui1R6)CcE#9* zsN&S0g*#iM-QVpUH3G(aM?b7XISZ)9T{Q_vxojy))Ebmf=ugxS&;Cn?;81u#MrLNmwZJ63)2uTF#ar;{73_2Al{R-O9 zk_aI|!a(p}(mqdlReA-Q2x_ZA@21MXYIv>2K|Z*Aq$aP4Si0O;bYC2T34@MKXRf9_ zXBFY#r2C=-dU%=rb$t+u-G^cxVSwNkFTXHD-lT&srz_@HO4v zpKbo!1=)pl7(wyyhg0@wL%_%=Y837>*26=U4=2SUo$4p+E1^mte#DQIOfD!&`>J1z z8G1MBJ-xe{Lm!Qqq-kI8_GS|4)1`M|$jq;6&5Y!QP-l!p=SYIFq^9c{{LUF4O#{_*m07+@(VId@M&=)F zq$cPjyu&iF@wHsIHMLN9AghUDn-#DSx1*8?`rW}6Bhr)SEiu-Z%GBbca|>ypnhHxmD*!Lu4oqgBsn7^O$)L=bJ)TV1)<$Jh ze}xWtR}G+cpzfP8Osy7_r4re*t%egidtejb*ewy zoAdFS=o4%S165kTv{;<>;BN)whM?>Ji2t7SRes17fKg_5+f*VaX!)q<@P!Mjwl?|g zHXMp^$3om4_m{Dx3dV&sGfQTtUVdeC0-tJO;L`40|0fDewQ=0X@sy*obgiNoe~kNOxp?f8OBzJd}o4J2)N*-ghN_mx8Ry)>T+luMI0IiUnFvm{DI@T2UTIKNe$NRzlxkGVszm zsw;3PXjYl$8$@v!z4Sv=1apF2`hVU+a{RooZ4ZD6G+Qg;XX`|;?ZheP>wl1$$ZKvT zfGjFQh1^=og7)eR>BKyHEas6?=EE>eaRjY7bT&xvV95b^_3U`l|HeX-iE>BTr)g6R zHp&d^vW%DJ*yo5c5p2!wVi8b2*BcUROzE+G$X(fw$EH>yOW(=~09#^I_Wko|>`ZkG z(Va+1P!)rpurb3!?oT#y0%;Hlo7@o%`cDhiLi*IP38Uw2_bYFbr?C5ftFq?-fMd(p z2kTft2xl#;Xt_zWs&1Sokmet9>Yr`iBrzYHPPjgR*6V2;_vv!j3MACwQiXo+(FGYM zS|uG$aBeeWH5<(BYWrFxH|joB8l95r5>-P`&w{^cSmKn$Y#u(w0y*f_^4 z_#<-yz<&D=JpxCU4+#a`=cOJQle+r5vci|$i{L4;m+X-1JR{+#w+d7*S%hx7{!Ptg z+he#y{^uIngF^Er?UsFW4R)hgvbB0cIDfM?&mY|klH4v=nXcjGMV&%WMyjtX(RHCi zR*cA~HeewL@y=WRf!Or;v!8NeMUeBF02t}G-s>f=^MxRs^YC<86^vH3@q!-3+B6t+-U0LfKqhrn=U^loBB;1(Q;(dg^okZJ|26$OgfR|%_PjYPR6hrP`hmcPfl5efU5 zBFX+0JD@LQn;hUq^oY);UWiaJQ+h&7V3N26=?J(87S2JVpvqc2)ME`L@ra_hZDALo z+h6jENYh(4;O~9SI?8Zx(HYm69tOh`tm!Jl@7D516-;P(DM3E zd3fcsP2c%Jt8%x(@e+9$=tpopf=_cw&D}1sffPp)gKZ0;l2jBS5@0e;~?Z!PL67J(kR{xCP zzC0pflLj@>Mwj@)#S%}?RAAM+lFrZKRdMdPTM+mYpjzinc4WWVTzJURac;S?MMK#= z7KI;@ao5}T_4bA+F~+!|1o%`W_qUoh|NB00tX0&#`t9-}nY~^zfJG8KkX_$LT ze>rj$+~4BAoWBvy2;jq2Wp=rX&VV-;{mUb3@iEh`osq9yz>*EdGwM9 z6P6+bE9CYvt2>0_=Ox2oy12^zD(YMsg$r$-soDS>RBykG6CJ9kE@5+xfqmslP9hh6 zo=F6`fij5TxRQlPb;cy=`%jJAh-B{_Zs-2Rr$3e`Fe4?vVZxuvQa2QKKgp}lDd)Q0 z=}C9f8&)xKe-55nmIj<&j6kj&Ay`0#ZYV}Bl|skO20FoIP(R`gm}|JpWU6VFOEb zc|}6ksg-9~>jOwzD0oq8t}FH7w@kEg*L{{yx?@zr2R>AQ__h)a-ZNE^?O%mlcJb0K z1U@=9xD7(^z9R3XST>eS>#QYQiIz-v4X}$@5H%y_*^FoR_1&Cwy@SNhZUn~OS@q5<@_%i8MTE%6rt$q3OK2+Vww)*@nh*5^0#;ldw}#BtLsZxIy+j+P0C?gn)BV zP)!shF&GnCY@`5C;P_nkOoZ#>ZRSr7jYYLB<+xEU41=&fC9JcS@QBbGyeiNN34x@x2^fy&Vs4;V%8IZC@nAVU)MsT@ zwnq~QxahV2V+M5Cbi=Rkqr^Q$mwcZs|636UV*ctAswa}IY^JV*O#4La=|@HOFarQw zjXmJ074Ms4k#{46?Rh~9ZbLE(Y){tZaFt}WPz9cR=*wnls9m0s2T_J1hoTuZ=fsJr z$I|MS(N^0!ztlP@2q!f+Q`>GyC5<3k#z==f!gCd zonHQaeuJrsJiSd8A85e7r&-}iEBj*h;XS#)!c&aAc)(aC&s?$4hFFp?H2)SR58r(> z?_$Jr`0<;MoJWWR;AB%=?a0{5dD0?TGa{dQbae<8EhjWphlg$%=`OLZFCxH^5f5E2 zl(rjppao6SItcF+YipFs2g?Me`xQoeeGSyq_LCCjBGT2|$yyBB3I&whAmV0{;oBWU^R4~Wj9F-x2*6N+Ib!CQ)@h_MRVKuWv!pl6E2)Op%m)5+Y_FO z-)kVDwobs%OjgI`4nvK*3)aW!CP|s-XW<$SAQTVdxsV-KeQ@FVB7g{hlD?Ci-DTE zZa$86Vjr2kw!kiwMamWECeFfkAufyKYm#uWP zbE^Yx5#F&++VTt>1VcXzP5G!O{E<<@XM0qxG})Fu8mF|d-D z$ndA~y^Z3nD_cxrTX3!FH|&2+KMFD$Gv95lJoSj>JDExEA&JzHwQ8I0&hl`i5JcIa ze7@z5wz&B|;?7F9Ky%rEhe5dem;Mwqu54tXc$2Uopg5@ahv;fRrePbc^U4MtJz~uc zI)5*ZY^t?7%Xh$`Kk zQw}l&9JS3GVnq-u0=0RD1!iH8BGmKfxu!6QruZ{o+^@nc52SC@k<-XEC)m1iR4?mJ zNP3AP=f;0uo3Rw%$xk&8p0$kt_IH+3d5?K)C=Oj|dQF!!H7om=9{d;l*@w2$u?;_KA%6Je zQp;%yG5>g7H7aR2Ni6nSG_@`pzC+ljSMez|f=eY9O5&3~&Y@Zlo5?>aUQ z0$eM;N4igl7PE5#48*v%Ctl3Xz-5ezfV~s3e3BsT z$^HIA%t{LHdc&^3&Rl*a$Q8s5l^}+{;3kLZOXM;LU0b?#S*H{Zx8mAJPeJgSpro~D znQBG4dIFItmY&=+VEk7WchggC$ap}`Yt+lZY&-^s! zo(Cpe#@2@XxE|gm-Tg2h-x6489v~O_R`3784Lie@2|0%9cu^B}85Wh5h{Rxia-8^C zmHB2W#h;UdO>}P$t8j=dtheJy_QkNVOcAufn0N&J>+x;inRR#aGh#|my^s22f9p@gp?8OBv2|$S z^7{{u>YtRMsfsTmtaTu~BW_1B-zTvDWuTL+_)>57>Xghx8y0`|SA?{-;F3F*83Z~c z>V4zSg;c$gj#n{j{q8s-(XS?L-}l<5J zYK$KXRE5i{H34_ROm+x1$)+VVNte)!+k;H^PEX}`m=SSpj?UrQFMjyHP-L2Tsh1Y+ literal 0 HcmV?d00001 diff --git a/PEMOQ/pemo1_SL.mexa64 b/PEMOQ/pemo1_SL.mexa64 new file mode 100644 index 0000000000000000000000000000000000000000..da1ff287aed89048b942f4f24d7038c0bc77d766 GIT binary patch literal 13080 zcmeHOZE#%Gc|NO^Y~c^B959%V#J9!_GBF$3SdMJhL2GHH6-cs?C2E4@^-6nJ+H|!m zcJCTl6vsifFcl1>67(@7a z&VV$DMzI*bpAa`o-Ksg!yxim0Du}5Nh6JdLMYz88amj)jGQ%`qq3{_^*^XGUJEiST zX*;Hd)@NF++=xn0y3XmQy6IZ2$CPv_?RwXWQojCUlAT9e6fBss-F>j5dj7hpRx5Vv z^2)_!ziuB>-5lAr>)^wStF5b7{MOgsI=}PMD_8!wb*lThz6W>Q^=8dRDuT>OkK`TABBJ0dwglb7oK`i0%d-@Vu@JAfz>)xf1`x{Q^0iSiloD3@V_d9Pn5yqW$^En z!5=PzQ#=e@?({5xQuY6*GWaWH@Ka^*Uf>2UcX|jwsdCo?H{8nLXKi#Ce;lbflGMup@X)7ES;r5=+aMVnjgYk@Irh7V@6RD)x z6X{QwN>*AD9vF>;WAS7p5g&u%JqN%UD;tX9UoB8 z?tw@m0{^4UX)|J(%0p+wO2_AqSVM3SCBB~c z0VP4%ZIR*O$gW5_GVDjiBAYjf&fp%=*3rJBIlR%o(cfT)p2n;ixi;iLDluj&gh&7V z?WKSA(fGXu*Qam~DG?!!llQ4HOc>o?N@RB%@;pAXH^uQ}uM{DNe5sa`zLI-wIX{o> z;l6~lggB|6yNplj_r+nukaEVEzqnfA^i1R~MT%1nT)QS%NP zy2pX*At|N%9r$vG{*VK&ao{5koZmGh9(CZ(=SOJ#C)J@7Ro}hA5TT=!R#pCeD2K*h zt3Km!ifp0bG*k+=y@p@k>L%ofA0nmcsRGh%FA?t{{(|KHnt190)02{ao_M-7(=SN= z8RF@dOn*!APZLkKV)_Zme}#Cu1=EKm|NF#K=%+^{{}}NU^6CAOA0wVZJ-tivY2vA? zO*cvYbHr1~r|Tu(PdtTs+K~K9*Bc`CMB%o7kTf*0{XQs!CaQ?^hStr7a_@%5XKN0u zM7WMYhr)GDUtucdTiqf=-x-QZ%mA?a*|L7Ma+FSLk)PJDr&o%v~t%{yzT=uQ#oWf6UKK~6=0lDDpWnXR8 zF>;a9$PAF+RDMje(68;Ajrsf+#}=QW{$<@u3eeh1#t#JZ`PH&s%GaaN-}~rH%vbk@ z5>`mz`d;X+pY)v0T^OJAjGwMJI_YaUeMy{`%?R}Drs}olKB5|4kTuM`pZ`Gnr7BeY z8~S4YU5#9%DqP6FE!kQVoxndeR$+Az4B6fwOsHl@WmxtF^AL2B;M~!FIn+oW2_okE z&&kWaJ8((yo>wmid@Cwn3!!>j&V2rdmu`Rw8UygR*S<{czEt%x+2t-`cp47|_NqI8 zjuyqJk}mUSKPnUe{EiKLM*%x*;8_J6wSliHz_5YG6yUXihZOL%4eVFIc^lZNfEgS3 zi~?@J@KzPOQvsW8;1&h+*}zf-Jc%cUGWqboWm__~#ub;bER_4+c){yCGK0`5hd;p0nqCI$(1~sN1ZiB&sm|YlOndIE%>R+hzb^BSDeHBTSt*(9Et0wTA(dhIf?KY% zIJBoD_kMfswZKHM&^eJfTMbX%{3D2u95ib~<7cYd@zUzpaw+?xOP>s#*ph$bzo;3w zMceW-GK)L6E&mp>M<#toK0$rqXBeD%ECu=k`})qH`SN|}c>Njj%g)57Z)PjKOpM-4 zWYFAdJlNIU913)G1v_@HH=2W8J;Cl^_YTg6KHFz(H`;_~ep`ZFyVuFjAfq{w7|15@X__%IsbSL?p)W~Uo?#6|EF+ROjihNt&8QJerHzO& zNZ*?lKEjPyJb`-kN0Pt0L|kJSb2Zv(7|HB#znL~tF=JpTl1!S34CSpM(}j>^ zk$Tw7SjIq6X$NIl(uVXYbgVxV9~}Cn2Dn$0U`N4bGqgULiCeN|6pXnb7--doPoZ7v z>(m%XjT}s)5NoT^n~Wd8{nutY6Y$ZWO%ECywrpv|D1udMh>5!kF;7TI~axTsC zd;h>YQTce)_#*GN3U&s)g|rZXg$OJ}U?Bnv5m<=8LIf5f@c%plYT6-qeVwDBKcmI0 z2AKEg4D+ zZ`C0LrFJV^k{uv6JMpob8GpAArWmR8nr&pFKC`Z=>|o_ zlrD&IL-SnENv+3g^IRU5BV_$g>WX8DK=NEJ)(s>-rR7+pERxS?{`a&T(*Y^>>hk`t z2=AM_%fDX-B&O*BO^<2%HBJ9S(-$>8tLfXC{#etWYr04W_IgbZ>-*Q%+`QGO>+R1b zt*p`D-{h~qb8}X58xC#Q?62SC-@H!Y#>V;$o9j2%Z*f<79*-%Qz2EN?da(n={jeDK z3hsZ!c%|TeT8vi-zMsYTBEkKn7+);7e-z`@g8N}HevRERNjl-MLr+llzOqlJnGz3(n2qGs(D|rSpQo4 z`B98tCtmot;e=$~D^}thNGbe!!TXI)BQHG#bbl+xSK$amDL*%g<0Wx^+Rb_KAX`~t zKb@Yd9;|k_>rc9ti23mx15SSI2P`Qg z9mC}*abKPTUaDQ?+WQ^pXN5TK+=K;oM&k410Jx<8H-S?m#5QCHiv(8bJQ*wI8{L1l=VX->q@Eb_D!-;2!vQ_mf*G5aLGhgkey| zPIVf0_n%F`sob}8yX2Y?thNI8l<3dR($Dz$;3e+c^Z>uxmY+ zc^Y_$p(~y%gTGh?|85!ld>QepwD0c00iLkoIC54V}gW-HyA zen?0_?rH3dHR~M>5jLX{DU5LH(pHWLw>;1l=xlF>2nXuI!H}+0sHGba4m@%7+B$ab2y}#ZwzlpL_Jn%^J34}> zG&S;Zg8mJ>kbF z8^*|Sj@|!7p1vb!j)TX7Oy}@Cq{Ni=C7el5!)%DYH?X}!k0F`(dzopS#(DnY{_Jdo zoL+Hk&)?Nd>3K_u%B3vZ^L(ob82M&<{{Cjl^FpY|q?^aA2q8yvFt+FKbf)XICfA?! zn9{jSn&UCg-}6kVeU#kwSE$;1qB$Yk^LIbfCaqY`KI8CBsyshnI!PcB*WVq#MO3!nd@;sUBX>RPc{~|DosW!9M1$ceH{RK>AIn!fcoc25~VLIhhMxM{4F{Sgq zPJ5oeFl9Yle~ufM&;8)buJU<4d2uxa zNOg+XziiJm4~et>_Po+68TUVCz+XKbTSiW-hqT3 z+i1W##kknIU3b^@(4MVJ_gEL(>dC6r0D|aJO||R6k1Ff>F@tzSqKF{O{+|2s=AAb) z-R?ho_K)72cb@zFp69vuxu5U7_r3X^dqsVL&Bj!)vu`ubO%O=66g+s9(qOEbRkE@8 z{Yy5L>lTfOM#T!J3P@4LG!9T6WAN;r##vB}M-(ly3%o*7Wmf>0>;|OWfV5LoBLxIS ziv`OHLFqXxi)y7arJkasOJVf~x}V}l$~ZfUh6T(NRdzSRj{N*ngHtLt$^7!w%dN6} zipt{fvXvqYs;=cLTiDa*k9~6X&_{ne_}rYWRj-sDe0}0u-IN8HUxr83zq1-luZACi zKW4vY!ooZM?nw^h`LVO{Y}dGKVISCwvLVhYxog>T;(0kgI|`XelE^~w;&x2k2TBzpjkoi z@8ZejPbd#w4}1n5Yv_VMCG4||6@MJSMQjo~Av1XzdOUc&?VZtxH?I3)y4TCR4b6>S ze;^iU55@IBthuo^9E}8;eXZeukmZ(m+q!+;U?}1XhqgiS`YqljnL=&Y7mo+xEWSCW zvv?@d7LBv6=vLMl=w2NQMf4yG#sUEr=s#8w-b8JG))+D31h;7!cPY!n*VhB$@`gGKvH|B!hzC{aIqkA1&zM^4St#_Vlo@=2IU%(qfOwxr+ya@{VY$WDe zWBQiOMyMULketStP5#xD>>f|^*I(e-DBGAw0pxpbWQ=4KsGm?EyDG%h{F1mVgc+%j z>4v<3t&wulFAaqzG5FUjF`;ia^woPuy)TC3`)`kYzbJf{#P@0%CCe19_@mPWPIJGx z?8i*Fn(N6vWx}TrAj3%$eu)A>&zSH_O?a8~XDl#PwVjk*6K;N0tTy4}RU)A4Ot`FKF1^)+Pc-Rwm~b?#5xPwHWCLZa+k`tz zc&`aBG2!=^aP$1O%YMT~9oWSbB)_%ZR5gksjpyV&bW5 zruTAw9`V#Q(z`i-CGpfH(z`f6gLvu+>0ZuXMm%+abQk9@BA%K)eJkh36HiT^Ud{Oe z;;E_A)to;&1w1uzx`Okkh^MAaYn(q$JT+0CZ&c8`KHDUVX1pu>NC7y;< z`WWY5BA%KuJ;eF_#8VTd2RVPdRAa$inOSE!_4F@!915QPA|frG*=Ib-4?KNmCU2R7 zW_S`h)C}V`W{w3N)9Wxd9iUndx*XHDBS|z(|B}tze&KDx{sP(WyI}VuPkQ#HmU?W5 zJioi3mm#q`WnvRWV#3~vpvFRYd7f(6X@557cC9u>iR3$O9xWX-xZsIiPC_jc2+()4>NJ?c= zNzfYygAUi=ccww@c}VFgI*tBt%5e*)sRo|5b5aR6b0&@k9cw>>_RG+m^r+B0j0_r& z!Y%g^vcHd`o7_X5{yS2h{shl|D7ET5dFgCK3gzUx8~BDZDCTzx6_Nrn1Y}Jym3NCs z9lFj*F~n2nzory6z~xU-0nUv>@cj+P-w!@T1I*#xzkN`@2+*`fC%S@EPQKdEx8L?GScW2%z|+z%k?X4i zgdpqphie?piHv>=ZI1Fk3_I4p;V8cVg2LSLE;4!z@t5(SrPl4ka}2sArDRbu)8^?L zEJT-fCI*oDd#JgYr#yX!Yzz*`1zjR3F-lqQ!8 z;A#VyB7i-HNudB54Pf{j&*eJ9qjVWt!&GX(GthDM11Xf72ROaSc$ z@ab7z2X`63I{;LR6wRh_ytC4?BZJpKDOW9;WX71GjdsOI+f?*2GEa>IOoQp?U}$jc zlk{;;Hw>z-C``!TQ;u~-M3Iaf=1QwCWBaa#4qNm9#Qn=LnSA=C(AD-OU+6ny zOB|<;?kv(Fm!9{Xa6X5cWvQnz`(OaLzkZHhH>pQCyWuc$aeqyYmGP`MV33e2>yh;) zj6CAk7nMWyA|Es6*$E@Qr15S1vTY}Id}F{GXdG4l^}6W8jE=9hdhP)SzmzD#`!z8Q!?PW0l+^hzGa1pcUt+=-pjr0gD?I&; zgAK`j4auV`lE>?kr&lCDNQN zdS}gARVR0zDDP#=TnC5eq8mD=ZN_j!m-xrQ?`1N5XA2XKhUDj-3$ul*M`hyrUc$@bx`6tv~%0V7N zZ8>N|@{niecZyR_!vwwV7+UCO&%gtK&eS6UM2{}{#ee|IQuhmxI@uGy2FNoxSL3Mt zL*F@DqH-XEMH0sAU~)f7v;9AKaZUZf0j z=|L*rTAs=a6ul17?Sqc(U!fbN_ukM>8V?8dh8_`=-N1zkD}#fFunO;r+#nAEwZK;ofxoQQA>=+)M9b?A#iXl*EMAg1B`V{;-9 zY3o73HFwjpIF2a-Ccl=gwb)UKwubzAhgKf+wdrW5SinV=RJNLgKLn^ptLHUk)USrh zR_VljsO(Ur-JEq$mO}WfvS(UDkO~S9QK@lFk7`%UyR#jHbMjEnOs6e4DF8JQ2yLcDj6}njC z33#Gme|&2+hOz|E4Z<#NUpafW?A#bJ=w4ktu~2)5zF2FCgtlPZgmkoEP+JxD2O_Nr zWJ-hNs(BR^3sHX42tcOD9JPj8^yyT)t?*V`sdZ~}Awxgg0{7wA>Lbxt&dqLoUh+edo|M#CE$qitd;W7<>0{rYQ9Yl+uPUph zeXpb-y|C%MZub+LmfEJ;uk8SWeFJ;{B<2~=1;d%lqo5Cep2_?kbP=YP&p{7k9xBBY zcN%m7=xj{E&7iA5=?v*E&<8+Y1^qSXb`-ON@l_Cuawrx$!*6p-SojA65 z7vwX6(HYh{UQI4~Z;D zegw}`kk1D-(vOVDPXPWcV}z7aW{!N2xL?B_uKA>^yg_60Za1gO4NKz<3@;X1Sb zR1fn;>VfpAU9LoXTx8a};XZqR;k`wDW9$W_JqLOZL0*cr2R|!Ts_?yUTzAwc1B8BzSz7QV{D*IIZwBd0J?5aYX^I#b;v zCpvZ3tj@{hXIlxHk*ybid0nUv#o!g-?{GT;AeqQ0A3^F1Uh;B`qql=QTuV`KweBI)Il zE|Rof(wig=O1eeTyCwaxq>o7Yq@;%=eO=NICEY9AdwFdw7CkMkiHM%i7P=O=Dy~|T z;N0Br&RyiHSm0VTTj1Kfin)s_7FASPi#Uq6Gn>(W%tHEJrTT3)Za4aCHeSe7zs|;s zm}>uQd<;|lDjOfmRDa6Gi<#=T+4wl7`gb-y`g=Xe?Q8pL6Ate^&> zNnmG_nVKg|Fw%FhK{XIf0y}n_>OO@DhJFdFR1=L!V8<6=HSd^U=wEEiH`#b8+dEA& zLp*M0Q*du27k&v--|fssc7|C)_P1=jjH!N@jbF+Ja_XzuaYWmOXkm`;-82u|u)4JR zPt{U@@8H(@r+Jw0YD>G&T90t+dxO@THa6lX`9Wa~9$Sw3+y@eux|{20jBn!wKC!P{w6`xc;J1v&bEF>tEqK~3Y9>izN>I7~@B`MWK7LE#8+jmNToQhD@;^Wd}x%$2^D2frc@K0gm$mj_>$ z2XD)RFVBNV^WZ&s@P6RA>g}OC`acJr>pssfxc=z-JpYy_omcbV@8rQh0G`Xwukzp* zV*!?{9jE5OuL3?JyXy$9EvCoyL@?-T!`+|jYn!|)8rC%9V$#T+B(L7-ZKGRBacKCX z-u7^`)fe{qabyzr`V!qtoc!So&~?qi%7wX==(1AC%TM^cIEIV$;2ILnKfV4$XJ-#g zOq>^*dX}nj9Ep29( z;YL#r*JIuQPFa17d&Y3Wi>0j ztLp35xSPGrHOp4Gy)qr-ABn~~8LuFIjjOuae7EV3-6f;jRq~D&u6|kUEO)*}qzgvw zXpP+b!fmn<4Z2e$Zkc6cl%H`u40poRT`=Po8d?dL&$5jycjK}yrqM;N?7L`rZ}9MgtTPu49_B*S)gDDe z58zUtiv}N*qP)Yq0v&R=b@-8p;)JOfcfvx&Z)<}Zvhj7|Y91v^$pfL$<3itOF3bS* zDIwSV^1sArIYnDrOnW@!yYDQ;SW{&OS|{O1j$RWysQr(k0`zb}$xd2pnz-u}yc)-yfT-mFAm!jk!1?9ieQ?vy# z+T$x;?avfd_E!ISqWJOSp)&zxul8|@s`X1g`#9{CjM7y5Jw*ozL{R=)>+cT8$bK0q zA*g+!L1l%JJ)MnM?Y|F*zZq5B(N< zwI5Yf=}C3Nf5rdUVsGr5*D{Ma^wDec~OYOha|0P)Ut^CWdrLBpw zSLcFzrU68->Rb6Yp+hbzd$k|l`zPAJ3tQ^D%3kgJ)&GrXvSL)dsr*$t+;7QWoevzH z4go^BAjYq4NDPGc8Y!m*le%P8}@41`hpu>pI<^iX&L3$G_E=z z2zIjBLFKE^uMwgAD-YmDi@8V{BkJ>Xxik8UKAEDOEyFK!n z;~jIiCy^@B%Ax#mkeSE`5+p%Tqy~~0{v$B*1;GenI*L`naZ{zW4c9eLqG{5jU8gdY z4>_@Q{k@ra$GzoU5%iz_!y$Lyyx+Vx^XBdB+|17Vn^<3OQy?HX1;w8W;x4r)q*E8n z->x$tonp6Gjr*U8jZ(IGsc2cwZBr0aAq)vn8`t3a+Mh}m)Q}mb^GymbYsz}WlHOUZ zcUJ2$HMBfa-3_q{I&@voUG>w=D5J!bWGM~)!c9{CiI8NM(QOJAOj++f=#ia&Y|?HL zhjn|6{PM8wA5+~O*|!_eSP<92fno9dBa6xQxu1XXz3cuivhJo6`Gfy>-@lM8q<%dv z_J3O^n3D!>3jbE{iM2Z)`SSNA(BKvn*NXD$+CsA5sD&i=72L0@gQpwd|JDGXX@GyO z0sedge4+tP{xER))7Jsiv;TYp{M`olFB{;)zztmfMB~ZugW$i7tDZf*z;CG|--(WG z5bss|MaN}5>%>pjs*(=K$z)k$6S=HaaP7QnS;C4B^;=0N?~J7ju9F|?@6P11&X7Ht zaa39T5-TxfTd8!`&ZJL3@ZMwAVO>La#x4|`f^eoL#e_38m``Wj6yDf zY=|_W{F(cNpfXDId?}IMZOHR`SiCvChGx7TPkE~@CwWyKtd{e<%k$+Nv?)EK=U2v0 zYJAEtWI5x^zp+W-H1~0zNO0bRUrT^aWe={WSXpw(gY(=-1y?+{_dRIW_RzFMb08Wm`PB_xCH}!Kn|Nc6)H|)Ys|3*%zc^dT{MBS=FEi*RP($@Au%X5>dB@J@`5g zKJLNSd+`F|pwT)*%&$-hiIxqRV_lL%~3v=em`ePq~d(K0a+_Pq+GM@@<>Jeh(9C;xX4sAMx zD&ex({c)+kvrFk;CH)t#2F=nX^TqN3GjQJg*;Th4l|7{^YgLsgeSdu?_rDJ@QT&)W zvww_o!>G?~HBaw<6I?llE|g8&M$R{11Q58q+@6K&Vf4%aQCz!l6q3qrw7&-O)R_gU zvT|`m`iJt-;pos242X8K^oMxq*U{2zmqQ!+-?|)n@WNu~fprhy;Vg#QniW0zV!7pN zr2^pAL0RyLc?gcY6?^nu_nzo%tYXggoNGyiw)_Ml*Mq*Z7_#05eGPR{eXlNt?tghP z8rsmhIgasW&YWv8_q#%+r6>QoaGqV(7t(7+J#? zQlbC;S7b(BiOzfwEM5nJw(@0E3IPnJ4bCg5ITq^r6)oCCPw7`Lpj)CGKo#-QYP0nB z<=c?G9BS)-1-ICRp3>W9>AZRR({ttR@Z~#qn5XxY3#cz%dWj4v%Pa);mH$DS9y&!+ z#X>u7)lnBc8ns5w;l;tc=$QtMj;+SHok<I=@XyrYGz%Ea~#o23BOQlxA}UerM6ey&x*(iY^Q@ zlD2E3*1SEfB~cn5SQ1nAFBcP#)BaJyvd8J{n32jmpDH@pgsu>BOYVaxt#G(8a57~Dr|e)^$&;tIwo==kp|l}XSKLB9mri?=|&w#f7PFL2_p2u!sFHny&A zIfK_`Gsf?t_-db*bt}JlI4Ou0g@n_pTzYTUNpIf;PU>@1wXrLw)x3x zW>yC`Nl{3xq!kaWcwof?D;`+!z={V}Jh0+{6%VX6@-FuBK!=uHlTQqit z?+Qn@?!Q-zMUn6*&*5b?0-Kjh%)`~Xm*IM~jF<%#7 zSB>wemB&d^^~GN*35t+ttBbc);~|!$$({}3Y{-kKd=MwsJb!pggNR7rbH-Zy2JxIe zf2_qng7^=^3(0&?+$i|`(OVc4c>DD{QH!?=9*4E~M-dWXcp;e&R%f}TBARYi)j40M z`80s9bbkA39&Lh?IlhK#qWP5Y2v0n^fW8|N?w|kZ+g(5`*{P>JaYW+F{dq#y!@nnB zWqJ(vK%M9Eb>Q{-)u_ME$$B=3N$+7Tupdf%xgP-6t^Z}<(m&b{^4lb^e|;LZ2}ujjW`;Pvb@8sK{y)YAq07Je9Bn!pK)i*uWB zLRi+l-G{Bd_>mz5+$_b~Snha$ zeLN*BhfaHhY>|v9BJo!+<{88t+KKn-Gz z5L3jlMAaZAUuaZK6cWY#)CMs%j6hn7Ig&;U)D4f4aaVO?t2o}ts)@0x?VSi0rYBr` z6x7WtI?lOl&UM0L*pZ1u3{l|fk2=8w>Zk8$Un;{|-wdGBp0pr!@ zbv{#h#_Lhw(lqgant#=H~G3F!(_be^S*}Z54_6A``iB!Fsx3iNdrEU0lG pE9>#SB2x0Kp%muy#6-h^qLsNn{QE)gXTA2PKBn}K`YQLg|36sR&G-NS literal 0 HcmV?d00001 diff --git a/PEMOQ/pemo4.mexa64 b/PEMOQ/pemo4.mexa64 new file mode 100644 index 0000000000000000000000000000000000000000..998db1ab9e918b0fd58c6fe298df50d800b150de GIT binary patch literal 13200 zcmeHOdvIITnLqL)F?n4(aljC0ZrOx&ildc79OAGcl5N=%Yda>fAPI34*^+F5EE(x4 zal%4V%p)Q=r0GmIWtrJb3%i-Nbh2HRO(_c|q;crdOn^?CmX^(Sfl6g}2uUf~0`~sC z^N_EurB3_T{;?d-z2EtM-{XAeanHTy==h87_SO=cjj2@1Zf3*{%@T-H3M%Vl2E@r$ zvzfSmjxFN4HR|=RZs_brs9mV;_rb?-}L6%pn zUAD{mDJrYO>sAdP)LdThYrLd)uf8GOxj<+P~5XW{-}W zz`s+vcXq>|1PjX`jH}d6#_1#|4|YBzb%3f z7Qz3a2>yBz{I84P)UG-%YuW*zQ2hT|1b?dtUS6afCW`3$fa|!di5_R&3WEP3u0nC% zT?B6hel0F*3M0-c_Q9kTe+a;QHkUQY0sn#=&&==dAB@NRNh6Rj{C?*5c6a!Lp+u-Z zk~Bhz?vCbYJQnH>^hQHMR#@Wi8xHuxkys!a*$u^8cKA0)56#g)G8sy;(C`o&itl8} zyAuY#h1pk2(`He$?L^8y)E%6iu z1PAU!q9fVgJ!}AP3mIJrK{xgVq5;Y`Cqe-uB=kE1Mj|roztQ0kN*#i9f_RgisliYp z(uejyQ}o^C?;E(w9}YyKEE?$@9JVLpJQK?0fY?_NHe?K-j3Dt2;`@Xo5DW|%c0?Fh zwUTwXH?g*MZ&R~>g?)v+A(vOP7xh*tpG(9{B6=@31q&7$<(NlA$JvDyR?l##vA zFXj5g=iVE{AUS6XdPKe-)O&zQ{IITb1%)gA>|%k_+-S2Jo{|5m@a!x9=gW2Gd2q}q&s=Xqar~3{4K;&S7pyn09e{YJatj_9Orj0)LHnzD(P&Fa3<`n{`Rc>9@JRs&pjh+e_f575B`ifo_%ZY2cB`G zMi~!oJtxZ#0I*%09#3RG;u%@TQnRz$pef=;`DwVL%4Esql|$t1Kt{g0}6Y3)&Xw zZti$hC{&gol!DCHKAoI|bPR1gfe_%(*>X{v`PRdljyc z6E@e#Nt;W1@fTk0ML4W=7U?z-k z(H!qDAfIodan8ewamaB|yVco6+OtrCvx}Fn)t-Xt38?Z4g}R`nG#XPzFZk+VZU{H8M{$O_~}yxWPggs1TNgNr7IVn!_#{~+xi-sscflqpjoS0%IiEMV`X@rs#96-sJqNN>Z=KB?~lQj zSEril^p{+x@k}ZUZ{6I{$cdUeTw7f|u5B)V&&v15zZ1Zz9I#&if0Hx0R{-D10Z{?`AP3wbfM4W*TLiEOQ@$vr zQ2<+VK&=4$xni#n0KFxI$qWJ9n$!63e|T;8=QQ3Hz|%S44FS-)TlhFGfbZut{!sw$ z;@u@Q4g=`nEo+|9Ho4N6lT$}r=^uO2{}%M@D_dGiV>RtQRX_H^CS+fvDI)XwrO8Qe z`UkGjmL(G-mnwEv<3Y2f<`Mtm0^V=uy|m zCnfjQ0k6zFD{wrKRUJ>z^MhfkJ$_t!_Vr^FyVi0%?LPEg^|<>Gu^5kQKSBRN`(D!i zZ)64!!(Hip^Wm;0(=Xi7Q1u;#cEjsJJM|=#%D9rAX%b2_YxEF4OMCDEIERDU%ynG0 zt^73X4h=w9z&Gtc&mTx(mLLH4lGN*B|7PO#a)h5whP$@#m^tHNUSgulbcGl^w8(ke;?e04Lj`l8pv#)!7YikdfZqEkyM$h{87GheO{2N@I zZEh}b`RXNIA?Xd0u9UQ)XBnbs4TXZeSRm@rczh_S8*x1u8B9g7$kapOaH!8n>hZA7 zSAPLo&)V_YMnZ)U;u@+LrFhC1N*Q`EU(&LOPX^>8GL?m!5- zXSWKE$~E`QjrPV>dhPO?bt5tu();2uES_W3Wh)$m$z^qVGQ|6jH`72oVW51yy=%E} zEqisUQ`YDkyV|M1V3Qs&^np+S+XMYNEwGaV@n{fHp!i(BnbrH^2`CN4V?l0#Ub#CJ zNQCsAfgQL-i~*gukdX?8^m^&1+vP)3cjrYOjwg7}>YWtFU{X~)uiiTX(NrjjS~k?_ z4URg!enn%QzM|m<;`pdUoXzo}U5QBlfU!pR#UeW}t|P`SYVGyWU?|p`O7!dXjg709 zuW&dTPzWEnI(1$I^+1(Te{ZMW8`1*^WmvC=ED+UqV1UIFc0QzV+q3SMPfa+nyV}%R z*V6s*=PyH^(z$hlo_|r$l}W6$dbTECEh?5!>NQr`cp)W$V0}TeR6b zEB3H;SKai58%f@=U}J%u+6bIq%KsHA!BjE=`j_m_>9KRRP%X6aW-;1x0? z(l>DZ5>v?=1bC%>+y2tgvPa8DW|aKP%T2E~{ix}*>m;lwh9`0DK^?YG9Lj&G&qhwk z-w9mLp>BLHOv%mVe|xHYsD73{Yk|*N;IkI^tOfobw1C>jR{Pef%JeFsbgj%Np0+0` zDW3Ms&HOqGPkZo`<_MDeT`#xCy+8=AlybGtzDV-&waBiLyxMQ4_Z}rB|8Qb5PW*XU z9=+!&?JW~j?N!&%I~YlQ@9dHKYEQXFh^O}Eb)6&9jPV5-5cOT4_PVc-c6%kD@Kds2 z73Y}LqYYC^syuAk@c5llAFB(_t8%gA;QX+ZWAh?`e#<*QPhtG};)mEEz@zG|7k@yTDm%CXAPC<6Qg_fkA1>KuSstrWk& zsBxN)mohaz^YOCWIL^n*ndqFVz`S@(+ zEQrrx+hu#_>(9-d1Ifp+{V(T*d|YFd1@X$X{SF*8!`f4_u(zYR6R^az= zns;rOt*r6Wd|QIw$<;X=<%#B9!jD?|i^bjsvfqr=kc%ga~tMAj?PZc|=CS>rGzK^R=JCLl9|8Wvvc0OWh|IjX~a4#^iMU)#cum*dmt2l@_} z&r8a`Z33?6L!2{;=zp;YPHXZ){=Z!W|ED52Z7UY?bEXLXW)b|?Meq-R7pgBBvtYTj zNnI!&9eAPhFH5-o^z$!{B7VF@a9n zr=u&kG;i>?dpCCDAkWkp9=|c@@1s*ZNoWM){{CpZHxTs)v5}PY2U5eVj~}uzLP7hD z4XYapE79Seh@Wp&`a>}zv5SRizt$g24G!*tiHY+=)5ud5XJIJ7)(x%>x8L2_LPw2w zxH2E2tx(2GVXk15bc4ybe;^PG(#fMxEC{4GnG^?=Sc-OiS!bd@%KR;xJ6#>#X2@{> z&+qogLOd-S0O0_WDz2@4eUq!*kDcI+?rwj#tEt`X=Pta6m+`3Sz>qk6^oPzt(QzSp z0tu&@oKEw}9_#U;De20oLq}8Rns6d&N`p=Yk?-6=Cmd)}2bpq5r7Y*Q=r~o_;=+6? zicTKk5R&(1BG|=rB%FK=jZx9btp~w19#pwHl)>fD%cv9cjiY}43B69hsSY}5KTdQ$pufD?- zt(BT;-K5rS%E4-6=#{7J)%UuhG!Iguaw)6q)jH1!jC?D5wLVpJ58e;vWMP>HLR$MN zd-eUV=t^m<;#Yc#`XHyZrQ+2(K~bt7C2RZwWh{sctwoi+T2ClS`|6a6*(Xgh(o^dW zMaKw4Qt?~c?_S8rcM~Zgsr8AXYJE=j^a*ISe;63Gsaz)K)-h^5WBszGvdrmzFlKwT z-ceMybQm3Wp`_?2@@9LrE>cwKaUYcLr^;9BJ4Js>N~UBj|2Smy zsjR|O>p*oL-l}iqPs3I-*{k)+%2(^uvx^}>sud)6SN4i#ATh_!bY195a#(P=yRuXC14zvFYF%_# z7bZ$iIk57y4kAp|e@vEtL>Csrg5>TJ&heihN98+mh{z~=>-ZyGN{Z4MP*pY;$Xmtr vB17lr&B?-YKWgln_0H*n8;~1^G#qZkrRuNVAM|_{D&JEp?8hw5t>ynOgP=P0 literal 0 HcmV?d00001 diff --git a/PEMOQ/toeplitzC.c b/PEMOQ/toeplitzC.c new file mode 100644 index 0000000..b3c3bd4 --- /dev/null +++ b/PEMOQ/toeplitzC.c @@ -0,0 +1,128 @@ +/*================================================================= + * + * toeplitzC.C Sample .MEX file corresponding to toeplitz.m + * Solves simple 3 body orbit problem + * + * The calling syntax is: + * + * [yp] = yprime(t, y) + * TOEPLITZ(C,R) is a non-symmetric Toeplitz matrix having C as its + * first column and R as its first row. + * + * TOEPLITZ(R) is a symmetric Toeplitz matrix for real R. + * For a complex vector R with a real first element, T = toeplitz(r) + * returns the Hermitian Toeplitz matrix formed from R. When the + * first element of R is not real, the resulting matrix is Hermitian + * off the main diagonal, i.e., T_{i,j} = conj(T_{j,i}) for i ~= j. + * + * You may also want to look at the corresponding M-code, yprime.m. + * + * This is a MEX-file for MATLAB. + * Copyright 1984-2006 The MathWorks, Inc. + * + *=================================================================*/ +/* $Revision: 1.10.6.4 $ */ +#include +#include "mex.h" + +/* Input Arguments */ + +#define C_IN prhs[0] +#define R_IN prhs[1] + + +/* Output Arguments */ + +#define T_OUT plhs[0] + + +static void toeplitzC( + double t[], + double c[], + double r[], + int m, int n + ) +{ + int i,j,m0; + + for (j=0;jm?m:j; + for (i=0;i<=m0;i++) + t[j*m+i] = r[j-i]; + for (i=j+1;i 2) { + mexErrMsgTxt("More than 2 input arguments."); + } else if (nrhs == 0) { + mexErrMsgTxt("1 or 2 input arguments required."); + } else if (nrhs == 1) { + mexErrMsgTxt("1 input argument: not implemented (yet), use toeplitz."); + } else if (nlhs > 1) { + mexErrMsgTxt("Too many output arguments."); + } + + if (nrhs == 1) { + mexErrMsgTxt("Not implemented (yet). Please use 2 input arguments or toeplitz.m with one input argument."); + } + + /* Check the dimensions of Y. Y can be 4 X 1 or 1 X 4. */ + + tm = mxGetNumberOfElements(C_IN); + tn = mxGetNumberOfElements(R_IN); + + /* Create a matrix for the return argument */ + if (mxIsComplex(C_IN) || mxIsComplex(R_IN)) + T_OUT = mxCreateDoubleMatrix(tm, tn, mxCOMPLEX); + else + T_OUT = mxCreateDoubleMatrix(tm, tn, mxREAL); + + /* Assign pointers to the various parameters */ + tr = mxGetPr(T_OUT); + + cr = mxGetPr(C_IN); + rr = mxGetPr(R_IN); + + /* Do the actual computations in a subroutine */ + toeplitzC(tr,cr,rr,tm,tn); + + /* Imaginary part */ + if (mxIsComplex(C_IN) || mxIsComplex(R_IN)){ + /*if (!mxIsComplex(C_IN)){ + mexErrMsgTxt("Not implemented (yet). A"); + } + else*/ + if (mxIsComplex(C_IN)) + ci = mxGetPi(C_IN); + else + ci = mxGetPr(mxCreateDoubleMatrix(tm, 1, mxREAL)); + /*if (!mxIsComplex(R_IN)){ + mexErrMsgTxt("Not implemented (yet). B"); + } + else*/ + if (mxIsComplex(R_IN)) + ri = mxGetPi(R_IN); + else + ri = mxGetPr(mxCreateDoubleMatrix(tn, 1, mxREAL)); + ti = mxGetPi(T_OUT); + toeplitzC(ti, ci, ri, tm, tn); + } + + return; + +} + + diff --git a/PQevalAudio/CB/PQCB.m b/PQevalAudio/CB/PQCB.m new file mode 100644 index 0000000..075cfa1 --- /dev/null +++ b/PQevalAudio/CB/PQCB.m @@ -0,0 +1,99 @@ +function [Nc, fc, fl, fu, dz] = PQCB (Version) +% Critical band parameters for the FFT model + +% P. Kabal $Revision: 1.1 $ $Date: 2003/12/07 13:32:57 $ + +B = inline ('7 * asinh (f / 650)'); +BI = inline ('650 * sinh (z / 7)'); + +fL = 80; +fU = 18000; + +% Critical bands - set up the tables +if (strcmp (Version, 'Basic')) + dz = 1/4; +elseif (strcmp (Version, 'Advanced')) + dz = 1/2; +else + error ('PQCB: Invalid version'); +end + +zL = B(fL); +zU = B(fU); +Nc = ceil((zU - zL) / dz); +zl = zL + (0:Nc-1) * dz; +zu = min (zL + (1:Nc) * dz, zU); +zc = 0.5 * (zl + zu); + +fl = BI (zl); +fc = BI (zc); +fu = BI (zu); + +if (strcmp (Version, 'Basic')) + fl = [ 80.000, 103.445, 127.023, 150.762, 174.694, ... + 198.849, 223.257, 247.950, 272.959, 298.317, ... + 324.055, 350.207, 376.805, 403.884, 431.478, ... + 459.622, 488.353, 517.707, 547.721, 578.434, ... + 609.885, 642.114, 675.161, 709.071, 743.884, ... + 779.647, 816.404, 854.203, 893.091, 933.119, ... + 974.336, 1016.797, 1060.555, 1105.666, 1152.187, ... + 1200.178, 1249.700, 1300.816, 1353.592, 1408.094, ... + 1464.392, 1522.559, 1582.668, 1644.795, 1709.021, ... + 1775.427, 1844.098, 1915.121, 1988.587, 2064.590, ... + 2143.227, 2224.597, 2308.806, 2395.959, 2486.169, ... + 2579.551, 2676.223, 2776.309, 2879.937, 2987.238, ... + 3098.350, 3213.415, 3332.579, 3455.993, 3583.817, ... + 3716.212, 3853.817, 3995.399, 4142.547, 4294.979, ... + 4452.890, 4616.482, 4785.962, 4961.548, 5143.463, ... + 5331.939, 5527.217, 5729.545, 5939.183, 6156.396, ... + 6381.463, 6614.671, 6856.316, 7106.708, 7366.166, ... + 7635.020, 7913.614, 8202.302, 8501.454, 8811.450, ... + 9132.688, 9465.574, 9810.536, 10168.013, 10538.460, ... + 10922.351, 11320.175, 11732.438, 12159.670, 12602.412, ... + 13061.229, 13536.710, 14029.458, 14540.103, 15069.295, ... + 15617.710, 16186.049, 16775.035, 17385.420 ]; + fc = [ 91.708, 115.216, 138.870, 162.702, 186.742, ... + 211.019, 235.566, 260.413, 285.593, 311.136, ... + 337.077, 363.448, 390.282, 417.614, 445.479, ... + 473.912, 502.950, 532.629, 562.988, 594.065, ... + 625.899, 658.533, 692.006, 726.362, 761.644, ... + 797.898, 835.170, 873.508, 912.959, 953.576, ... + 995.408, 1038.511, 1082.938, 1128.746, 1175.995, ... + 1224.744, 1275.055, 1326.992, 1380.623, 1436.014, ... + 1493.237, 1552.366, 1613.474, 1676.641, 1741.946, ... + 1809.474, 1879.310, 1951.543, 2026.266, 2103.573, ... + 2183.564, 2266.340, 2352.008, 2440.675, 2532.456, ... + 2627.468, 2725.832, 2827.672, 2933.120, 3042.309, ... + 3155.379, 3272.475, 3393.745, 3519.344, 3649.432, ... + 3784.176, 3923.748, 4068.324, 4218.090, 4373.237, ... + 4533.963, 4700.473, 4872.978, 5051.700, 5236.866, ... + 5428.712, 5627.484, 5833.434, 6046.825, 6267.931, ... + 6497.031, 6734.420, 6980.399, 7235.284, 7499.397, ... + 7773.077, 8056.673, 8350.547, 8655.072, 8970.639, ... + 9297.648, 9636.520, 9987.683, 10351.586, 10728.695, ... + 11119.490, 11524.470, 11944.149, 12379.066, 12829.775, ... + 13294.850, 13780.887, 14282.503, 14802.338, 15341.057, ... + 15899.345, 16477.914, 17077.504, 17690.045 ]; + fu = [ 103.445, 127.023, 150.762, 174.694, 198.849, ... + 223.257, 247.950, 272.959, 298.317, 324.055, ... + 350.207, 376.805, 403.884, 431.478, 459.622, ... + 488.353, 517.707, 547.721, 578.434, 609.885, ... + 642.114, 675.161, 709.071, 743.884, 779.647, ... + 816.404, 854.203, 893.091, 933.113, 974.336, ... + 1016.797, 1060.555, 1105.666, 1152.187, 1200.178, ... + 1249.700, 1300.816, 1353.592, 1408.094, 1464.392, ... + 1522.559, 1582.668, 1644.795, 1709.021, 1775.427, ... + 1844.098, 1915.121, 1988.587, 2064.590, 2143.227, ... + 2224.597, 2308.806, 2395.959, 2486.169, 2579.551, ... + 2676.223, 2776.309, 2879.937, 2987.238, 3098.350, ... + 3213.415, 3332.579, 3455.993, 3583.817, 3716.212, ... + 3853.348, 3995.399, 4142.547, 4294.979, 4452.890, ... + 4643.482, 4785.962, 4961.548, 5143.463, 5331.939, ... + 5527.217, 5729.545, 5939.183, 6156.396, 6381.463, ... + 6614.671, 6856.316, 7106.708, 7366.166, 7635.020, ... + 7913.614, 8202.302, 8501.454, 8811.450, 9132.688, ... + 9465.574, 9810.536, 10168.013, 10538.460, 10922.351, ... + 11320.175, 11732.438, 12159.670, 12602.412, 13061.229, ... + 13536.710, 14029.458, 14540.103, 15069.295, 15617.710, ... + 16186.049, 16775.035, 17385.420, 18000.000 ]; +end diff --git a/PQevalAudio/CB/PQDFTFrame.m b/PQevalAudio/CB/PQDFTFrame.m new file mode 100644 index 0000000..4719334 --- /dev/null +++ b/PQevalAudio/CB/PQDFTFrame.m @@ -0,0 +1,60 @@ +function X2 = PQDFTFrame (x) +% Calculate the DFT of a frame of data (NF values), returning the +% squared-magnitude DFT vector (NF/2 + 1 values) + +% P. Kabal $Revision: 1.1 $ $Date: 2003/12/07 13:32:57 $ + +persistent hw + +NF = 2048; % Frame size (samples) + +if (isempty (hw)) + Amax = 32768; + fc = 1019.5; + Fs = 48000; + Lp = 92; + % Set up the window (including all gains) + GL = PQ_GL (NF, Amax, fc/Fs, Lp); + hw = GL * PQHannWin (NF); +end + +% Window the data +xw = hw .* x; + +% DFT (output is real followed by imaginary) +X = PQRFFT (xw, NF, 1); + +% Squared magnitude +X2 = PQRFFTMSq (X, NF); + +%---------------------------------------- +function GL = PQ_GL (NF, Amax, fcN, Lp) +% Scaled Hann window, including loudness scaling + +% Calculate the gain for the Hann Window +% - level Lp (SPL) corresponds to a sine with normalized frequency +% fcN and a peak value of Amax + +W = NF - 1; +gp = PQ_gp (fcN, NF, W); +GL = 10^(Lp / 20) / (gp * Amax/4 * W); + +%---------- +function gp = PQ_gp (fcN, NF, W) +% Calculate the peak factor. The signal is a sinusoid windowed with +% a Hann window. The sinusoid frequency falls between DFT bins. The +% peak of the frequency response (on a continuous frequency scale) falls +% between DFT bins. The largest DFT bin value is the peak factor times +% the peak of the continuous response. +% fcN - Normalized sinusoid frequency (0-1) +% NF - Frame (DFT) length samples +% NW - Window length samples + +% Distance to the nearest DFT bin +df = 1 / NF; +k = floor (fcN / df); +dfN = min ((k+1) * df - fcN, fcN - k * df); + +dfW = dfN * W; +gp = sin(pi * dfW) / (pi * dfW * (1 - dfW^2)); + diff --git a/PQevalAudio/CB/PQeval.m b/PQevalAudio/CB/PQeval.m new file mode 100644 index 0000000..ec98373 --- /dev/null +++ b/PQevalAudio/CB/PQeval.m @@ -0,0 +1,112 @@ +function [MOVI, Fmem] = PQeval (xR, xT, Fmem) +% PEAQ - Process one frame with the FFT model + +% P. Kabal $Revision: 1.1 $ $Date: 2003/12/07 13:32:58 $ + +NF = 2048; +Version = 'Basic'; + +% Windowed DFT +X2(1,:) = PQDFTFrame (xR); +X2(2,:) = PQDFTFrame (xT); + +% Critical band grouping and frequency spreading +[EbN, Es] = PQ_excitCB (X2); + +% Time domain smoothing => "Excitation patterns" +[Ehs(1,:), Fmem.TDS.Ef(1,:)] = PQ_timeSpread (Es(1,:), Fmem.TDS.Ef(1,:)); +[Ehs(2,:), Fmem.TDS.Ef(2,:)] = PQ_timeSpread (Es(2,:), Fmem.TDS.Ef(2,:)); + +% Level and pattern adaptation => "Spectrally adapted patterns" +[EP, Fmem.Adap] = PQadapt (Ehs, Fmem.Adap, Version, 'FFT'); + +% Modulation patterns +[M, ERavg, Fmem.Env] = PQmodPatt (Es, Fmem.Env); + +% Loudness +MOVI.Loud.NRef = PQloud (Ehs(1,:), Version, 'FFT'); +MOVI.Loud.NTest = PQloud (Ehs(2,:), Version, 'FFT'); + +% Modulation differences +MOVI.MDiff = PQmovModDiffB (M, ERavg); + +% Noise Loudness +MOVI.NLoud.NL = PQmovNLoudB (M, EP); + +% Bandwidth +MOVI.BW = PQmovBW (X2); + +% Noise-to-mask ratios +MOVI.NMR = PQmovNMRB (EbN, Ehs(1,:)); + +% Probability of detection +MOVI.PD = PQmovPD (Ehs); + +% Error harmonic structure +MOVI.EHS.EHS = PQmovEHS (xR, xT, X2); + +%-------------------- +function [EbN, Es] = PQ_excitCB (X2) + +persistent W2 EIN + +NF = 2048; +Version = 'Basic'; +if (isempty (W2)) + Fs = 48000; + f = linspace (0, Fs/2, NF/2+1); + W2 = PQWOME (f); + [Nc, fc] = PQCB (Version); + EIN = PQIntNoise (fc); +end + +% Allocate storage +XwN2 = zeros (1, NF/2+1); + +% Outer and middle ear filtering +Xw2(1,:) = W2 .* X2(1,1:NF/2+1); +Xw2(2,:) = W2 .* X2(2,1:NF/2+1); + +% Form the difference magnitude signal +for (k = 0:NF/2) + XwN2(k+1) = (Xw2(1,k+1) - 2 * sqrt (Xw2(1,k+1) * Xw2(2,k+1)) ... + + Xw2(2,k+1)); +end + +% Group into partial critical bands +Eb(1,:) = PQgroupCB (Xw2(1,:), Version); +Eb(2,:) = PQgroupCB (Xw2(2,:), Version); +EbN = PQgroupCB (XwN2, Version); + +% Add the internal noise term => "Pitch patterns" +E(1,:) = Eb(1,:) + EIN; +E(2,:) = Eb(2,:) + EIN; + +% Critical band spreading => "Unsmeared excitation patterns" +Es(1,:) = PQspreadCB (E(1,:), Version); +Es(2,:) = PQspreadCB (E(2,:), Version); + +%-------------------- +function [Ehs, Ef] = PQ_timeSpread (Es, Ef) + +persistent Nc a b + +if (isempty (Nc)) + [Nc, fc] = PQCB ('Basic'); + Fs = 48000; + NF = 2048; + Nadv = NF / 2; + Fss = Fs / Nadv; + t100 = 0.030; + tmin = 0.008; + [a, b] = PQtConst (t100, tmin, fc, Fss); +end + +% Allocate storage +Ehs = zeros (1, Nc); + +% Time domain smoothing +for (m = 0:Nc-1) + Ef(m+1) = a(m+1) * Ef(m+1) + b(m+1) * Es(m+1); + Ehs(m+1) = max(Ef(m+1), Es(m+1)); +end diff --git a/PQevalAudio/CB/PQgroupCB.m b/PQevalAudio/CB/PQgroupCB.m new file mode 100644 index 0000000..ae6c79d --- /dev/null +++ b/PQevalAudio/CB/PQgroupCB.m @@ -0,0 +1,63 @@ +function Eb = PQgroupCB (X2, Ver) +% Group a DFT energy vector into critical bands +% X2 - Squared-magnitude vector (DFT bins) +% Eb - Excitation vector (fractional critical bands) + +% P. Kabal $Revision: 1.2 $ $Date: 2004/02/05 04:25:46 $ + +persistent Nc kl ku Ul Uu Version + +Emin = 1e-12; + +if (~ strcmp (Ver, Version)) + Version = Ver; + % Set up the DFT bin to critical band mapping + NF = 2048; + Fs = 48000; + [Nc, kl, ku, Ul, Uu] = PQ_CBMapping (NF, Fs, Version); +end + +% Allocate storage +Eb = zeros (1, Nc); + +% Compute the excitation in each band +for (i = 0:Nc-1) + Ea = Ul(i+1) * X2(kl(i+1)+1); % First bin + for (k = (kl(i+1)+1):(ku(i+1)-1)) + Ea = Ea + X2(k+1); % Middle bins + end + Ea = Ea + Uu(i+1) * X2(ku(i+1)+1); % Last bin + Eb(i+1) = max(Ea, Emin); +end + +%--------------------------------------- +function [Nc, kl, ku, Ul, Uu] = PQ_CBMapping (NF, Fs, Version) + +[Nc, fc, fl, fu] = PQCB (Version); + +% Fill in the DFT bin to critical band mappings +df = Fs / NF; +for (i = 0:Nc-1) + fli = fl(i+1); + fui = fu(i+1); + for (k = 0:NF/2) + if ((k+0.5)*df > fli) + kl(i+1) = k; % First bin in band i + Ul(i+1) = (min(fui, (k+0.5)*df) ... + - max(fli, (k-0.5)*df)) / df; + break; + end + end + for (k = NF/2:-1:0) + if ((k-0.5)*df < fui) + ku(i+1) = k; % Last bin in band i + if (kl(i+1) == ku(i+1)) + Uu(i+1) = 0; % Single bin in band + else + Uu(i+1) = (min(fui, (k+0.5)*df) ... + - max(fli, (k-0.5)*df)) / df; + end + break; + end + end +end diff --git a/PQevalAudio/CB/PQspreadCB.m b/PQevalAudio/CB/PQspreadCB.m new file mode 100644 index 0000000..d0efe16 --- /dev/null +++ b/PQevalAudio/CB/PQspreadCB.m @@ -0,0 +1,66 @@ +function Es = PQspreadCB (E, Ver) +% Spread an excitation vector (pitch pattern) - FFT model +% Both E and Es are powers + +% P. Kabal $Revision: 1.1 $ $Date: 2003/12/07 13:32:58 $ + +persistent Bs Version + +if (~ strcmp (Ver, Version)) + Version = Ver; + Nc = length (E); + Bs = PQ_SpreadCB (ones(1,Nc), ones(1,Nc), Version); +end + +Es = PQ_SpreadCB (E, Bs, Version); + +%------------------------- +function Es = PQ_SpreadCB (E, Bs, Ver); + +persistent Nc dz fc aL aUC Version + +% Power law for addition of spreading +e = 0.4; + +if (~ strcmp (Ver, Version)) + Version = Ver; + [Nc, fc, fl, fu, dz] = PQCB (Version); +end + +% Allocate storage +aUCEe = zeros (1, Nc); +Ene = zeros (1, Nc); +Es = zeros (1, Nc); + +% Calculate energy dependent terms +aL = 10^(-2.7 * dz); +for (m = 0:Nc-1) + aUC = 10^((-2.4 - 23 / fc(m+1)) * dz); + aUCE = aUC * E(m+1)^(0.2 * dz); + gIL = (1 - aL^(m+1)) / (1 - aL); + gIU = (1 - aUCE^(Nc-m)) / (1 - aUCE); + En = E(m+1) / (gIL + gIU - 1); + aUCEe(m+1) = aUCE^e; + Ene(m+1) = En^e; +end + +% Lower spreading +Es(Nc-1+1) = Ene(Nc-1+1); +aLe = aL^e; +for (m = Nc-2:-1:0) + Es(m+1) = aLe * Es(m+1+1) + Ene(m+1); +end + +% Upper spreading i > m +for (m = 0:Nc-2) + r = Ene(m+1); + a = aUCEe(m+1); + for (i = m+1:Nc-1) + r = r * a; + Es(i+1) = Es(i+1) + r; + end +end + +for (i = 0:Nc-1) + Es(i+1) = (Es(i+1))^(1/e) / Bs(i+1); +end diff --git a/PQevalAudio/MOV/PQavgMOVB.m b/PQevalAudio/MOV/PQavgMOVB.m new file mode 100644 index 0000000..a78e7ea --- /dev/null +++ b/PQevalAudio/MOV/PQavgMOVB.m @@ -0,0 +1,263 @@ +function MOV = PQavgMOVB (MOVC, Nchan, Nwup) +% Time average MOV precursors + +% P. Kabal $Revision: 1.1 $ $Date: 2003/12/07 13:34:46 $ + +Fs = 48000; +NF = 2048; +Nadv = NF / 2; +Fss = Fs / Nadv; +tdel = 0.5; +tex = 0.050; + +% BandwidthRefB, BandwidthTestB +[MOV(0+1), MOV(1+1)] = PQ_avgBW (MOVC.BW); + +% Total NMRB, RelDistFramesB +[MOV(2+1), MOV(10+1)] = PQ_avgNMRB (MOVC.NMR); + +% WinModDiff1B, AvgModDiff1B, AvgModDiff2B +N500ms = ceil (tdel * Fss); +Ndel = max (0, N500ms - Nwup); +[MOV(3+1), MOV(6+1), MOV(7+1)] = PQ_avgModDiffB (Ndel, MOVC.MDiff); + +% RmsNoiseLoudB +N50ms = ceil (tex * Fss); +Nloud = PQloudTest (MOVC.Loud); +Ndel = max (Nloud + N50ms, Ndel); +MOV(8+1) = PQ_avgNLoudB (Ndel, MOVC.NLoud); + +% ADBB, MFPDB +[MOV(4+1), MOV(9+1)] = PQ_avgPD (MOVC.PD); + +% EHSB +MOV(5+1) = PQ_avgEHS (MOVC.EHS); + +%----------------------------------------- +function EHSB = PQ_avgEHS (EHS) + +[Nchan, Np] = size (EHS.EHS); + +s = 0; +for (j = 0:Nchan-1) + s = s + PQ_LinPosAvg (EHS.EHS(j+1,:)); +end +EHSB = 1000 * s / Nchan; + + +%----------------------------------------- +function [ADBB, MFPDB] = PQ_avgPD (PD) + +global PQopt + +c0 = 0.9; +if (isempty (PQopt)) + c1 = 1; +else + c1 = PQopt.PDfactor; +end + +N = length (PD.Pc); +Phc = 0; +Pcmax = 0; +Qsum = 0; +nd = 0; +for (i = 0:N-1) + Phc = c0 * Phc + (1 - c0) * PD.Pc(i+1); + Pcmax = max (Pcmax * c1, Phc); + + if (PD.Pc(i+1) > 0.5) + nd = nd + 1; + Qsum = Qsum + PD.Qc(i+1); + end +end + +if (nd == 0) + ADBB = 0; +elseif (Qsum > 0) + ADBB = log10 (Qsum / nd); +else + ADBB = -0.5; +end + +MFPDB = Pcmax; + +%----------------------------------------- +function [TotalNMRB, RelDistFramesB] = PQ_avgNMRB (NMR) + +[Nchan, Np] = size (NMR.NMRavg); +Thr = 10^(1.5 / 10); + +s = 0; +for (j = 0:Nchan-1) + s = s + 10 * log10 (PQ_LinAvg (NMR.NMRavg(j+1,:))); +end +TotalNMRB = s / Nchan; + +s = 0; +for (j = 0:Nchan-1) + s = s + PQ_FractThr (Thr, NMR.NMRmax(j+1,:)); +end +RelDistFramesB = s / Nchan; + +%----------------------------------------- +function [BandwidthRefB, BandwidthTestB] = PQ_avgBW (BW) + +[Nchan, Np] = size (BW.BWRef); + +sR = 0; +sT = 0; +for (j = 0:Nchan-1) + sR = sR + PQ_LinPosAvg (BW.BWRef(j+1,:)); + sT = sT + PQ_LinPosAvg (BW.BWTest(j+1,:)); +end +BandwidthRefB = sR / Nchan; +BandwidthTestB = sT / Nchan; + +%----------------------------------------- +function [WinModDiff1B, AvgModDiff1B, AvgModDiff2B] = PQ_avgModDiffB (Ndel, MDiff) + +NF = 2048; +Nadv = NF / 2; +Fs = 48000; + +Fss = Fs / Nadv; +tavg = 0.1; + +[Nchan, Np] = size (MDiff.Mt1B); + +% Sliding window average - delayed average +L = floor (tavg * Fss); % 100 ms sliding window length +s = 0; +for (j = 0:Nchan-1) + s = s + PQ_WinAvg (L, MDiff.Mt1B(j+1,Ndel+1:Np-1+1)); +end +WinModDiff1B = s / Nchan; + +% Weighted linear average - delayed average +s = 0; +for (j = 0:Nchan-1) + s = s + PQ_WtAvg (MDiff.Mt1B(j+1,Ndel+1:Np-1+1), MDiff.Wt(j+1,Ndel+1:Np-1+1)); +end +AvgModDiff1B = s / Nchan; + +% Weighted linear average - delayed average +s = 0; +for (j = 0:Nchan-1) + s = s + PQ_WtAvg (MDiff.Mt2B(j+1,Ndel+1:Np-1+1), MDiff.Wt(j+1,Ndel+1:Np-1+1)); +end +AvgModDiff2B = s / Nchan; + +%----------------------------------------- +function RmsNoiseLoudB = PQ_avgNLoudB (Ndel, NLoud) + +[Nchan, Np] = size (NLoud.NL); + +% RMS average - delayed average and loudness threshold +s = 0; +for (j = 0:Nchan-1) + s = s + PQ_RMSAvg (NLoud.NL(j+1,Ndel+1:Np-1+1)); +end +RmsNoiseLoudB = s / Nchan; + +%----------------------------------- +% Average values values, omitting values which are negative +function s = PQ_LinPosAvg (x) + +N = length(x); + +Nv = 0; +s = 0; +for (i = 0:N-1) + if (x(i+1) >= 0) + s = s + x(i+1); + Nv = Nv + 1; + end +end + +if (Nv > 0) + s = s / Nv; +end + +%---------- +% Fraction of values above a threshold +function Fd = PQ_FractThr (Thr, x) + +N = length (x); + +Nv = 0; +for (i = 0:N-1) + if (x(i+1) > Thr) + Nv = Nv + 1; + end +end + +if (N > 0) + Fd = Nv / N; +else + Fd = 0; +end + +%----------- +% Sliding window (L samples) average +function s = PQ_WinAvg (L, x) + +N = length (x); + +s = 0; +for (i = L-1:N-1) + t = 0; + for (m = 0:L-1) + t = t + sqrt (x(i-m+1)); + end + s = s + (t / L)^4; +end + +if (N >= L) + s = sqrt (s / (N - L + 1)); +end + +%---------- +% Weighted average +function s = PQ_WtAvg (x, W) + +N = length (x); + +s = 0; +sW = 0; +for (i = 0:N-1) + s = s + W(i+1) * x(i+1); + sW = sW + W(i+1); +end + +if (N > 0) + s = s / sW; +end + +%---------- +% Linear average +function LinAvg = PQ_LinAvg (x) + +N = length (x); +s = 0; +for (i = 0:N-1) + s = s + x(i+1); +end + +LinAvg = s / N; + +%---------- +% Square root of average of squared values +function RMSAvg = PQ_RMSAvg (x) + +N = length (x); +s = 0; +for (i = 0:N-1) + s = s + x(i+1)^2; +end + +if (N > 0) + RMSAvg = sqrt(s / N); +else + RMSAvg = 0; +end diff --git a/PQevalAudio/MOV/PQframeMOV.m b/PQevalAudio/MOV/PQframeMOV.m new file mode 100644 index 0000000..af3f9bf --- /dev/null +++ b/PQevalAudio/MOV/PQframeMOV.m @@ -0,0 +1,65 @@ +function PQframeMOV (i, MOVI) +% Copy instantaneous MOV values to a new structure +% The output struct MOVC is a global. +% For most MOV's, they are just copied to the output structure. +% The exception is for the probability of detection, where the +% MOV's measure the maximum frequency-by-frequecy between channels. + +% P. Kabal $Revision: 1.1 $ $Date: 2003/12/07 13:34:46 $ + +global MOVC + +[Nchan,Nc] = size (MOVC.MDiff.Mt1B); + +for (j = 1:Nchan) + + % Modulation differences + MOVC.MDiff.Mt1B(j,i+1) = MOVI(j).MDiff.Mt1B; + MOVC.MDiff.Mt2B(j,i+1) = MOVI(j).MDiff.Mt2B; + MOVC.MDiff.Wt(j,i+1) = MOVI(j).MDiff.Wt; + + % Noise loudness + MOVC.NLoud.NL(j,i+1) = MOVI(j).NLoud.NL; + + % Total loudness + MOVC.Loud.NRef(j,i+1) = MOVI(j).Loud.NRef; + MOVC.Loud.NTest(j,i+1) = MOVI(j).Loud.NTest; + + % Bandwidth + MOVC.BW.BWRef(j,i+1) = MOVI(j).BW.BWRef; + MOVC.BW.BWTest(j,i+1) = MOVI(j).BW.BWTest; + + % Noise-to-mask ratio + MOVC.NMR.NMRavg(j,i+1) = MOVI(j).NMR.NMRavg; + MOVC.NMR.NMRmax(j,i+1) = MOVI(j).NMR.NMRmax; + + % Error harmonic structure + MOVC.EHS.EHS(j,i+1) = MOVI(j).EHS.EHS; +end + +% Probability of detection (collapse frequency bands) +[MOVC.PD.Pc(i+1), MOVC.PD.Qc(i+1)] = PQ_ChanPD (MOVI); + +%---------------------------------------- +function [Pc, Qc] = PQ_ChanPD (MOVI) + +Nc = length (MOVI(1).PD.p); +Nchan = length (MOVI); + +Pr = 1; +Qc = 0; +if (Nchan > 1) + for (m = 0:Nc-1) + pbin = max (MOVI(1).PD.p(m+1), MOVI(2).PD.p(m+1)); + qbin = max (MOVI(1).PD.q(m+1), MOVI(2).PD.q(m+1)); + Pr = Pr * (1 - pbin); + Qc = Qc + qbin; + end +else + for (m = 0:Nc-1) + Pr = Pr * (1 - MOVI.PD.p(m+1)); + Qc = Qc + MOVI.PD.q(m+1); + end +end + +Pc = 1 - Pr; diff --git a/PQevalAudio/MOV/PQloudTest.m b/PQevalAudio/MOV/PQloudTest.m new file mode 100644 index 0000000..69fab55 --- /dev/null +++ b/PQevalAudio/MOV/PQloudTest.m @@ -0,0 +1,29 @@ +function Ndel = PQloudTest (Loud) +% Loudness threshold + +% P. Kabal $Revision: 1.1 $ $Date: 2003/12/07 13:34:46 $ + +[Nchan, Np] = size (Loud.NRef); + +Thr = 0.1; + +% Loudness threshold +Ndel = Np; +for (j = 0:Nchan-1) + Ndel = min (Ndel, PQ_LThresh (Thr, Loud.NRef(j+1,:), Loud.NTest(j+1,:))); +end + +%----------- +function it = PQ_LThresh (Thr, NRef, NTest) +% Loudness check: Look for the first time, the loudness exceeds a threshold +% for both the test and reference signals. + +Np = length (NRef); + +it = Np; +for (i = 0:Np-1) + if (NRef(i+1) > Thr & NTest(i+1) > Thr) + it = i; + break; + end +end diff --git a/PQevalAudio/MOV/PQmovBW.m b/PQevalAudio/MOV/PQmovBW.m new file mode 100644 index 0000000..2e4d45e --- /dev/null +++ b/PQevalAudio/MOV/PQmovBW.m @@ -0,0 +1,45 @@ +function BW = PQmovBW (X2) +% Bandwidth tests + +% P. Kabal $Revision: 1.1 $ $Date: 2003/12/07 13:34:46 $ + +persistent kx kl FR FT N + +if (isempty (kx)) + NF = 2048; + Fs = 48000; + fx = 21586; + kx = round (fx / Fs * NF); % 921 + fl = 8109; + kl = round (fl / Fs * NF); % 346 + FRdB = 10; + FR = 10^(FRdB / 10); + FTdB = 5; + FT = 10^(FTdB / 10); + N = NF / 2; % Limit from pseudo-code +end + +Xth = X2(2,kx+1); +for (k = kx+1:N-1) + Xth = max (Xth, X2(2,k+1)); +end + +% BWRef and BWTest remain negative if the BW of the test signal +% does not exceed FR * Xth for kx-1 <= k <= kl+1 +BW.BWRef = -1; +XthR = FR * Xth; +for (k = kx-1:-1:kl+1) + if (X2(1,k+1) >= XthR) + BW.BWRef = k + 1; + break; + end +end + +BW.BWTest = -1; +XthT = FT * Xth; +for (k = BW.BWRef-1:-1:0) + if (X2(2,k+1) >= XthT) + BW.BWTest = k + 1; + break; + end +end diff --git a/PQevalAudio/MOV/PQmovEHS.m b/PQevalAudio/MOV/PQmovEHS.m new file mode 100644 index 0000000..36e6ca5 --- /dev/null +++ b/PQevalAudio/MOV/PQmovEHS.m @@ -0,0 +1,135 @@ +function EHS = PQmovEHS (xR, xT, X2) +% Calculate the EHS MOV values + +% P. Kabal $Revision: 1.2 $ $Date: 2004/02/05 04:26:19 $ + +persistent NF Nadv NL M Hw + +if (isempty (NL)) + NF = 2048; + Nadv = NF / 2; + Fs = 48000; + Fmax = 9000; + NL = 2^(PQ_log2(NF * Fmax / Fs)); + M = NL; + Hw = (1 / M) * sqrt(8 / 3) * PQHannWin (M); +end + +EnThr = 8000; +kmax = NL + M - 1; + +EnRef = xR(Nadv+1:NF-1+1) * xR(Nadv+1:NF-1+1)'; +EnTest = xT(Nadv+1:NF-1+1) * xT(Nadv+1:NF-1+1)'; + +% Set the return value to be negative for small energy frames +if (EnRef < EnThr & EnTest < EnThr) + EHS = -1; + return; +end + +% Allocate storage +D = zeros (1, kmax); + +% Differences of log values +for (k = 0:kmax-1) + D(k+1) = log (X2(2,k+1) / X2(1,k+1)); +end + +% Correlation computation +C = PQ_Corr (D, NL, M); + +% Normalize the correlations +Cn = PQ_NCorr (C, D, NL, M); +Cnm = (1 / NL) * sum (Cn(1:NL)); + +% Window the correlation +Cw = Hw .* (Cn - Cnm); + +% DFT +cp = PQRFFT (Cw, NL, 1); + +% Squared magnitude +c2 = PQRFFTMSq (cp, NL); + +% Search for a peak after a valley +EHS = PQ_FindPeak (c2, NL/2+1); + +%---------------------------------------- +function log2 = PQ_log2 (a) + +log2 = 0; +m = 1; +while (m < a) + log2 = log2 + 1; + m = 2 * m; +end +log2 = log2 - 1; + +%---------- +function C = PQ_Corr (D, NL, M) +% Correlation calculation + +% Direct computation of the correlation +% for (i = 0:NL-1) +% s = 0; +% for (j = 0:M-1) +% s = s + D(j+1) * D(i+j+1); +% end +% C(i+1) = s; +% end + +% Calculate the correlation indirectly +NFFT = 2 * NL; +D0 = [D(1:M) zeros(1,NFFT-M)]; +D1 = [D(1:M+NL-1) zeros(1,NFFT-(M+NL-1))]; + +% DFTs of the zero-padded sequences +d0 = PQRFFT (D0, NFFT, 1); +d1 = PQRFFT (D1, NFFT, 1); + +% Multiply (complex) sequences +dx(0+1) = d0(0+1) * d1(0+1); +for (n = 1:NFFT/2-1) + m = NFFT/2 + n; + dx(n+1) = d0(n+1) * d1(n+1) + d0(m+1) * d1(m+1); + dx(m+1) = d0(n+1) * d1(m+1) - d0(m+1) * d1(n+1); +end +dx(NFFT/2+1) = d0(NFFT/2+1) * d1(NFFT/2+1); + +% Inverse DFT +Cx = PQRFFT (dx, NFFT, -1); +C = Cx(1:NL); + +%---------- +function Cn = PQ_NCorr (C, D, NL, M) +% Normalize the correlation + +Cn = zeros (1, NL); + +s0 = C(0+1); +sj = s0; +Cn(0+1) = 1; +for (i = 1:NL-1) + sj = sj + (D(i+M-1+1)^2 - D(i-1+1)^2); + d = s0 * sj; + if (d <= 0) + Cn(i+1) = 1; + else + Cn(i+1) = C(i+1) / sqrt (d); + end +end + +%---------- +function EHS = PQ_FindPeak (c2, N) +% Search for a peak after a valley + +cprev = c2(0+1); +cmax = 0; +for (n = 1:N-1) + if (c2(n+1) > cprev) % Rising from a valley + if (c2(n+1) > cmax) + cmax = c2(n+1); + end + end +end +EHS = cmax; diff --git a/PQevalAudio/MOV/PQmovModDiffB.m b/PQevalAudio/MOV/PQmovModDiffB.m new file mode 100644 index 0000000..925fff8 --- /dev/null +++ b/PQevalAudio/MOV/PQmovModDiffB.m @@ -0,0 +1,43 @@ +function MDiff = PQmovModDiffB (M, ERavg) +% Modulation difference related MOV precursors (Basic version) + +% P. Kabal $Revision: 1.1 $ $Date: 2003/12/07 13:34:46 $ + +persistent Nc Ete + +if (isempty (Nc)) + e = 0.3; + [Nc, fc] = PQCB ('Basic'); + Et = PQIntNoise (fc); + for (m = 0:Nc-1) + Ete(m+1) = Et(m+1)^e; + end +end + +% Parameters +negWt2B = 0.1; +offset1B = 1.0; +offset2B = 0.01; +levWt = 100; + +s1B = 0; +s2B = 0; +Wt = 0; +for (m = 0:Nc-1) + if (M(1,m+1) > M(2,m+1)) + num1B = M(1,m+1) - M(2,m+1); + num2B = negWt2B * num1B; + else + num1B = M(2,m+1) - M(1,m+1); + num2B = num1B; + end + MD1B = num1B / (offset1B + M(1,m+1)); + MD2B = num2B / (offset2B + M(1,m+1)); + s1B = s1B + MD1B; + s2B = s2B + MD2B; + Wt = Wt + ERavg(m+1) / (ERavg(m+1) + levWt * Ete(m+1)); +end + +MDiff.Mt1B = (100 / Nc) * s1B; +MDiff.Mt2B = (100 / Nc) * s2B; +MDiff.Wt = Wt; diff --git a/PQevalAudio/MOV/PQmovNLoudB.m b/PQevalAudio/MOV/PQmovNLoudB.m new file mode 100644 index 0000000..5d028bf --- /dev/null +++ b/PQevalAudio/MOV/PQmovNLoudB.m @@ -0,0 +1,33 @@ +function NL = PQmovNLoudB (M, EP) +% Noise Loudness + +% P. Kabal $Revision: 1.1 $ $Date: 2003/12/07 13:34:47 $ + +persistent Nc Et + +if (isempty (Nc)) + [Nc, fc] = PQCB ('Basic'); + Et = PQIntNoise (fc); +end + +% Parameters +alpha = 1.5; +TF0 = 0.15; +S0 = 0.5; +NLmin = 0; +e = 0.23; + +s = 0; +for (m = 0:Nc-1) + sref = TF0 * M(1,m+1) + S0; + stest = TF0 * M(2,m+1) + S0; + beta = exp (-alpha * (EP(2,m+1) - EP(1,m+1)) / EP(1,m+1)); + a = max (stest * EP(2,m+1) - sref * EP(1,m+1), 0); + b = Et(m+1) + sref * EP(1,m+1) * beta; + s = s + (Et(m+1) / stest)^e * ((1 + a / b)^e - 1); +end + +NL = (24 / Nc) * s; +if (NL < NLmin) + NL = 0; +end diff --git a/PQevalAudio/MOV/PQmovNMRB.m b/PQevalAudio/MOV/PQmovNMRB.m new file mode 100644 index 0000000..78c29bc --- /dev/null +++ b/PQevalAudio/MOV/PQmovNMRB.m @@ -0,0 +1,36 @@ +function NMR = PQmovNMRB (EbN, Ehs) +% Noise-to-mask ratio - Basic version +% NMR(1) average NMR +% NMR(2) max NMR + +% P. Kabal $Revision: 1.1 $ $Date: 2003/12/07 13:34:47 $ + +persistent Nc gm + +if (isempty (Nc)) + [Nc, fc, fl, fu, dz] = PQCB ('Basic'); + gm = PQ_MaskOffset (dz, Nc); +end + +NMR.NMRmax = 0; +s = 0; +for (m = 0:Nc-1) + NMRm = EbN(m+1) / (gm(m+1) * Ehs(m+1)); + s = s + NMRm; + if (NMRm > NMR.NMRmax) + NMR.NMRmax = NMRm; + end +end +NMR.NMRavg = s / Nc; + +%---------------------------------------- +function gm = PQ_MaskOffset (dz, Nc) + +for (m = 0:Nc-1) + if (m <= 12 / dz) + mdB = 3; + else + mdB = 0.25 * m * dz; + end + gm(m+1) = 10^(-mdB / 10); +end diff --git a/PQevalAudio/MOV/PQmovPD.m b/PQevalAudio/MOV/PQmovPD.m new file mode 100644 index 0000000..7c62a8c --- /dev/null +++ b/PQevalAudio/MOV/PQmovPD.m @@ -0,0 +1,42 @@ +function PD = PQmovPD (Ehs) +% Probability of detection + +% P. Kabal $Revision: 1.1 $ $Date: 2003/12/07 13:34:47 $ + +Nc = length (Ehs); + +% Allocate storage +PD.p = zeros (1, Nc); +PD.q = zeros (1, Nc); + +persistent c g d1 d2 bP bM + +if (isempty (c)) + c = [-0.198719 0.0550197 -0.00102438 5.05622e-6 9.01033e-11]; + d1 = 5.95072; + d2 = 6.39468; + g = 1.71332; + bP = 4; + bM = 6; +end + +for (m = 0:Nc-1) + EdBR = 10 * log10 (Ehs(1,m+1)); + EdBT = 10 * log10 (Ehs(2,m+1)); + edB = EdBR - EdBT; + if (edB > 0) + L = 0.3 * EdBR + 0.7 * EdBT; + b = bP; + else + L = EdBT; + b = bM; + end + if (L > 0) + s = d1 * (d2 / L)^g ... + + c(1) + L * (c(2) + L * (c(3) + L * (c(4) + L * c(5)))); + else + s = 1e30; + end + PD.p(m+1) = 1 - 0.5^((edB / s)^b); % Detection probability + PD.q(m+1) = abs (fix(edB)) / s; % Steps above threshold +end diff --git a/PQevalAudio/MOV/PQprtMOV.m b/PQevalAudio/MOV/PQprtMOV.m new file mode 100644 index 0000000..5e7b064 --- /dev/null +++ b/PQevalAudio/MOV/PQprtMOV.m @@ -0,0 +1,35 @@ +function PQprtMOV (MOV, ODG) +% Print MOV values (PEAQ Basic version) + +% P. Kabal $Revision: 1.1 $ $Date: 2003/12/07 13:34:47 $ + +N = length (MOV); +PQ_NMOV_B = 11; +PQ_NMOV_A = 5; + +fprintf ('Model Output Variables:\n'); +if (N == PQ_NMOV_B) + fprintf (' BandwidthRefB: %g\n', MOV(1)); + fprintf (' BandwidthTestB: %g\n', MOV(2)); + fprintf (' Total NMRB: %g\n', MOV(3)); + fprintf (' WinModDiff1B: %g\n', MOV(4)); + fprintf (' ADBB: %g\n', MOV(5)); + fprintf (' EHSB: %g\n', MOV(6)); + fprintf (' AvgModDiff1B: %g\n', MOV(7)); + fprintf (' AvgModDiff2B: %g\n', MOV(8)); + fprintf (' RmsNoiseLoudB: %g\n', MOV(9)); + fprintf (' MFPDB: %g\n', MOV(10)); + fprintf (' RelDistFramesB: %g\n', MOV(11)); +elseif (N == NMOV_A) + fprintf (' RmsModDiffA: %g\n', MOV(1)); + fprintf (' RmsNoiseLoudAsymA: %g\n', MOV(2)); + fprintf (' Segmental NMRB: %g\n', MOV(3)); + fprintf (' EHSB: %g\n', MOV(4)); + fprintf (' AvgLinDistA: %g\n', MOV(5)); +else + error ('Invalid number of MOVs'); +end + +fprintf ('Objective Difference Grade: %.3f\n', ODG); + +return; diff --git a/PQevalAudio/MOV/PQprtMOVCi.m b/PQevalAudio/MOV/PQprtMOVCi.m new file mode 100644 index 0000000..d080c48 --- /dev/null +++ b/PQevalAudio/MOV/PQprtMOVCi.m @@ -0,0 +1,40 @@ +function PQprtMOVCi (Nchan, i, MOVC) +% Print MOV precursors + +% P. Kabal $Revision: 1.1 $ $Date: 2003/12/07 13:34:47 $ + +fprintf ('Frame: %d\n', i); + +if (Nchan == 1) + fprintf (' Ntot : %g %g\n', ... + MOVC.Loud.NRef(1,i+1), MOVC.Loud.NTest(1,i+1)); + fprintf (' ModDiff: %g %g %g\n', ... + MOVC.MDiff.Mt1B(1,i+1), MOVC.MDiff.Mt2B(1,i+1), MOVC.MDiff.Wt(1,i+1)); + fprintf (' NL : %g\n', MOVC.NLoud.NL(1,i+1)); + fprintf (' BW : %g %g\n', ... + MOVC.BW.BWRef(1,i+1), MOVC.BW.BWTest(1,i+1)); + fprintf (' NMR : %g %g\n', ... + MOVC.NMR.NMRavg(1,i+1), MOVC.NMR.NMRmax(1,i+1)); + fprintf (' PD : %g %g\n', MOVC.PD.Pc(i+1), MOVC.PD.Qc(i+1)); + fprintf (' EHS : %g\n', 1000 * MOVC.EHS.EHS(1,i+1)); +else + fprintf (' Ntot : %g %g // %g %g\n', ... + MOVC.Loud.NRef(1,i+1), MOVC.Loud.NTest(1,i+1), ... + MOVC.Loud.NRef(2,i+1), MOVC.Loud.NTest(2,i+1)); + fprintf (' ModDiff: %g %g %g // %g %g %g\n', ... + MOVC.MDiff.Mt1B(1,i+1), MOVC.MDiff.Mt2B(1,i+1), MOVC.MDiff.Wt(1,i+1), ... + MOVC.MDiff.Mt1B(2,i+1), MOVC.MDiff.Mt2B(2,i+1), MOVC.MDiff.Wt(2,i+1)); + fprintf (' NL : %g // %g\n', ... + MOVC.NLoud.NL(1,i+1), ... + MOVC.NLoud.NL(2,i+1)); + fprintf (' BW : %g %g // %g %g\n', ... + MOVC.BW.BWRef(1,i+1), MOVC.BW.BWTest(1,i+1), ... + MOVC.BW.BWRef(2,i+1), MOVC.BW.BWTest(2,i+1)); + fprintf (' NMR : %g %g // %g %g\n', ... + MOVC.NMR.NMRavg(1,i+1), MOVC.NMR.NMRmax(1,i+1), ... + MOVC.NMR.NMRavg(2,i+1), MOVC.NMR.NMRmax(2,i+1)); + fprintf (' PD : %g %g\n', MOVC.PD.Pc(i+1), MOVC.PD.Qc(i+1)); + fprintf (' EHS : %g // %g\n', ... + 1000 * MOVC.EHS.EHS(1,i+1), ... + 1000 * MOVC.EHS.EHS(2,i+1)); +end diff --git a/PQevalAudio/Misc/PQHannWin.m b/PQevalAudio/Misc/PQHannWin.m new file mode 100644 index 0000000..f028dd4 --- /dev/null +++ b/PQevalAudio/Misc/PQHannWin.m @@ -0,0 +1,10 @@ +function hw = PQHannWin (NF) +% Hann window + +% P. Kabal $Revision: 1.1 $ $Date: 2003/12/07 13:34:10 $ + +hw = zeros (1, NF); + +for (n = 0:NF-1) + hw(n+1) = 0.5 * (1 - cos(2 * pi * n / (NF-1))); +end diff --git a/PQevalAudio/Misc/PQIntNoise.m b/PQevalAudio/Misc/PQIntNoise.m new file mode 100644 index 0000000..7dcb462 --- /dev/null +++ b/PQevalAudio/Misc/PQIntNoise.m @@ -0,0 +1,10 @@ +function EIN = PQIntNoise (f) +% Generate the internal noise energy vector + +% P. Kabal $Revision: 1.1 $ $Date: 2003/12/07 13:34:10 $ + +N = length (f); +for (m = 0:N-1) + INdB = 1.456 * (f(m+1) / 1000)^(-0.8); + EIN(m+1) = 10^(INdB / 10); +end diff --git a/PQevalAudio/Misc/PQRFFT.m b/PQevalAudio/Misc/PQRFFT.m new file mode 100644 index 0000000..9ac3c8c --- /dev/null +++ b/PQevalAudio/Misc/PQRFFT.m @@ -0,0 +1,33 @@ +function X = PQRFFT (x, N, ifn) +% Calculate the DFT of a real N-point sequence or the inverse +% DFT corresponding to a real N-point sequence. +% ifn > 0, forward transform +% input x(n) - N real values +% output X(k) - The first N/2+1 points are the real +% parts of the transform, the next N/2-1 points +% are the imaginary parts of the transform. However +% the imaginary part for the first point and the +% middle point which are known to be zero are not +% stored. +% ifn < 0, inverse transform +% input X(k) - The first N/2+1 points are the real +% parts of the transform, the next N/2-1 points +% are the imaginary parts of the transform. However +% the imaginary part for the first point and the +% middle point which are known to be zero are not +% stored. +% output x(n) - N real values + +% P. Kabal $Revision: 1.1 $ $Date: 2003/12/07 13:34:11 $ + +if (ifn > 0) + X = fft (x, N); + XR = real(X(0+1:N/2+1)); + XI = imag(X(1+1:N/2-1+1)); + X = [XR XI]; +else + xR = [x(0+1:N/2+1)]; + xI = [0 x(N/2+1+1:N-1+1) 0]; + x = complex ([xR xR(N/2-1+1:-1:1+1)], [xI -xI(N/2-1+1:-1:1+1)]); + X = real (ifft (x, N)); +end diff --git a/PQevalAudio/Misc/PQRFFTMSq.m b/PQevalAudio/Misc/PQRFFTMSq.m new file mode 100644 index 0000000..a032b76 --- /dev/null +++ b/PQevalAudio/Misc/PQRFFTMSq.m @@ -0,0 +1,13 @@ +function X2 = PQRFFTMSq (X, N) +% Calculate the magnitude squared frequency response from the +% DFT values corresponding to a real signal (assumes N is even) + +% P. Kabal $Revision: 1.1 $ $Date: 2003/12/07 13:34:11 $ + +X2 = zeros (1, N/2+1); + +X2(0+1) = X(0+1)^2; +for (k = 1:N/2-1) + X2(k+1) = X(k+1)^2 + X(N/2+k+1)^2; +end +X2(N/2+1) = X(N/2+1)^2; diff --git a/PQevalAudio/Misc/PQWOME.m b/PQevalAudio/Misc/PQWOME.m new file mode 100644 index 0000000..e83f147 --- /dev/null +++ b/PQevalAudio/Misc/PQWOME.m @@ -0,0 +1,13 @@ +function W2 = PQWOME (f) +% Generate the weighting for the outer & middle ear filtering +% Note: The output is a magnitude-squared vector + +% P. Kabal $Revision: 1.1 $ $Date: 2003/12/07 13:34:11 $ + +N = length (f); +for (k = 0:N-1) + fkHz = f(k+1) / 1000; + AdB = -2.184 * fkHz^(-0.8) + 6.5 * exp(-0.6 * (fkHz - 3.3)^2) ... + - 0.001 * fkHz^(3.6); + W2(k+1) = 10^(AdB / 10); +end diff --git a/PQevalAudio/Misc/PQdataBoundary.m b/PQevalAudio/Misc/PQdataBoundary.m new file mode 100644 index 0000000..6825f81 --- /dev/null +++ b/PQevalAudio/Misc/PQdataBoundary.m @@ -0,0 +1,106 @@ +function Lim = PQdataBoundary (WAV, Nchan, StartS, Ns) +% Search for the data boundaries in a file +% StartS - starting sample frame +% Ns - Number of sample frames + +% P. Kabal $Revision: 1.1 $ $Date: 2003/12/07 13:34:10 $ + +PQ_L = 5; +Amax = 32768; +NBUFF = 2048; +PQ_ATHR = 200 * (Amax / 32768); + +% Search from the beginning of the file +Lim(1) = -1; +is = StartS; +EndS = StartS + Ns - 1; +while (is <= EndS) + Nf = min (EndS - is + 1, NBUFF); + x = PQgetData (WAV, is, Nf); + for (k = 0:Nchan-1) + Lim(1) = max (Lim(1), PQ_DataStart (x(k+1,:), Nf, PQ_L, PQ_ATHR)); + end + if (Lim(1) >= 0) + Lim(1) = Lim(1) + is; + break + end + is = is + NBUFF - (PQ_L-1); +end + +% Search from the end of the file +% This loop is written as if it is going in a forward direction +% - When the "forward" position is i, the "backward" position is +% EndS - (i - StartS + 1) + 1 +Lim(2) = -1; +is = StartS; +while (is <= EndS) + Nf = min (EndS - is + 1, NBUFF); + ie = is + Nf - 1; % Forward limits [is, ie] + js = EndS - (ie - StartS + 1) + 1; % Backward limits [js, js+Nf-1] + x = PQgetData (WAV, js, Nf); + for (k = 0:Nchan-1) + Lim(2) = max (Lim(2), PQ_DataEnd (x(k+1,:), Nf, PQ_L, PQ_ATHR)); + end + if (Lim(2) >= 0) + Lim(2) = Lim(2) + js; + break + end + is = is + NBUFF - (PQ_L-1); +end + +% Sanity checks +if (~ ((Lim(1) >= 0 & Lim(2) >= 0) | (Lim(1) < 0 & Lim(2) < 0))) + error ('>>> PQdataBoundary: limits have difference signs'); +end +if (~(Lim(1) <= Lim(2))) + error ('>>> PQdataBoundary: Lim(1) > Lim(2)'); +end + +if (Lim(1) < 0) + Lim(1) = 0; + Lim(2) = 0; +end + +%---------- +function ib = PQ_DataStart (x, N, L, Thr) + +ib = -1; +s = 0; +M = min (N, L); +for (i = 0:M-1) + s = s + abs (x(i+1)); +end +if (s > Thr) + ib = 0; + return +end + +for (i = 1:N-L) % i is the first sample + s = s + (abs (x(i+L-1+1)) - abs (x(i-1+1))); % L samples apart + if (s > Thr) + ib = i; + return + end +end + +%---------- +function ie = PQ_DataEnd (x, N, L, Thr) + +ie = -1; +s = 0; +M = min (N, L); +for (i = N-M:N-1) + s = s + abs (x(i+1)); +end +if (s > Thr) + ie = N-1; + return +end + +for (i = N-2:-1:L-1) % i is the last sample + s = s + (abs (x(i-L+1+1)) - abs (x(i+1+1))); % L samples apart + if (s > Thr) + ie = i; + return + end +end diff --git a/PQevalAudio/Misc/PQgetData.m b/PQevalAudio/Misc/PQgetData.m new file mode 100644 index 0000000..8921e23 --- /dev/null +++ b/PQevalAudio/Misc/PQgetData.m @@ -0,0 +1,59 @@ +function x = PQgetData (WAV, i, N) +% Get data from internal buffer or file +% i - file position +% N - number of samples +% x - output data (scaled to the range -32768 to +32767) + +% Only two files can be "active" at a time. +% N = 0 resets the buffer + + +% P. Kabal $Revision: 1.1 $ $Date: 2003/12/07 13:34:10 $ + +persistent Buff + +iB = WAV.iB + 1; +if (N == 0) + Buff(iB).N = 20 * 1024; % Fixed size + Buff(iB).x = PQ_ReadWAV (WAV, i, Buff(iB).N); + Buff(iB).i = i; +end + +if (N > Buff(iB).N) + error ('>>> PQgetData: Request exceeds buffer size'); +end + +% Check if requested data is not already in the buffer +is = i - Buff(iB).i; +if (is < 0 | is + N - 1 > Buff(iB).N - 1) + Buff(iB).x = PQ_ReadWAV (WAV, i, Buff(iB).N); + Buff(iB).i = i; +end + +% Copy the data +Nchan = WAV.Nchan; +is = i - Buff(iB).i; +x = Buff(iB).x(1:Nchan,is+1:is+N-1+1); + +%------ +function x = PQ_ReadWAV (WAV, i, N) +% This function considers the data to extended with zeros before and +% after the data in the file. If the starting offset i is negative, +% zeros are filled in before the data starts at offset 0. If the request +% extends beyond the end of data in the file, zeros are appended. + +Amax = 32768; +Nchan = WAV.Nchan; + +x = zeros (Nchan, N); + +Nz = 0; +if (i < 0) + Nz = min (-i, N); + i = i + Nz; +end + +Ns = min (N - Nz, WAV.Nframe - i); +if (i >= 0 & Ns > 0) + x(1:Nchan,Nz+1:Nz+Ns-1+1) = Amax * (audioread (WAV.Fname, [i+1 i+Ns-1+1]))'; +end diff --git a/PQevalAudio/Misc/PQinitFMem.m b/PQevalAudio/Misc/PQinitFMem.m new file mode 100644 index 0000000..71d077d --- /dev/null +++ b/PQevalAudio/Misc/PQinitFMem.m @@ -0,0 +1,13 @@ +function Fmem = PQinitFMem (Nc, PCinit) +% Initialize the filter memories + +% P. Kabal $Revision: 1.1 $ $Date: 2003/12/07 13:34:10 $ + +Fmem.TDS.Ef(1:2,1:Nc) = 0; +Fmem.Adap.P(1:2,1:Nc) = 0; +Fmem.Adap.Rn(1:Nc) = 0; +Fmem.Adap.Rd(1:Nc) = 0; +Fmem.Adap.PC(1:2,1:Nc) = PCinit; +Fmem.Env.Ese(1:2,1:Nc) = 0; +Fmem.Env.DE(1:2,1:Nc) = 0; +Fmem.Env.Eavg(1:2,1:Nc) = 0; diff --git a/PQevalAudio/Misc/PQtConst.m b/PQevalAudio/Misc/PQtConst.m new file mode 100644 index 0000000..ddc78d8 --- /dev/null +++ b/PQevalAudio/Misc/PQtConst.m @@ -0,0 +1,13 @@ +function [a, b] = PQtConst (t100, tmin, f , Fs) +% Calculate the difference equation parameters. The time +% constant of the difference equation depends on the center +% frequencies. + +% P. Kabal $Revision: 1.1 $ $Date: 2003/12/07 13:34:11 $ + +N = length (f); +for (m = 0:N-1) + t = tmin + (100 / f(m+1)) * (t100 - tmin); + a(m+1) = exp (-1 / (Fs * t)); + b(m+1) = (1 - a(m+1)); +end diff --git a/PQevalAudio/Misc/PQwavFilePar.m b/PQevalAudio/Misc/PQwavFilePar.m new file mode 100644 index 0000000..e49d817 --- /dev/null +++ b/PQevalAudio/Misc/PQwavFilePar.m @@ -0,0 +1,32 @@ +function WAV = PQwavFilePar (File) +% Print a WAVE file header, pick up the file parameters + +% P. Kabal $Revision: 1.1 $ $Date: 2003/12/07 13:34:11 $ + +persistent iB + +if (isempty (iB)) + iB = 0; +else + iB = mod (iB + 1, 2); % Only two files can be "active" at a time +end + +%[size WAV.Fs Nbit] = wavread (File, 'size'); +[sound, WAV.Fs] = audioread(File); +Nbit = 16; +WAV.Fname = File; +WAV.Nframe = size(sound, 1); +WAV.Nchan = size(sound, 2); +WAV.iB = iB; % Buffer number + +% Initialize the buffer +PQgetData (WAV, 0, 0); + +fprintf (' WAVE file: %s\n', File); +if (WAV.Nchan == 1) + fprintf (' Number of samples : %d (%.4g s)\n', WAV.Nframe, WAV.Nframe / WAV.Fs); +else + fprintf (' Number of frames : %d (%.4g s)\n', WAV.Nframe, WAV.Nframe / WAV.Fs); +end +fprintf (' Sampling frequency: %g\n', WAV.Fs); +fprintf (' Number of channels: %d (%d-bit integer)\n', WAV.Nchan, Nbit); diff --git a/PQevalAudio/PEAQTest.m b/PQevalAudio/PEAQTest.m new file mode 100644 index 0000000..17b2ef4 --- /dev/null +++ b/PQevalAudio/PEAQTest.m @@ -0,0 +1,9 @@ +clear; clc; +%addpath('PQevalAudio', 'PQevalAudio/CB','PQevalAudio/Misc','PQevalAudio/MOV', 'PQevalAudio/Patt') + +ref = '2_clean_48.wav' % 16 times oversampling +test = '2_noise_48.wav' % 4 times oversampling +%test = '../2_output_48.wav' % 4 times oversampling + + +[odg, movb] = PQevalAudio(ref, test) diff --git a/PQevalAudio/PQevalAudio.m b/PQevalAudio/PQevalAudio.m new file mode 100644 index 0000000..3eafa52 --- /dev/null +++ b/PQevalAudio/PQevalAudio.m @@ -0,0 +1,177 @@ +function [ODG, MOVB]= PQevalAudio (Fref, Ftest, StartS, EndS) +% Perceptual evaluation of audio quality. + +% - StartS shifts the frames, so that the first frame starts at that sample. +% This is a two element array, one element for each input file. If StartS is +% a scalar, it applies to both files. +% - EndS marks the end of data. The processing stops with the last frame that +% contains that sample. This is a two element array, one element for each +% input file. If EndS is as scalar, it applies to both files. + +% P. Kabal $Revision: 1.2 $ $Date: 2004/02/05 04:25:24 $ + +% Globals (to save on copying in/out of functions) +global MOVC PQopt + +% Analysis parameters +NF = 2048; +Nadv = NF / 2; +Version = 'Basic'; + +% Options +PQopt.ClipMOV = 0; +PQopt.PCinit = 0; +PQopt.PDfactor = 1; +PQopt.Ni = 1; +PQopt.DelayOverlap = 1; +PQopt.DataBounds = 1; +PQopt.EndMin = NF / 2; + +%addpath ('CB', 'MOV', 'Misc', 'Patt'); + +if (nargin < 3) + StartS = [0, 0]; +end +if (nargin < 4) + EndS = []; +end + +% Get the number of samples and channels for each file +WAV(1) = PQwavFilePar (Fref); +WAV(2) = PQwavFilePar (Ftest); + +% Reconcile file differences +PQ_CheckWAV (WAV); +if (WAV(1).Nframe ~= WAV(2).Nframe) + disp ('>>> Number of samples differ: using the minimum'); +end + +% Data boundaries +Nchan = WAV(1).Nchan; +[StartS, Fstart, Fend] = PQ_Bounds (WAV, Nchan, StartS, EndS, PQopt); + +% Number of PEAQ frames +Np = Fend - Fstart + 1; +if (PQopt.Ni < 0) + PQopt.Ni = ceil (Np / abs(PQopt.Ni)); +end + +% Initialize the MOV structure +MOVC = PQ_InitMOVC (Nchan, Np); + +% Initialize the filter memory +Nc = PQCB (Version); +for (j = 0:Nchan-1) + Fmem(j+1) = PQinitFMem (Nc, PQopt.PCinit); +end + +is = 0; +for (i = -Fstart:Np-1) + + % Read a frame of data + xR = PQgetData (WAV(1), StartS(1) + is, NF); % Reference file + xT = PQgetData (WAV(2), StartS(2) + is, NF); % Test file + is = is + Nadv; + + % Process a frame + for (j = 0:Nchan-1) + [MOVI(j+1), Fmem(j+1)] = PQeval (xR(j+1,:), xT(j+1,:), Fmem(j+1)); + end + + if (i >= 0) + % Move the MOV precursors into a new structure + PQframeMOV (i, MOVI); % Output is in global MOVC + + % Print the MOV precursors + %if (PQopt.Ni ~= 0 & mod (i, PQopt.Ni) == 0) + % PQprtMOVCi (Nchan, i, MOVC); + %nd + end +end + +% Time average of the MOV values +if (PQopt.DelayOverlap) + Nwup = Fstart; +else + Nwup = 0; +end +MOVB = PQavgMOVB (MOVC, Nchan, Nwup); + +% Neural net +ODG = PQnNet (MOVB); + +% Summary printout +%PQprtMOV (MOVB, ODG); + +%---------- +function PQ_CheckWAV (WAV) +% Check the file parameters + +Fs = 48000; + +if (WAV(1).Nchan ~= WAV(2).Nchan) + error ('>>> Number of channels differ'); +end +if (WAV(1).Nchan > 2) + error ('>>> Too many input channels'); +end +if (WAV(1).Nframe ~= WAV(2).Nframe) + disp ('>>> Number of samples differ'); +end +if (WAV(1).Fs ~= WAV(2).Fs) + error ('>>> Sampling frequencies differ'); +end +if (WAV(1).Fs ~= Fs) + error ('>>> Invalid Sampling frequency: only 48 kHz supported'); +end + +%---------- +function [StartS, Fstart, Fend] = PQ_Bounds (WAV, Nchan, StartS, EndS, PQopt) + +PQ_NF = 2048; +PQ_NADV = (PQ_NF / 2); + +if (isempty (StartS)) + StartS(1) = 0; + StartS(2) = 0; +elseif (length (StartS) == 1) + StartS(2) = StartS(1); +end +Ns = WAV(1).Nframe; + +% Data boundaries (determined from the reference file) +if (PQopt.DataBounds) + Lim = PQdataBoundary (WAV(1), Nchan, StartS(1), Ns); + %fprintf ('PEAQ Data Boundaries: %ld (%.3f s) - %ld (%.3f s)\n', ... + % Lim(1), Lim(1)/WAV(1).Fs, Lim(2), Lim(2)/WAV(1).Fs); +else + Lim = [Starts(1), StartS(1) + Ns - 1]; +end + +% Start frame number +Fstart = floor ((Lim(1) - StartS(1)) / PQ_NADV); + +% End frame number +Fend = floor ((Lim(2) - StartS(1) + 1 - PQopt.EndMin) / PQ_NADV); + +%---------- +function MOVC = PQ_InitMOVC (Nchan, Np) +MOVC.MDiff.Mt1B = zeros (Nchan, Np); +MOVC.MDiff.Mt2B = zeros (Nchan, Np); +MOVC.MDiff.Wt = zeros (Nchan, Np); + +MOVC.NLoud.NL = zeros (Nchan, Np); + +MOVC.Loud.NRef = zeros (Nchan, Np); +MOVC.Loud.NTest = zeros (Nchan, Np); + +MOVC.BW.BWRef = zeros (Nchan, Np); +MOVC.BW.BWTest = zeros (Nchan, Np); + +MOVC.NMR.NMRavg = zeros (Nchan, Np); +MOVC.NMR.NMRmax = zeros (Nchan, Np); + +MOVC.PD.Pc = zeros (1, Np); +MOVC.PD.Qc = zeros (1, Np); + +MOVC.EHS.EHS = zeros (Nchan, Np); diff --git a/PQevalAudio/PQnNet.m b/PQevalAudio/PQnNet.m new file mode 100644 index 0000000..641e93e --- /dev/null +++ b/PQevalAudio/PQnNet.m @@ -0,0 +1,100 @@ +function ODG = PQnNetB (MOV) +% Neural net to get the final ODG + +% P. Kabal $Revision: 1.1 $ $Date: 2003/12/07 13:27:44 $ + +persistent amin amax wx wxb wy wyb bmin bmax I J CLIPMOV +global PQopt + +if (isempty (amin)) + I = length (MOV); + if (I == 11) + [amin, amax, wx, wxb, wy, wyb, bmin, bmax] = NNetPar ('Basic'); + else + [amin, amax, wx, wxb, wy, wyb, bmin, bmax] = NNetPar ('Advanced'); + end + [I, J] = size (wx); +end + +sigmoid = inline ('1 / (1 + exp(-x))'); + +% Scale the MOV's +Nclip = 0; +MOVx = zeros (1, I); +for (i = 0:I-1) + MOVx(i+1) = (MOV(i+1) - amin(i+1)) / (amax(i+1) - amin(i+1)); + if (~ isempty (PQopt) & PQopt.ClipMOV ~= 0) + if (MOVx(i+1) < 0) + MOVx(i+1) = 0; + Nclip = Nclip + 1; + elseif (MOVx(i+1) > 1) + MOVx(i+1) = 1; + Nclip = Nclip + 1; + end + end +end +if (Nclip > 0) + fprintf ('>>> %d MOVs clipped\n', Nclip); +end + +% Neural network +DI = wyb; +for (j = 0:J-1) + arg = wxb(j+1); + for (i = 0:I-1) + arg = arg + wx(i+1,j+1) * MOVx(i+1); + end + DI = DI + wy(j+1) * sigmoid (arg); +end + +ODG = bmin + (bmax - bmin) * sigmoid (DI); + +function [amin, amax, wx, wxb, wy, wyb, bmin, bmax] = NNetPar (Version) + +if (strcmp (Version, 'Basic')) + amin = ... + [393.916656, 361.965332, -24.045116, 1.110661, -0.206623, ... + 0.074318, 1.113683, 0.950345, 0.029985, 0.000101, ... + 0]; + amax = ... + [921, 881.131226, 16.212030, 107.137772, 2.886017, ... + 13.933351, 63.257874, 1145.018555, 14.819740, 1, ... + 1]; + wx = ... + [ [ -0.502657, 0.436333, 1.219602 ]; + [ 4.307481, 3.246017, 1.123743 ]; + [ 4.984241, -2.211189, -0.192096 ]; + [ 0.051056, -1.762424, 4.331315 ]; + [ 2.321580, 1.789971, -0.754560 ]; + [ -5.303901, -3.452257, -10.814982 ]; + [ 2.730991, -6.111805, 1.519223 ]; + [ 0.624950, -1.331523, -5.955151 ]; + [ 3.102889, 0.871260, -5.922878 ]; + [ -1.051468, -0.939882, -0.142913 ]; + [ -1.804679, -0.503610, -0.620456 ] ]; + wxb = ... + [ -2.518254, 0.654841, -2.207228 ]; + wy = ... + [ -3.817048, 4.107138, 4.629582 ]; + wyb = -0.307594; + bmin = -3.98; + bmax = 0.22; +else + amin = ... + [ 13.298751, 0.041073, -25.018791, 0.061560, 0.024523 ]; + amax = ... + [ 2166.5, 13.24326, 13.46708, 10.226771, 14.224874 ]; + wx = ... + [ [ 21.211773, -39.913052, -1.382553, -14.545348, -0.320899 ]; + [ -8.981803, 19.956049, 0.935389, -1.686586, -3.238586 ]; + [ 1.633830, -2.877505, -7.442935, 5.606502, -1.783120 ]; + [ 6.103821, 19.587435, -0.240284, 1.088213, -0.511314 ]; + [ 11.556344, 3.892028, 9.720441, -3.287205, -11.031250 ] ]; + wxb = ... + [ 1.330890, 2.686103, 2.096598, -1.327851, 3.087055 ]; + wy = ... + [ -4.696996, -3.289959, 7.004782, 6.651897, 4.009144 ]; + wyb = -1.360308; + bmin = -3.98; + bmax = 0.22; +end diff --git a/PQevalAudio/Patt/PQadapt.m b/PQevalAudio/Patt/PQadapt.m new file mode 100644 index 0000000..5beb8b3 --- /dev/null +++ b/PQevalAudio/Patt/PQadapt.m @@ -0,0 +1,107 @@ +function [EP, Fmem] = PQadapt (Ehs, Fmem, Ver, Mod) +% Level and pattern adaptation + +% P. Kabal $Revision: 1.1 $ $Date: 2003/12/07 13:35:08 $ + +persistent a b Nc M1 M2 Version Model + +if (~strcmp (Ver, Version) | ~strcmp (Mod, Model)) + Version = Ver; + Model = Mod; + if (strcmp (Model, 'FFT')) + [Nc, fc] = PQCB (Version); + NF = 2048; + Nadv = NF / 2; + else + [Nc, fc] = PQFB; + Nadv = 192; + end + Version = Ver; + Model = Mod; + Fs = 48000; + Fss = Fs / Nadv; + t100 = 0.050; + tmin = 0.008; + [a b] = PQtConst (t100, tmin, fc, Fss); + [M1, M2] = PQ_M1M2 (Version, Model); +end + +% Allocate memory +EP = zeros (2, Nc); +R = zeros (2, Nc); + +% Smooth the excitation patterns +% Calculate the correlation terms +sn = 0; +sd = 0; +for (m = 0:Nc-1) + Fmem.P(1,m+1) = a(m+1) * Fmem.P(1,m+1) + b(m+1) * Ehs(1,m+1); + Fmem.P(2,m+1) = a(m+1) * Fmem.P(2,m+1) + b(m+1) * Ehs(2,m+1); + sn = sn + sqrt (Fmem.P(2,m+1) * Fmem.P(1,m+1)); + sd = sd + Fmem.P(2,m+1); +end + +% Level correlation +CL = (sn / sd)^2; + +for (m = 0:Nc-1) + +% Scale one of the signals to match levels + if (CL > 1) + EP(1,m+1) = Ehs(1,m+1) / CL; + EP(2,m+1) = Ehs(2,m+1); + else + EP(1,m+1) = Ehs(1,m+1); + EP(2,m+1) = Ehs(2,m+1) * CL; + end + +% Calculate a pattern match correction factor + Fmem.Rn(m+1) = a(m+1) * Fmem.Rn(m+1) + EP(2,m+1) * EP(1,m+1); + Fmem.Rd(m+1) = a(m+1) * Fmem.Rd(m+1) + EP(1,m+1) * EP(1,m+1); + if (Fmem.Rd(m+1) <= 0 | Fmem.Rn(m+1) <= 0) + error ('>>> PQadap: Rd or Rn is zero'); + end + if (Fmem.Rn(m+1) >= Fmem.Rd(m+1)) + R(1,m+1) = 1; + R(2,m+1) = Fmem.Rd(m+1) / Fmem.Rn(m+1); + else + R(1,m+1) = Fmem.Rn(m+1) / Fmem.Rd(m+1); + R(2,m+1) = 1; + end +end + +% Average the correction factors over M channels and smooth with time +% Create spectrally adapted patterns +for (m = 0:Nc-1) + iL = max (m - M1, 0); + iU = min (m + M2, Nc-1); + s1 = 0; + s2 = 0; + for (i = iL:iU) + s1 = s1 + R(1,i+1); + s2 = s2 + R(2,i+1); + end + Fmem.PC(1,m+1) = a(m+1) * Fmem.PC(1,m+1) + b(m+1) * s1 / (iU-iL+1); + Fmem.PC(2,m+1) = a(m+1) * Fmem.PC(2,m+1) + b(m+1) * s2 / (iU-iL+1); + + % Final correction factor => spectrally adapted patterns + EP(1,m+1) = EP(1,m+1) * Fmem.PC(1,m+1); + EP(2,m+1) = EP(2,m+1) * Fmem.PC(2,m+1); +end + +%-------------------------------------- +function [M1, M2] = PQ_M1M2 (Version, Model) +% Return band averaging parameters + +if (strcmp (Version, 'Basic')) + M1 = 3; + M2 = 4; +elseif (strcmp (Version, 'Advanced')) + if (strcmp (Model, 'FFT')) + M1 = 1; + M2 = 2; + else + M1 = 1; + M2 = 1; + end +end diff --git a/PQevalAudio/Patt/PQloud.m b/PQevalAudio/Patt/PQloud.m new file mode 100644 index 0000000..d46f40d --- /dev/null +++ b/PQevalAudio/Patt/PQloud.m @@ -0,0 +1,53 @@ +function Ntot = PQloud (Ehs, Ver, Mod) +% Calculate the loudness + +% P. Kabal $Revision: 1.1 $ $Date: 2003/12/07 13:35:09 $ + +e = 0.23; + +persistent Nc s Et Ets Version Model + +if (~strcmp (Ver, Version) | ~strcmp (Mod, Model)) + Version = Ver; + Model = Mod; + if (strcmp (Model, 'FFT')) + [Nc, fc] = PQCB (Version); + c = 1.07664; + else + [Nc, fc] = PQFB; + c = 1.26539; + end + E0 = 1e4; + Et = PQ_enThresh (fc); + s = PQ_exIndex (fc); + for (m = 0:Nc-1) + Ets(m+1) = c * (Et(m+1) / (s(m+1) * E0))^e; + end +end + +sN = 0; +for (m = 0:Nc-1) + Nm = Ets(m+1) * ((1 - s(m+1) + s(m+1) * Ehs(m+1) / Et(m+1))^e - 1); + sN = sN + max(Nm, 0); +end +Ntot = (24 / Nc) * sN; + +%==================== +function s = PQ_exIndex (f) +% Excitation index + +N = length (f); +for (m = 0:N-1) + sdB = -2 - 2.05 * atan(f(m+1) / 4000) - 0.75 * atan((f(m+1) / 1600)^2); + s(m+1) = 10^(sdB / 10); +end + +%-------------------- +function Et = PQ_enThresh (f) +% Excitation threshold + +N = length (f); +for (m = 0:N-1) + EtdB = 3.64 * (f(m+1) / 1000)^(-0.8); + Et(m+1) = 10^(EtdB / 10); +end diff --git a/PQevalAudio/Patt/PQmodPatt.m b/PQevalAudio/Patt/PQmodPatt.m new file mode 100644 index 0000000..d7a5ee4 --- /dev/null +++ b/PQevalAudio/Patt/PQmodPatt.m @@ -0,0 +1,34 @@ +function [M, ERavg, Fmem] = PQmodPatt (Es, Fmem) +% Modulation pattern processing + +% P. Kabal $Revision: 1.1 $ $Date: 2003/12/07 13:35:09 $ + +persistent Nc a b Fss + +if (isempty (Nc)) + Fs = 48000; + NF = 2048; + Fss = Fs / (NF/2); + [Nc, fc] = PQCB ('Basic'); + t100 = 0.050; + t0 = 0.008; + [a, b] = PQtConst (t100, t0, fc, Fss); +end + +% Allocate memory +M = zeros (2, Nc); + +e = 0.3; +for (i = 1:2) + for (m = 0:Nc-1) + Ee = Es(i,m+1)^e; + Fmem.DE(i,m+1) = a(m+1) * Fmem.DE(i,m+1) ... + + b(m+1) * Fss * abs (Ee - Fmem.Ese(i,m+1)); + Fmem.Eavg(i,m+1) = a(m+1) * Fmem.Eavg(i,m+1) + b(m+1) * Ee; + Fmem.Ese(i,m+1) = Ee; + + M(i,m+1) = Fmem.DE(i,m+1) / (1 + Fmem.Eavg(i,m+1)/0.3); + end +end + +ERavg = Fmem.Eavg(1,:); diff --git a/conf/conf.yaml b/conf/conf.yaml new file mode 100644 index 0000000..9a33608 --- /dev/null +++ b/conf/conf.yaml @@ -0,0 +1,112 @@ +defaults: + - dset: dataset + - hydra/job_logging: colorlog + - hydra/hydra_logging: colorlog + +path_experiment: "non-specified_name" #there should be a better way to do this +tensorboard_logs: "/scratch/work/molinee2/tensorboard_logs/unet_historical" #path with tensorboard +# Dataset related +fs: 44100 #default is 44100, better NEVER change +seg_len_s_train: 5 #length of the train (and val) segments in seconds +freq_inference: 10 #we do inference after * epochs +seg_len_s_test: 15 #lenum_test_segments: 10 #number of test segments (inferenced every epoch) +num_real_test_segments: 5 #number of real recordings inferenced every epoch +num_test_segments: 15 +buffer_size: 1000 # buffer size for shuffling datasets (train and val) +# Dataset Augmentation +overlap: 0 #overlap when extracting audio segments, default is 0, augment if more data is needed + +# Logging and printing, and does not impact training +#device: cuda +verbose: 0 +use_tensorboard: True +use_soft_denoising: False + + +num_workers: 10 + + +# Checkpointing, by default automatically load last checkpoint +checkpoint: true +continue_from: '' # Path the a checkpoint.th file to start from. + # this is not used in the name of the experiment! + # so use a dummy=something not to mixup experiments. +continue_best: false # continue from best, not last state if continue_from is set. +only_inference: false + +# Optimization related +optim: adam +lr: 1e-4 #used +variable_lr: True +beta1: 0.5 #used +beta2: 0.9 #used + +loss: "mae" #choose loss: + +epochs: 73 +batch_size: 16 +val_take: -1 +steps_per_epoch: 1000 + + +sp: + method: "wiener" + +#STFT parameteres +stft: + win_size: 2048 #STFT window size + hop_size: 512 + +#inference param +inference: + audio: None +# Models +model: unet # either demucs or dwave +unet: + activation: "elu" + use_csff: False + use_SAM: True + use_cam: False + use_fam: False + use_fencoding: True + use_tdf: False + use_alttdfs: False + num_tfc: 3 + num_stages: 3 + depth: 6 + f_dim: 1025 #hardcoded, depends on the stft window + + + +# Hydra config +hydra: + job: + config: + # configuration for the ${hydra.job.override_dirname} runtime variable + override_dirname: + kv_sep: '=' + item_sep: ',' + # Remove all paths, as the / in them would mess up things + # Remove params that would not impact the training itself + # Remove all slurm and submit params. + # This is ugly I know... + exclude_keys: ['path_experiment', + 'hydra.job_logging.handles.file.filename'] + job_logging: + handlers: + file: + class: logging.FileHandler + mode: w + formatter: colorlog + filename: trainer.log + console: + class: logging.StreamHandler + formatter: colorlog + stream: ext://sys.stderr + + hydra_logging: + handlers: + console: + class: logging.StreamHandler + formatter: colorlog + stream: ext://sys.stderr diff --git a/conf/dset/dataset.yaml b/conf/dset/dataset.yaml new file mode 100644 index 0000000..2dc6d94 --- /dev/null +++ b/conf/dset/dataset.yaml @@ -0,0 +1,13 @@ +dset: + path_music_train: [ "/scratch/work/molinee2/datasets/MusicNet_with_opera/train"] + path_music_validation: ["/scratch/work/molinee2/datasets/MusicNet_with_opera/validation"] + path_music_test: ["/scratch/work/molinee2/datasets/MusicNet_extended/test"] + + path_piano_test: ["/scratch/work/molinee2/datasets/MusicNet_with_opera/test/piano"] + path_strings_test: ["/scratch/work/molinee2/datasets/MusicNet_with_opera/test/strings"] + path_orchestra_test: ["/scratch/work/molinee2/datasets/MusicNet_with_opera/test/orchestra"] + path_opera_test: ["/scratch/work/molinee2/datasets/MusicNet_with_opera/test/opera"] + + path_noise: "/scratch/work/molinee2/datasets/noise_set_refined_last" #path with noise samples. There must be an info.csv file there + path_recordings: "/scratch/work/molinee2/datasets/real_noisy_data_test" #path with real noisy recordings. There must be an audio_files.txt file there + diff --git a/experiments/trained_model/noti b/experiments/trained_model/noti new file mode 100644 index 0000000..e69de29