libcm is a C development framework with an emphasis on audio signal processing applications.
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

cmProc2.c 165KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095109610971098109911001101110211031104110511061107110811091110111111121113111411151116111711181119112011211122112311241125112611271128112911301131113211331134113511361137113811391140114111421143114411451146114711481149115011511152115311541155115611571158115911601161116211631164116511661167116811691170117111721173117411751176117711781179118011811182118311841185118611871188118911901191119211931194119511961197119811991200120112021203120412051206120712081209121012111212121312141215121612171218121912201221122212231224122512261227122812291230123112321233123412351236123712381239124012411242124312441245124612471248124912501251125212531254125512561257125812591260126112621263126412651266126712681269127012711272127312741275127612771278127912801281128212831284128512861287128812891290129112921293129412951296129712981299130013011302130313041305130613071308130913101311131213131314131513161317131813191320132113221323132413251326132713281329133013311332133313341335133613371338133913401341134213431344134513461347134813491350135113521353135413551356135713581359136013611362136313641365136613671368136913701371137213731374137513761377137813791380138113821383138413851386138713881389139013911392139313941395139613971398139914001401140214031404140514061407140814091410141114121413141414151416141714181419142014211422142314241425142614271428142914301431143214331434143514361437143814391440144114421443144414451446144714481449145014511452145314541455145614571458145914601461146214631464146514661467146814691470147114721473147414751476147714781479148014811482148314841485148614871488148914901491149214931494149514961497149814991500150115021503150415051506150715081509151015111512151315141515151615171518151915201521152215231524152515261527152815291530153115321533153415351536153715381539154015411542154315441545154615471548154915501551155215531554155515561557155815591560156115621563156415651566156715681569157015711572157315741575157615771578157915801581158215831584158515861587158815891590159115921593159415951596159715981599160016011602160316041605160616071608160916101611161216131614161516161617161816191620162116221623162416251626162716281629163016311632163316341635163616371638163916401641164216431644164516461647164816491650165116521653165416551656165716581659166016611662166316641665166616671668166916701671167216731674167516761677167816791680168116821683168416851686168716881689169016911692169316941695169616971698169917001701170217031704170517061707170817091710171117121713171417151716171717181719172017211722172317241725172617271728172917301731173217331734173517361737173817391740174117421743174417451746174717481749175017511752175317541755175617571758175917601761176217631764176517661767176817691770177117721773177417751776177717781779178017811782178317841785178617871788178917901791179217931794179517961797179817991800180118021803180418051806180718081809181018111812181318141815181618171818181918201821182218231824182518261827182818291830183118321833183418351836183718381839184018411842184318441845184618471848184918501851185218531854185518561857185818591860186118621863186418651866186718681869187018711872187318741875187618771878187918801881188218831884188518861887188818891890189118921893189418951896189718981899190019011902190319041905190619071908190919101911191219131914191519161917191819191920192119221923192419251926192719281929193019311932193319341935193619371938193919401941194219431944194519461947194819491950195119521953195419551956195719581959196019611962196319641965196619671968196919701971197219731974197519761977197819791980198119821983198419851986198719881989199019911992199319941995199619971998199920002001200220032004200520062007200820092010201120122013201420152016201720182019202020212022202320242025202620272028202920302031203220332034203520362037203820392040204120422043204420452046204720482049205020512052205320542055205620572058205920602061206220632064206520662067206820692070207120722073207420752076207720782079208020812082208320842085208620872088208920902091209220932094209520962097209820992100210121022103210421052106210721082109211021112112211321142115211621172118211921202121212221232124212521262127212821292130213121322133213421352136213721382139214021412142214321442145214621472148214921502151215221532154215521562157215821592160216121622163216421652166216721682169217021712172217321742175217621772178217921802181218221832184218521862187218821892190219121922193219421952196219721982199220022012202220322042205220622072208220922102211221222132214221522162217221822192220222122222223222422252226222722282229223022312232223322342235223622372238223922402241224222432244224522462247224822492250225122522253225422552256225722582259226022612262226322642265226622672268226922702271227222732274227522762277227822792280228122822283228422852286228722882289229022912292229322942295229622972298229923002301230223032304230523062307230823092310231123122313231423152316231723182319232023212322232323242325232623272328232923302331233223332334233523362337233823392340234123422343234423452346234723482349235023512352235323542355235623572358235923602361236223632364236523662367236823692370237123722373237423752376237723782379238023812382238323842385238623872388238923902391239223932394239523962397239823992400240124022403240424052406240724082409241024112412241324142415241624172418241924202421242224232424242524262427242824292430243124322433243424352436243724382439244024412442244324442445244624472448244924502451245224532454245524562457245824592460246124622463246424652466246724682469247024712472247324742475247624772478247924802481248224832484248524862487248824892490249124922493249424952496249724982499250025012502250325042505250625072508250925102511251225132514251525162517251825192520252125222523252425252526252725282529253025312532253325342535253625372538253925402541254225432544254525462547254825492550255125522553255425552556255725582559256025612562256325642565256625672568256925702571257225732574257525762577257825792580258125822583258425852586258725882589259025912592259325942595259625972598259926002601260226032604260526062607260826092610261126122613261426152616261726182619262026212622262326242625262626272628262926302631263226332634263526362637263826392640264126422643264426452646264726482649265026512652265326542655265626572658265926602661266226632664266526662667266826692670267126722673267426752676267726782679268026812682268326842685268626872688268926902691269226932694269526962697269826992700270127022703270427052706270727082709271027112712271327142715271627172718271927202721272227232724272527262727272827292730273127322733273427352736273727382739274027412742274327442745274627472748274927502751275227532754275527562757275827592760276127622763276427652766276727682769277027712772277327742775277627772778277927802781278227832784278527862787278827892790279127922793279427952796279727982799280028012802280328042805280628072808280928102811281228132814281528162817281828192820282128222823282428252826282728282829283028312832283328342835283628372838283928402841284228432844284528462847284828492850285128522853285428552856285728582859286028612862286328642865286628672868286928702871287228732874287528762877287828792880288128822883288428852886288728882889289028912892289328942895289628972898289929002901290229032904290529062907290829092910291129122913291429152916291729182919292029212922292329242925292629272928292929302931293229332934293529362937293829392940294129422943294429452946294729482949295029512952295329542955295629572958295929602961296229632964296529662967296829692970297129722973297429752976297729782979298029812982298329842985298629872988298929902991299229932994299529962997299829993000300130023003300430053006300730083009301030113012301330143015301630173018301930203021302230233024302530263027302830293030303130323033303430353036303730383039304030413042304330443045304630473048304930503051305230533054305530563057305830593060306130623063306430653066306730683069307030713072307330743075307630773078307930803081308230833084308530863087308830893090309130923093309430953096309730983099310031013102310331043105310631073108310931103111311231133114311531163117311831193120312131223123312431253126312731283129313031313132313331343135313631373138313931403141314231433144314531463147314831493150315131523153315431553156315731583159316031613162316331643165316631673168316931703171317231733174317531763177317831793180318131823183318431853186318731883189319031913192319331943195319631973198319932003201320232033204320532063207320832093210321132123213321432153216321732183219322032213222322332243225322632273228322932303231323232333234323532363237323832393240324132423243324432453246324732483249325032513252325332543255325632573258325932603261326232633264326532663267326832693270327132723273327432753276327732783279328032813282328332843285328632873288328932903291329232933294329532963297329832993300330133023303330433053306330733083309331033113312331333143315331633173318331933203321332233233324332533263327332833293330333133323333333433353336333733383339334033413342334333443345334633473348334933503351335233533354335533563357335833593360336133623363336433653366336733683369337033713372337333743375337633773378337933803381338233833384338533863387338833893390339133923393339433953396339733983399340034013402340334043405340634073408340934103411341234133414341534163417341834193420342134223423342434253426342734283429343034313432343334343435343634373438343934403441344234433444344534463447344834493450345134523453345434553456345734583459346034613462346334643465346634673468346934703471347234733474347534763477347834793480348134823483348434853486348734883489349034913492349334943495349634973498349935003501350235033504350535063507350835093510351135123513351435153516351735183519352035213522352335243525352635273528352935303531353235333534353535363537353835393540354135423543354435453546354735483549355035513552355335543555355635573558355935603561356235633564356535663567356835693570357135723573357435753576357735783579358035813582358335843585358635873588358935903591359235933594359535963597359835993600360136023603360436053606360736083609361036113612361336143615361636173618361936203621362236233624362536263627362836293630363136323633363436353636363736383639364036413642364336443645364636473648364936503651365236533654365536563657365836593660366136623663366436653666366736683669367036713672367336743675367636773678367936803681368236833684368536863687368836893690369136923693369436953696369736983699370037013702370337043705370637073708370937103711371237133714371537163717371837193720372137223723372437253726372737283729373037313732373337343735373637373738373937403741374237433744374537463747374837493750375137523753375437553756375737583759376037613762376337643765376637673768376937703771377237733774377537763777377837793780378137823783378437853786378737883789379037913792379337943795379637973798379938003801380238033804380538063807380838093810381138123813381438153816381738183819382038213822382338243825382638273828382938303831383238333834383538363837383838393840384138423843384438453846384738483849385038513852385338543855385638573858385938603861386238633864386538663867386838693870387138723873387438753876387738783879388038813882388338843885388638873888388938903891389238933894389538963897389838993900390139023903390439053906390739083909391039113912391339143915391639173918391939203921392239233924392539263927392839293930393139323933393439353936393739383939394039413942394339443945394639473948394939503951395239533954395539563957395839593960396139623963396439653966396739683969397039713972397339743975397639773978397939803981398239833984398539863987398839893990399139923993399439953996399739983999400040014002400340044005400640074008400940104011401240134014401540164017401840194020402140224023402440254026402740284029403040314032403340344035403640374038403940404041404240434044404540464047404840494050405140524053405440554056405740584059406040614062406340644065406640674068406940704071407240734074407540764077407840794080408140824083408440854086408740884089409040914092409340944095409640974098409941004101410241034104410541064107410841094110411141124113411441154116411741184119412041214122412341244125412641274128412941304131413241334134413541364137413841394140414141424143414441454146414741484149415041514152415341544155415641574158415941604161416241634164416541664167416841694170417141724173417441754176417741784179418041814182418341844185418641874188418941904191419241934194419541964197419841994200420142024203420442054206420742084209421042114212421342144215421642174218421942204221422242234224422542264227422842294230423142324233423442354236423742384239424042414242424342444245424642474248424942504251425242534254425542564257425842594260426142624263426442654266426742684269427042714272427342744275427642774278427942804281428242834284428542864287428842894290429142924293429442954296429742984299430043014302430343044305430643074308430943104311431243134314431543164317431843194320432143224323432443254326432743284329433043314332433343344335433643374338433943404341434243434344434543464347434843494350435143524353435443554356435743584359436043614362436343644365436643674368436943704371437243734374437543764377437843794380438143824383438443854386438743884389439043914392439343944395439643974398439944004401440244034404440544064407440844094410441144124413441444154416441744184419442044214422442344244425442644274428442944304431443244334434443544364437443844394440444144424443444444454446444744484449445044514452445344544455445644574458445944604461446244634464446544664467446844694470447144724473447444754476447744784479448044814482448344844485448644874488448944904491449244934494449544964497449844994500450145024503450445054506450745084509451045114512451345144515451645174518451945204521452245234524452545264527452845294530453145324533453445354536453745384539454045414542454345444545454645474548454945504551455245534554455545564557455845594560456145624563456445654566456745684569457045714572457345744575457645774578457945804581458245834584458545864587458845894590459145924593459445954596459745984599460046014602460346044605460646074608460946104611461246134614461546164617461846194620462146224623462446254626462746284629463046314632463346344635463646374638463946404641464246434644464546464647464846494650465146524653465446554656465746584659466046614662466346644665466646674668466946704671467246734674467546764677467846794680468146824683468446854686468746884689469046914692469346944695469646974698469947004701470247034704470547064707470847094710471147124713471447154716471747184719472047214722472347244725472647274728472947304731473247334734473547364737473847394740474147424743474447454746474747484749475047514752475347544755475647574758475947604761476247634764476547664767476847694770477147724773477447754776477747784779478047814782478347844785478647874788478947904791479247934794479547964797479847994800480148024803480448054806480748084809481048114812481348144815481648174818481948204821482248234824482548264827482848294830483148324833483448354836483748384839484048414842484348444845484648474848484948504851485248534854485548564857485848594860486148624863486448654866486748684869487048714872487348744875487648774878487948804881488248834884488548864887488848894890489148924893489448954896489748984899490049014902490349044905490649074908490949104911491249134914491549164917491849194920492149224923492449254926492749284929493049314932493349344935493649374938493949404941494249434944494549464947494849494950495149524953495449554956495749584959496049614962496349644965496649674968496949704971497249734974497549764977497849794980498149824983498449854986498749884989499049914992499349944995499649974998499950005001500250035004500550065007500850095010501150125013501450155016501750185019502050215022502350245025502650275028502950305031503250335034503550365037503850395040504150425043504450455046504750485049505050515052505350545055505650575058505950605061506250635064506550665067506850695070507150725073507450755076507750785079508050815082508350845085508650875088508950905091509250935094509550965097509850995100510151025103510451055106510751085109511051115112511351145115511651175118511951205121512251235124512551265127512851295130513151325133513451355136513751385139514051415142514351445145514651475148514951505151515251535154515551565157515851595160516151625163516451655166516751685169517051715172517351745175517651775178517951805181518251835184518551865187518851895190519151925193519451955196519751985199520052015202520352045205520652075208520952105211521252135214521552165217521852195220522152225223522452255226522752285229523052315232523352345235523652375238523952405241524252435244524552465247524852495250525152525253525452555256525752585259526052615262526352645265526652675268526952705271527252735274527552765277527852795280528152825283528452855286528752885289529052915292529352945295529652975298529953005301530253035304530553065307530853095310531153125313531453155316531753185319532053215322532353245325532653275328532953305331533253335334533553365337533853395340534153425343534453455346534753485349535053515352535353545355535653575358535953605361536253635364536553665367536853695370537153725373537453755376537753785379538053815382538353845385538653875388538953905391539253935394539553965397539853995400540154025403540454055406540754085409541054115412541354145415541654175418541954205421542254235424542554265427542854295430543154325433543454355436543754385439544054415442544354445445544654475448544954505451545254535454545554565457545854595460546154625463546454655466546754685469547054715472547354745475547654775478547954805481548254835484548554865487548854895490549154925493549454955496549754985499550055015502550355045505550655075508550955105511551255135514551555165517551855195520552155225523552455255526552755285529553055315532553355345535553655375538553955405541554255435544554555465547554855495550555155525553555455555556555755585559556055615562556355645565556655675568556955705571557255735574557555765577557855795580558155825583558455855586558755885589559055915592559355945595559655975598559956005601560256035604560556065607560856095610561156125613561456155616561756185619562056215622562356245625562656275628562956305631563256335634563556365637563856395640564156425643564456455646564756485649565056515652565356545655565656575658565956605661566256635664566556665667566856695670567156725673567456755676567756785679568056815682568356845685568656875688568956905691569256935694569556965697569856995700570157025703570457055706570757085709571057115712571357145715571657175718571957205721572257235724572557265727572857295730573157325733573457355736573757385739574057415742574357445745574657475748574957505751575257535754575557565757575857595760576157625763576457655766576757685769577057715772577357745775577657775778577957805781578257835784578557865787578857895790579157925793579457955796579757985799580058015802580358045805580658075808580958105811581258135814581558165817581858195820582158225823582458255826582758285829583058315832583358345835583658375838583958405841584258435844584558465847584858495850585158525853585458555856585758585859586058615862586358645865586658675868586958705871587258735874587558765877587858795880588158825883588458855886588758885889589058915892589358945895589658975898589959005901590259035904590559065907590859095910591159125913591459155916591759185919592059215922592359245925592659275928592959305931593259335934593559365937593859395940594159425943594459455946594759485949595059515952595359545955595659575958595959605961596259635964596559665967596859695970597159725973597459755976597759785979598059815982598359845985598659875988598959905991599259935994599559965997599859996000600160026003600460056006600760086009601060116012601360146015601660176018601960206021602260236024602560266027602860296030603160326033603460356036603760386039604060416042604360446045604660476048604960506051605260536054605560566057605860596060606160626063606460656066606760686069607060716072607360746075607660776078607960806081608260836084608560866087608860896090609160926093609460956096609760986099610061016102610361046105610661076108610961106111611261136114611561166117611861196120612161226123612461256126612761286129613061316132613361346135613661376138613961406141614261436144614561466147614861496150615161526153615461556156615761586159616061616162616361646165616661676168616961706171617261736174617561766177617861796180618161826183618461856186618761886189619061916192619361946195619661976198619962006201620262036204620562066207620862096210621162126213621462156216621762186219622062216222622362246225622662276228622962306231623262336234623562366237623862396240624162426243624462456246624762486249625062516252625362546255625662576258625962606261626262636264626562666267626862696270627162726273627462756276627762786279628062816282628362846285628662876288628962906291629262936294629562966297629862996300630163026303630463056306630763086309631063116312631363146315631663176318631963206321632263236324632563266327632863296330633163326333633463356336633763386339634063416342634363446345634663476348634963506351635263536354635563566357635863596360636163626363636463656366636763686369637063716372637363746375637663776378637963806381638263836384638563866387638863896390639163926393639463956396639763986399640064016402640364046405640664076408640964106411641264136414641564166417641864196420642164226423642464256426642764286429643064316432643364346435643664376438643964406441644264436444644564466447644864496450645164526453645464556456645764586459646064616462646364646465646664676468646964706471647264736474647564766477647864796480648164826483648464856486648764886489649064916492649364946495649664976498649965006501650265036504650565066507650865096510651165126513651465156516651765186519652065216522652365246525652665276528652965306531653265336534653565366537653865396540654165426543654465456546654765486549655065516552655365546555655665576558655965606561656265636564656565666567656865696570657165726573657465756576657765786579658065816582658365846585658665876588658965906591659265936594659565966597659865996600660166026603660466056606660766086609661066116612661366146615661666176618661966206621662266236624662566266627662866296630
  1. #include "cmPrefix.h"
  2. #include "cmGlobal.h"
  3. #include "cmRpt.h"
  4. #include "cmErr.h"
  5. #include "cmCtx.h"
  6. #include "cmMem.h"
  7. #include "cmMallocDebug.h"
  8. #include "cmLinkedHeap.h"
  9. #include "cmSymTbl.h"
  10. #include "cmFloatTypes.h"
  11. #include "cmComplexTypes.h"
  12. #include "cmFile.h"
  13. #include "cmFileSys.h"
  14. #include "cmProcObj.h"
  15. #include "cmProcTemplate.h"
  16. #include "cmAudioFile.h"
  17. #include "cmMath.h"
  18. #include "cmProc.h"
  19. #include "cmVectOps.h"
  20. #include "cmKeyboard.h"
  21. #include "cmGnuPlot.h"
  22. #include "cmTime.h"
  23. #include "cmMidi.h"
  24. #include "cmProc2.h"
  25. //------------------------------------------------------------------------------------------------------------
  26. cmArray* cmArrayAllocate( cmCtx* c, cmArray* ap, unsigned eleCnt, unsigned eleByteCnt, unsigned flags )
  27. {
  28. cmArray* p = cmObjAlloc( cmArray, c, ap );
  29. if( eleCnt > 0 && eleByteCnt > 0 )
  30. if( cmArrayInit(p, eleCnt, eleByteCnt, flags ) != cmOkRC )
  31. cmArrayFree(&p);
  32. return cmOkRC;
  33. }
  34. cmRC_t cmArrayFree( cmArray** pp )
  35. {
  36. cmRC_t rc = cmOkRC;
  37. cmArray* p = *pp;
  38. if( pp == NULL || *pp == NULL )
  39. return cmOkRC;
  40. if((rc = cmArrayFinal(p)) != cmOkRC )
  41. return rc;
  42. cmMemPtrFree(&p->ptr);
  43. cmObjFree(pp);
  44. return rc;
  45. }
  46. cmRC_t cmArrayInit( cmArray* p, unsigned eleCnt, unsigned eleByteCnt, unsigned flags )
  47. {
  48. cmRC_t rc = cmOkRC;
  49. if((rc = cmArrayFinal(p)) != cmOkRC )
  50. return rc;
  51. p->allocByteCnt = eleCnt * eleByteCnt;
  52. p->ptr = cmIsFlag(flags,kZeroArrayFl) ? cmMemResizeZ( char, p->ptr, p->allocByteCnt ) : cmMemResize( char, p->ptr, p->allocByteCnt );
  53. p->eleCnt = eleCnt;
  54. p->eleByteCnt = eleByteCnt;
  55. return rc;
  56. }
  57. cmRC_t cmArrayFinal( cmArray* p )
  58. { return cmOkRC; }
  59. char* cmArrayReallocDestroy(cmArray* p, unsigned newEleCnt, unsigned newEleByteCnt, unsigned flags )
  60. {
  61. // if memory is expanding
  62. if( newEleCnt * newEleByteCnt > p->allocByteCnt )
  63. cmArrayInit( p, newEleCnt, newEleByteCnt, flags );
  64. else
  65. {
  66. // ... otherwise memory is contrcmting
  67. p->eleCnt = newEleCnt;
  68. p->eleByteCnt = newEleByteCnt;
  69. if( cmIsFlag( flags, kZeroArrayFl ))
  70. memset(p->ptr,0,p->eleByteCnt);
  71. }
  72. return p->ptr;
  73. }
  74. void cmArrayReallocDestroyV(cmArray* p, int eleByteCnt,unsigned flags, ... )
  75. {
  76. unsigned i;
  77. unsigned eleCnt = 0;
  78. unsigned argCnt = 0;
  79. va_list vl0,vl1;
  80. assert(eleByteCnt>0);
  81. va_start(vl0,flags);
  82. va_copy(vl1,vl0);
  83. while( va_arg(vl0,void**) != NULL )
  84. {
  85. int argEleCnt = va_arg(vl0,int);
  86. assert(argEleCnt>0);
  87. eleCnt += argEleCnt;
  88. ++argCnt;
  89. }
  90. va_end(vl0);
  91. char* a = cmArrayReallocDestroy(p,eleCnt,eleByteCnt,flags);
  92. for(i=0; i<argCnt; ++i)
  93. {
  94. void** vp = va_arg(vl1,void**);
  95. unsigned n = va_arg(vl1,unsigned);
  96. *vp = a;
  97. a += n*eleByteCnt;
  98. }
  99. va_end(vl1);
  100. }
  101. char* cmArrayReallocPreserve(cmArray* p, unsigned newEleCnt, unsigned newEleByteCnt, unsigned flags )
  102. {
  103. unsigned cn = p->eleCnt * p->eleByteCnt;
  104. unsigned dn = newEleCnt * newEleByteCnt;
  105. if( dn > p->allocByteCnt )
  106. p->allocByteCnt = dn;
  107. p->ptr = cmIsFlag(flags,kZeroArrayFl ) ? cmMemResizePZ( char, p->ptr, cn ) : cmMemResizeP( char, p->ptr, cn);
  108. p->eleCnt = newEleCnt;
  109. p->eleByteCnt= newEleByteCnt;
  110. return p->ptr;
  111. }
  112. //------------------------------------------------------------------------------------------------------------
  113. cmAudioFileWr* cmAudioFileWrAlloc( cmCtx* c, cmAudioFileWr* ap, unsigned procSmpCnt, const char* fn, double srate, unsigned chCnt, unsigned bitsPerSample )
  114. {
  115. cmAudioFileWr* p = cmObjAlloc( cmAudioFileWr, c, ap );
  116. if( cmAudioFileWrInit( p, procSmpCnt, fn, srate, chCnt, bitsPerSample ) != cmOkRC )
  117. cmObjFree(&p);
  118. return p;
  119. }
  120. cmRC_t cmAudioFileWrFree( cmAudioFileWr** pp )
  121. {
  122. cmRC_t rc = cmOkRC;
  123. if( pp != NULL && *pp != NULL )
  124. {
  125. cmAudioFileWr* p = *pp;
  126. if((rc = cmAudioFileWrFinal(p)) == cmOkRC )
  127. {
  128. cmMemPtrFree(&p->bufV);
  129. cmMemPtrFree(&p->fn );
  130. cmObjFree(pp);
  131. }
  132. }
  133. return rc;
  134. }
  135. cmRC_t cmAudioFileWrInit( cmAudioFileWr* p, unsigned procSmpCnt, const char* fn, double srate, unsigned chCnt, unsigned bitsPerSample )
  136. {
  137. cmRC_t rc;
  138. cmRC_t afRC;
  139. if((rc = cmAudioFileWrFinal( p)) != cmOkRC )
  140. return rc;
  141. p->h = cmAudioFileNewCreate( fn, srate, bitsPerSample, chCnt, &afRC, p->obj.err.rpt );
  142. if( afRC != kOkAfRC )
  143. return cmCtxRtCondition( &p->obj, afRC, "Unable to open the audio file:'%s'", fn );
  144. p->bufV = cmMemResize( cmSample_t, p->bufV, procSmpCnt * chCnt );
  145. p->procSmpCnt = procSmpCnt;
  146. p->chCnt = chCnt;
  147. p->curChCnt = 0;
  148. p->fn = cmMemResizeZ( cmChar_t, p->fn, strlen(fn)+1 );
  149. strcpy(p->fn,fn);
  150. return rc;
  151. }
  152. cmRC_t cmAudioFileWrFinal( cmAudioFileWr* p )
  153. {
  154. cmRC_t afRC;
  155. if( p != NULL )
  156. {
  157. if( cmAudioFileIsValid( p->h ) )
  158. if(( afRC = cmAudioFileDelete( &p->h )) != kOkAfRC )
  159. return cmCtxRtCondition( &p->obj, afRC, "Unable to close the audio file:'%s'", p->fn );
  160. }
  161. return cmOkRC;
  162. }
  163. cmRC_t cmAudioFileWrExec( cmAudioFileWr* p, unsigned chIdx, const cmSample_t* sp, unsigned sn )
  164. {
  165. cmRC_t afRC;
  166. assert( sn <= p->procSmpCnt && chIdx < p->chCnt );
  167. cmSample_t* buf = p->bufV + (chIdx * p->procSmpCnt);
  168. cmVOS_Copy( buf, sn, sp);
  169. if( sn < p->procSmpCnt )
  170. cmVOS_Fill( buf+sn, p->procSmpCnt-sn, 0 );
  171. p->curChCnt++;
  172. if( p->curChCnt == p->chCnt )
  173. {
  174. p->curChCnt = 0;
  175. cmSample_t* bufPtrPtr[ p->chCnt ];
  176. unsigned i = 0;
  177. for(i=0; i<p->chCnt; ++i)
  178. bufPtrPtr[i] = p->bufV + (i*p->procSmpCnt);
  179. if((afRC = cmAudioFileWriteSample( p->h, p->procSmpCnt, p->chCnt, bufPtrPtr )) != kOkAfRC )
  180. return cmCtxRtCondition( &p->obj, afRC, "Write failed on audio file:'%s'", p->fn );
  181. }
  182. return cmOkRC;
  183. }
  184. void cmAudioFileWrTest( cmRpt_t* rpt, cmLHeapH_t lhH, cmSymTblH_t stH )
  185. {
  186. const char* fn = "/home/kevin/src/cm/test0.aif";
  187. double durSecs = 10;
  188. double srate = 44100;
  189. unsigned chCnt = 2;
  190. unsigned bitsPerSmp = 16;
  191. unsigned procSmpCnt = 64;
  192. double hz = 1000;
  193. unsigned overToneCnt= 1;
  194. unsigned smpCnt = durSecs * srate;
  195. unsigned i;
  196. cmCtx* c = cmCtxAlloc( NULL, rpt, lhH, stH );
  197. cmSigGen* sgp = cmSigGenAlloc( c, NULL, procSmpCnt, srate, kWhiteWfId, hz, overToneCnt );
  198. cmAudioFileWr* awp = cmAudioFileWrAlloc( c, NULL, procSmpCnt, fn, srate, chCnt, bitsPerSmp );
  199. for(i=0; i<smpCnt; ++i)
  200. {
  201. cmSigGenExec( sgp );
  202. cmAudioFileWrExec( awp, 0, sgp->outV, sgp->outN );
  203. cmAudioFileWrExec( awp, 1, sgp->outV, sgp->outN );
  204. i += sgp->outN;
  205. }
  206. printf("Frames:%i\n",smpCnt);
  207. cmAudioFileWrFree(&awp);
  208. cmSigGenFree( &sgp );
  209. cmCtxFree(&c);
  210. cmAudioFileReportFn( fn, 0, 20, rpt );
  211. }
  212. //------------------------------------------------------------------------------------------------------------
  213. cmMatrixBuf* cmMatrixBufAllocFile( cmCtx* c, cmMatrixBuf* p, const char* fn )
  214. {
  215. cmRC_t rc;
  216. cmMatrixBuf* op = cmObjAlloc( cmMatrixBuf, c, p );
  217. if( fn != NULL )
  218. if((rc = cmMatrixBufInitFile(op,fn)) != cmOkRC )
  219. cmObjFree(&op);
  220. return op;
  221. }
  222. cmMatrixBuf* cmMatrixBufAllocCopy(cmCtx* c, cmMatrixBuf* p, unsigned rn, unsigned cn, const cmSample_t* sp )
  223. {
  224. cmRC_t rc;
  225. cmMatrixBuf* op = cmObjAlloc( cmMatrixBuf, c, p );
  226. if( sp != NULL && rn > 0 && cn > 0 )
  227. if((rc = cmMatrixBufInitCopy(op,rn,cn,sp)) != cmOkRC )
  228. cmObjFree(&op);
  229. return op;
  230. }
  231. cmMatrixBuf* cmMatrixBufAlloc( cmCtx* c, cmMatrixBuf* p, unsigned rn, unsigned cn )
  232. {
  233. cmRC_t rc;
  234. cmMatrixBuf* op = cmObjAlloc( cmMatrixBuf, c, p );
  235. if( rn > 0 && cn > 0 )
  236. if((rc = cmMatrixBufInit(op,rn,cn)) != cmOkRC )
  237. cmObjFree(&op);
  238. return op;
  239. }
  240. cmRC_t cmMatrixBufFree( cmMatrixBuf** pp )
  241. {
  242. cmRC_t rc = cmOkRC;
  243. if( pp != NULL && *pp != NULL )
  244. {
  245. cmMatrixBuf* p = *pp;
  246. if((rc = cmMatrixBufFinal(p)) == cmOkRC )
  247. {
  248. cmMemPtrFree(&p->bufPtr);
  249. cmObjFree(pp);
  250. }
  251. }
  252. return rc;
  253. }
  254. void _cmMatrixBufGetFileSize( FILE* fp, unsigned* lineCharCntPtr, unsigned* lineCntPtr )
  255. {
  256. *lineCharCntPtr = 0;
  257. *lineCntPtr = 0;
  258. while( !feof(fp) )
  259. {
  260. char ch;
  261. unsigned charCnt = 0;
  262. while( (ch = getc(fp)) != EOF )
  263. {
  264. charCnt++;
  265. if( ch == '\n' )
  266. break;
  267. }
  268. *lineCntPtr += 1;
  269. if(charCnt > *lineCharCntPtr )
  270. *lineCharCntPtr = charCnt;
  271. }
  272. *lineCharCntPtr += 5; // add a safety margin
  273. }
  274. cmRC_t _cmMatrixBufGetMatrixSize( cmObj* op, FILE* fp, unsigned lineCharCnt, unsigned lineCnt, unsigned* rowCntPtr, unsigned * colCntPtr, const char* fn )
  275. {
  276. unsigned i;
  277. char lineBuf[ lineCharCnt + 1 ];
  278. *rowCntPtr = 0;
  279. *colCntPtr = 0;
  280. for(i=0; i<lineCnt; ++i)
  281. {
  282. if(fgets(lineBuf,lineCharCnt,fp)==NULL)
  283. {
  284. // if the last line is blank then return from here
  285. if( feof(fp) )
  286. return cmOkRC;
  287. return cmCtxRtCondition( op, cmSystemErrorRC, "A read error occured on the matrix file:'%s'.",fn);
  288. }
  289. assert( strlen(lineBuf) < lineCharCnt );
  290. char* lp = lineBuf;
  291. char* tp;
  292. // eat any leading white space
  293. while( (*lp) && isspace(*lp) )
  294. ++lp;
  295. // if the line was blank then skip it
  296. if( strlen(lp) == 0 || *lp == '#' )
  297. continue;
  298. (*rowCntPtr) += 1;
  299. unsigned colCnt;
  300. for(colCnt=0; (tp = strtok(lp," ")) != NULL; ++colCnt )
  301. lp = NULL;
  302. if( colCnt > *colCntPtr )
  303. *colCntPtr = colCnt;
  304. }
  305. return cmOkRC;
  306. }
  307. double _cmMatrixBufStrToNum( cmObj* op, const char* cp )
  308. {
  309. double v;
  310. if( sscanf(cp,"%le ",&v) != 1 )
  311. cmCtxRtCondition( op, cmArgAssertRC, "Parse error reading matrix file.");
  312. return v;
  313. }
  314. cmRC_t _cmMatrixBufReadFile(cmObj* op, FILE* fp, cmSample_t* p, unsigned lineCharCnt, unsigned rn, unsigned cn)
  315. {
  316. char lineBuf[ lineCharCnt+1 ];
  317. unsigned ci = 0;
  318. unsigned ri = 0;
  319. while( fgets(lineBuf,lineCharCnt,fp) != NULL )
  320. {
  321. char* lp = lineBuf;
  322. char* tp;
  323. while( (*lp) && isspace(*lp) )
  324. lp++;
  325. if( strlen(lp) == 0 || *lp == '#' )
  326. continue;
  327. for(ci=0; (tp = strtok(lp," ")) != NULL; ++ci )
  328. {
  329. p[ (ci*rn) + ri ] = _cmMatrixBufStrToNum(op,tp); //atof(tp);
  330. lp = NULL;
  331. }
  332. ++ri;
  333. }
  334. return cmOkRC;
  335. }
  336. cmRC_t cmMatrixBufInitFile( cmMatrixBuf* p, const char* fn )
  337. {
  338. cmRC_t rc;
  339. FILE* fp;
  340. unsigned lineCharCnt;
  341. unsigned lineCnt;
  342. unsigned rn;
  343. unsigned cn;
  344. if((fp = fopen(fn,"rt")) == NULL )
  345. return cmCtxRtCondition( &p->obj, cmSystemErrorRC, "Unable to open the matrix file:'%s'", fn );
  346. // get the length of the longest line in the file
  347. _cmMatrixBufGetFileSize(fp,&lineCharCnt,&lineCnt);
  348. rewind(fp);
  349. // get the count of matrix rows and columns
  350. if((rc=_cmMatrixBufGetMatrixSize( &p->obj, fp, lineCharCnt, lineCnt, &rn, &cn, fn )) != cmOkRC )
  351. goto errLabel;
  352. rewind(fp);
  353. // allocate the matrix memory
  354. cmMatrixBufInit(p,rn,cn);
  355. // fill the matrix from the file
  356. rc = _cmMatrixBufReadFile(&p->obj,fp,p->bufPtr,lineCharCnt,rn,cn);
  357. errLabel:
  358. if( rc != cmOkRC )
  359. cmMatrixBufFinal(p);
  360. fclose(fp);
  361. return rc;
  362. }
  363. cmRC_t cmMatrixBufInitCopy( cmMatrixBuf* p, unsigned rn, unsigned cn, const cmSample_t* sp )
  364. {
  365. cmRC_t rc;
  366. if((rc = cmMatrixBufInit(p,rn,cn)) != cmOkRC )
  367. return rc;
  368. cmVOS_Copy(p->bufPtr,(rn*cn),sp);
  369. return rc;
  370. }
  371. cmRC_t cmMatrixBufInit( cmMatrixBuf* p, unsigned rn, unsigned cn )
  372. {
  373. cmRC_t rc;
  374. if((rc = cmMatrixBufFinal(p)) != cmOkRC )
  375. return rc;
  376. p->rn = rn;
  377. p->cn = cn;
  378. p->bufPtr = cmMemResize( cmSample_t, p->bufPtr, rn*cn );
  379. return cmOkRC;
  380. }
  381. cmRC_t cmMatrixBufFinal( cmMatrixBuf* p )
  382. { return cmOkRC; }
  383. cmSample_t* cmMatrixBufColPtr( cmMatrixBuf* p, unsigned ci )
  384. { assert(ci<p->cn); return p->bufPtr + (ci * p->rn); }
  385. cmSample_t* cmMatrixBufRowPtr( cmMatrixBuf* p, unsigned ri )
  386. { assert(ri<p->rn); return p->bufPtr + ri; }
  387. void cmMatrixBufTest( cmRpt_t* rpt, cmLHeapH_t lhH, cmSymTblH_t stH )
  388. {
  389. cmSample_t v[] = {1,2,2,3};
  390. cmCtx* c = cmCtxAlloc(NULL,rpt,lhH,stH);
  391. cmMatrixBuf* mbp = cmMatrixBufAllocFile(c, NULL, "temp.mat" );
  392. cmMatrixBuf* vbp = cmMatrixBufAllocCopy(c, NULL, 4,1,v);
  393. unsigned i;
  394. printf("rn:%i cn:%i\n",mbp->rn,mbp->cn);
  395. //cmVOS_Print( stdout, 10, 10, mbp->bufPtr );
  396. printf("%.1f\n ",cmVOS_Median( cmMatrixBufColPtr(vbp,0),vbp->rn));
  397. for(i=0; i<mbp->cn; ++i)
  398. {
  399. //cmVOS_Print( stdout, 1, mbp->cn, cmMatrixBufColPtr(c,mbp,i) );
  400. printf("%.1f, ",cmVOS_Median( cmMatrixBufColPtr(mbp,i),mbp->rn));
  401. }
  402. printf("\n");
  403. cmMatrixBufFree(&mbp);
  404. cmMatrixBufFree(&vbp);
  405. cmCtxFree(&c);
  406. }
  407. //------------------------------------------------------------------------------------------------------------
  408. cmSigGen* cmSigGenAlloc( cmCtx* c, cmSigGen* p, unsigned procSmpCnt, double srate, unsigned wfId, double fundFrqHz, unsigned overToneCnt )
  409. {
  410. cmSigGen* op = cmObjAlloc( cmSigGen, c, p );
  411. if( procSmpCnt > 0 && srate > 0 && wfId != kInvalidWfId )
  412. if( cmSigGenInit( op, procSmpCnt, srate, wfId, fundFrqHz, overToneCnt ) != cmOkRC )
  413. cmObjFree(&op);
  414. return op;
  415. }
  416. cmRC_t cmSigGenFree( cmSigGen** pp )
  417. {
  418. cmRC_t rc = cmOkRC;
  419. if( pp != NULL && *pp != NULL )
  420. {
  421. cmSigGen* p = *pp;
  422. if((rc = cmSigGenFinal(p)) == cmOkRC )
  423. {
  424. cmMemPtrFree(&p->outV);
  425. cmObjFree(pp);
  426. }
  427. }
  428. return rc;
  429. }
  430. cmRC_t cmSigGenInit( cmSigGen* p, unsigned procSmpCnt, double srate, unsigned wfId, double fundFrqHz, unsigned overToneCnt )
  431. {
  432. assert( srate > 0 && procSmpCnt > 0 );
  433. p->outV = cmMemResize( cmSample_t, p->outV, procSmpCnt );
  434. p->outN = procSmpCnt;
  435. p->wfId = wfId;
  436. p->overToneCnt = overToneCnt;
  437. p->fundFrqHz = fundFrqHz;
  438. p->phase = 0;
  439. p->delaySmp = 0;
  440. p->srate = srate;
  441. return cmOkRC;
  442. }
  443. cmRC_t cmSigGenFinal( cmSigGen* p )
  444. { return cmOkRC; }
  445. cmRC_t cmSigGenExec( cmSigGen* p )
  446. {
  447. switch( p->wfId )
  448. {
  449. case kSineWfId: p->phase = cmVOS_SynthSine( p->outV, p->outN, p->phase, p->srate, p->fundFrqHz ); break;
  450. case kCosWfId: p->phase = cmVOS_SynthCosine( p->outV, p->outN, p->phase, p->srate, p->fundFrqHz ); break;
  451. case kSquareWfId: p->phase = cmVOS_SynthSquare( p->outV, p->outN, p->phase, p->srate, p->fundFrqHz, p->overToneCnt ); break;
  452. case kTriangleWfId: p->phase = cmVOS_SynthTriangle( p->outV, p->outN, p->phase, p->srate, p->fundFrqHz, p->overToneCnt ); break;
  453. case kSawtoothWfId: p->phase = cmVOS_SynthSawtooth( p->outV, p->outN, p->phase, p->srate, p->fundFrqHz, p->overToneCnt ); break;
  454. case kWhiteWfId: cmVOS_Random( p->outV, p->outN, -1.0, 1.0 ); break;
  455. case kPinkWfId: p->delaySmp = cmVOS_SynthPinkNoise(p->outV, p->outN, p->delaySmp ); break;
  456. case kPulseWfId: p->phase = cmVOS_SynthPulseCos( p->outV, p->outN, p->phase, p->srate, p->fundFrqHz, p->overToneCnt ); break;
  457. case kImpulseWfId: p->phase = cmVOS_SynthImpulse( p->outV, p->outN, p->phase, p->srate, p->fundFrqHz ); break;
  458. case kSilenceWfId: cmVOS_Fill( p->outV, p->outN, 0 ); break;
  459. case kPhasorWfId: p->phase = cmVOS_SynthPhasor( p->outV, p->outN, p->phase, p->srate, p->fundFrqHz ); break;
  460. case kSeqWfId: p->phase = cmVOS_Seq( p->outV, p->outN, p->phase, 1 ); break;
  461. case kInvalidWfId:
  462. default:
  463. return cmCtxRtAssertFailed( &p->obj, 0, "Invalid waveform shape.");
  464. }
  465. return cmOkRC;
  466. }
  467. //------------------------------------------------------------------------------------------------------------
  468. cmDelay* cmDelayAlloc( cmCtx* c, cmDelay* ap, unsigned procSmpCnt, unsigned delaySmpCnt )
  469. {
  470. cmDelay* p = cmObjAlloc( cmDelay, c, ap );
  471. if( procSmpCnt > 0 && delaySmpCnt > 0 )
  472. if( cmDelayInit( p,procSmpCnt,delaySmpCnt) != cmOkRC && ap == NULL )
  473. cmObjFree(&p);
  474. return p;
  475. }
  476. cmRC_t cmDelayFree( cmDelay** pp )
  477. {
  478. cmRC_t rc = cmOkRC;
  479. if( pp != NULL && *pp != NULL )
  480. {
  481. cmDelay* p = *pp;
  482. if((rc = cmDelayFinal(*pp)) == cmOkRC )
  483. {
  484. cmMemPtrFree(&p->bufPtr);
  485. cmObjFree(pp);
  486. }
  487. }
  488. return rc;
  489. }
  490. cmRC_t cmDelayInit( cmDelay* p, unsigned procSmpCnt, unsigned delaySmpCnt )
  491. {
  492. p->procSmpCnt = procSmpCnt;
  493. p->delaySmpCnt = delaySmpCnt;
  494. p->bufSmpCnt = delaySmpCnt + procSmpCnt;
  495. p->bufPtr = cmMemResizeZ( cmSample_t, p->bufPtr, p->bufSmpCnt);
  496. p->delayInIdx = 0;
  497. p->outCnt = 1;
  498. p->outV[0] = p->bufPtr;
  499. p->outN[0] = p->procSmpCnt;
  500. p->outV[1] = NULL;
  501. p->outN[1] = 0;
  502. return cmOkRC;
  503. }
  504. cmRC_t cmDelayFinal( cmDelay* p )
  505. { return cmOkRC; }
  506. cmRC_t cmDelayCopyIn( cmDelay* p, const cmSample_t* sp, unsigned sn )
  507. {
  508. assert(sn<=p->procSmpCnt);
  509. unsigned n0 = cmMin(sn,p->bufSmpCnt - p->delayInIdx);
  510. // copy as many samples as possible from the input to the delayInIdx
  511. cmVOS_Copy(p->bufPtr + p->delayInIdx, n0, sp);
  512. p->delayInIdx = (p->delayInIdx + n0) % p->bufSmpCnt;
  513. // if there was not enough room to copy all the samples into the end of the buffer ....
  514. if( n0 < sn )
  515. {
  516. assert( p->delayInIdx == 0 );
  517. // ... then copy the rest to the beginning of the buffer
  518. unsigned n1 = sn - n0;
  519. cmVOS_Copy(p->bufPtr,n1, sp + n0 );
  520. p->delayInIdx = (p->delayInIdx + n1) % p->bufSmpCnt;
  521. }
  522. return cmOkRC;
  523. }
  524. cmRC_t cmDelayAdvance( cmDelay* p, unsigned sn )
  525. {
  526. // advance the output by sn and make sn samples available
  527. int delayOutIdx = ((p->outV[0] - p->bufPtr) + sn) % p->bufSmpCnt;
  528. p->outV[0] = p->bufPtr + delayOutIdx;
  529. p->outN[0] = cmMin(p->bufSmpCnt - delayOutIdx , sn );
  530. p->outCnt = p->outN[0] == sn ? 1 : 2 ;
  531. p->outV[1] = p->outCnt == 1 ? NULL : p->bufPtr;
  532. p->outN[1] = p->outCnt == 1 ? 0 : sn - p->outN[0];
  533. return cmOkRC;
  534. }
  535. cmRC_t cmDelayExec( cmDelay* p, const cmSample_t* sp, unsigned sn, bool bypassFl )
  536. {
  537. cmRC_t rc = cmOkRC;
  538. if( bypassFl )
  539. memcpy(p->outV,sp,sn*sizeof(cmSample_t));
  540. else
  541. {
  542. cmDelayCopyIn(p,sp,sn);
  543. rc = cmDelayAdvance(p,sn);
  544. }
  545. return rc;
  546. }
  547. void cmDelayTest( cmRpt_t* rpt, cmLHeapH_t lhH, cmSymTblH_t stH )
  548. {
  549. cmCtx ctx;
  550. cmDelay delay;
  551. cmSigGen sigGen;
  552. unsigned procCnt = 4;
  553. unsigned procSmpCnt = 5;
  554. unsigned delaySmpCnt = 3;
  555. unsigned i;
  556. cmCtx* c = cmCtxAlloc( &ctx, rpt, lhH, stH );
  557. cmDelay* dlp = cmDelayAlloc( c, &delay, procSmpCnt, delaySmpCnt );
  558. cmSigGen* sgp = cmSigGenAlloc( c, &sigGen, procSmpCnt, 0, kSeqWfId,0, 0);
  559. for(i=0; i<procCnt; ++i)
  560. {
  561. cmSigGenExec(sgp);
  562. cmDelayExec(dlp,sgp->outV,sgp->outN,false);
  563. //cmVOS_Print( c->outFp, 1, sgp->outN, sgp->outV, 5, 0 );
  564. cmCtxPrint(c,"%i %i : ",i,0);
  565. cmVOS_PrintE( rpt, 1, dlp->outN[0], dlp->outV[0] );
  566. if( dlp->outN[1] > 0 )
  567. {
  568. cmCtxPrint(c,"%i %i : ",i,1);
  569. cmVOS_PrintE( rpt, 1, dlp->outN[1], dlp->outV[1] );
  570. }
  571. }
  572. cmSigGenFinal(sgp);
  573. cmDelayFinal(dlp);
  574. cmCtxFinal(c);
  575. }
  576. //------------------------------------------------------------------------------------------------------------
  577. cmFIR* cmFIRAllocKaiser(cmCtx* c, cmFIR* p, unsigned procSmpCnt, double srate, double passHz, double stopHz, double passDb, double stopDb )
  578. {
  579. cmFIR* op = cmObjAlloc( cmFIR, c, p );
  580. if( procSmpCnt > 0 && srate > 0 )
  581. if( cmFIRInitKaiser(op,procSmpCnt,srate,passHz,stopHz,passDb,stopDb) != cmOkRC )
  582. cmObjFree(&op);
  583. return op;
  584. }
  585. cmFIR* cmFIRAllocSinc( cmCtx* c, cmFIR* p, unsigned procSmpCnt, double srate, unsigned sincSmpCnt, double fcHz, unsigned flags )
  586. {
  587. cmFIR* op = cmObjAlloc( cmFIR, c, p );
  588. if( srate > 0 && sincSmpCnt > 0 )
  589. if( cmFIRInitSinc(op,procSmpCnt,srate,sincSmpCnt,fcHz,flags) != cmOkRC )
  590. cmObjFree(&op);
  591. return op;
  592. }
  593. cmRC_t cmFIRFree( cmFIR** pp )
  594. {
  595. cmRC_t rc = cmOkRC;
  596. if( pp != NULL && *pp != NULL)
  597. {
  598. cmFIR* p = *pp;
  599. if((rc = cmFIRFinal(*pp)) != cmOkRC )
  600. {
  601. cmMemPtrFree(&p->coeffV);
  602. cmMemPtrFree(&p->outV);
  603. cmMemPtrFree(&p->delayV);
  604. cmObjFree(pp);
  605. }
  606. }
  607. return rc;
  608. }
  609. cmRC_t cmFIRInitKaiser( cmFIR* p, unsigned procSmpCnt, double srate, double passHz, double stopHz, double passDb, double stopDb )
  610. {
  611. // pass/stop frequencies above nyquist produce incorrect results
  612. assert( passHz <= srate/2 && stopHz<=srate/2);
  613. // based on Orfandis, Introduction to Signal Processing, p.551 Prentice Hall, 1996
  614. double fcHz = (passHz + stopHz) / 2; // fc is half way between passHz and stopHz
  615. double dHz = fabs(stopHz-passHz);
  616. //double signFcmt = stopHz > passHz ? 1 : -1;
  617. // convert ripple spec from db to linear
  618. double dPass = (pow(10,passDb/20)-1) / (pow(10,passDb/20)+1);
  619. double dStop = pow(10,-stopDb/20);
  620. // in prcmtice the ripple must be equal in the stop and pass band - so take the minimum between the two
  621. double d = cmMin(dPass,dStop);
  622. // convert the ripple bcmk to db
  623. double A = -20 * log10(d);
  624. // compute the kaiser alpha coeff
  625. double alpha = 0;
  626. if( A >= 50.0 ) // for ripple > 50
  627. alpha = 0.1102 * (A-8.7);
  628. else // for ripple <= 21
  629. {
  630. if( A > 21 )
  631. alpha = (0.5842 * (pow(A-21.0,0.4))) + (0.07886*(A-21));
  632. }
  633. double D = 0.922;
  634. if( A > 21 )
  635. D = (A - 7.95) / 14.36;
  636. // compute the filter order
  637. unsigned N = (unsigned)(floor(D * srate / dHz) + 1);
  638. //if N is even
  639. if( cmIsEvenU(N) )
  640. N = N + 1;
  641. //printf("fc=%f df=%f dPass=%f dStop=%f d=%f alpha=%f A=%f D=%f N=%i\n",fcHz,dHz,dPass,dStop,d,alpha,A,D,N);
  642. // form an ideal FIR LP impulse response based on a sinc function
  643. cmFIRInitSinc(p,procSmpCnt,srate,N,fcHz,0);
  644. // compute a kaiser function to truncate the sinc
  645. double wnd[ N ];
  646. cmVOD_Kaiser( wnd, N, alpha );
  647. // apply the kaiser window to the sinc function
  648. cmVOD_MultVV( p->coeffV, p->coeffCnt, wnd );
  649. double sum = cmVOD_Sum(p->coeffV,p->coeffCnt);
  650. cmVOD_DivVS(p->coeffV,p->coeffCnt,sum );
  651. //cmVOD_Print(stdout,1,p->coeffCnt,p->coeffV);
  652. return cmOkRC;
  653. }
  654. cmRC_t cmFIRInitSinc( cmFIR* p, unsigned procSmpCnt, double srate, unsigned sincSmpCnt, double fcHz, unsigned flags )
  655. {
  656. cmRC_t rc;
  657. if((rc = cmFIRFinal(p)) != cmOkRC )
  658. return rc;
  659. p->coeffCnt = sincSmpCnt;
  660. p->outV = cmMemResizeZ( cmSample_t, p->outV, procSmpCnt );
  661. p->outN = procSmpCnt;
  662. p->coeffV = cmMemResizeZ( double, p->coeffV, p->coeffCnt );
  663. p->delayV = cmMemResizeZ( double, p->delayV, p->coeffCnt-1 ); // there is always one less delay than coeff
  664. p->delayIdx = 0;
  665. cmVOD_LP_Sinc(p->coeffV, p->coeffCnt, srate, fcHz, cmIsFlag(flags,kHighPassFIRFl) ? kHighPass_LPSincFl : 0 );
  666. return cmOkRC;
  667. }
  668. cmRC_t cmFIRFinal( cmFIR* p )
  669. { return cmOkRC; }
  670. cmRC_t cmFIRExec( cmFIR* p, const cmSample_t* sbp, unsigned sn )
  671. {
  672. unsigned delayCnt = p->coeffCnt-1;
  673. int di = p->delayIdx;
  674. const cmSample_t* sep = sbp + sn;
  675. cmSample_t* op = p->outV;
  676. assert( di < delayCnt );
  677. assert( sn <= p->outN );
  678. // for each input sample
  679. while( sbp < sep )
  680. {
  681. // advance the delay line
  682. p->delayIdx = (p->delayIdx + 1) % delayCnt;
  683. const double* cbp = p->coeffV;
  684. const double* cep = cbp + p->coeffCnt;
  685. // mult input sample by coeff[0]
  686. double v = *sbp * *cbp++;
  687. // calc the output sample
  688. while( cbp<cep)
  689. {
  690. // note that the delay is being iterated bcmkwards
  691. if( di == -1 )
  692. di=delayCnt-1;
  693. v += *cbp++ * p->delayV[di];
  694. --di;
  695. }
  696. // store the output sample
  697. *op++ = v;
  698. // insert the input sample
  699. p->delayV[ p->delayIdx ] = *sbp++;
  700. // store the position of the newest ele in the delay line
  701. di = p->delayIdx;
  702. }
  703. return cmOkRC;
  704. }
  705. void cmFIRTest( cmRpt_t* rpt, cmLHeapH_t lhH, cmSymTblH_t stH )
  706. {
  707. unsigned N = 512;
  708. cmKbRecd kb;
  709. cmCtx c;
  710. cmCtxInit(&c,rpt,lhH,stH);
  711. double srate = N;
  712. unsigned procSmpCnt = N;
  713. cmPlotSetup("Test Proc Impl",2,1);
  714. cmSample_t in[ procSmpCnt ];
  715. cmVOS_Fill(in,procSmpCnt,0);
  716. in[0] = 1;
  717. cmVOS_Random(in,procSmpCnt, -1.0, 1.0 );
  718. cmFIR* ffp = cmFIRAllocKaiser( &c, NULL, procSmpCnt,srate, srate*0.025, srate/2, 10, 60 );
  719. //cmFIR* ffp = cmFIRAllocSinc( &c, NULL, 32, 1000, 0 );
  720. cmFftSR* ftp = cmFftAllocSR( &c, NULL, ffp->outV, ffp->outN, kToPolarFftFl );
  721. cmFIRExec( ffp, in, procSmpCnt );
  722. cmFftExecSR( ftp, NULL, 0 );
  723. cmVOR_AmplitudeToDb(ftp->magV,ftp->binCnt,ftp->magV);
  724. printf("coeff cnt:%i\n",ffp->coeffCnt );
  725. cmPlotClear();
  726. cmPlotLineR( "test", NULL, ftp->magV, NULL, ftp->binCnt, NULL, kSolidPlotLineId );
  727. cmPlotDraw();
  728. cmKeyPress(&kb);
  729. cmFftFreeSR(&ftp);
  730. cmFIRFree(&ffp);
  731. }
  732. //------------------------------------------------------------------------------------------------------------
  733. cmFuncFilter* cmFuncFilterAlloc( cmCtx* c, cmFuncFilter* p, unsigned procSmpCnt, cmFuncFiltPtr_t funcPtr, void* userPtr, unsigned wndSmpCnt )
  734. {
  735. cmRC_t rc;
  736. cmFuncFilter* op = cmObjAlloc( cmFuncFilter,c, p );
  737. if( procSmpCnt > 0 && funcPtr != NULL && wndSmpCnt > 0 )
  738. {
  739. if( cmShiftBufAlloc(c,&p->shiftBuf,0,0,0) != NULL )
  740. if((rc = cmFuncFilterInit(op,procSmpCnt,funcPtr,userPtr,wndSmpCnt)) != cmOkRC )
  741. {
  742. cmShiftBuf* sbp = &p->shiftBuf;
  743. cmShiftBufFree(&sbp);
  744. cmObjFree(&op);
  745. }
  746. }
  747. return op;
  748. }
  749. cmRC_t cmFuncFilterFree( cmFuncFilter** pp )
  750. {
  751. cmRC_t rc = cmOkRC;
  752. if( pp!=NULL && *pp != NULL )
  753. {
  754. cmFuncFilter* p = *pp;
  755. if((rc = cmFuncFilterFinal(*pp)) == cmOkRC )
  756. {
  757. cmShiftBuf* sbp = &p->shiftBuf;
  758. cmShiftBufFree(&sbp);
  759. cmMemPtrFree(&p->outV);
  760. cmObjFree(pp);
  761. }
  762. }
  763. return rc;
  764. }
  765. cmRC_t cmFuncFilterInit( cmFuncFilter* p, unsigned procSmpCnt, cmFuncFiltPtr_t funcPtr, void* userPtr, unsigned wndSmpCnt )
  766. {
  767. cmRC_t rc;
  768. if(( rc = cmFuncFilterFinal(p)) != cmOkRC )
  769. return rc;
  770. // The shift buffer always consits of the p->wndSmpCnt-1 samples from the previous
  771. // exec followed by the latest procSmpCnt samples at the end of the buffer
  772. cmShiftBufInit( &p->shiftBuf, procSmpCnt, wndSmpCnt + procSmpCnt - 1, procSmpCnt );
  773. p->outV = cmMemResizeZ( cmSample_t, p->outV, procSmpCnt);
  774. p->outN = procSmpCnt;
  775. p->funcPtr = funcPtr;
  776. p->curWndSmpCnt = 1;
  777. p->wndSmpCnt = wndSmpCnt;
  778. return rc;
  779. }
  780. cmRC_t cmFuncFilterFinal( cmFuncFilter* p )
  781. { return cmOkRC; }
  782. cmRC_t cmFuncFilterExec( cmFuncFilter* p, const cmSample_t* sp, unsigned sn )
  783. {
  784. assert( sn <= p->outN);
  785. // The window used by this function is always causal. At the very beginning of the signal
  786. // the window length begins at 1 and increases until is has the length p->wndSmpCnt.
  787. // Note that this approach ignores any zeros automatically prepended to the beginning of the
  788. // signal by the shift buffer. The first window processed always has a length of 1 and
  789. // begins with the first actual sample given to the shift buffer. Successive windows increase
  790. // by one and start at the first actual sample until the full window length is available
  791. // from the shift buffer. At this point the window length remains constant and it is hopped
  792. // by one sample for each window.
  793. while(cmShiftBufExec(&p->shiftBuf,sp,sn))
  794. {
  795. const cmSample_t* fsp = p->shiftBuf.outV;
  796. cmSample_t* dp = p->outV;
  797. cmSample_t* ep = p->outV + sn; // produce as many output values as there are input samples
  798. // for each output sample
  799. while( dp < ep )
  800. {
  801. // the source range should never extend outside the shift buffer
  802. assert( fsp + p->curWndSmpCnt <= p->shiftBuf.outV + p->shiftBuf.wndSmpCnt );
  803. // calc the next output value
  804. *dp++ = p->funcPtr( fsp, p->curWndSmpCnt, p->userPtr );
  805. // if the window has not yet achieved its full length ...
  806. if( p->curWndSmpCnt < p->wndSmpCnt )
  807. ++p->curWndSmpCnt; // ... then increase its length by 1
  808. else
  809. ++fsp; // ... otherwise shift it ahead by 1
  810. }
  811. }
  812. return cmOkRC;
  813. }
  814. cmSample_t cmFuncFiltTestFunc( const cmSample_t* sp, unsigned sn, void* vp )
  815. {
  816. //printf("% f % f %p % i\n",*sp,*sp+(sn-1),sp,sn);
  817. cmSample_t v = cmVOS_Median(sp,sn);
  818. printf("%f ",v);
  819. return v;
  820. //return *sp;
  821. }
  822. void cmFuncFilterTest( cmRpt_t* rpt, cmLHeapH_t lhH, cmSymTblH_t stH )
  823. {
  824. unsigned procSmpCnt = 1;
  825. unsigned N = 32;
  826. cmSample_t v[N];
  827. cmCtx c;
  828. unsigned i;
  829. cmCtxAlloc(&c,rpt,lhH,stH);
  830. cmVOS_Seq(v,N,0,1);
  831. cmVOS_Print(rpt,1,32,v);
  832. cmFuncFilter* ffp = NULL;
  833. ffp = cmFuncFilterAlloc( &c, NULL, procSmpCnt, cmFuncFiltTestFunc, NULL, 5 );
  834. for(i=0; i<N; ++i)
  835. cmFuncFilterExec(ffp,v+(i*procSmpCnt),procSmpCnt);
  836. cmFuncFilterFree( &ffp );
  837. cmCtxFinal(&c);
  838. //unsigned v1n = 9;
  839. //cmSample_t v1[9] = { 1, 75, 91, 35, 6, 80, 40, 91, 79};
  840. //cmSample_t v1[10] = {53, 64, 48, 78, 30, 59, 7, 50, 71, 53 };
  841. //printf("Median: %f \n",cmVOS_Median(v1,v1+v1n));
  842. }
  843. //------------------------------------------------------------------------------------------------------------
  844. cmDhmm* cmDhmmAlloc( cmCtx* c, cmDhmm* ap, unsigned stateN, unsigned symN, cmReal_t* initV, cmReal_t* transM, cmReal_t* stsM )
  845. {
  846. cmDhmm* p = cmObjAlloc( cmDhmm, c, ap );
  847. if( stateN > 0 && symN > 0 )
  848. if( cmDhmmInit(p, stateN, symN, initV, transM, stsM ) != cmOkRC )
  849. cmObjFree(&p);
  850. return p;
  851. }
  852. cmRC_t cmDhmmFree( cmDhmm** pp )
  853. {
  854. cmRC_t rc = cmOkRC;
  855. cmDhmm* p = *pp;
  856. if( pp==NULL || *pp==NULL )
  857. return cmOkRC;
  858. if((rc = cmDhmmFinal(p)) != cmOkRC )
  859. return cmOkRC;
  860. cmObjFree(pp);
  861. return rc;
  862. }
  863. cmRC_t cmDhmmInit( cmDhmm* p, unsigned stateN, unsigned symN, cmReal_t* initV, cmReal_t* transM, cmReal_t* stsM )
  864. {
  865. cmRC_t rc;
  866. if((rc = cmDhmmFinal(p)) != cmOkRC )
  867. return rc;
  868. p->stateN = stateN;
  869. p->symN = symN;
  870. p->initV = initV;
  871. p->transM = transM;
  872. p->stsM = stsM;
  873. return cmOkRC;
  874. }
  875. cmRC_t cmDhmmFinal( cmDhmm* p )
  876. { return cmOkRC; }
  877. cmRC_t cmDhmmExec( cmDhmm* p )
  878. {
  879. return cmOkRC;
  880. }
  881. // Generate a random matrix with rows that sum to 1.0.
  882. void _cmDhmmGenRandMatrix( cmReal_t* dp, unsigned rn, unsigned cn )
  883. {
  884. cmReal_t v[ cn ];
  885. unsigned i,j;
  886. for(i=0; i<rn; ++i)
  887. {
  888. cmVOR_Random( v, cn, 0.0, 1.0 );
  889. cmVOR_NormalizeProbability( v, cn);
  890. for(j=0; j<cn; ++j)
  891. dp[ (j * rn) + i ] = v[j];
  892. }
  893. }
  894. enum { kEqualProbHmmFl=0x01, kRandProbHmmFl=0x02, kManualProbHmmFl=0x04 };
  895. void _cmDhmmGenProb( cmReal_t* dp, unsigned rn, unsigned cn, unsigned flags, const cmReal_t* sp )
  896. {
  897. switch( flags )
  898. {
  899. case kRandProbHmmFl:
  900. _cmDhmmGenRandMatrix( dp, rn, cn );
  901. break;
  902. case kEqualProbHmmFl:
  903. {
  904. // equal prob
  905. cmReal_t pr = 1.0/cn;
  906. unsigned i,j;
  907. for(i=0; i<rn; ++i)
  908. for(j=0; j<cn; ++j)
  909. dp[ (j*rn) + i ] = pr;
  910. }
  911. break;
  912. case kManualProbHmmFl:
  913. cmVOR_Copy( dp, (rn*cn), sp );
  914. break;
  915. default:
  916. assert(0);
  917. }
  918. }
  919. // generate a random integer in the range 0 to probN-1 where probV[ probN ] contains
  920. // the probability of generating each of the possible values.
  921. unsigned _cmDhmmGenRandInt( const cmReal_t* probV, unsigned probN, unsigned stride )
  922. {
  923. cmReal_t tmp[ probN ];
  924. cmReal_t cumSumV[ probN+1 ];
  925. const cmReal_t* sp = probV;
  926. cumSumV[0] = 0;
  927. if( stride > 1 )
  928. {
  929. cmVOR_CopyStride( tmp, probN, probV, stride );
  930. sp = tmp;
  931. }
  932. cmVOR_CumSum( cumSumV+1, probN, sp );
  933. return cmVOR_BinIndex( cumSumV, probN+1, (cmReal_t)rand()/RAND_MAX );
  934. }
  935. cmRC_t cmDhmmGenObsSequence( cmDhmm* p, unsigned* dbp, unsigned dn )
  936. {
  937. const unsigned* dep = dbp + dn;
  938. // generate the first state based on the init state prob. vector
  939. unsigned state = _cmDhmmGenRandInt( p->initV, p->stateN, 1 );
  940. // generate an observation from the state based on the symbol prob. vector
  941. *dbp++ = _cmDhmmGenRandInt( p->stsM + state, p->symN, p->stateN );
  942. while( dbp < dep )
  943. {
  944. // get the next state based on the previous state
  945. state = _cmDhmmGenRandInt( p->transM + state, p->stateN, p->stateN );
  946. // given a state generate an observation
  947. *dbp++ = _cmDhmmGenRandInt( p->stsM + state, p->symN, p->stateN );
  948. }
  949. return cmOkRC;
  950. }
  951. /// Perform a forward evaluation of the model given a set of observations.
  952. /// initPrV[ stateN ] is the probability of the model being in each state at the start of the evaluation.
  953. /// alphaM[ stateN, obsN ] is a return value and represents the probability of seeing each symbol at each time step.
  954. enum { kNoLogScaleHmmFl = 0x00, kLogScaleHmmFl = 0x01 };
  955. cmRC_t cmDHmmForwardEval( cmDhmm* p, const cmReal_t* initPrV, const unsigned* obsV, unsigned obsN, cmReal_t* alphaM, unsigned flags, cmReal_t* logProbPtr )
  956. {
  957. bool scaleFl = cmIsFlag(flags,kLogScaleHmmFl);
  958. cmReal_t logProb = 0;
  959. cmReal_t* abp = alphaM; // define first dest. column
  960. const cmReal_t* aep = abp + p->stateN;
  961. const cmReal_t* sts = p->stsM + (obsV[0] * p->stateN); // stsM[] column for assoc'd with first obs. symbol
  962. unsigned i;
  963. // calc the prob of begining in each state given the obs. symbol
  964. for(i=0; abp < aep; ++i )
  965. *abp++ = *initPrV++ * *sts++;
  966. // scale to prevent underflow
  967. if( scaleFl )
  968. {
  969. cmReal_t sum = cmVOR_Sum(abp-p->stateN,p->stateN);
  970. if( sum > 0 )
  971. {
  972. cmVOR_DivVS( abp-p->stateN,p->stateN,sum);
  973. logProb += log(sum);
  974. }
  975. }
  976. // for each time step
  977. for(i=1; i<obsN; ++i)
  978. {
  979. // next state 0 (first col, first row) is calc'd first
  980. const cmReal_t* tm = p->transM;
  981. // pick the stsM[] column assoc'd with ith observation symbol
  982. const cmReal_t* sts = p->stsM + (obsV[i] * p->stateN);
  983. // store a pointer to the alpha column assoc'd with obsV[i-1]
  984. const cmReal_t* app0 = abp - p->stateN;
  985. aep = abp + p->stateN;
  986. // for each dest state
  987. while( abp < aep )
  988. {
  989. // prob of being in each state on the previous time step
  990. const cmReal_t* app = app0;
  991. const cmReal_t* ape = app + p->stateN;
  992. *abp = 0;
  993. // for each src state - calc prob. of trans from src to dest
  994. while( app<ape )
  995. *abp += *app++ * *tm++;
  996. // calc prob of obs symbol in dest state
  997. *abp++ *= *sts++;
  998. }
  999. // scale to prevent underflow
  1000. if( scaleFl )
  1001. {
  1002. cmReal_t sum = cmVOR_Sum(abp-p->stateN,p->stateN);
  1003. if( sum > 0 )
  1004. {
  1005. cmVOR_DivVS( abp-p->stateN,p->stateN,sum);
  1006. logProb += log(sum);
  1007. }
  1008. }
  1009. }
  1010. if( logProbPtr != NULL )
  1011. *logProbPtr = logProb;
  1012. return cmOkRC;
  1013. }
  1014. cmRC_t cmDHmmBcmkwardEval( cmDhmm* p, const unsigned* obsV, unsigned obsN, cmReal_t* betaM, unsigned flags )
  1015. {
  1016. bool scaleFl = cmIsFlag(flags,kLogScaleHmmFl);
  1017. int i,j,t;
  1018. cmVOR_Fill(betaM+((obsN-1)*p->stateN),p->stateN,1.0);
  1019. // for each time step
  1020. for(t=obsN-2; t>=0; --t)
  1021. {
  1022. // for each state at t
  1023. for(i=0; i<p->stateN; ++i)
  1024. {
  1025. double Bt = 0;
  1026. // for each state at t+1
  1027. for(j=0; j<p->stateN; ++j)
  1028. {
  1029. double aij = p->transM[ (j * p->stateN) + i ];
  1030. double bj = p->stsM[ (obsV[t+1] * p->stateN) + j ];
  1031. double Bt1 = betaM[ ((t+1) * p->stateN) + j ];
  1032. Bt += aij * bj * Bt1;
  1033. }
  1034. betaM[ (t * p->stateN) + i ] = Bt;
  1035. }
  1036. if( scaleFl )
  1037. {
  1038. double* bp = betaM + (t * p->stateN);
  1039. double sum = cmVOR_Sum(bp, p->stateN );
  1040. if( sum > 0 )
  1041. cmVOR_DivVS( bp, p->stateN, sum );
  1042. }
  1043. }
  1044. return cmOkRC;
  1045. }
  1046. void _cmDhmmNormRow( cmReal_t* p, unsigned pn, unsigned stride, const cmReal_t* sp )
  1047. {
  1048. if( sp == NULL )
  1049. sp = p;
  1050. cmReal_t sum = 0;
  1051. unsigned n = pn * stride;
  1052. const cmReal_t* bp = sp;
  1053. const cmReal_t* ep = bp + n;
  1054. for(; bp<ep; bp+=stride)
  1055. sum += *bp;
  1056. for(ep = p+n; p<ep; p+=stride,sp+=stride)
  1057. *p = *sp / sum;
  1058. }
  1059. void _cmDhmmNormMtxRows( cmReal_t* dp, unsigned rn, unsigned cn, const cmReal_t* sp )
  1060. {
  1061. const cmReal_t* erp = sp + rn;
  1062. while( sp < erp )
  1063. _cmDhmmNormRow( dp++, cn, rn, sp++ );
  1064. }
  1065. cmRC_t cmDhmmTrainEM( cmDhmm* p, const unsigned* obsV, unsigned obsN, unsigned iterCnt, unsigned flags )
  1066. {
  1067. unsigned i,j,k,t;
  1068. cmReal_t alphaM[ p->stateN * obsN ];
  1069. cmReal_t betaM[ p->stateN * obsN ];
  1070. cmReal_t g[ p->stateN * obsN ];
  1071. cmReal_t q[ p->stateN * p->symN ];
  1072. cmReal_t E[ p->stateN * p->stateN ];
  1073. cmReal_t logProb = 0;
  1074. //cmDhmmReport(p->obj.ctx,p);
  1075. for(k=0; k<iterCnt; ++k)
  1076. {
  1077. cmVOR_Fill( q, (p->stateN * p->symN), 0 );
  1078. cmVOR_Fill( E, (p->stateN * p->stateN), 0 );
  1079. // calculate alpha and beta
  1080. cmDHmmForwardEval( p, p->initV, obsV, obsN, alphaM, flags, &logProb);
  1081. cmDHmmBcmkwardEval( p, obsV, obsN, betaM, flags );
  1082. // gamma[ stateN, obsN ] = alphaM .* betaM - gamma is the probability of being in each state at each time step
  1083. cmVOR_MultVVV( g,(p->stateN*obsN), alphaM, betaM );
  1084. // normalize gamma
  1085. for(i=0; i<obsN; ++i)
  1086. cmVOR_NormalizeProbability( g + (i*p->stateN), p->stateN );
  1087. //printf("ITER:%i logProb:%f\n",k,logProb);
  1088. // count the number of times state i emits obsV[0] in the starting location
  1089. cmVOR_Copy( q + (obsV[0] * p->stateN), p->stateN, g );
  1090. for(t=0; t<obsN-1; ++t)
  1091. {
  1092. // point to alpha[:,t] and beta[:,t+1]
  1093. const cmReal_t* alpha_t0 = alphaM + (t*p->stateN);
  1094. const cmReal_t* beta_t1 = betaM + ((t+1)*p->stateN);
  1095. cmReal_t Et[ p->stateN * p->stateN ];
  1096. cmReal_t Esum = 0;
  1097. // for each source state
  1098. for(i=0; i<p->stateN; ++i)
  1099. {
  1100. // for each dest state
  1101. for(j=0; j<p->stateN; ++j)
  1102. {
  1103. // prob of transitioning from state i to j and emitting obs[t] at time t
  1104. cmReal_t Eps = alpha_t0[i] * p->transM[ (j*p->stateN) + i ] * p->stsM[ (obsV[t+1]*p->stateN) + j ] * beta_t1[j];
  1105. // count the number of transitions from i to j
  1106. Et[ (j*p->stateN) + i ] = Eps;
  1107. Esum += Eps;
  1108. }
  1109. // count the number of times state i emits obsV[t]
  1110. q[ (obsV[t+1] * p->stateN) + i ] += g[ ((t+1)*p->stateN) + i ];
  1111. }
  1112. // normalize Et and sum it into E
  1113. cmVOR_DivVS( Et, (p->stateN*p->stateN), Esum );
  1114. cmVOR_AddVV( E, (p->stateN*p->stateN), Et );
  1115. }
  1116. // update the model
  1117. _cmDhmmNormMtxRows( p->initV, 1, p->stateN, g );
  1118. _cmDhmmNormMtxRows( p->transM, p->stateN, p->stateN, E );
  1119. _cmDhmmNormMtxRows( p->stsM, p->stateN, p->symN, q );
  1120. }
  1121. return cmOkRC;
  1122. }
  1123. cmRC_t cmDhmmReport( cmDhmm* p )
  1124. {
  1125. cmVOR_PrintL("initV:\n", p->obj.err.rpt, 1, p->stateN, p->initV );
  1126. cmVOR_PrintL("transM:\n", p->obj.err.rpt, p->stateN, p->stateN, p->transM );
  1127. cmVOR_PrintL("symM:\n", p->obj.err.rpt, p->stateN, p->symN, p->stsM );
  1128. return cmOkRC;
  1129. }
  1130. void cmDhmmTest( cmRpt_t* rpt, cmLHeapH_t lhH, cmSymTblH_t stH )
  1131. {
  1132. unsigned stateN = 2;
  1133. unsigned symN = 3;
  1134. unsigned obsN = 4;
  1135. unsigned iterN = 10;
  1136. cmReal_t initV0[ stateN ];
  1137. cmReal_t transM0[ stateN * stateN ];
  1138. cmReal_t stsM0[ symN * stateN ];
  1139. cmReal_t initV1[ stateN ];
  1140. cmReal_t transM1[ stateN * stateN ];
  1141. cmReal_t stsM1[ symN * stateN ];
  1142. unsigned obsV[ obsN ];
  1143. unsigned hist[ symN ];
  1144. unsigned genFl = kManualProbHmmFl;
  1145. cmReal_t initV[] =
  1146. {
  1147. 0.44094,
  1148. 0.55906
  1149. };
  1150. cmReal_t transM[] =
  1151. {
  1152. 0.48336,
  1153. 0.81353,
  1154. 0.51664,
  1155. 0.18647,
  1156. };
  1157. cmReal_t stsM[] =
  1158. {
  1159. 0.20784,
  1160. 0.18698,
  1161. 0.43437,
  1162. 0.24102,
  1163. 0.35779,
  1164. 0.57199
  1165. };
  1166. unsigned obsM[] = { 2, 2, 2, 0 };
  1167. cmReal_t initV2[] = { 0.79060, 0.20940 };
  1168. cmReal_t transM2[] = { 0.508841, 0.011438, 0.491159, 0.988562 };
  1169. cmReal_t stsM2[] = { 0.25789, 0.35825, 0.61981, 0.42207, 0.12230, 0.21969 };
  1170. //srand( time(NULL) );
  1171. // generate a random HMM
  1172. _cmDhmmGenProb( initV0, 1, stateN, genFl, initV );
  1173. _cmDhmmGenProb( transM0, stateN, stateN, genFl, transM );
  1174. _cmDhmmGenProb( stsM0, stateN, symN, genFl, stsM );
  1175. cmCtx* c = cmCtxAlloc( NULL, rpt, lhH, stH);
  1176. cmDhmm* h0p = cmDhmmAlloc( c, NULL, stateN, symN, initV0, transM0, stsM0 );
  1177. // generate an observation sequence based on the random HMM
  1178. //cmDhmmGenObsSequence(c, h0p, obsV, obsN );
  1179. memcpy(obsV,obsM,obsN*sizeof(unsigned));
  1180. if( 0 )
  1181. {
  1182. // print the HMM
  1183. cmDhmmReport( h0p);
  1184. // print the observation symbols
  1185. cmVOU_PrintL("obs:\n", rpt, 1, obsN, obsV );
  1186. // print the histogram of the obs. symbols
  1187. cmVOU_Hist( hist, symN, obsV, obsN );
  1188. cmVOU_PrintL("hist:\n", rpt, 1, symN, hist );
  1189. // calc alpha (the forward probabilities)
  1190. cmReal_t alphaM[ h0p->stateN*obsN ];
  1191. cmReal_t logProb=0;
  1192. cmDHmmForwardEval( h0p, h0p->initV, obsV, obsN, alphaM, kLogScaleHmmFl, &logProb);
  1193. printf("log prob:%f\n alpha:\n", logProb );
  1194. cmVOR_Print( rpt, h0p->stateN, obsN, alphaM );
  1195. // calc beta (the bcmkward probabilities)
  1196. cmReal_t betaM[ h0p->stateN*obsN ];
  1197. logProb=0;
  1198. cmDHmmBcmkwardEval( h0p, obsV, obsN, betaM, kLogScaleHmmFl);
  1199. printf("log prob:%f\n beta:\n", logProb );
  1200. cmVOR_Print( h0p->obj.err.rpt, h0p->stateN, obsN, betaM );
  1201. }
  1202. // initialize a second HMM with random probabilities
  1203. _cmDhmmGenProb( initV1, 1, stateN, kManualProbHmmFl, initV2 );
  1204. _cmDhmmGenProb( transM1, stateN, stateN, kManualProbHmmFl, transM2 );
  1205. _cmDhmmGenProb( stsM1, stateN, symN, kManualProbHmmFl, stsM2 );
  1206. cmDhmm* h1p = cmDhmmAlloc( c, NULL, stateN, symN, initV1, transM1, stsM1 );
  1207. cmDhmmTrainEM( h1p, obsV, obsN, iterN, kLogScaleHmmFl );
  1208. cmDhmmFree(&h1p);
  1209. cmDhmmFree(&h0p);
  1210. cmCtxFree(&c);
  1211. }
  1212. //------------------------------------------------------------------------------------------------------------
  1213. cmConvolve* cmConvolveAlloc( cmCtx* c, cmConvolve* ap, const cmSample_t* h, unsigned hn, unsigned procSmpCnt )
  1214. {
  1215. cmConvolve* p = cmObjAlloc( cmConvolve, c, ap);
  1216. p->fft = cmFftAllocSR( c,NULL,NULL,0,kNoConvertFftFl);
  1217. p->ifft= cmIFftAllocRS(c,NULL,p->fft->binCnt);
  1218. if( hn > 0 && procSmpCnt > 0 )
  1219. if( cmConvolveInit(p,h,hn,procSmpCnt) != cmOkRC )
  1220. cmObjFree(&p);
  1221. return p;
  1222. }
  1223. cmRC_t cmConvolveFree( cmConvolve** pp )
  1224. {
  1225. cmRC_t rc;
  1226. cmConvolve* p = *pp;
  1227. if( pp == NULL || *pp == NULL )
  1228. return cmOkRC;
  1229. if((rc = cmConvolveFinal(p)) != cmOkRC )
  1230. return cmOkRC;
  1231. cmFftFreeSR(&p->fft);
  1232. cmIFftFreeRS(&p->ifft);
  1233. cmMemPtrFree(&p->H);
  1234. cmMemPtrFree(&p->outV);
  1235. cmObjFree(pp);
  1236. return cmOkRC;
  1237. }
  1238. cmRC_t cmConvolveInit( cmConvolve* p, const cmSample_t* h, unsigned hn, unsigned procSmpCnt )
  1239. {
  1240. cmRC_t rc;
  1241. unsigned i;
  1242. unsigned cn = cmNextPowerOfTwo( hn + procSmpCnt - 1 );
  1243. if((rc = cmConvolveFinal(p)) != cmOkRC )
  1244. return rc;
  1245. cmFftInitSR( p->fft, NULL, cn, kNoConvertFftFl );
  1246. cmIFftInitRS( p->ifft, p->fft->binCnt);
  1247. p->H = cmMemResizeZ( cmComplexR_t,p->H, p->fft->binCnt );
  1248. p->outV = cmMemResizeZ( cmSample_t,p->outV, cn );
  1249. p->olaV = p->outV + procSmpCnt;
  1250. p->outN = procSmpCnt;
  1251. p->hn = hn;
  1252. // take the FFT of the impulse response
  1253. cmFftExecSR( p->fft, h, hn );
  1254. // copy the FFT of the impulse response to p->H[]
  1255. for(i=0; i<p->fft->binCnt; ++i)
  1256. p->H[i] = p->fft->complexV[i] / p->fft->wndSmpCnt;
  1257. return cmOkRC;
  1258. }
  1259. cmRC_t cmConvolveFinal( cmConvolve* p )
  1260. { return cmOkRC; }
  1261. cmRC_t cmConvolveExec( cmConvolve* p, const cmSample_t* x, unsigned xn )
  1262. {
  1263. unsigned i;
  1264. // take FT of input signal
  1265. cmFftExecSR( p->fft, x, xn );
  1266. // multiply the signal spectra of the input signal and impulse response
  1267. for(i=0; i<p->fft->binCnt; ++i)
  1268. p->ifft->complexV[i] = p->H[i] * p->fft->complexV[i];
  1269. // take the IFFT of the convolved spectrum
  1270. cmIFftExecRS(p->ifft,NULL);
  1271. // sum with previous impulse response tail
  1272. cmVOS_AddVVV( p->outV, p->outN-1, p->olaV, p->ifft->outV );
  1273. // first sample of the impulse response tail is complete
  1274. p->outV[p->outN-1] = p->ifft->outV[p->outN-1];
  1275. // store the new impulse response tail
  1276. cmVOS_Copy(p->olaV,p->hn-1,p->ifft->outV + p->outN );
  1277. return cmOkRC;
  1278. }
  1279. cmRC_t cmConvolveSignal( cmCtx* c, const cmSample_t* h, unsigned hn, const cmSample_t* x, unsigned xn, cmSample_t* y, unsigned yn )
  1280. {
  1281. cmConvolve* p = cmConvolveAlloc(c,NULL,h,hn,xn);
  1282. cmConvolveExec(p,x,xn);
  1283. unsigned n = cmMin(p->outN,yn);
  1284. cmVOS_Copy(y,n,p->outV);
  1285. if( yn > p->outN )
  1286. {
  1287. unsigned m = cmMin(yn-p->outN,p->hn-1);
  1288. cmVOS_Copy(y+n,m,p->olaV);
  1289. }
  1290. cmConvolveFree(&p);
  1291. return cmOkRC;
  1292. }
  1293. cmRC_t cmConvolveTest(cmRpt_t* rpt, cmLHeapH_t lhH, cmSymTblH_t stH )
  1294. {
  1295. cmCtx *c = cmCtxAlloc(NULL,rpt,lhH,stH);
  1296. cmSample_t h[] = { 1, .5, .25, 0, 0 };
  1297. unsigned hn = sizeof(h) / sizeof(h[0]);
  1298. cmSample_t x[] = { 1, 0, 0, 0, 1, 0, 0, 0 };
  1299. unsigned xn = sizeof(x) / sizeof(x[0]);
  1300. unsigned yn = xn+hn-1;
  1301. cmSample_t y[yn];
  1302. //cmVOS_Hann(h,5);
  1303. cmConvolve* p = cmConvolveAlloc(c,NULL,h,hn,xn);
  1304. cmConvolveExec(p,x,xn);
  1305. cmVOS_Print( rpt, 1, p->outN, p->outV );
  1306. cmVOS_Print( rpt, 1, p->hn-1, p->olaV );
  1307. cmConvolveFree(&p);
  1308. cmConvolveSignal(c,h,hn,x,xn,y,yn);
  1309. cmVOS_Print( rpt, 1, hn+xn-1, y );
  1310. cmCtxFree(&c);
  1311. return cmOkRC;
  1312. }
  1313. //------------------------------------------------------------------------------------------------------------
  1314. cmBfcc* cmBfccAlloc( cmCtx* ctx, cmBfcc* ap, unsigned bandCnt, unsigned binCnt, double binHz )
  1315. {
  1316. cmBfcc* p = cmObjAlloc( cmBfcc, ctx, ap );
  1317. if( bandCnt > 0 )
  1318. if( cmBfccInit( p, bandCnt, binCnt, binHz ) != cmOkRC )
  1319. cmBfccFree(&p);
  1320. return p;
  1321. }
  1322. cmRC_t cmBfccFree( cmBfcc** pp )
  1323. {
  1324. cmRC_t rc;
  1325. if( pp== NULL || *pp==NULL)
  1326. return cmOkRC;
  1327. cmBfcc* p = *pp;
  1328. if((rc = cmBfccFinal(p)) != cmOkRC )
  1329. return rc;
  1330. cmMemPtrFree(&p->dctMtx);
  1331. cmMemPtrFree(&p->filtMask);
  1332. cmMemPtrFree(&p->outV);
  1333. cmObjFree(pp);
  1334. return rc;
  1335. }
  1336. cmRC_t cmBfccInit( cmBfcc* p, unsigned bandCnt, unsigned binCnt, double binHz )
  1337. {
  1338. cmRC_t rc;
  1339. if((rc = cmBfccFinal(p)) != cmOkRC )
  1340. return rc;
  1341. p->dctMtx = cmMemResizeZ( cmReal_t, p->dctMtx, bandCnt*bandCnt);
  1342. p->filtMask = cmMemResizeZ( cmReal_t, p->filtMask, bandCnt*binCnt);
  1343. p->outV = cmMemResizeZ( cmReal_t, p->outV, bandCnt );
  1344. p->binCnt = binCnt;
  1345. p->bandCnt = bandCnt;
  1346. cmVOR_BarkMask( p->filtMask, bandCnt, binCnt, binHz );
  1347. cmVOR_DctMatrix(p->dctMtx, bandCnt, bandCnt );
  1348. return rc;
  1349. }
  1350. cmRC_t cmBfccFinal( cmBfcc* p )
  1351. { return cmOkRC; }
  1352. cmRC_t cmBfccExec( cmBfcc* p, const cmReal_t* magV, unsigned binCnt )
  1353. {
  1354. assert( binCnt <= p->binCnt );
  1355. cmReal_t t[ p->bandCnt ];
  1356. cmReal_t v[ binCnt ];
  1357. // convert magnitude to power
  1358. cmVOR_PowVVS(v,binCnt,magV,2.0);
  1359. // apply the filter mask to the power spectrum
  1360. cmVOR_MultVMV( t, p->bandCnt, p->filtMask, binCnt, v );
  1361. //cmVOR_PrintL("\t:\n", p->obj.ctx->outFuncPtr, 1, p->bandCnt, t);
  1362. cmVOR_ReplaceLte( t, p->bandCnt, t, 0, 0.1e-5 );
  1363. cmVOR_LogV( t, p->bandCnt, t );
  1364. // decorellate the bands with a DCT
  1365. cmVOR_MultVMV( p->outV, p->bandCnt, p->dctMtx, p->bandCnt, t );
  1366. return cmOkRC;
  1367. }
  1368. void cmBfccTest( cmRpt_t* rpt, cmLHeapH_t lhH, cmSymTblH_t stH )
  1369. {
  1370. double srate = 11025;
  1371. unsigned binCnt = 129;
  1372. double binHz = srate/(binCnt-1);
  1373. unsigned bandCnt = kDefaultBarkBandCnt;
  1374. cmCtx* c = cmCtxAlloc( NULL, rpt, lhH, stH );
  1375. cmBfcc* b = cmBfccAlloc( c, NULL, bandCnt, binCnt, binHz );
  1376. cmReal_t powV[] = {
  1377. 0.8402, 0.3944, 0.7831, 0.7984, 0.9116, 0.1976, 0.3352, 0.7682, 0.2778, 0.5540,
  1378. 0.4774, 0.6289, 0.3648, 0.5134, 0.9522, 0.9162, 0.6357, 0.7173, 0.1416, 0.6070,
  1379. 0.0163, 0.2429, 0.1372, 0.8042, 0.1567, 0.4009, 0.1298, 0.1088, 0.9989, 0.2183,
  1380. 0.5129, 0.8391, 0.6126, 0.2960, 0.6376, 0.5243, 0.4936, 0.9728, 0.2925, 0.7714,
  1381. 0.5267, 0.7699, 0.4002, 0.8915, 0.2833, 0.3525, 0.8077, 0.9190, 0.0698, 0.9493,
  1382. 0.5260, 0.0861, 0.1922, 0.6632, 0.8902, 0.3489, 0.0642, 0.0200, 0.4577, 0.0631,
  1383. 0.2383, 0.9706, 0.9022, 0.8509, 0.2667, 0.5398, 0.3752, 0.7602, 0.5125, 0.6677,
  1384. 0.5316, 0.0393, 0.4376, 0.9318, 0.9308, 0.7210, 0.2843, 0.7385, 0.6400, 0.3540,
  1385. 0.6879, 0.1660, 0.4401, 0.8801, 0.8292, 0.3303, 0.2290, 0.8934, 0.3504, 0.6867,
  1386. 0.9565, 0.5886, 0.6573, 0.8587, 0.4396, 0.9240, 0.3984, 0.8148, 0.6842, 0.9110,
  1387. 0.4825, 0.2158, 0.9503, 0.9201, 0.1477, 0.8811, 0.6411, 0.4320, 0.6196, 0.2811,
  1388. 0.7860, 0.3075, 0.4470, 0.2261, 0.1875, 0.2762, 0.5564, 0.4165, 0.1696, 0.9068,
  1389. 0.1032, 0.1261, 0.4954, 0.7605, 0.9848, 0.9350, 0.6844, 0.3832, 0.7498 };
  1390. //cmVOR_Random(powV, binCnt, 0.0, 1.0 );
  1391. cmBfccExec(b,powV,binCnt);
  1392. //cmVOR_PrintL("\nin:\n", rpt, 1, binCnt, powV);
  1393. //cmVOR_PrintL("\nfilt:\n", rpt, bandCnt, binCnt, b->filtMask);
  1394. //cmVOR_PrintL("\ndct:\n", rpt, bandCnt, bandCnt,b->dctMtx);
  1395. cmVOR_PrintL("\nbfcc:\n", rpt, 1, bandCnt, b->outV);
  1396. cmBfccFree(&b);
  1397. cmCtxFree(&c);
  1398. }
  1399. //------------------------------------------------------------------------------------------------------------
  1400. cmCeps* cmCepsAlloc( cmCtx* ctx, cmCeps* ap, unsigned binCnt, unsigned outN )
  1401. {
  1402. cmCeps* p = cmObjAlloc( cmCeps, ctx, ap );
  1403. //cmIFftAllocRR( ctx, &p->ft, 0 );
  1404. if( binCnt > 0 )
  1405. if( cmCepsInit( p, binCnt, outN ) != cmOkRC )
  1406. cmCepsFree(&p);
  1407. return p;
  1408. }
  1409. cmRC_t cmCepsFree( cmCeps** pp )
  1410. {
  1411. cmRC_t rc;
  1412. if( pp== NULL || *pp==NULL)
  1413. return cmOkRC;
  1414. cmCeps* p = *pp;
  1415. if((rc = cmCepsFinal(p)) != cmOkRC )
  1416. return rc;
  1417. //cmObjFreeStatic( cmIFftFreeRR, cmIFftRR, p->ft );
  1418. cmMemPtrFree(&p->dctM);
  1419. cmMemPtrFree(&p->outV);
  1420. cmObjFree(pp);
  1421. return rc;
  1422. }
  1423. cmRC_t cmCepsInit( cmCeps* p, unsigned binCnt, unsigned outN )
  1424. {
  1425. cmRC_t rc;
  1426. if((rc = cmCepsFinal(p)) != cmOkRC )
  1427. return rc;
  1428. //cmIFftInitRR( &p->ft, binCnt );
  1429. p->dct_cn = (binCnt-1)*2;
  1430. p->dctM = cmMemResize( cmReal_t, p->dctM, outN*p->dct_cn );
  1431. p->outN = outN; //p->ft.outN;
  1432. p->outV = cmMemResizeZ( cmReal_t, p->outV, outN ); //p->ft.outV;
  1433. p->binCnt = binCnt;
  1434. assert( outN <= p->dct_cn );
  1435. cmVOR_DctMatrix( p->dctM, outN, p->dct_cn );
  1436. return rc;
  1437. }
  1438. cmRC_t cmCepsFinal( cmCeps* p )
  1439. { return cmOkRC; }
  1440. cmRC_t cmCepsExec( cmCeps* p, const cmReal_t* magV, const cmReal_t* phsV, unsigned binCnt )
  1441. {
  1442. assert( binCnt == p->binCnt );
  1443. cmReal_t v[ p->dct_cn ];
  1444. // guard against zeros in the magn spectrum
  1445. cmVOR_ReplaceLte(v,binCnt,magV,0.0,0.00001);
  1446. // take the log of the spectrum
  1447. cmVOR_LogV(v,binCnt,v);
  1448. // reconstruct the negative frequencies
  1449. int i,j;
  1450. for(i=1,j=p->dct_cn-1; i<binCnt; ++i,--j)
  1451. v[j] = v[i];
  1452. // take the DCT
  1453. cmVOR_MultVMV( p->outV, p->outN, p->dctM, p->dct_cn, v );
  1454. //cmIFftExecPolarRR( &p->ft, v, phsV );
  1455. return cmOkRC;
  1456. }
  1457. //------------------------------------------------------------------------------------------------------------
  1458. cmOla* cmOlaAlloc( cmCtx* ctx, cmOla* ap, unsigned wndSmpCnt, unsigned hopSmpCnt, unsigned procSmpCnt, unsigned wndTypeId )
  1459. {
  1460. cmOla* p = cmObjAlloc( cmOla, ctx, ap );
  1461. cmWndFuncAlloc( ctx, &p->wf, kHannWndId, wndSmpCnt, 0);
  1462. if( wndSmpCnt > 0 )
  1463. if( cmOlaInit(p,wndSmpCnt,hopSmpCnt,procSmpCnt,wndTypeId) != cmOkRC )
  1464. cmOlaFree(&p);
  1465. return p;
  1466. }
  1467. cmRC_t cmOlaFree( cmOla** pp )
  1468. {
  1469. cmRC_t rc;
  1470. if( pp==NULL || *pp==NULL )
  1471. return cmOkRC;
  1472. cmOla* p = *pp;
  1473. if(( rc = cmOlaFinal(p)) != cmOkRC )
  1474. return rc;
  1475. cmMemPtrFree(&p->bufV);
  1476. cmMemPtrFree(&p->outV);
  1477. cmObjFreeStatic( cmWndFuncFree, cmWndFunc, p->wf );
  1478. cmObjFree(pp);
  1479. return rc;
  1480. }
  1481. cmRC_t cmOlaInit( cmOla* p, unsigned wndSmpCnt, unsigned hopSmpCnt, unsigned procSmpCnt, unsigned wndTypeId )
  1482. {
  1483. cmRC_t rc;
  1484. if((rc = cmOlaFinal(p)) != cmOkRC )
  1485. return rc;
  1486. if((rc = cmWndFuncInit( &p->wf, wndTypeId, wndSmpCnt, 0)) != cmOkRC )
  1487. return rc;
  1488. p->bufV = cmMemResizeZ( cmSample_t, p->bufV, wndSmpCnt );
  1489. p->outV = cmMemResizeZ( cmSample_t, p->outV, hopSmpCnt );
  1490. p->outPtr = p->outV + hopSmpCnt;
  1491. // hopSmpCnt must be an even multiple of procSmpCnt
  1492. assert( hopSmpCnt % procSmpCnt == 0 );
  1493. assert( wndSmpCnt >= hopSmpCnt );
  1494. p->wndSmpCnt = wndSmpCnt;
  1495. p->hopSmpCnt = hopSmpCnt;
  1496. p->procSmpCnt= procSmpCnt;
  1497. p->idx = 0;
  1498. return rc;
  1499. }
  1500. cmRC_t cmOlaFinal( cmOla* p )
  1501. { return cmOkRC; }
  1502. // The incoming buffer and the ola buf (bufV)
  1503. // can be divided into three logical parts:
  1504. //
  1505. // [head][middle][tail]
  1506. //
  1507. // head = hopSmpCnt values
  1508. // tail = hopSmpCnt values
  1509. // middle = wndSmpCnt - (2*hopSmpCnt) values
  1510. //
  1511. // Each exec can be broken into three phases:
  1512. //
  1513. // outV = bufV[tail] + in[head]
  1514. // bufV[middle] += in[middle]
  1515. // bufV[tail] = in[tail]
  1516. //
  1517. cmRC_t cmOlaExecS( cmOla* p, const cmSample_t* sp, unsigned sN )
  1518. {
  1519. assert( sN == p->wndSmpCnt );
  1520. cmRC_t rc = cmOkRC;
  1521. const cmSample_t* ep = sp + sN;
  1522. const cmSample_t* wp = p->wf.wndV;
  1523. int i,j,k,n;
  1524. // [Sum head of incoming samples with tail of ola buf]
  1525. // fill outV with the bufV[idx:idx+hopSmpCnt] + sp[hopSmpCnt]
  1526. for(i=0; i<p->hopSmpCnt; ++i)
  1527. {
  1528. p->outV[i] = p->bufV[p->idx++] + (*sp++ * *wp++);
  1529. if( p->idx == p->wndSmpCnt )
  1530. p->idx = 0;
  1531. }
  1532. // [Sum middle of incoming samples with middle of ola buf]
  1533. // sum next wndSmpCnt - hopSmpCnt samples of sp[] into bufV[]
  1534. n = p->wndSmpCnt - (2*p->hopSmpCnt);
  1535. k = p->idx;
  1536. for(j=0; j<n; ++j)
  1537. {
  1538. p->bufV[k++] += (*sp++ * *wp++);
  1539. if( k == p->wndSmpCnt )
  1540. k = 0;
  1541. }
  1542. // [Assign tail of incoming to tail of ola buf]
  1543. // assign ending samples from sp[] into bufV[]
  1544. while( sp < ep )
  1545. {
  1546. p->bufV[k++] = (*sp++ * *wp++);
  1547. if( k == p->wndSmpCnt )
  1548. k = 0;
  1549. }
  1550. p->outPtr = p->outV;
  1551. return rc;
  1552. }
  1553. cmRC_t cmOlaExecR( cmOla* p, const cmReal_t* sp, unsigned sN )
  1554. {
  1555. assert( sN == p->wndSmpCnt );
  1556. cmRC_t rc = cmOkRC;
  1557. const cmReal_t* ep = sp + sN;
  1558. const cmSample_t* wp = p->wf.wndV;
  1559. int i,j,k,n;
  1560. // fill outV with the bufV[idx:idx+hopSmpCnt] + sp[hopSmpCnt]
  1561. for(i=0; i<p->hopSmpCnt; ++i)
  1562. {
  1563. p->outV[i] = p->bufV[p->idx++] + (*sp++ * *wp++);
  1564. if( p->idx == p->wndSmpCnt )
  1565. p->idx = 0;
  1566. }
  1567. // sum next wndSmpCnt - hopSmpCnt samples of sp[] into bufV[]
  1568. n = p->wndSmpCnt - (2*p->hopSmpCnt);
  1569. k = p->idx;
  1570. for(j=0; j<n; ++j)
  1571. {
  1572. p->bufV[k++] += (*sp++ * *wp++);
  1573. if( k == p->wndSmpCnt )
  1574. k = 0;
  1575. }
  1576. // assign ending samples from sp[] into bufV[]
  1577. while( sp < ep )
  1578. {
  1579. p->bufV[k++] = (*sp++ * *wp++);
  1580. if( k == p->wndSmpCnt )
  1581. k = 0;
  1582. }
  1583. p->outPtr = p->outV;
  1584. return rc;
  1585. }
  1586. const cmSample_t* cmOlaExecOut(cmOla* p)
  1587. {
  1588. const cmSample_t* sp = p->outPtr;
  1589. if( sp >= p->outV + p->hopSmpCnt )
  1590. return NULL;
  1591. p->outPtr += p->procSmpCnt;
  1592. return sp;
  1593. }
  1594. //------------------------------------------------------------------------------------------------------------
  1595. cmPhsToFrq* cmPhsToFrqAlloc( cmCtx* c, cmPhsToFrq* ap, double srate, unsigned binCnt, unsigned hopSmpCnt )
  1596. {
  1597. cmPhsToFrq* p = cmObjAlloc( cmPhsToFrq, c, ap );
  1598. if( srate != 0 )
  1599. if( cmPhsToFrqInit( p, srate, binCnt, hopSmpCnt ) != cmOkRC )
  1600. cmPhsToFrqFree(&p);
  1601. return p;
  1602. }
  1603. cmRC_t cmPhsToFrqFree( cmPhsToFrq** pp )
  1604. {
  1605. cmRC_t rc = cmOkRC;
  1606. cmPhsToFrq* p = *pp;
  1607. if( pp==NULL || *pp== NULL )
  1608. return rc;
  1609. if((rc = cmPhsToFrqFinal(p)) != cmOkRC )
  1610. return rc;
  1611. cmMemPtrFree(&p->hzV);
  1612. cmMemPtrFree(&p->phsV);
  1613. cmMemPtrFree(&p->wV);
  1614. cmObjFree(pp);
  1615. return rc;
  1616. }
  1617. cmRC_t cmPhsToFrqInit( cmPhsToFrq* p, double srate, unsigned binCnt, unsigned hopSmpCnt )
  1618. {
  1619. cmRC_t rc;
  1620. unsigned i;
  1621. if((rc = cmPhsToFrqFinal(p)) != cmOkRC )
  1622. return rc;
  1623. p->hzV = cmMemResizeZ( cmReal_t, p->hzV, binCnt );
  1624. p->phsV = cmMemResizeZ( cmReal_t, p->phsV, binCnt );
  1625. p->wV = cmMemResizeZ( cmReal_t, p->wV, binCnt );
  1626. p->srate = srate;
  1627. p->binCnt = binCnt;
  1628. p->hopSmpCnt = hopSmpCnt;
  1629. for(i=0; i<binCnt; ++i)
  1630. p->wV[i] = M_PI * i * hopSmpCnt / (binCnt-1);
  1631. return rc;
  1632. }
  1633. cmRC_t cmPhsToFrqFinal(cmPhsToFrq* p )
  1634. { return cmOkRC; }
  1635. cmRC_t cmPhsToFrqExec( cmPhsToFrq* p, const cmReal_t* phsV )
  1636. {
  1637. cmRC_t rc = cmOkRC;
  1638. unsigned i;
  1639. double twoPi = 2.0 * M_PI;
  1640. double den = twoPi * p->hopSmpCnt;
  1641. for(i=0; i<p->binCnt; ++i)
  1642. {
  1643. cmReal_t dPhs = phsV[i] - p->phsV[i];
  1644. // unwrap phase - see phase_study.m for explanation
  1645. cmReal_t k = round( (p->wV[i] - dPhs) / twoPi);
  1646. // convert phase change to Hz
  1647. p->hzV[i] = (k * twoPi + dPhs) * p->srate / den;
  1648. // store phase for next iteration
  1649. p->phsV[i] = phsV[i];
  1650. }
  1651. return rc;
  1652. }
  1653. //------------------------------------------------------------------------------------------------------------
  1654. cmPvAnl* cmPvAnlAlloc( cmCtx* ctx, cmPvAnl* ap, unsigned procSmpCnt, double srate, unsigned wndSmpCnt, unsigned hopSmpCnt, unsigned flags )
  1655. {
  1656. cmRC_t rc;
  1657. cmPvAnl* p = cmObjAlloc( cmPvAnl, ctx, ap );
  1658. cmShiftBufAlloc(ctx, &p->sb, procSmpCnt, wndSmpCnt, hopSmpCnt );
  1659. cmWndFuncAlloc( ctx, &p->wf, kHannWndId, wndSmpCnt, 0);
  1660. cmFftAllocSR( ctx, &p->ft, p->wf.outV, wndSmpCnt, kToPolarFftFl);
  1661. cmPhsToFrqAlloc(ctx, &p->pf, srate, p->ft.binCnt, hopSmpCnt );
  1662. if( procSmpCnt > 0 )
  1663. if((rc = cmPvAnlInit(p,procSmpCnt,srate,wndSmpCnt,hopSmpCnt,flags)) != cmOkRC )
  1664. cmPvAnlFree(&p);
  1665. return p;
  1666. }
  1667. cmRC_t cmPvAnlFree( cmPvAnl** pp )
  1668. {
  1669. cmRC_t rc;
  1670. if( pp==NULL || *pp==NULL )
  1671. return cmOkRC;
  1672. cmPvAnl* p = *pp;
  1673. if((rc = cmPvAnlFinal(p) ) != cmOkRC )
  1674. return rc;
  1675. cmObjFreeStatic( cmPhsToFrqFree, cmPhsToFrq, p->pf );
  1676. cmObjFreeStatic( cmFftFreeSR, cmFftSR, p->ft );
  1677. cmObjFreeStatic( cmWndFuncFree, cmWndFunc, p->wf );
  1678. cmObjFreeStatic( cmShiftBufFree, cmShiftBuf, p->sb );
  1679. cmObjFree(pp);
  1680. return rc;
  1681. }
  1682. cmRC_t cmPvAnlInit( cmPvAnl* p, unsigned procSmpCnt, double srate, unsigned wndSmpCnt, unsigned hopSmpCnt, unsigned flags )
  1683. {
  1684. cmRC_t rc;
  1685. if((rc = cmPvAnlFinal(p)) != cmOkRC )
  1686. return rc;
  1687. if((rc = cmShiftBufInit( &p->sb, procSmpCnt, wndSmpCnt, hopSmpCnt )) != cmOkRC )
  1688. return rc;
  1689. if((rc = cmWndFuncInit( &p->wf, kHannWndId | kNormByLengthWndFl, wndSmpCnt, 0)) != cmOkRC )
  1690. return rc;
  1691. if((rc = cmFftInitSR( &p->ft, p->wf.outV, wndSmpCnt, kToPolarFftFl)) != cmOkRC )
  1692. return rc;
  1693. if((rc = cmPhsToFrqInit( &p->pf, srate, p->ft.binCnt, hopSmpCnt)) != cmOkRC )
  1694. return rc;
  1695. // if the window was just initialized
  1696. // divide the window to indirectly apply the magnitude normalization
  1697. //if( p->wndSmpCnt != wndSmpCnt )
  1698. // cmVOS_DivVS( p->wf.wndV, p->wf.outN, wndSmpCnt );
  1699. p->flags = flags;
  1700. p->procSmpCnt = procSmpCnt;
  1701. p->wndSmpCnt = wndSmpCnt;
  1702. p->hopSmpCnt = hopSmpCnt;
  1703. p->binCnt = p->ft.binCnt;
  1704. p->magV = p->ft.magV;
  1705. p->phsV = p->ft.phsV;
  1706. p->hzV = p->pf.hzV;
  1707. return rc;
  1708. }
  1709. cmRC_t cmPvAnlFinal(cmPvAnl* p )
  1710. { return cmOkRC; }
  1711. bool cmPvAnlExec( cmPvAnl* p, const cmSample_t* x, unsigned xN )
  1712. {
  1713. bool fl = false;
  1714. while( cmShiftBufExec(&p->sb,x,xN) )
  1715. {
  1716. cmWndFuncExec(&p->wf, p->sb.outV, p->sb.wndSmpCnt );
  1717. cmFftExecSR(&p->ft,NULL,0);
  1718. if( cmIsFlag(p->flags,kCalcHzPvaFl) )
  1719. cmPhsToFrqExec(&p->pf,p->phsV);
  1720. fl = true;
  1721. }
  1722. return fl;
  1723. }
  1724. //------------------------------------------------------------------------------------------------------------
  1725. cmPvSyn* cmPvSynAlloc( cmCtx* ctx, cmPvSyn* ap, unsigned procSmpCnt, double outSrate, unsigned wndSmpCnt, unsigned hopSmpCnt, unsigned wndTypeId )
  1726. {
  1727. cmRC_t rc;
  1728. cmPvSyn* p = cmObjAlloc( cmPvSyn, ctx, ap );
  1729. cmWndFuncAlloc( ctx, &p->wf, kHannWndId, wndSmpCnt, 0);
  1730. cmIFftAllocRS( ctx, &p->ft, wndSmpCnt/2+1 );
  1731. cmOlaAlloc( ctx, &p->ola, wndSmpCnt, hopSmpCnt, procSmpCnt, wndTypeId );
  1732. if( procSmpCnt )
  1733. if((rc = cmPvSynInit(p,procSmpCnt,outSrate,wndSmpCnt,hopSmpCnt,wndTypeId)) != cmOkRC )
  1734. cmPvSynFree(&p);
  1735. return p;
  1736. }
  1737. cmRC_t cmPvSynFree( cmPvSyn** pp )
  1738. {
  1739. cmRC_t rc;
  1740. if( pp==NULL || *pp==NULL )
  1741. return cmOkRC;
  1742. cmPvSyn* p = *pp;
  1743. if((rc = cmPvSynFinal(p)) != cmOkRC )
  1744. return rc;
  1745. cmMemPtrFree(&p->minRphV);
  1746. cmMemPtrFree(&p->maxRphV);
  1747. cmMemPtrFree(&p->itrV);
  1748. cmMemPtrFree(&p->phs0V);
  1749. cmMemPtrFree(&p->phsV);
  1750. cmMemPtrFree(&p->mag0V);
  1751. cmMemPtrFree(&p->magV);
  1752. cmObjFreeStatic( cmOlaFree, cmOla, p->ola);
  1753. cmObjFreeStatic( cmIFftFreeRS, cmIFftRS, p->ft );
  1754. cmObjFreeStatic( cmWndFuncFree, cmWndFunc, p->wf );
  1755. cmObjFree(pp);
  1756. return cmOkRC;
  1757. }
  1758. cmRC_t cmPvSynInit( cmPvSyn* p, unsigned procSmpCnt, double outSrate, unsigned wndSmpCnt, unsigned hopSmpCnt, unsigned wndTypeId )
  1759. {
  1760. cmRC_t rc;
  1761. int k;
  1762. double twoPi = 2.0 * M_PI;
  1763. bool useHannFl = true;
  1764. int m = useHannFl ? 2 : 1;
  1765. if((rc = cmPvSynFinal(p)) != cmOkRC )
  1766. return rc;
  1767. p->outSrate = outSrate;
  1768. p->procSmpCnt = procSmpCnt;
  1769. p->wndSmpCnt = wndSmpCnt;
  1770. p->hopSmpCnt = hopSmpCnt;
  1771. p->binCnt = wndSmpCnt / 2 + 1;
  1772. p->minRphV = cmMemResizeZ( cmReal_t, p->minRphV, p->binCnt );
  1773. p->maxRphV = cmMemResizeZ( cmReal_t, p->maxRphV, p->binCnt );
  1774. p->itrV = cmMemResizeZ( cmReal_t, p->itrV, p->binCnt );
  1775. p->phs0V = cmMemResizeZ( cmReal_t, p->phs0V, p->binCnt );
  1776. p->phsV = cmMemResizeZ( cmReal_t, p->phsV, p->binCnt );
  1777. p->mag0V = cmMemResizeZ( cmReal_t, p->mag0V, p->binCnt );
  1778. p->magV = cmMemResizeZ( cmReal_t, p->magV, p->binCnt );
  1779. if((rc = cmWndFuncInit( &p->wf, wndTypeId, wndSmpCnt, 0)) != cmOkRC )
  1780. return rc;
  1781. if((rc = cmIFftInitRS( &p->ft, p->binCnt )) != cmOkRC )
  1782. return rc;
  1783. if((rc = cmOlaInit( &p->ola, wndSmpCnt, hopSmpCnt, procSmpCnt, wndTypeId )) != cmOkRC )
  1784. return rc;
  1785. for(k=0; k<p->binCnt; ++k)
  1786. {
  1787. // complete revolutions per hop in radians
  1788. p->itrV[k] = twoPi * floor((double)k * hopSmpCnt / wndSmpCnt );
  1789. p->minRphV[k] = ((cmReal_t)(k-m)) * hopSmpCnt * twoPi / wndSmpCnt;
  1790. p->maxRphV[k] = ((cmReal_t)(k+m)) * hopSmpCnt * twoPi / wndSmpCnt;
  1791. //printf("%f %f %f\n",p->itrV[k],p->minRphV[k],p->maxRphV[k]);
  1792. }
  1793. return rc;
  1794. }
  1795. cmRC_t cmPvSynFinal(cmPvSyn* p )
  1796. { return cmOkRC; }
  1797. cmRC_t cmPvSynExec( cmPvSyn* p, const cmReal_t* magV, const cmReal_t* phsV )
  1798. {
  1799. double twoPi = 2.0 * M_PI;
  1800. unsigned k;
  1801. for(k=0; k<p->binCnt; ++k)
  1802. {
  1803. // phase dist between cur and prv frame
  1804. cmReal_t dp = phsV[k] - p->phs0V[k];
  1805. // dist must be positive (accum phase always increases)
  1806. if( dp < -0.00001 )
  1807. dp += twoPi;
  1808. // add in complete revolutions based on the bin frequency
  1809. // (these would have been lost from 'dp' due to phase wrap)
  1810. dp += p->itrV[k];
  1811. // constrain the phase change to lie within the range of the kth bin
  1812. if( dp < p->minRphV[k] )
  1813. dp += twoPi;
  1814. if( dp > p->maxRphV[k] )
  1815. dp -= twoPi;
  1816. p->phsV[k] = p->phs0V[k] + dp;
  1817. p->magV[k] = p->mag0V[k];
  1818. p->phs0V[k] = phsV[k];
  1819. p->mag0V[k] = magV[k];
  1820. }
  1821. cmIFftExecPolarRS( &p->ft, magV, phsV );
  1822. cmOlaExecS( &p->ola, p->ft.outV, p->ft.outN );
  1823. //printf("%i %i\n",p->binCnt,p->ft.binCnt );
  1824. //cmVOR_Print( p->obj.ctx->outFuncPtr, 1, p->binCnt, magV );
  1825. //cmVOR_Print( p->obj.ctx->outFuncPtr, 1, p->binCnt, p->phsV );
  1826. //cmVOS_Print( p->obj.ctx->outFuncPtr, 1, 10, p->ft.outV );
  1827. return cmOkRC;
  1828. }
  1829. cmRC_t cmPvSynDoIt( cmPvSyn* p, const cmSample_t* v )
  1830. {
  1831. cmOlaExecS( &p->ola, v, p->wndSmpCnt );
  1832. //printf("%f\n",cmVOS_RMS(s,p->wndSmpCnt,p->wndSmpCnt));
  1833. return cmOkRC;
  1834. }
  1835. const cmSample_t* cmPvSynExecOut(cmPvSyn* p )
  1836. { return cmOlaExecOut(&p->ola); }
  1837. //------------------------------------------------------------------------------------------------------------
  1838. cmMidiSynth* cmMidiSynthAlloc( cmCtx* ctx, cmMidiSynth* ap, const cmMidiSynthPgm* pgmArray, unsigned pgmCnt, unsigned voiceCnt, unsigned procSmpCnt, unsigned outChCnt, cmReal_t srate )
  1839. {
  1840. cmMidiSynth* p = cmObjAlloc( cmMidiSynth, ctx, ap );
  1841. if( pgmArray != NULL )
  1842. if( cmMidiSynthInit( p, pgmArray, pgmCnt, voiceCnt, procSmpCnt, outChCnt, srate ) != cmOkRC )
  1843. cmMidiSynthFree(&p);
  1844. return p;
  1845. }
  1846. cmRC_t cmMidiSynthFree( cmMidiSynth** pp )
  1847. {
  1848. cmRC_t rc;
  1849. if( pp==NULL || *pp==NULL)
  1850. return cmOkRC;
  1851. cmMidiSynth* p = *pp;
  1852. if((rc = cmMidiSynthFinal(p)) != cmOkRC )
  1853. return rc;
  1854. cmMemPtrFree(&p->voiceArray);
  1855. cmMemPtrFree(&p->outM);
  1856. cmMemPtrFree(&p->outChArray);
  1857. cmObjFree(pp);
  1858. return cmOkRC;
  1859. }
  1860. cmRC_t cmMidiSynthInit( cmMidiSynth* p, const cmMidiSynthPgm* pgmArray, unsigned pgmCnt, unsigned voiceCnt, unsigned procSmpCnt, unsigned outChCnt, cmReal_t srate )
  1861. {
  1862. // at least one pgm must be given
  1863. assert( pgmCnt > 0 );
  1864. unsigned i;
  1865. cmRC_t rc;
  1866. if((rc = cmMidiSynthFinal(p)) != cmOkRC )
  1867. return rc;
  1868. p->voiceArray = cmMemResizeZ( cmMidiVoice, p->voiceArray, voiceCnt );
  1869. p->outM = cmMemResizeZ( cmSample_t, p->outM, outChCnt * procSmpCnt );
  1870. p->outChArray = cmMemResizeZ( cmSample_t*, p->outChArray, outChCnt );
  1871. p->avail = p->voiceArray;
  1872. p->voiceCnt = voiceCnt;
  1873. p->activeVoiceCnt = 0;
  1874. p->voiceStealCnt = 0;
  1875. p->procSmpCnt = procSmpCnt;
  1876. p->outChCnt = outChCnt;
  1877. p->srate = srate;
  1878. for(i=0; i<outChCnt; ++i)
  1879. p->outChArray[i] = p->outM + (i*procSmpCnt);
  1880. for(i=0; i<kMidiChCnt; ++i)
  1881. {
  1882. p->chArray[i].pgm = 0;
  1883. p->chArray[i].active = NULL;
  1884. p->chArray[i].pitchBend = 0;
  1885. p->chArray[i].synthPtr = p;
  1886. memset(p->chArray[i].midiCtl, 0, kMidiCtlCnt * sizeof(cmMidiByte_t));
  1887. }
  1888. for(i=0; i<voiceCnt; ++i)
  1889. {
  1890. p->voiceArray[i].index = i;
  1891. p->voiceArray[i].flags = 0;
  1892. p->voiceArray[i].pitch = kInvalidMidiPitch;
  1893. p->voiceArray[i].velocity = kInvalidMidiVelocity;
  1894. p->voiceArray[i].pgm.pgm = kInvalidMidiPgm;
  1895. p->voiceArray[i].pgm.cbPtr = NULL;
  1896. p->voiceArray[i].pgm.cbDataPtr = NULL;
  1897. p->voiceArray[i].chPtr = NULL;
  1898. p->voiceArray[i].link = i<voiceCnt-1 ? p->voiceArray + i + 1 : NULL;
  1899. }
  1900. for(i=0; i<pgmCnt; ++i)
  1901. {
  1902. unsigned idx = pgmArray[i].pgm;
  1903. if( idx >= kMidiPgmCnt )
  1904. rc = cmCtxRtCondition( &p->obj, cmArgAssertRC, "MIDI program change values must be less than %i.",kMidiPgmCnt);
  1905. else
  1906. {
  1907. p->pgmArray[ idx ].cbPtr = pgmArray[i].cbPtr;
  1908. p->pgmArray[ idx ].cbDataPtr = pgmArray[i].cbDataPtr;
  1909. p->pgmArray[ idx ].pgm = idx;
  1910. }
  1911. }
  1912. return rc;
  1913. }
  1914. cmRC_t cmMidiSynthFinal( cmMidiSynth* p )
  1915. { return cmOkRC; }
  1916. cmRC_t _cmMidiSynthOnNoteOn( cmMidiSynth* p, cmMidiByte_t ch, cmMidiByte_t pitch, cmMidiByte_t vel )
  1917. {
  1918. assert( ch < kMidiChCnt );
  1919. if( p->activeVoiceCnt == p->voiceCnt )
  1920. {
  1921. ++p->voiceStealCnt;
  1922. return cmOkRC;
  1923. }
  1924. assert( p->avail != NULL );
  1925. cmMidiSynthCh* chPtr = p->chArray + ch;
  1926. cmMidiVoice* vp = p->avail;
  1927. ++p->activeVoiceCnt;
  1928. // update avail
  1929. p->avail = p->avail->link;
  1930. // update active
  1931. vp->flags |= kActiveMsFl | kKeyGateMsFl;
  1932. vp->pitch = pitch;
  1933. vp->velocity = vel;
  1934. vp->pgm = p->pgmArray[ chPtr->pgm ];
  1935. vp->chPtr = chPtr;
  1936. vp->link = chPtr->active;
  1937. chPtr->active = vp;
  1938. vp->pgm.cbPtr( vp, kAttackMsId, NULL, 0 );
  1939. return cmOkRC;
  1940. }
  1941. cmRC_t _cmMidiSynthOnNoteOff( cmMidiSynth* p, cmMidiByte_t ch, cmMidiByte_t pitch, cmMidiByte_t vel )
  1942. {
  1943. assert( ch < kMidiChCnt );
  1944. cmMidiSynthCh* cp = p->chArray + ch;
  1945. cmMidiVoice* vp = cp->active;
  1946. // find the voice for the given pitch
  1947. while( vp != NULL )
  1948. {
  1949. if( (vp->pitch == pitch) && (cmIsFlag(vp->flags,kKeyGateMsFl)==true) )
  1950. break;
  1951. vp = vp->link;
  1952. }
  1953. // if no voice had a key down on this pitch
  1954. if( vp == NULL )
  1955. {
  1956. return cmOkRC;
  1957. }
  1958. // mark the key as 'up'
  1959. vp->flags = cmClrFlag(vp->flags,kKeyGateMsFl);
  1960. // if the sustain pedal is up
  1961. if( cp->midiCtl[ kSustainCtlMdId ] == 0 )
  1962. vp->pgm.cbPtr( vp, kReleaseMsId, NULL, 0 );
  1963. return cmOkRC;
  1964. }
  1965. cmRC_t _cmMidiSynthOnCtl( cmMidiSynth* p, cmMidiByte_t ch, cmMidiByte_t ctlId, cmMidiByte_t ctlValue )
  1966. {
  1967. assert( ch < kMidiChCnt && ctlId < kMidiCtlCnt );
  1968. cmMidiSynthCh* cp = p->chArray + ch;
  1969. cp->midiCtl[ ctlId ] = ctlValue;
  1970. // if the sustain pedal is going up
  1971. if( ctlId == kSustainCtlMdId && ctlValue == 0 )
  1972. {
  1973. cmMidiVoice* vp = cp->active;
  1974. while(vp != NULL)
  1975. {
  1976. if( cmIsFlag(vp->flags,kKeyGateMsFl)==false )
  1977. vp->pgm.cbPtr(vp, kReleaseMsId, NULL, 0 );
  1978. vp = vp->link;
  1979. }
  1980. }
  1981. //printf("%i %i %f\n",ctlId,ctlValue,ctlValue/127.0);
  1982. return cmOkRC;
  1983. }
  1984. cmRC_t cmMidiSynthOnMidi( cmMidiSynth* p, const cmMidiPacket_t* pktArray, unsigned pktCnt )
  1985. {
  1986. unsigned i=0;
  1987. for(i=0; i<pktCnt; ++i)
  1988. if( pktArray[i].msgArray != NULL )
  1989. {
  1990. unsigned j;
  1991. for(j=0; j<pktArray[i].msgCnt; ++j)
  1992. {
  1993. const cmMidiMsg* mp = pktArray[i].msgArray + j;
  1994. cmMidiByte_t ch = mp->status & 0x0f;
  1995. cmMidiByte_t status = mp->status & 0xf0;
  1996. switch( status )
  1997. {
  1998. case kNoteOnMdId:
  1999. if( mp->d1 != 0 )
  2000. {
  2001. _cmMidiSynthOnNoteOn(p,ch,mp->d0,mp->d1);
  2002. break;
  2003. }
  2004. // fall through
  2005. case kNoteOffMdId:
  2006. _cmMidiSynthOnNoteOff(p,ch,mp->d0,mp->d1);
  2007. break;
  2008. case kPolyPresMdId:
  2009. break;
  2010. case kCtlMdId:
  2011. _cmMidiSynthOnCtl( p, ch, mp->d0, mp->d1 );
  2012. break;
  2013. case kPgmMdId:
  2014. break;
  2015. case kChPresMdId:
  2016. break;
  2017. case kPbendMdId:
  2018. break;
  2019. default:
  2020. printf("Unknown MIDI status:%i %i\n",(int)status,(int)mp->status);
  2021. break;
  2022. }
  2023. }
  2024. }
  2025. return cmOkRC;
  2026. }
  2027. cmRC_t cmMidiSynthExec( cmMidiSynth* p, cmSample_t* outChArray[], unsigned outChCnt )
  2028. {
  2029. unsigned i;
  2030. cmSample_t** chArray = outChArray == NULL ? p->outChArray : outChArray;
  2031. unsigned chCnt = outChArray == NULL ? p->outChCnt : outChCnt;
  2032. // FIX: make one active chain attached to cmMidiSynth rather than many
  2033. // active chains attached to each channel - this will avoid the extra
  2034. // iterations below.
  2035. // for each channel
  2036. for(i=0; i<kMidiChCnt; ++i)
  2037. {
  2038. cmMidiVoice* vp = p->chArray[i].active;
  2039. cmMidiVoice* prv = NULL;
  2040. // for each voice assigned to this channel
  2041. while(vp != NULL)
  2042. {
  2043. // tell the voice to perform its DSP function - returns 0 if the voice is no longer active
  2044. if( vp->pgm.cbPtr( vp, kDspMsId, chArray, chCnt ) )
  2045. {
  2046. prv = vp;
  2047. vp = vp->link;
  2048. }
  2049. else
  2050. {
  2051. cmMidiVoice* nvp = vp->link;
  2052. // remove vp from the active chain
  2053. if( prv != NULL )
  2054. prv->link = vp->link;
  2055. else
  2056. {
  2057. assert( vp == p->chArray[i].active );
  2058. // vp is first recd on active chain, nvp becomes first ...
  2059. p->chArray[i].active = vp->link;
  2060. prv = NULL; // ... so prv must be NULL
  2061. }
  2062. // insert this voice on the available chain
  2063. vp->link = p->avail;
  2064. p->avail = vp;
  2065. --p->activeVoiceCnt;
  2066. vp = nvp;
  2067. }
  2068. }
  2069. }
  2070. return cmOkRC;
  2071. }
  2072. //------------------------------------------------------------------------------------------------------------
  2073. cmWtVoice* cmWtVoiceAlloc( cmCtx* ctx, cmWtVoice* ap, unsigned procSmpCnt, cmReal_t hz )
  2074. {
  2075. cmWtVoice* p = cmObjAlloc( cmWtVoice, ctx, ap );
  2076. if( procSmpCnt != 0 )
  2077. if( cmWtVoiceInit( p, procSmpCnt, hz ) != cmOkRC )
  2078. cmWtVoiceFree(&p);
  2079. return p;
  2080. }
  2081. cmRC_t cmWtVoiceFree( cmWtVoice** pp )
  2082. {
  2083. cmRC_t rc = cmOkRC;
  2084. if( pp==NULL || *pp==NULL )
  2085. return cmOkRC;
  2086. cmWtVoice* p = *pp;
  2087. if((rc = cmWtVoiceFinal(p)) != cmOkRC )
  2088. return rc;
  2089. cmMemPtrFree(&p->outV);
  2090. cmObjFree(pp);
  2091. return rc;
  2092. }
  2093. cmRC_t cmWtVoiceInit( cmWtVoice* p, unsigned procSmpCnt, cmReal_t hz )
  2094. {
  2095. p->outV = cmMemResizeZ( cmSample_t, p->outV, procSmpCnt );
  2096. p->outN = procSmpCnt;
  2097. p->hz = hz;
  2098. p->level = 0;
  2099. p->durSmpCnt = 0;
  2100. p->phase = 0;
  2101. p->state = kOffWtId;
  2102. return cmOkRC;
  2103. }
  2104. cmRC_t cmWtVoiceFinal( cmWtVoice* p )
  2105. { return cmOkRC; }
  2106. int cmWtVoiceExec( cmWtVoice* p, struct cmMidiVoice_str* mvp, unsigned sel, cmSample_t* outChArray[], unsigned outChCnt )
  2107. {
  2108. switch( sel )
  2109. {
  2110. case kAttackMsId:
  2111. p->state = kAtkWtId;
  2112. p->hz = (13.75 * pow(2,(-9.0/12.0))) * pow(2,((double)mvp->pitch / 12));
  2113. //printf("%fhz\n",p->hz);
  2114. break;
  2115. case kReleaseMsId:
  2116. p->state = kRlsWtId;
  2117. //printf("rls:%f\n",p->phase);
  2118. break;
  2119. case kDspMsId:
  2120. {
  2121. if( p->state == kRlsWtId )
  2122. {
  2123. p->state = kOffWtId;
  2124. return 0;
  2125. }
  2126. cmMidiSynth* sp = mvp->chPtr->synthPtr;
  2127. cmSample_t* dp = outChCnt == 0 ? p->outV : outChArray[0];
  2128. cmSample_t* ep = dp + p->outN;
  2129. cmReal_t rps = (2.0 * M_PI * p->hz) / sp->srate;
  2130. cmReal_t sum=0;
  2131. unsigned i=0;
  2132. for(; dp < ep; ++dp)
  2133. {
  2134. *dp += (cmSample_t)(0.5 * sin( p->phase ));
  2135. sum += *dp;
  2136. ++i;
  2137. p->phase += rps;
  2138. }
  2139. //printf("(%f %f %i %i %p) ",p->phase,sum,i,p->outN,outChArray[0] );
  2140. }
  2141. break;
  2142. default:
  2143. assert(0);
  2144. break;
  2145. }
  2146. return 1;
  2147. }
  2148. //------------------------------------------------------------------------------------------------------------
  2149. cmWtVoiceBank* cmWtVoiceBankAlloc( cmCtx* ctx, cmWtVoiceBank* ap, double srate, unsigned procSmpCnt, unsigned voiceCnt, unsigned chCnt )
  2150. {
  2151. cmWtVoiceBank* p = cmObjAlloc( cmWtVoiceBank, ctx, ap );
  2152. if( srate != 0 )
  2153. if( cmWtVoiceBankInit( p, srate, procSmpCnt, voiceCnt, chCnt ) != cmOkRC )
  2154. cmWtVoiceBankFree(&p);
  2155. return p;
  2156. }
  2157. cmRC_t cmWtVoiceBankFree( cmWtVoiceBank** pp )
  2158. {
  2159. cmRC_t rc;
  2160. if( pp==NULL || *pp==NULL)
  2161. return cmOkRC;
  2162. cmWtVoiceBank* p = *pp;
  2163. if((rc = cmWtVoiceBankFinal(p)) != cmOkRC )
  2164. return rc;
  2165. cmMemPtrFree(&p->voiceArray);
  2166. cmMemPtrFree(&p->chArray);
  2167. cmMemPtrFree(&p->buf);
  2168. cmObjFree(pp);
  2169. return rc;
  2170. }
  2171. cmRC_t cmWtVoiceBankInit( cmWtVoiceBank* p, double srate, unsigned procSmpCnt, unsigned voiceCnt, unsigned chCnt )
  2172. {
  2173. cmRC_t rc;
  2174. unsigned i;
  2175. if((rc = cmWtVoiceBankFinal(p)) != cmOkRC )
  2176. return rc;
  2177. p->voiceArray = cmMemResizeZ( cmWtVoice*, p->voiceArray, voiceCnt );
  2178. for(i=0; i<voiceCnt; ++i)
  2179. p->voiceArray[i] = cmWtVoiceAlloc(p->obj.ctx,NULL,procSmpCnt,0);
  2180. p->voiceCnt = voiceCnt;
  2181. p->buf = cmMemResizeZ( cmSample_t, p->buf, chCnt * procSmpCnt );
  2182. p->chArray = cmMemResizeZ( cmSample_t*, p->chArray, chCnt );
  2183. for(i=0; i<chCnt; ++i)
  2184. p->chArray[i] = p->buf + (i*procSmpCnt);
  2185. p->chCnt = chCnt;
  2186. p->procSmpCnt = procSmpCnt;
  2187. return cmOkRC;
  2188. }
  2189. cmRC_t cmWtVoiceBankFinal( cmWtVoiceBank* p )
  2190. {
  2191. unsigned i;
  2192. for(i=0; i<p->voiceCnt; ++i)
  2193. cmWtVoiceFree(&p->voiceArray[i]);
  2194. return cmOkRC;
  2195. }
  2196. int cmWtVoiceBankExec( cmWtVoiceBank* p, struct cmMidiVoice_str* voicePtr, unsigned sel, cmSample_t* outChArray[], unsigned outChCnt )
  2197. {
  2198. cmWtVoice* vp = p->voiceArray[ voicePtr->index ];
  2199. bool fl = outChArray==NULL || outChCnt==0;
  2200. cmSample_t** chArray = fl ? p->chArray : outChArray;
  2201. unsigned chCnt = fl ? p->chCnt : outChCnt;
  2202. return cmWtVoiceExec( vp, voicePtr, sel, chArray, chCnt );
  2203. }
  2204. //------------------------------------------------------------------------------------------------------------
  2205. cmAudioFileBuf* cmAudioFileBufAlloc( cmCtx* ctx, cmAudioFileBuf* ap, unsigned procSmpCnt, const char* fn, unsigned audioChIdx, unsigned begSmpIdx, unsigned durSmpCnt )
  2206. {
  2207. cmAudioFileBuf* p = cmObjAlloc( cmAudioFileBuf, ctx, ap );
  2208. if( procSmpCnt != 0 )
  2209. if( cmAudioFileBufInit( p, procSmpCnt, fn, audioChIdx, begSmpIdx, durSmpCnt ) != cmOkRC )
  2210. cmAudioFileBufFree(&p);
  2211. return p;
  2212. }
  2213. cmRC_t cmAudioFileBufFree( cmAudioFileBuf** pp )
  2214. {
  2215. cmRC_t rc;
  2216. if( pp==NULL || *pp==NULL)
  2217. return cmOkRC;
  2218. cmAudioFileBuf* p = *pp;
  2219. if((rc = cmAudioFileBufFinal(p)) != cmOkRC )
  2220. return rc;
  2221. cmMemPtrFree(&p->bufV);
  2222. cmMemPtrFree(&p->fn);
  2223. cmObjFree(pp);
  2224. return rc;
  2225. }
  2226. cmRC_t cmAudioFileBufInit( cmAudioFileBuf* p, unsigned procSmpCnt, const char* fn, unsigned audioChIdx, unsigned begSmpIdx, unsigned durSmpCnt )
  2227. {
  2228. cmAudioFileH_t afH;
  2229. cmRC_t rc;
  2230. if((rc = cmAudioFileBufFinal(p)) != cmOkRC )
  2231. return rc;
  2232. // open the audio file for reading
  2233. if( cmAudioFileIsValid( afH = cmAudioFileNewOpen( fn, &p->info, &rc, p->obj.err.rpt ))==false || rc != kOkAfRC )
  2234. return cmCtxRtCondition(&p->obj, cmArgAssertRC,"The audio file '%s' could not be opend.",fn );
  2235. // validate the audio channel
  2236. if( audioChIdx >= p->info.chCnt )
  2237. return cmCtxRtCondition(&p->obj, cmArgAssertRC,"The audio file channel index %i is out of range for the audio file '%s'.",audioChIdx,fn);
  2238. // validate the start sample index
  2239. if( begSmpIdx > p->info.frameCnt )
  2240. return cmCtxRtCondition(&p->obj, cmOkRC, "The start sample index %i is past the end of the audio file '%s'.",begSmpIdx,fn);
  2241. if( durSmpCnt == cmInvalidCnt )
  2242. durSmpCnt = p->info.frameCnt - begSmpIdx;
  2243. // validate the duration
  2244. if( begSmpIdx + durSmpCnt > p->info.frameCnt )
  2245. {
  2246. unsigned newDurSmpCnt = p->info.frameCnt - begSmpIdx;
  2247. cmCtxRtCondition(&p->obj, cmOkRC, "The selected sample duration %i is past the end of the audio file '%s' and has been shorted to %i samples.",durSmpCnt,fn,newDurSmpCnt);
  2248. durSmpCnt = newDurSmpCnt;
  2249. }
  2250. // seek to the starting sample
  2251. if( cmAudioFileSeek( afH, begSmpIdx ) != kOkAfRC )
  2252. return cmCtxRtCondition(&p->obj, cmArgAssertRC,"Seek to sample index %i failed on the audio file '%s'.",begSmpIdx,fn);
  2253. // allocate the buffer memory
  2254. p->bufV = cmMemResizeZ( cmSample_t, p->bufV, durSmpCnt );
  2255. p->fn = cmMemResize( char, p->fn, strlen(fn)+1 );
  2256. p->bufN = durSmpCnt;
  2257. p->begSmpIdx = begSmpIdx;
  2258. p->chIdx = audioChIdx;
  2259. strcpy(p->fn,fn);
  2260. cmSample_t* outV = p->bufV;
  2261. // read the file into the buffer
  2262. unsigned rdSmpCnt = cmMin(4096,durSmpCnt);
  2263. unsigned cmc = 0;
  2264. while( cmc < durSmpCnt )
  2265. {
  2266. unsigned actualReadCnt = 0;
  2267. unsigned n = rdSmpCnt;
  2268. cmSample_t* chArray[] = {outV};
  2269. if( cmc + n > durSmpCnt )
  2270. n = durSmpCnt - cmc;
  2271. if((rc=cmAudioFileReadSample( afH, n, audioChIdx, 1, chArray, &actualReadCnt)) != kOkAfRC )
  2272. break;
  2273. cmc += actualReadCnt;
  2274. outV += actualReadCnt;
  2275. }
  2276. if( rc==kOkAfRC || (rc != kOkAfRC && cmAudioFileIsEOF(afH)))
  2277. rc = cmOkRC;
  2278. return rc;
  2279. }
  2280. cmRC_t cmAudioFileBufFinal(cmAudioFileBuf* p )
  2281. { return cmOkRC; }
  2282. unsigned cmAudioFileBufExec( cmAudioFileBuf* p, unsigned smpIdx, cmSample_t* outV, unsigned outN, bool sumIntoOutFl )
  2283. {
  2284. if( outV == NULL || outN == 0 || smpIdx >= p->bufN )
  2285. return 0;
  2286. unsigned n = cmMin(outN,p->bufN-smpIdx);
  2287. if( sumIntoOutFl )
  2288. cmVOS_AddVV(outV,n,p->bufV + smpIdx);
  2289. else
  2290. cmVOS_Copy(outV,n,p->bufV + smpIdx );
  2291. if( n < outN )
  2292. memset(outV+n,0,(outN-n)*sizeof(cmSample_t));
  2293. return n;
  2294. }
  2295. //------------------------------------------------------------------------------------------------------------
  2296. cmMDelay* cmMDelayAlloc( cmCtx* ctx, cmMDelay* ap, unsigned procSmpCnt, cmReal_t srate, cmReal_t fbCoeff, unsigned delayCnt, const cmReal_t* delayMsArray, const cmReal_t* delayGainArray )
  2297. {
  2298. cmMDelay* p = cmObjAlloc( cmMDelay, ctx, ap );
  2299. if( procSmpCnt != 0 )
  2300. if( cmMDelayInit( p, procSmpCnt, srate, fbCoeff, delayCnt, delayMsArray, delayGainArray ) != cmOkRC )
  2301. cmMDelayFree(&p);
  2302. return p;
  2303. }
  2304. cmRC_t cmMDelayFree( cmMDelay** pp )
  2305. {
  2306. cmRC_t rc;
  2307. if( pp == NULL || *pp==NULL)
  2308. return cmOkRC;
  2309. cmMDelay* p = *pp;
  2310. if((rc = cmMDelayFinal(p)) != cmOkRC )
  2311. return rc;
  2312. unsigned i;
  2313. for(i=0; i<p->delayCnt; ++i)
  2314. cmMemPtrFree(&p->delayArray[i].delayBuf);
  2315. cmMemPtrFree(&p->delayArray);
  2316. cmMemPtrFree(&p->outV);
  2317. cmObjFree(pp);
  2318. return cmOkRC;
  2319. }
  2320. cmRC_t cmMDelayInit( cmMDelay* p, unsigned procSmpCnt, cmReal_t srate, cmReal_t fbCoeff, unsigned delayCnt, const cmReal_t* delayMsArray, const cmReal_t* delayGainArray )
  2321. {
  2322. cmRC_t rc;
  2323. if((rc = cmMDelayFinal(p)) != cmOkRC )
  2324. return rc;
  2325. if( delayCnt <= 0 )
  2326. return rc;
  2327. p->delayArray = cmMemResizeZ( cmMDelayHead, p->delayArray, delayCnt );
  2328. unsigned i;
  2329. for(i=0; i<delayCnt; ++i)
  2330. {
  2331. p->delayArray[i].delayGain = delayGainArray == NULL ? 1.0 : delayGainArray[i];
  2332. p->delayArray[i].delayMs = delayMsArray[i];
  2333. p->delayArray[i].delaySmpFrac = delayMsArray[i] * srate / 1000.0;
  2334. p->delayArray[i].delayBufSmpCnt = ceil(delayMsArray[i] * srate / 1000)+2;
  2335. p->delayArray[i].delayBuf = cmMemResizeZ( cmSample_t, p->delayArray[i].delayBuf, p->delayArray[i].delayBufSmpCnt );
  2336. p->delayArray[i].inIdx = 0;
  2337. }
  2338. p->delayCnt= delayCnt;
  2339. p->outV = cmMemResizeZ( cmSample_t, p->outV, procSmpCnt );
  2340. p->outN = procSmpCnt;
  2341. p->fbCoeff = fbCoeff;
  2342. p->srate = srate;
  2343. return cmOkRC;
  2344. }
  2345. cmRC_t cmMDelayFinal( cmMDelay* p )
  2346. { return cmOkRC; }
  2347. void _cmMDelayExec( cmMDelay* p, cmMDelayHead* hp, const cmSample_t inV[], cmSample_t outV[], unsigned sigN )
  2348. {
  2349. cmSample_t* dl = hp->delayBuf; // ptr to the base of the delay line
  2350. cmReal_t dfi = (cmReal_t)(hp->inIdx - hp->delaySmpFrac) + hp->delayBufSmpCnt; // fractional delay in samples
  2351. int dii0 = ((int)dfi) % hp->delayBufSmpCnt; // index to the sample just before the delay position
  2352. int dii1 = (dii0 + 1) % hp->delayBufSmpCnt; // index to the sample just after the delay position
  2353. //cmReal_t frac = 0; //dfi - dii0; // interpolation coeff.
  2354. unsigned i;
  2355. for(i=0; i<sigN; i++)
  2356. {
  2357. /*
  2358. outPtr[i] = -(((f+0)*(f-1)*(f-2)/6) * _wtPtr[iPhs0])
  2359. +(((f+1)*(f-1)*(f-2)/2) * _wtPtr[iPhs0+1])
  2360. -(((f+1)*(f-0)*(f-2)/2) * _wtPtr[iPhs0+2])
  2361. +(((f+1)*(f-0)*(f-1)/6) * _wtPtr[iPhs0+3]);
  2362. */
  2363. cmSample_t outSmp = dl[dii0]; // + (frac * (dl[dii1]-dl[dii0]));
  2364. outV[i] += outSmp/p->delayCnt;
  2365. dl[hp->inIdx] = (p->fbCoeff * outSmp) + inV[i];
  2366. hp->inIdx = (hp->inIdx+1) % hp->delayBufSmpCnt;
  2367. dii0 = (dii0+1) % hp->delayBufSmpCnt;
  2368. dii1 = (dii1+1) % hp->delayBufSmpCnt;
  2369. }
  2370. }
  2371. cmRC_t cmMDelayExec( cmMDelay* p, const cmSample_t* inV, cmSample_t* outV, unsigned sigN, bool bypassFl )
  2372. {
  2373. assert( sigN <= p->outN);
  2374. if( outV == NULL )
  2375. {
  2376. outV = p->outV;
  2377. sigN = cmMin(sigN,p->outN);
  2378. cmVOS_Fill(outV,sigN,0);
  2379. }
  2380. else
  2381. {
  2382. cmVOS_Zero(outV,sigN);
  2383. }
  2384. if( inV == NULL )
  2385. return cmOkRC;
  2386. if( bypassFl )
  2387. {
  2388. memcpy(outV,inV,sigN*sizeof(cmSample_t));
  2389. return cmOkRC;
  2390. }
  2391. unsigned di;
  2392. for( di=0; di<p->delayCnt; ++di)
  2393. {
  2394. cmMDelayHead* hp = p->delayArray + di;
  2395. hp->delaySmpFrac = hp->delayMs * p->srate / 1000.0;
  2396. _cmMDelayExec(p,hp,inV,outV,sigN);
  2397. }
  2398. return cmOkRC;
  2399. }
  2400. void cmMDelaySetTapMs( cmMDelay* p, unsigned tapIdx, cmReal_t ms )
  2401. {
  2402. assert( tapIdx < p->delayCnt );
  2403. p->delayArray[tapIdx].delayMs = ms;
  2404. }
  2405. void cmMDelaySetTapGain(cmMDelay* p, unsigned tapIdx, cmReal_t gain )
  2406. {
  2407. assert( tapIdx < p->delayCnt );
  2408. p->delayArray[tapIdx].delayGain = gain;
  2409. }
  2410. void cmMDelayReport( cmMDelay* p, cmRpt_t* rpt )
  2411. {
  2412. cmRptPrintf(rpt,"tap cnt:%i fb:%f sr:%f\n",p->delayCnt,p->fbCoeff,p->srate);
  2413. }
  2414. //------------------------------------------------------------------------------------------------------------
  2415. cmAudioSegPlayer* cmAudioSegPlayerAlloc( cmCtx* ctx, cmAudioSegPlayer* ap, unsigned procSmpCnt, unsigned outChCnt )
  2416. {
  2417. cmAudioSegPlayer* p = cmObjAlloc( cmAudioSegPlayer, ctx, ap );
  2418. if( procSmpCnt != 0 )
  2419. if( cmAudioSegPlayerInit( p, procSmpCnt, outChCnt ) != cmOkRC )
  2420. cmAudioSegPlayerFree(&p);
  2421. return p;
  2422. }
  2423. cmRC_t cmAudioSegPlayerFree( cmAudioSegPlayer** pp )
  2424. {
  2425. if( pp == NULL || *pp == NULL )
  2426. return cmOkRC;
  2427. cmAudioSegPlayer* p = *pp;
  2428. cmMemPtrFree(&p->segArray);
  2429. cmMemPtrFree(&p->outM);
  2430. cmObjFree(pp);
  2431. return cmOkRC;
  2432. }
  2433. cmRC_t cmAudioSegPlayerInit( cmAudioSegPlayer* p, unsigned procSmpCnt, unsigned outChCnt )
  2434. {
  2435. cmRC_t rc = cmOkRC;
  2436. if((rc = cmAudioSegPlayerFinal(p)) != cmOkRC )
  2437. return rc;
  2438. p->procSmpCnt = procSmpCnt;
  2439. p->outChCnt = outChCnt;
  2440. p->segCnt = 0;
  2441. if( outChCnt )
  2442. {
  2443. unsigned i;
  2444. p->outM = cmMemResizeZ( cmSample_t, p->outM, procSmpCnt * outChCnt );
  2445. p->outChArray = cmMemResizeZ( cmSample_t*, p->outChArray, outChCnt );
  2446. for(i=0; i<outChCnt; ++i)
  2447. p->outChArray[i] = p->outM + (i*procSmpCnt);
  2448. }
  2449. return rc;
  2450. }
  2451. cmRC_t cmAudioSegPlayerFinal( cmAudioSegPlayer* p )
  2452. { return cmOkRC; }
  2453. cmRC_t _cmAudioSegPlayerSegSetup( cmAudioSeg* sp, unsigned id, cmAudioFileBuf* bufPtr, unsigned smpIdx, unsigned smpCnt, unsigned outChIdx )
  2454. {
  2455. sp->bufPtr = bufPtr;
  2456. sp->id = id;
  2457. sp->smpIdx = smpIdx;
  2458. sp->smpCnt = smpCnt;
  2459. sp->outChIdx = outChIdx;
  2460. sp->outSmpIdx = 0;
  2461. sp->flags = 0;
  2462. return cmOkRC;
  2463. }
  2464. cmAudioSeg* _cmAudioSegPlayerIdToSegPtr( cmAudioSegPlayer* p, unsigned id, bool ignoreErrFl )
  2465. {
  2466. unsigned i = 0;
  2467. for(i=0; i<p->segCnt; ++i)
  2468. if( p->segArray[i].id == id )
  2469. return p->segArray + i;
  2470. if( !ignoreErrFl )
  2471. cmCtxRtCondition(&p->obj, cmArgAssertRC,"Unable to locate an audio segment with id=%i.",id);
  2472. return NULL;
  2473. }
  2474. cmRC_t cmAudioSegPlayerInsert( cmAudioSegPlayer* p, unsigned id, cmAudioFileBuf* bufPtr, unsigned smpIdx, unsigned smpCnt, unsigned outChIdx )
  2475. {
  2476. cmRC_t rc;
  2477. assert( _cmAudioSegPlayerIdToSegPtr( p, id, true ) == NULL );
  2478. p->segArray = cmMemResizePZ( cmAudioSeg, p->segArray, p->segCnt + 1 );
  2479. cmAudioSeg* sp = p->segArray + p->segCnt;
  2480. if((rc = _cmAudioSegPlayerSegSetup( sp, id, bufPtr, smpIdx, smpCnt, outChIdx )) == cmOkRC )
  2481. ++p->segCnt;
  2482. return rc;
  2483. }
  2484. cmRC_t cmAudioSegPlayerEdit( cmAudioSegPlayer* p, unsigned id, cmAudioFileBuf* bufPtr, unsigned smpIdx, unsigned smpCnt, unsigned outChIdx )
  2485. {
  2486. cmAudioSeg* sp = _cmAudioSegPlayerIdToSegPtr(p,id,false);
  2487. return _cmAudioSegPlayerSegSetup( sp, id, bufPtr, smpIdx, smpCnt, outChIdx );
  2488. }
  2489. cmRC_t cmAudioSegPlayerRemove( cmAudioSegPlayer* p, unsigned id, bool delFl )
  2490. {
  2491. cmAudioSeg* sp = _cmAudioSegPlayerIdToSegPtr(p,id,false);
  2492. if( sp == NULL )
  2493. return cmArgAssertRC;
  2494. sp->flags = cmEnaFlag( sp->flags, kDelAspFl, delFl );
  2495. return cmOkRC;
  2496. }
  2497. cmRC_t cmAudioSegPlayerEnable( cmAudioSegPlayer* p, unsigned id, bool enableFl, unsigned outSmpIdx )
  2498. {
  2499. cmAudioSeg* sp = _cmAudioSegPlayerIdToSegPtr(p,id,false);
  2500. if( sp == NULL )
  2501. return cmArgAssertRC;
  2502. if( outSmpIdx != cmInvalidIdx )
  2503. sp->outSmpIdx = outSmpIdx;
  2504. sp->flags = cmEnaFlag( sp->flags, kEnableAspFl, enableFl );
  2505. return cmOkRC;
  2506. }
  2507. void _cmAudioSegPlayerResetSeg( cmAudioSeg* sp )
  2508. {
  2509. sp->outSmpIdx = 0;
  2510. sp->flags = cmClrFlag(sp->flags, kEnableAspFl );
  2511. }
  2512. cmRC_t cmAudioSegPlayerReset( cmAudioSegPlayer* p )
  2513. {
  2514. unsigned i;
  2515. for(i=0; i<p->segCnt; ++i)
  2516. {
  2517. cmAudioSeg* sp = p->segArray + i;
  2518. _cmAudioSegPlayerResetSeg(sp);
  2519. }
  2520. return cmOkRC;
  2521. }
  2522. cmRC_t cmAudioSegPlayerExec( cmAudioSegPlayer* p, cmSample_t** outChPtr, unsigned outChCnt, unsigned procSmpCnt )
  2523. {
  2524. unsigned i;
  2525. if( outChPtr == NULL || outChCnt == 0 )
  2526. {
  2527. assert( p->outChCnt > 0 );
  2528. outChPtr = p->outChArray;
  2529. outChCnt = p->outChCnt;
  2530. assert( p->procSmpCnt <= procSmpCnt );
  2531. }
  2532. for(i=0; i<p->segCnt; ++i)
  2533. {
  2534. cmAudioSeg* sp = p->segArray + i;
  2535. // if the output channel is valid and the segment is enabled and not deleted
  2536. if( sp->outChIdx < outChCnt && (sp->flags & (kEnableAspFl | kDelAspFl)) == kEnableAspFl )
  2537. {
  2538. unsigned bufSmpIdx = sp->smpIdx + sp->outSmpIdx;
  2539. unsigned bufSmpCnt = 0;
  2540. // if all the samples have been played
  2541. if( sp->bufPtr->bufN <= bufSmpIdx )
  2542. _cmAudioSegPlayerResetSeg(sp);
  2543. else
  2544. {
  2545. // prevent playing past the end of the buffer
  2546. bufSmpCnt = cmMin( procSmpCnt, sp->bufPtr->bufN - bufSmpIdx );
  2547. // limit the number of samples to the segment length
  2548. bufSmpCnt = cmMin( bufSmpCnt, sp->smpCnt - sp->outSmpIdx );
  2549. // sum the samples into the output channel
  2550. cmVOS_AddVV( outChPtr[ sp->outChIdx ], bufSmpCnt, sp->bufPtr->bufV + bufSmpIdx );
  2551. // incr the next output sample index
  2552. sp->outSmpIdx += bufSmpCnt;
  2553. }
  2554. if( bufSmpCnt < procSmpCnt )
  2555. cmVOS_Zero( outChPtr[ sp->outChIdx ] + bufSmpCnt, procSmpCnt - bufSmpCnt );
  2556. }
  2557. }
  2558. return cmOkRC;
  2559. }
  2560. //------------------------------------------------------------------------------------------------------------
  2561. /*
  2562. cmCluster0* cmCluster0Alloc( cmCtx* ctx, cmCluster0* ap, unsigned stateCnt, unsigned binCnt, unsigned flags, cmCluster0DistFunc_t distFunc, void* dstUserPtr )
  2563. {
  2564. cmCluster0* p = cmObjAlloc( cmCluster0, ctx, ap );
  2565. if( stateCnt != 0 )
  2566. if( cmCluster0Init( p, stateCnt, binCnt, flags, distFunc, distUserPtr ) != cmOkRC )
  2567. cmCluster0Free(&p);
  2568. return p;
  2569. }
  2570. cmRC_t cmCluster0Free( cmCluster0** pp )
  2571. {
  2572. if( pp == NULL || *pp == NULL )
  2573. return cmOkRC;
  2574. cmCluster0* p = *pp;
  2575. cmMemPtrFree(&p->oM);
  2576. cmMemPtrFree(&p->tM);
  2577. cmMemPtrFree(&p->dV);
  2578. cmObjFree(pp);
  2579. return cmOkRC;
  2580. }
  2581. cmRC_t cmCluster0Init( cmCluster0* p, unsigned stateCnt, unsigned binCnt, unsigned flags, cmCluster0DistFunc_t distFunc, void* distUserPtr )
  2582. {
  2583. cmRC_t rc;
  2584. if((rc = cmCluster0Final(p)) != cmOkRC )
  2585. return rc;
  2586. p->oM = cmMemResizeZ( cmReal_t, p->oM, binCnt * stateCnt );
  2587. p->tM = cmMemResizeZ( cmReal_t, p->tM, stateCnt * stateCnt );
  2588. p->stateCnt = stateCnt;
  2589. p->binCnt = binCnt;
  2590. p->flags = flags;
  2591. p->distFunc = distFunc;
  2592. p->distUserPtr = distUserPtr;
  2593. p->cnt = 0;
  2594. }
  2595. cmRC_t cmCluster0Final( cmCluster0* p )
  2596. { return cmOkRC; }
  2597. cmRC_t cmCluster0Exec( cmCluster0* p, const cmReal_t* v, unsigned vn )
  2598. {
  2599. assert( vn <= p->binCnt );
  2600. ++cnt;
  2601. if( cnt <= stateCnt )
  2602. {
  2603. cmVOR_Copy( p->oM + ((cnt-1)*binCnt), vn, v );
  2604. return cmOkRC;
  2605. }
  2606. return cmOkRC;
  2607. }
  2608. */
  2609. cmNmf_t* cmNmfAlloc( cmCtx* ctx, cmNmf_t* ap, unsigned n, unsigned m, unsigned r, unsigned maxIterCnt, unsigned convergeCnt )
  2610. {
  2611. cmNmf_t* p = cmObjAlloc( cmNmf_t, ctx, ap );
  2612. if( n != 0 )
  2613. if( cmNmfInit( p, n, m, r, maxIterCnt, convergeCnt ) != cmOkRC )
  2614. cmNmfFree(&p);
  2615. return p;
  2616. }
  2617. cmRC_t cmNmfFree( cmNmf_t** pp )
  2618. {
  2619. if( pp== NULL || *pp == NULL )
  2620. return cmOkRC;
  2621. cmNmf_t* p = *pp;
  2622. cmMemPtrFree(&p->V);
  2623. cmMemPtrFree(&p->W);
  2624. cmMemPtrFree(&p->H);
  2625. cmMemPtrFree(&p->tr);
  2626. cmMemPtrFree(&p->x);
  2627. cmMemPtrFree(&p->t0nm);
  2628. cmMemPtrFree(&p->t1nm);
  2629. cmMemPtrFree(&p->Wt);
  2630. cmMemPtrFree(&p->trm);
  2631. cmMemPtrFree(&p->crm);
  2632. cmMemPtrFree(&p->c0);
  2633. cmMemPtrFree(&p->c1);
  2634. cmMemPtrFree(&p->idxV);
  2635. cmObjFree(pp);
  2636. return cmOkRC;
  2637. }
  2638. cmRC_t cmNmfInit( cmNmf_t* p, unsigned n, unsigned m, unsigned r, unsigned maxIterCnt, unsigned convergeCnt )
  2639. {
  2640. cmRC_t rc;
  2641. if((rc = cmNmfFinal(p)) != cmOkRC )
  2642. return rc;
  2643. p->n = n;
  2644. p->m = m;
  2645. p->r = r;
  2646. p->maxIterCnt = maxIterCnt;
  2647. p->convergeCnt= convergeCnt;
  2648. p->V = cmMemResizeZ(cmReal_t, p->V, n*m );
  2649. p->W = cmMemResize( cmReal_t, p->W, n*r );
  2650. p->H = cmMemResize( cmReal_t, p->H, r*m );
  2651. p->tr = cmMemResize( cmReal_t, p->tr, r );
  2652. p->x = cmMemResize( cmReal_t, p->x, r*cmMax(m,n) );
  2653. p->t0nm = cmMemResize( cmReal_t, p->t0nm, cmMax(r,n)*m );
  2654. p->Ht = p->t0nm;
  2655. p->t1nm = cmMemResize( cmReal_t, p->t1nm, n*m );
  2656. p->Wt = cmMemResize( cmReal_t, p->Wt, r*n );
  2657. p->trm = cmMemResize( cmReal_t, p->trm, r*cmMax(m,n) );
  2658. p->crm = cmMemResizeZ(unsigned, p->crm, r*m);
  2659. p->tnr = p->trm;
  2660. p->c0 = cmMemResizeZ(unsigned, p->c0, m*m);
  2661. p->c1 = cmMemResizeZ(unsigned, p->c1, m*m);
  2662. p->idxV = cmMemResizeZ(unsigned, p->idxV, m );
  2663. p->c0m = p->c0;
  2664. p->c1m = p->c1;
  2665. cmVOR_Random(p->W,n*r,0.0,1.0);
  2666. cmVOR_Random(p->H,r*m,0.0,1.0);
  2667. return rc;
  2668. }
  2669. cmRC_t cmNmfFinal(cmNmf_t* p )
  2670. { return cmOkRC; }
  2671. // NMF base on: Lee and Seung, 2001, Algo's for Non-negative Matrix Fcmtorization
  2672. // Connectivity stopping technique based on: http://www.broadinstitute.org/mpr/publications/projects/NMF/nmf.m
  2673. cmRC_t cmNmfExec( cmNmf_t* p, const cmReal_t* vM, unsigned cn )
  2674. {
  2675. cmRC_t rc = cmOkRC;
  2676. unsigned i,j,k;
  2677. unsigned n = p->n;
  2678. unsigned m = p->m;
  2679. unsigned r = p->r;
  2680. unsigned stopIter = 0;
  2681. assert(cn <= m );
  2682. // shift in the incoming columns of V[]
  2683. if( cn < m )
  2684. cmVOR_Shift(p->V, n*m, n*cn,0);
  2685. cmVOR_Copy( p->V, n*cn,vM );
  2686. // shift H[] by the same amount as V[]
  2687. if( cn < m )
  2688. cmVOR_Shift( p->H, r*m, r*cn,0);
  2689. cmVOR_Random(p->H, r*cn, 0.0, 1.0 );
  2690. cmVOU_Zero( p->c1m, m*m );
  2691. for(i=0,j=0; i<p->maxIterCnt && stopIter<p->convergeCnt; ++i)
  2692. {
  2693. // x[r,m] =repmat(sum(W,1)',1,m);
  2694. cmVOR_SumM( p->W, n, r, p->tr );
  2695. for(j=0; j<m; ++j)
  2696. cmVOR_Copy( p->x + (j*r), r, p->tr );
  2697. cmVOR_Transpose(p->Wt,p->W,n,r);
  2698. //H=H.*(W'*(V./(W*H)))./x;
  2699. cmVOR_MultMMM(p->t0nm,n,m,p->W,p->H,r); // t0nm[n,m] = W*H
  2700. cmVOR_DivVVV( p->t1nm,n*m,p->V,p->t0nm); // t1nm[n,m] = V./(W*H)
  2701. cmVOR_MultMMM(p->trm,r,m,p->Wt,p->t1nm,n); // trm[r,m] = W'*(V./(W*H))
  2702. cmVOR_MultVV(p->H,r*m,p->trm); // H[r,m] = H .* (W'*(V./(W*H)))
  2703. cmVOR_DivVV(p->H,r*m, p->x ); // H[r,m] = (H .* (W'*(V./(W*H)))) ./ x
  2704. // x[n,r]=repmat(sum(H,2)',n,1);
  2705. cmVOR_SumMN(p->H, r, m, p->tr );
  2706. for(j=0; j<n; ++j)
  2707. cmVOR_CopyN(p->x + j, r, n, p->tr, 1 );
  2708. cmVOR_Transpose(p->Ht,p->H,r,m);
  2709. // W=W.*((V./(W*H))*H')./x;
  2710. cmVOR_MultMMM(p->tnr,n,r,p->t1nm,p->Ht,m); // tnr[n,r] = (V./(W*H))*Ht
  2711. cmVOR_MultVV(p->W,n*r,p->tnr); // W[n,r] = W.*(V./(W*H))*Ht
  2712. cmVOR_DivVV(p->W,n*r,p->x); // W[n,r] = W.*(V./(W*H))*Ht ./x
  2713. if( i % 10 == 0 )
  2714. {
  2715. cmVOR_ReplaceLte( p->H, r*m, p->H, 2.2204e-16, 2.2204e-16 );
  2716. cmVOR_ReplaceLte( p->W, n*r, p->W, 2.2204e-16, 2.2204e-16 );
  2717. cmVOR_MaxIndexM( p->idxV, p->H, r, m );
  2718. unsigned mismatchCnt = 0;
  2719. for(j=0; j<m; ++j)
  2720. for(k=0; k<m; ++k)
  2721. {
  2722. unsigned c_idx = (j*m)+k;
  2723. p->c0m[ c_idx ] = p->idxV[j] == p->idxV[k];
  2724. mismatchCnt += p->c0m[ c_idx ] != p->c1m[ c_idx ];
  2725. }
  2726. if( mismatchCnt == 0 )
  2727. ++stopIter;
  2728. else
  2729. stopIter = 0;
  2730. printf("%i %i %i\n",i,stopIter,mismatchCnt);
  2731. fflush(stdout);
  2732. unsigned* tcm = p->c0m;
  2733. p->c0m = p->c1m;
  2734. p->c1m = tcm;
  2735. }
  2736. }
  2737. return rc;
  2738. }
  2739. //------------------------------------------------------------------------------------------------------------
  2740. cmRC_t _cmVectArrayAppend( cmVectArray_t* p, const void* v, unsigned typeByteCnt, unsigned valCnt )
  2741. {
  2742. cmRC_t rc = cmSubSysFailRC;
  2743. cmVectArrayVect_t* ep = NULL;
  2744. unsigned byteCnt = typeByteCnt * valCnt;
  2745. if( byteCnt == 0 || v == NULL )
  2746. return rc;
  2747. // verify that all vectors written to this vector array contain the same data type.
  2748. if( (cmIsFlag(p->flags,kFloatVaFl) && typeByteCnt!=sizeof(float)) || (cmIsFlag(p->flags,kDoubleVaFl) && typeByteCnt!=sizeof(double)))
  2749. return cmCtxRtCondition(&p->obj,cmInvalidArgRC,"All data stored to a cmVectArray_t must be a consistent type.");
  2750. // allocate space for the link record
  2751. if((ep = cmMemAllocZ(cmVectArrayVect_t,1)) == NULL )
  2752. goto errLabel;
  2753. // allocate space for the vector data
  2754. if((ep->u.v = cmMemAlloc(char,typeByteCnt*valCnt)) == NULL )
  2755. goto errLabel;
  2756. // append the link recd to the end of the element list
  2757. if( p->ep != NULL )
  2758. p->ep->link = ep;
  2759. else
  2760. {
  2761. p->bp = ep;
  2762. p->cur = p->bp;
  2763. }
  2764. p->ep = ep;
  2765. // store the length of the vector
  2766. ep->n = valCnt;
  2767. // copy in the vector data
  2768. memcpy(ep->u.v,v,byteCnt);
  2769. // track the number of vectors stored
  2770. p->vectCnt += 1;
  2771. // track the longest data vector
  2772. if( valCnt > p->maxEleCnt )
  2773. p->maxEleCnt = valCnt;
  2774. rc = cmOkRC;
  2775. errLabel:
  2776. if(rc != cmOkRC )
  2777. {
  2778. cmMemFree(ep->u.v);
  2779. cmMemFree(ep);
  2780. }
  2781. return rc;
  2782. }
  2783. cmVectArray_t* cmVectArrayAlloc( cmCtx* ctx, unsigned flags )
  2784. {
  2785. cmRC_t rc = cmOkRC;
  2786. cmVectArray_t* p = cmObjAlloc(cmVectArray_t,ctx,NULL);
  2787. assert(p != NULL);
  2788. switch( flags & (kFloatVaFl | kDoubleVaFl | kSampleVaFl | kRealVaFl ) )
  2789. {
  2790. case kFloatVaFl:
  2791. p->flags |= kFloatVaFl;
  2792. p->typeByteCnt = sizeof(float);
  2793. break;
  2794. case kDoubleVaFl:
  2795. p->flags |= kDoubleVaFl;
  2796. p->typeByteCnt = sizeof(double);
  2797. break;
  2798. case kSampleVaFl:
  2799. if( sizeof(cmSample_t) == sizeof(float) )
  2800. p->flags |= kFloatVaFl;
  2801. else
  2802. {
  2803. if( sizeof(cmSample_t) == sizeof(double))
  2804. p->flags |= kDoubleVaFl;
  2805. else
  2806. rc = cmCtxRtCondition(&p->obj,cmInvalidArgRC,"The size of the cmSample_t is not consistent with float or double.");
  2807. }
  2808. p->typeByteCnt = sizeof(cmSample_t);
  2809. break;
  2810. case kRealVaFl:
  2811. if( sizeof(cmReal_t) == sizeof(float) )
  2812. p->flags |= kFloatVaFl;
  2813. else
  2814. {
  2815. if( sizeof(cmReal_t) == sizeof(double))
  2816. p->flags |= kDoubleVaFl;
  2817. else
  2818. rc = cmCtxRtCondition(&p->obj,cmInvalidArgRC,"The size of the cmReal_t is not consistent with float or double.");
  2819. }
  2820. p->typeByteCnt = sizeof(cmReal_t);
  2821. break;
  2822. default:
  2823. rc = cmCtxRtCondition(&p->obj,cmInvalidArgRC,"The vector array value type flag was not recognized.");
  2824. }
  2825. if(rc != cmOkRC)
  2826. cmVectArrayFree(&p);
  2827. return p;
  2828. }
  2829. cmVectArray_t* cmVectArrayAllocFromFile(cmCtx* ctx, const char* fn )
  2830. {
  2831. cmRC_t rc = cmOkRC;
  2832. FILE* fp = NULL;
  2833. char* buf = NULL;
  2834. cmVectArray_t* p = NULL;
  2835. unsigned hn = 4;
  2836. unsigned hdr[hn];
  2837. // create the file
  2838. if((fp = fopen(fn,"rb")) == NULL )
  2839. {
  2840. rc = cmCtxRtCondition(&ctx->obj,cmSystemErrorRC,"The vector array file '%s' could not be opened.",cmStringNullGuard(fn));
  2841. goto errLabel;
  2842. }
  2843. if( fread(hdr,sizeof(unsigned),hn,fp) != hn )
  2844. {
  2845. rc = cmCtxRtCondition(&ctx->obj,cmSystemErrorRC,"The vector array file header could not be read from '%s'.",cmStringNullGuard(fn));
  2846. goto errLabel;
  2847. }
  2848. unsigned flags = hdr[0];
  2849. unsigned typeByteCnt = hdr[1];
  2850. unsigned vectCnt = hdr[2];
  2851. unsigned maxEleCnt = hdr[3];
  2852. unsigned i;
  2853. buf = cmMemAlloc(char,maxEleCnt*typeByteCnt);
  2854. if((p = cmVectArrayAlloc(ctx, flags )) == NULL )
  2855. goto errLabel;
  2856. for(i=0; i<vectCnt; ++i)
  2857. {
  2858. unsigned vn;
  2859. if( fread(&vn,sizeof(unsigned),1,fp) != 1 )
  2860. {
  2861. rc = cmCtxRtCondition(&p->obj,cmSystemErrorRC,"The vector array file element count read failed on vector index:%i in '%s'.",i,cmStringNullGuard(fn));
  2862. goto errLabel;
  2863. }
  2864. assert( vn <= maxEleCnt );
  2865. if( fread(buf,typeByteCnt,vn,fp) != vn )
  2866. {
  2867. rc = cmCtxRtCondition(&p->obj,cmSystemErrorRC,"The vector array data read failed on vector index:%i in '%s'.",i,cmStringNullGuard(fn));
  2868. goto errLabel;
  2869. }
  2870. if((rc = _cmVectArrayAppend(p,buf, typeByteCnt, vn )) != cmOkRC )
  2871. {
  2872. rc = cmCtxRtCondition(&p->obj,rc,"The vector array data store failed on vector index:%i in '%s'.",i,cmStringNullGuard(fn));
  2873. goto errLabel;
  2874. }
  2875. }
  2876. errLabel:
  2877. if( fp != NULL )
  2878. fclose(fp);
  2879. cmMemFree(buf);
  2880. if(rc != cmOkRC && p != NULL)
  2881. cmVectArrayFree(&p);
  2882. return p;
  2883. }
  2884. cmRC_t cmVectArrayFree( cmVectArray_t** pp )
  2885. {
  2886. cmRC_t rc = cmOkRC;
  2887. if( pp == NULL || *pp == NULL )
  2888. return rc;
  2889. cmVectArray_t* p = *pp;
  2890. if((rc = cmVectArrayClear(p)) != cmOkRC )
  2891. return rc;
  2892. cmMemFree(p->tempV);
  2893. cmObjFree(pp);
  2894. return rc;
  2895. }
  2896. cmRC_t cmVectArrayClear( cmVectArray_t* p )
  2897. {
  2898. cmVectArrayVect_t* ep = p->bp;
  2899. while( ep!=NULL )
  2900. {
  2901. cmVectArrayVect_t* np = ep->link;
  2902. cmMemFree(ep->u.v);
  2903. cmMemFree(ep);
  2904. ep = np;
  2905. }
  2906. p->bp = NULL;
  2907. p->ep = NULL;
  2908. p->maxEleCnt = 0;
  2909. p->vectCnt = 0;
  2910. return cmOkRC;
  2911. }
  2912. unsigned cmVectArrayCount( const cmVectArray_t* p )
  2913. { return p->vectCnt; }
  2914. /*
  2915. unsigned cmVectArrayEleCount( const cmVectArray_t* p )
  2916. {
  2917. cmVectArrayVect_t* ep = p->bp;
  2918. unsigned n = 0;
  2919. for(; ep!=NULL; ep=ep->link )
  2920. n += ep->n;
  2921. return n;
  2922. }
  2923. */
  2924. cmRC_t cmVectArrayAppendS( cmVectArray_t* p, const cmSample_t* v, unsigned vn )
  2925. { return _cmVectArrayAppend(p,v,sizeof(v[0]),vn); }
  2926. cmRC_t cmVectArrayAppendR( cmVectArray_t* p, const cmReal_t* v, unsigned vn )
  2927. { return _cmVectArrayAppend(p,v,sizeof(v[0]),vn); }
  2928. cmRC_t cmVectArrayAppendF( cmVectArray_t* p, const float* v, unsigned vn )
  2929. { return _cmVectArrayAppend(p,v,sizeof(v[0]),vn); }
  2930. cmRC_t cmVectArrayAppendD( cmVectArray_t* p, const double* v, unsigned vn )
  2931. { return _cmVectArrayAppend(p,v,sizeof(v[0]),vn); }
  2932. cmRC_t cmVectArrayAppendI( cmVectArray_t* p, const int* v, unsigned vn )
  2933. {
  2934. p->tempV = cmMemResize(double,p->tempV,vn);
  2935. unsigned i;
  2936. for(i=0; i<vn; ++i)
  2937. p->tempV[i] = v[i];
  2938. cmRC_t rc = cmVectArrayAppendD(p,p->tempV,vn);
  2939. return rc;
  2940. }
  2941. cmRC_t cmVectArrayAppendU( cmVectArray_t* p, const unsigned* v, unsigned vn )
  2942. {
  2943. p->tempV = cmMemResize(double,p->tempV,vn);
  2944. unsigned i;
  2945. for(i=0; i<vn; ++i)
  2946. p->tempV[i] = v[i];
  2947. cmRC_t rc = cmVectArrayAppendD(p,p->tempV,vn);
  2948. return rc;
  2949. }
  2950. cmRC_t cmVectArrayWrite( cmVectArray_t* p, const char* fn )
  2951. {
  2952. cmRC_t rc = cmOkRC;
  2953. FILE* fp = NULL;
  2954. cmVectArrayVect_t* ep;
  2955. unsigned i;
  2956. unsigned hn = 4;
  2957. unsigned hdr[hn];
  2958. hdr[0] = p->flags;
  2959. hdr[1] = p->typeByteCnt;
  2960. hdr[2] = p->vectCnt;
  2961. hdr[3] = p->maxEleCnt;
  2962. // create the file
  2963. if((fp = fopen(fn,"wb")) == NULL )
  2964. return cmCtxRtCondition(&p->obj,cmSystemErrorRC,"The vector array file '%s' could not be created.",cmStringNullGuard(fn));
  2965. // write the header
  2966. if( fwrite(hdr,sizeof(unsigned),hn,fp) != hn )
  2967. {
  2968. rc = cmCtxRtCondition(&p->obj,cmSystemErrorRC,"Vector array file header write failed in '%s'.",cmStringNullGuard(fn));
  2969. goto errLabel;
  2970. }
  2971. // write each vector element
  2972. for(ep=p->bp,i=0; ep!=NULL; ep=ep->link,++i)
  2973. {
  2974. // write the count of data values in the vector
  2975. if( fwrite(&ep->n,sizeof(ep->n),1,fp) != 1 )
  2976. {
  2977. rc = cmCtxRtCondition(&p->obj,cmSystemErrorRC,"Vector array file write failed on element header %i in '%s'.",i,cmStringNullGuard(fn));
  2978. goto errLabel;
  2979. }
  2980. // write the vector
  2981. if(fwrite(ep->u.v,p->typeByteCnt,ep->n,fp) != ep->n )
  2982. {
  2983. rc = cmCtxRtCondition(&p->obj,cmSystemErrorRC,"Vector array file write failed on data vector %i in '%s'.",i,cmStringNullGuard(fn));
  2984. goto errLabel;
  2985. }
  2986. }
  2987. errLabel:
  2988. if( fp != NULL )
  2989. fclose(fp);
  2990. return rc;
  2991. }
  2992. unsigned cmVectArrayForEachS( cmVectArray_t* p, unsigned idx, unsigned cnt, cmVectArrayForEachFuncS_t func, void* arg )
  2993. {
  2994. cmVectArrayVect_t* ep = p->bp;
  2995. unsigned i = 0;
  2996. unsigned n = 0;
  2997. // for each sub-array
  2998. for(; ep!=NULL && n<cnt; ep=ep->link )
  2999. {
  3000. // if the cur sub-array is in the range of idx:idx+cnt
  3001. if( i <= idx && idx < i + ep->n )
  3002. {
  3003. unsigned j = idx - i; // starting idx into cur sub-array
  3004. assert(j<ep->n);
  3005. unsigned m = cmMin(ep->n - j,cnt-n); // cnt of ele's to send from cur sub-array
  3006. // do callback
  3007. if( func(arg, idx, ep->u.sV + j, m ) != cmOkRC )
  3008. break;
  3009. idx += m;
  3010. n += m;
  3011. }
  3012. i += ep->n;
  3013. }
  3014. return n;
  3015. }
  3016. cmRC_t cmVectArrayWriteVectorS( cmCtx* ctx, const char* fn, const cmSample_t* v, unsigned vn )
  3017. {
  3018. cmRC_t rc = cmOkRC;
  3019. cmVectArray_t* p;
  3020. if((p = cmVectArrayAlloc(ctx,kSampleVaFl)) == NULL )
  3021. return cmCtxRtCondition(&ctx->obj,cmSubSysFailRC,"Unable to allocate an cmVectArray_t in cmVectArrayWriteVectorS().");
  3022. if((rc = cmVectArrayAppendS(p,v,vn)) != cmOkRC )
  3023. {
  3024. rc = cmCtxRtCondition(&p->obj,rc,"Vector append failed in cmVectArrayWriteVectorS().");
  3025. goto errLabel;
  3026. }
  3027. if((rc = cmVectArrayWrite(p,fn)) != cmOkRC )
  3028. rc = cmCtxRtCondition(&p->obj,rc,"Vector array write failed in cmVectArrayWriteVectorS().");
  3029. errLabel:
  3030. if((rc = cmVectArrayFree(&p)) != cmOkRC )
  3031. rc = cmCtxRtCondition(&ctx->obj,rc,"Free failed on cmVectArrayFree() in cmVectArrayWriteVectorS().");
  3032. return rc;
  3033. }
  3034. cmRC_t cmVectArrayWriteMatrixS( cmCtx* ctx, const char* fn, const cmSample_t* m, unsigned rn, unsigned cn )
  3035. {
  3036. cmRC_t rc = cmOkRC;
  3037. cmVectArray_t* p;
  3038. unsigned i;
  3039. if((p = cmVectArrayAlloc(ctx,kSampleVaFl)) == NULL )
  3040. return cmCtxRtCondition(&ctx->obj,cmSubSysFailRC,"Unable to allocate an cmVectArray_t in cmVectArrayWriteVectorS().");
  3041. for(i=0; i<cn; ++i)
  3042. {
  3043. if((rc = cmVectArrayAppendS(p,m + (i*rn), rn)) != cmOkRC )
  3044. {
  3045. rc = cmCtxRtCondition(&p->obj,rc,"Vector append failed in cmVectArrayWriteVectorS().");
  3046. goto errLabel;
  3047. }
  3048. }
  3049. if((rc = cmVectArrayWrite(p,fn)) != cmOkRC )
  3050. rc = cmCtxRtCondition(&p->obj,rc,"Vector array write failed in cmVectArrayWriteVectorS().");
  3051. errLabel:
  3052. if((rc = cmVectArrayFree(&p)) != cmOkRC )
  3053. rc = cmCtxRtCondition(&ctx->obj,rc,"Free failed on cmVectArrayFree() in cmVectArrayWriteVectorS().");
  3054. return rc;
  3055. }
  3056. cmRC_t cmVectArrayWriteMatrixR( cmCtx* ctx, const char* fn, const cmReal_t* m, unsigned rn, unsigned cn )
  3057. {
  3058. cmRC_t rc = cmOkRC;
  3059. cmVectArray_t* p;
  3060. unsigned i;
  3061. if((p = cmVectArrayAlloc(ctx,kRealVaFl)) == NULL )
  3062. return cmCtxRtCondition(&ctx->obj,cmSubSysFailRC,"Unable to allocate an cmVectArray_t in cmVectArrayWriteVectorR().");
  3063. for(i=0; i<cn; ++i)
  3064. {
  3065. if((rc = cmVectArrayAppendR(p,m + (i*rn), rn)) != cmOkRC )
  3066. {
  3067. rc = cmCtxRtCondition(&p->obj,rc,"Vector append failed in cmVectArrayWriteVectorR().");
  3068. goto errLabel;
  3069. }
  3070. }
  3071. if((rc = cmVectArrayWrite(p,fn)) != cmOkRC )
  3072. rc = cmCtxRtCondition(&p->obj,rc,"Vector array write failed in cmVectArrayWriteVectorR().");
  3073. errLabel:
  3074. if((rc = cmVectArrayFree(&p)) != cmOkRC )
  3075. rc = cmCtxRtCondition(&ctx->obj,rc,"Free failed on cmVectArrayFree() in cmVectArrayWriteVectorR().");
  3076. return rc;
  3077. }
  3078. cmRC_t cmVectArrayWriteVectorI( cmCtx* ctx, const char* fn, const int* v, unsigned vn )
  3079. {
  3080. cmRC_t rc = cmOkRC;
  3081. cmVectArray_t* p;
  3082. if((p = cmVectArrayAlloc(ctx,kIntVaFl)) == NULL )
  3083. return cmCtxRtCondition(&ctx->obj,cmSubSysFailRC,"Unable to allocate an cmVectArray_t in cmVectArrayWriteVectorS().");
  3084. if((rc = cmVectArrayAppendI(p,v,vn)) != cmOkRC )
  3085. {
  3086. rc = cmCtxRtCondition(&p->obj,rc,"Vector append failed in cmVectArrayWriteVectorS().");
  3087. goto errLabel;
  3088. }
  3089. if((rc = cmVectArrayWrite(p,fn)) != cmOkRC )
  3090. rc = cmCtxRtCondition(&p->obj,rc,"Vector array write failed in cmVectArrayWriteVectorS().");
  3091. errLabel:
  3092. if((rc = cmVectArrayFree(&p)) != cmOkRC )
  3093. rc = cmCtxRtCondition(&ctx->obj,rc,"Free failed on cmVectArrayFree() in cmVectArrayWriteVectorS().");
  3094. return rc;
  3095. }
  3096. cmRC_t cmVectArrayWriteMatrixI( cmCtx* ctx, const char* fn, const int* m, unsigned rn, unsigned cn )
  3097. {
  3098. cmRC_t rc = cmOkRC;
  3099. cmVectArray_t* p;
  3100. unsigned i;
  3101. if((p = cmVectArrayAlloc(ctx,kIntVaFl)) == NULL )
  3102. return cmCtxRtCondition(&ctx->obj,cmSubSysFailRC,"Unable to allocate an cmVectArray_t in cmVectArrayWriteVectorI().");
  3103. for(i=0; i<cn; ++i)
  3104. {
  3105. if((rc = cmVectArrayAppendI(p,m + (i*rn), rn)) != cmOkRC )
  3106. {
  3107. rc = cmCtxRtCondition(&p->obj,rc,"Vector append failed in cmVectArrayWriteVectorI().");
  3108. goto errLabel;
  3109. }
  3110. }
  3111. if((rc = cmVectArrayWrite(p,fn)) != cmOkRC )
  3112. rc = cmCtxRtCondition(&p->obj,rc,"Vector array write failed in cmVectArrayWriteVectorI().");
  3113. errLabel:
  3114. if((rc = cmVectArrayFree(&p)) != cmOkRC )
  3115. rc = cmCtxRtCondition(&ctx->obj,rc,"Free failed on cmVectArrayFree() in cmVectArrayWriteVectorI().");
  3116. return rc;
  3117. }
  3118. cmRC_t cmVectArrayForEachTextFuncS( void* arg, unsigned idx, const cmSample_t* xV, unsigned xN )
  3119. {
  3120. assert(0);
  3121. return cmOkRC;
  3122. }
  3123. cmRC_t cmVectArrayRewind( cmVectArray_t* p )
  3124. {
  3125. p->cur = p->bp;
  3126. return cmOkRC;
  3127. }
  3128. cmRC_t cmVectArrayAdvance( cmVectArray_t* p, unsigned n )
  3129. {
  3130. unsigned i;
  3131. for(i=0; i<n; ++i)
  3132. {
  3133. if( p->cur == NULL )
  3134. break;
  3135. p->cur = p->cur->link;
  3136. }
  3137. return cmOkRC;
  3138. }
  3139. bool cmVectArrayIsEOL( const cmVectArray_t* p )
  3140. { return p->cur == NULL; }
  3141. unsigned cmVectArrayEleCount( const cmVectArray_t* p )
  3142. {
  3143. if( p->cur == NULL )
  3144. return 0;
  3145. return p->cur->n;
  3146. }
  3147. cmRC_t cmVectArrayGetF( cmVectArray_t* p, float* v, unsigned* vnRef )
  3148. {
  3149. assert( cmIsFlag(p->flags,kFloatVaFl) );
  3150. if( cmVectArrayIsEOL(p) )
  3151. return cmCtxRtCondition(&p->obj,cmSubSysFailRC,"VectArray get failed because the state is EOL.");
  3152. unsigned n = cmMin(*vnRef,p->cur->n);
  3153. cmVOF_Copy(v,n,p->cur->u.fV);
  3154. *vnRef = n;
  3155. return cmOkRC;
  3156. }
  3157. cmRC_t cmVectArrayGetD( cmVectArray_t* p, double* v, unsigned* vnRef )
  3158. {
  3159. assert( cmIsFlag(p->flags,kDoubleVaFl) );
  3160. if( cmVectArrayIsEOL(p) )
  3161. return cmCtxRtCondition(&p->obj,cmSubSysFailRC,"VectArray get failed because the state is EOL.");
  3162. unsigned n = cmMin(*vnRef,p->cur->n);
  3163. cmVOD_Copy(v,n,p->cur->u.dV);
  3164. *vnRef = n;
  3165. return cmOkRC;
  3166. }
  3167. cmRC_t cmVectArrayGetI( cmVectArray_t* p, int* v, unsigned* vnRef )
  3168. {
  3169. if( cmVectArrayIsEOL(p) )
  3170. return cmCtxRtCondition(&p->obj,cmSubSysFailRC,"VectArray get failed because the state is EOL.");
  3171. unsigned n = cmMin(*vnRef,p->cur->n);
  3172. unsigned i;
  3173. if( cmIsFlag(p->flags,kDoubleVaFl) )
  3174. {
  3175. for(i=0; i<n; ++i)
  3176. v[i] = (int)(p->cur->u.dV[i]);
  3177. }
  3178. else
  3179. {
  3180. if( cmIsFlag(p->flags,kFloatVaFl) )
  3181. {
  3182. for(i=0; i<n; ++i)
  3183. v[i] = (int)(p->cur->u.fV[i]);
  3184. }
  3185. else
  3186. {
  3187. assert(0);
  3188. }
  3189. }
  3190. *vnRef = n;
  3191. return cmOkRC;
  3192. }
  3193. cmRC_t cmVectArrayGetU( cmVectArray_t* p, unsigned* v, unsigned* vnRef )
  3194. {
  3195. if( cmVectArrayIsEOL(p) )
  3196. return cmCtxRtCondition(&p->obj,cmSubSysFailRC,"VectArray get failed because the state is EOL.");
  3197. unsigned n = cmMin(*vnRef,p->cur->n);
  3198. unsigned i;
  3199. if( cmIsFlag(p->flags,kDoubleVaFl) )
  3200. {
  3201. for(i=0; i<n; ++i)
  3202. v[i] = (unsigned)(p->cur->u.dV[i]);
  3203. }
  3204. else
  3205. {
  3206. if( cmIsFlag(p->flags,kFloatVaFl) )
  3207. {
  3208. for(i=0; i<n; ++i)
  3209. v[i] = (unsigned)(p->cur->u.fV[i]);
  3210. }
  3211. else
  3212. {
  3213. assert(0);
  3214. }
  3215. }
  3216. *vnRef = n;
  3217. return cmOkRC;
  3218. }
  3219. unsigned cmVectArrayVectEleCount( cmVectArray_t* p, unsigned groupIdx, unsigned groupCnt )
  3220. {
  3221. unsigned n = 0;
  3222. cmVectArrayVect_t* pos = p->cur;
  3223. if( cmVectArrayRewind(p) != cmOkRC )
  3224. goto errLabel;
  3225. if( cmVectArrayAdvance(p,groupIdx) != cmOkRC )
  3226. goto errLabel;
  3227. while( !cmVectArrayIsEOL(p) )
  3228. {
  3229. n += cmVectArrayEleCount(p);
  3230. if(cmVectArrayAdvance(p,groupCnt) != cmOkRC )
  3231. goto errLabel;
  3232. }
  3233. errLabel:
  3234. p->cur = pos;
  3235. return n;
  3236. }
  3237. cmRC_t cmVectArrayFormVectF( cmVectArray_t* p, unsigned groupIdx, unsigned groupCnt, float** vRef, unsigned* vnRef )
  3238. {
  3239. cmRC_t rc = cmOkRC;
  3240. *vRef = NULL;
  3241. *vnRef = 0;
  3242. unsigned N = cmVectArrayVectEleCount(p,groupIdx,groupCnt);
  3243. if( N == 0 )
  3244. return rc;
  3245. float* v = cmMemAllocZ(float,N);
  3246. unsigned i = 0;
  3247. cmVectArrayVect_t* pos = p->cur;
  3248. if( cmVectArrayRewind(p) != cmOkRC )
  3249. goto errLabel;
  3250. if( cmVectArrayAdvance(p,groupIdx) != cmOkRC )
  3251. goto errLabel;
  3252. while( !cmVectArrayIsEOL(p) )
  3253. {
  3254. unsigned n = cmVectArrayEleCount(p);
  3255. assert(i+n <= N);
  3256. cmVectArrayGetF(p,v+i,&n);
  3257. i += n;
  3258. cmVectArrayAdvance(p,groupCnt);
  3259. }
  3260. *vRef = v;
  3261. *vnRef = i;
  3262. errLabel:
  3263. p->cur = pos;
  3264. return rc;
  3265. }
  3266. cmRC_t cmVectArrayFormVectColF( cmVectArray_t* p, unsigned groupIdx, unsigned groupCnt, unsigned colIdx, float** vRef, unsigned* vnRef )
  3267. {
  3268. cmRC_t rc = cmOkRC;
  3269. *vRef = NULL;
  3270. *vnRef = 0;
  3271. // assume there will be one output element for each group
  3272. unsigned N = cmVectArrayCount(p)/groupCnt + 1;
  3273. if( N == 0 )
  3274. return rc;
  3275. float* v = cmMemAllocZ(float,N);
  3276. unsigned i = 0;
  3277. cmVectArrayVect_t* pos = p->cur;
  3278. if( cmVectArrayRewind(p) != cmOkRC )
  3279. goto errLabel;
  3280. if( cmVectArrayAdvance(p,groupIdx) != cmOkRC )
  3281. goto errLabel;
  3282. while( i<N && !cmVectArrayIsEOL(p) )
  3283. {
  3284. unsigned tn = cmVectArrayEleCount(p);
  3285. float tv[tn];
  3286. // read the sub-vector
  3287. cmVectArrayGetF(p,tv,&tn);
  3288. // store the output value
  3289. if( colIdx < tn )
  3290. {
  3291. v[i] = tv[colIdx];
  3292. i += 1;
  3293. }
  3294. cmVectArrayAdvance(p,groupCnt);
  3295. }
  3296. *vRef = v;
  3297. *vnRef = i;
  3298. errLabel:
  3299. p->cur = pos;
  3300. return rc;
  3301. }
  3302. cmRC_t cmVectArrayFormVectColU( cmVectArray_t* p, unsigned groupIdx, unsigned groupCnt, unsigned colIdx, unsigned** vRef, unsigned* vnRef )
  3303. {
  3304. cmRC_t rc = cmOkRC;
  3305. *vRef = NULL;
  3306. *vnRef = 0;
  3307. // assume there will be one output element for each group
  3308. unsigned N = cmVectArrayCount(p)/groupCnt + 1;
  3309. if( N == 0 )
  3310. return rc;
  3311. unsigned* v = cmMemAllocZ(unsigned,N);
  3312. unsigned i = 0;
  3313. cmVectArrayVect_t* pos = p->cur;
  3314. if( cmVectArrayRewind(p) != cmOkRC )
  3315. goto errLabel;
  3316. if( cmVectArrayAdvance(p,groupIdx) != cmOkRC )
  3317. goto errLabel;
  3318. while( i<N && !cmVectArrayIsEOL(p) )
  3319. {
  3320. unsigned tn = cmVectArrayEleCount(p);
  3321. unsigned tv[tn];
  3322. // read the sub-vector
  3323. cmVectArrayGetU(p,tv,&tn);
  3324. assert( colIdx < tn );
  3325. // store the output value
  3326. if( colIdx < tn )
  3327. v[i++] = tv[colIdx];
  3328. cmVectArrayAdvance(p,groupCnt);
  3329. }
  3330. *vRef = v;
  3331. *vnRef = i;
  3332. errLabel:
  3333. p->cur = pos;
  3334. return rc;
  3335. }
  3336. cmRC_t cmVectArrayTest( cmCtx* ctx, const char* fn )
  3337. {
  3338. cmRC_t rc = cmOkRC;
  3339. cmVectArray_t* p = NULL;
  3340. unsigned flags = kSampleVaFl;
  3341. cmSample_t v[] = { 0, 1, 2, 3, 4, 5 };
  3342. if( fn == NULL || strlen(fn)==0 )
  3343. return cmCtxRtCondition(&p->obj,cmSubSysFailRC,"Invalid test output file name.");
  3344. if( (p = cmVectArrayAlloc(ctx,flags)) == NULL )
  3345. return cmCtxRtCondition(&p->obj,cmSubSysFailRC,"The vectory array object allocation failed.");
  3346. if( cmVectArrayAppendS(p,v,1) != cmOkRC )
  3347. {
  3348. rc = cmCtxRtCondition(&p->obj,cmSubSysFailRC,"Vector append 1 failed.");
  3349. goto errLabel;
  3350. }
  3351. if( cmVectArrayAppendS(p,v+1,2) != cmOkRC )
  3352. {
  3353. rc = cmCtxRtCondition(&p->obj,cmSubSysFailRC,"Vector append 2 failed.");
  3354. goto errLabel;
  3355. }
  3356. if( cmVectArrayAppendS(p,v+3,3) != cmOkRC )
  3357. {
  3358. rc = cmCtxRtCondition(&p->obj,cmSubSysFailRC,"Vector append 3 failed.");
  3359. goto errLabel;
  3360. }
  3361. if( cmVectArrayWrite(p,fn) != cmOkRC )
  3362. {
  3363. rc = cmCtxRtCondition(&p->obj,cmSubSysFailRC,"Vector array write failed.");
  3364. goto errLabel;
  3365. }
  3366. //cmVectArrayForEachS(p,0,cmVectArrayEleCount(p),cmVectArrayForEachTextFuncS,&ctx->printRpt);
  3367. if( cmVectArrayFree(&p) != cmOkRC )
  3368. {
  3369. rc = cmCtxRtCondition(&p->obj,cmSubSysFailRC,"The vectory array release failed.");
  3370. goto errLabel;
  3371. }
  3372. if((p = cmVectArrayAllocFromFile(ctx, fn )) == NULL )
  3373. {
  3374. rc = cmCtxRtCondition(&p->obj,cmSubSysFailRC,"VectArray alloc from file failed.");
  3375. goto errLabel;
  3376. }
  3377. while(!cmVectArrayIsEOL(p))
  3378. {
  3379. unsigned n = cmVectArrayEleCount(p);
  3380. cmSample_t v[n];
  3381. if( cmVectArrayGetS(p,v,&n) != cmOkRC )
  3382. {
  3383. rc = cmCtxRtCondition(&p->obj,cmSubSysFailRC,"VectArrayGetS() failed.");
  3384. goto errLabel;
  3385. }
  3386. cmVOS_PrintL("v:",NULL,1,n,v);
  3387. cmVectArrayAdvance(p,1);
  3388. }
  3389. errLabel:
  3390. if( cmVectArrayFree(&p) != cmOkRC )
  3391. rc = cmCtxRtCondition(&p->obj,cmSubSysFailRC,"The vectory array release failed.");
  3392. return rc;
  3393. }
  3394. //-----------------------------------------------------------------------------------------------------------------------
  3395. cmWhFilt* cmWhFiltAlloc( cmCtx* c, cmWhFilt* p, unsigned binCnt, cmReal_t binHz, cmReal_t coeff, cmReal_t maxHz )
  3396. {
  3397. cmWhFilt* op = cmObjAlloc(cmWhFilt,c,p);
  3398. if( binCnt > 0 )
  3399. if( cmWhFiltInit(op,binCnt,binHz,coeff,maxHz) != cmOkRC )
  3400. cmWhFiltFree(&op);
  3401. return op;
  3402. }
  3403. cmRC_t cmWhFiltFree( cmWhFilt** pp )
  3404. {
  3405. cmRC_t rc = cmOkRC;
  3406. if( pp==NULL || *pp==NULL )
  3407. return rc;
  3408. cmWhFilt* p = *pp;
  3409. if((rc = cmWhFiltFinal(p)) != cmOkRC )
  3410. return rc;
  3411. cmMemFree(p->whM);
  3412. cmMemFree(p->whiV);
  3413. cmMemFree(p->iV);
  3414. cmObjFree(pp);
  3415. return rc;
  3416. }
  3417. cmRC_t cmWhFiltInit( cmWhFilt* p, unsigned binCnt, cmReal_t binHz, cmReal_t coeff, cmReal_t maxHz )
  3418. {
  3419. cmRC_t rc;
  3420. if((rc = cmWhFiltFinal(p)) != cmOkRC )
  3421. return rc;
  3422. p->binCnt = binCnt;
  3423. p->binHz = binHz;
  3424. p->bandCnt = maxHz == 0 ? 34 : ceil(log10(maxHz/229.0 + 1) * 21.4 - 1)-1;
  3425. if( p->bandCnt <= 0 )
  3426. return cmCtxRtCondition(&p->obj, cmInvalidArgRC, "Max. Hz too low to form any frequency bands.");
  3427. cmReal_t flV[ p->bandCnt ];
  3428. cmReal_t fcV[ p->bandCnt ];
  3429. cmReal_t fhV[ p->bandCnt ];
  3430. int i;
  3431. for(i=0; i<p->bandCnt; ++i)
  3432. {
  3433. fcV[i] = 229.0 * (pow(10.0,(i+2)/21.4) - 1.0);
  3434. flV[i] = i==0 ? 0 : fcV[i-1];
  3435. }
  3436. for(i=0; i<p->bandCnt-1; ++i)
  3437. fhV[i] = fcV[i+1];
  3438. fhV[p->bandCnt-1] = fcV[p->bandCnt-1] + (fcV[p->bandCnt-1] - fcV[p->bandCnt-2]);
  3439. //cmVOR_PrintL("flV",NULL,1,p->bandCnt,flV);
  3440. //cmVOR_PrintL("fcV",NULL,1,p->bandCnt,fcV);
  3441. //cmVOR_PrintL("fhV",NULL,1,p->bandCnt,fhV);
  3442. cmReal_t* tM = cmMemAlloc(cmReal_t,p->bandCnt * p->binCnt);
  3443. p->whM = cmMemResizeZ(cmReal_t,p->whM,p->binCnt * p->bandCnt);
  3444. p->iV = cmMemResizeZ(cmReal_t,p->iV,p->binCnt);
  3445. // generate the bin index values
  3446. for(i=0; i<p->binCnt; ++i)
  3447. p->iV[i] = i;
  3448. cmReal_t stSpread = 0; // set stSpread to 0 to use flV/fhV[]
  3449. cmVOR_TriangleMask(tM, p->bandCnt, p->binCnt, fcV, p->binHz, stSpread, flV, fhV );
  3450. cmVOR_Transpose(p->whM, tM, p->bandCnt, p->binCnt );
  3451. cmMemFree(tM);
  3452. //cmVOR_PrintL("whM",NULL,p->bandCnt,p->binCnt,p->whM);
  3453. //cmVectArrayWriteMatrixR(p->obj.ctx, "/home/kevin/temp/frqtrk/whM.va", p->whM, p->binCnt, p->bandCnt );
  3454. unsigned whiN = p->bandCnt+2;
  3455. p->whiV = cmMemResizeZ(cmReal_t,p->whiV,whiN);
  3456. for(i=0; i<whiN; ++i)
  3457. {
  3458. if( i == 0 )
  3459. p->whiV[i] = 0;
  3460. else
  3461. if( i == whiN-1 )
  3462. p->whiV[i] = fhV[p->bandCnt-1]/binHz;
  3463. else
  3464. p->whiV[i] = fcV[i-1]/binHz;
  3465. }
  3466. //cmVOR_PrintL("whiV",NULL,1,whiN,p->whiV);
  3467. //cmVectArrayWriteMatrixR(p->obj.ctx, "/home/kevin/temp/frqtrk/whiV.va", p->whiV, whiN, 1 );
  3468. return rc;
  3469. }
  3470. cmRC_t cmWhFiltFinal( cmWhFilt* p )
  3471. { return cmOkRC; }
  3472. cmRC_t cmWhFiltExec( cmWhFilt* p, const cmReal_t* xV, cmReal_t* yV, unsigned xyN )
  3473. {
  3474. assert( xyN == p->binCnt);
  3475. cmRC_t rc = cmOkRC;
  3476. unsigned whiN = p->bandCnt + 2;
  3477. unsigned mbi = cmMin(xyN, floor(p->whiV[whiN-1]));
  3478. // calculate the level in each band to form a composite filter
  3479. cmReal_t y0V[ whiN ];
  3480. cmReal_t* b0V = y0V + 1;
  3481. cmVOR_MultVVM(b0V, p->bandCnt, xV, p->binCnt, p->whM );
  3482. //cmVOR_PrintL("b0V",NULL,1,p->bandCnt,b0V);
  3483. // BEWARE: zeros in b0V will generate Inf's when sent
  3484. // through the cmVOR_PowVS() function.
  3485. int i;
  3486. for(i=0; i<p->bandCnt; ++i)
  3487. if( b0V[i] < 0.000001 )
  3488. b0V[i] = 0.000001;
  3489. // apply a non-linear expansion function to each band
  3490. cmVOR_PowVS(b0V,p->bandCnt,p->coeff-1);
  3491. //cmVOR_PrintL("b0V",NULL,1,p->bandCnt,b0V);
  3492. // add edge values to the filter
  3493. y0V[0] = b0V[0];
  3494. y0V[whiN-1] = b0V[p->bandCnt-1];
  3495. //cmVOR_PrintL("y0V",NULL,1,whiN,y0V);
  3496. cmVOR_Interp1(yV,p->iV,p->binCnt,p->whiV,y0V,whiN);
  3497. cmVOR_Fill(yV+mbi,xyN-mbi,1.0);
  3498. //cmVOR_PrintL("yV",NULL,1,p->binCnt,yV);
  3499. cmVOR_MultVV(yV,xyN,xV);
  3500. return rc;
  3501. }
  3502. //-----------------------------------------------------------------------------------------------------------------------
  3503. cmFrqTrk* cmFrqTrkAlloc( cmCtx* c, cmFrqTrk* p, const cmFrqTrkArgs_t* a )
  3504. {
  3505. cmFrqTrk* op = cmObjAlloc(cmFrqTrk,c,p);
  3506. op->logVa = cmVectArrayAlloc(c,kRealVaFl);
  3507. op->levelVa = cmVectArrayAlloc(c,kRealVaFl);
  3508. op->specVa = cmVectArrayAlloc(c,kRealVaFl);
  3509. op->attenVa = cmVectArrayAlloc(c,kRealVaFl);
  3510. op->wf = cmWhFiltAlloc(c,NULL,0,0,0,0);
  3511. if( a != NULL )
  3512. if( cmFrqTrkInit(op,a) != cmOkRC )
  3513. cmFrqTrkFree(&op);
  3514. return op;
  3515. }
  3516. cmRC_t cmFrqTrkFree( cmFrqTrk** pp )
  3517. {
  3518. cmRC_t rc = cmOkRC;
  3519. if( pp==NULL || *pp==NULL )
  3520. return rc;
  3521. cmFrqTrk* p = *pp;
  3522. if((rc = cmFrqTrkFinal(p)) != cmOkRC )
  3523. return rc;
  3524. unsigned i;
  3525. for(i=0; i<p->a.chCnt; ++i)
  3526. {
  3527. cmMemFree(p->ch[i].dbV);
  3528. cmMemFree(p->ch[i].hzV);
  3529. }
  3530. cmMemFree(p->ch);
  3531. cmMemFree(p->dbM);
  3532. cmMemFree(p->pkiV);
  3533. cmMemFree(p->dbV);
  3534. cmMemFree(p->aV);
  3535. cmVectArrayFree(&p->logVa);
  3536. cmVectArrayFree(&p->levelVa);
  3537. cmVectArrayFree(&p->specVa);
  3538. cmVectArrayFree(&p->attenVa);
  3539. cmWhFiltFree(&p->wf);
  3540. cmMemFree(p->logFn);
  3541. cmMemFree(p->levelFn);
  3542. cmMemFree(p->specFn);
  3543. cmObjFree(pp);
  3544. return rc;
  3545. }
  3546. cmRC_t cmFrqTrkInit( cmFrqTrk* p, const cmFrqTrkArgs_t* a )
  3547. {
  3548. cmRC_t rc;
  3549. if((rc = cmFrqTrkFinal(p)) != cmOkRC )
  3550. return rc;
  3551. p->a = *a;
  3552. p->ch = cmMemResizeZ(cmFrqTrkCh_t,p->ch,a->chCnt );
  3553. p->hN = cmMax(1,a->wndSecs * a->srate / a->hopSmpCnt );
  3554. p->sN = 4*p->hN;
  3555. p->binHz = a->srate / ((p->a.binCnt-1)*2);
  3556. p->bN = cmMin( p->a.binCnt, ceil(p->a.pkMaxHz / p->binHz ));
  3557. p->dbM = cmMemResizeZ(cmReal_t,p->dbM,p->hN*p->bN);
  3558. p->hi = 0;
  3559. p->fN = 0;
  3560. p->dbV = cmMemResizeZ(cmReal_t,p->dbV,p->bN);
  3561. p->pkiV = cmMemResizeZ(unsigned,p->pkiV,p->bN);
  3562. p->deadN_max = a->maxTrkDeadSec * a->srate / a->hopSmpCnt;
  3563. p->minTrkN = a->minTrkSec * a->srate / a->hopSmpCnt;
  3564. p->nextTrkId = 1;
  3565. p->aV = cmMemResizeZ(cmReal_t,p->aV,p->a.binCnt);
  3566. p->attenDlyPhsMax = cmMax(3,a->attenDlySec * a->srate / a->hopSmpCnt );
  3567. p->attenPhsMax = cmMax(3,a->attenAtkSec * a->srate / a->hopSmpCnt );
  3568. if( a->logFn != NULL )
  3569. p->logFn = cmMemResizeStr(p->logFn,a->logFn);
  3570. if( a->levelFn != NULL )
  3571. p->levelFn = cmMemResizeStr(p->levelFn,a->levelFn);
  3572. if( a->specFn != NULL )
  3573. p->specFn = cmMemResizeStr(p->specFn,a->specFn);
  3574. if( a->attenFn != NULL )
  3575. p->attenFn = cmMemResizeStr(p->attenFn,a->attenFn);
  3576. if(cmWhFiltInit(p->wf,p->bN,p->binHz,p->a.whFiltCoeff,p->a.pkMaxHz) != cmOkRC )
  3577. cmCtxRtCondition(&p->obj, cmSubSysFailRC, "Whitening filter intitialization failed.");
  3578. unsigned i;
  3579. for(i=0; i<p->a.chCnt; ++i)
  3580. {
  3581. p->ch[i].dbV = cmMemResizeZ(cmReal_t,p->ch[i].dbV,p->sN);
  3582. p->ch[i].hzV = cmMemResizeZ(cmReal_t,p->ch[i].hzV,p->sN);
  3583. }
  3584. return rc;
  3585. }
  3586. cmRC_t cmFrqTrkFinal( cmFrqTrk* p )
  3587. {
  3588. cmRC_t rc = cmOkRC;
  3589. if( p->logFn != NULL )
  3590. cmVectArrayWrite(p->logVa,p->logFn);
  3591. if( p->levelFn != NULL )
  3592. cmVectArrayWrite(p->levelVa,p->levelFn);
  3593. if( p->specFn != NULL )
  3594. cmVectArrayWrite(p->specVa,p->specFn);
  3595. if( p->attenFn != NULL )
  3596. cmVectArrayWrite(p->attenVa,p->attenFn);
  3597. cmWhFiltFinal(p->wf);
  3598. return rc;
  3599. }
  3600. // Return an available channel record or NULL if all channel records are in use.
  3601. cmFrqTrkCh_t* _cmFrqTrkFindAvailCh( cmFrqTrk* p )
  3602. {
  3603. unsigned i;
  3604. for(i=0; i<p->a.chCnt; ++i)
  3605. if( p->ch[i].activeFl == false )
  3606. return p->ch + i;
  3607. return NULL;
  3608. }
  3609. // Estimate the peak frequency by parabolic interpolotion into hzV[p->bN]
  3610. void _cmFrqTrkMagnToHz( cmFrqTrk* p, const cmReal_t* dbV, unsigned* pkiV, unsigned pkN, cmReal_t* hzV )
  3611. {
  3612. unsigned i;
  3613. for(i=0; i<pkN; ++i)
  3614. if( pkiV[i] != cmInvalidIdx )
  3615. {
  3616. unsigned pki = pkiV[i];
  3617. cmReal_t y0 = pki>0 ? dbV[ pki-1 ] : dbV[pki];
  3618. cmReal_t y1 = dbV[ pki ];
  3619. cmReal_t y2 = pki<p->bN-1 ? dbV[ pki+1 ] : dbV[pki];
  3620. cmReal_t den = y0 - (2.*y1) + y2;
  3621. cmReal_t offs = den==0 ? 0 : 0.5 * ((y0 - y2) / den);
  3622. hzV[pki] = p->binHz * (pki+offs);
  3623. //if( hzV[pki] < 0 )
  3624. //{
  3625. // printf("%f : %f %f %f : %f %f\n",hzV[pki],y0,y1,y2,den,offs);
  3626. //}
  3627. }
  3628. }
  3629. unsigned _cmFrqTrkActiveChCount( cmFrqTrk* p )
  3630. {
  3631. unsigned n = 0;
  3632. unsigned i;
  3633. for(i=0; i<p->a.chCnt; ++i)
  3634. if( p->ch[i].activeFl )
  3635. ++n;
  3636. return n;
  3637. }
  3638. void _cmFrqTrkWriteLevel( cmFrqTrk* p, const cmReal_t* dbV, const cmReal_t* hzV, unsigned bN )
  3639. {
  3640. if( p->levelFn != NULL )
  3641. {
  3642. double maxHz = 5000.0;
  3643. unsigned maxBinIdx = cmMin(bN,maxHz / p->binHz);
  3644. unsigned vn = 3;
  3645. cmReal_t v[vn];
  3646. unsigned idx = cmVOR_MaxIndex(dbV,maxBinIdx,1);
  3647. v[0] = cmVOR_Mean(dbV,maxBinIdx);
  3648. v[1] = dbV[idx];
  3649. v[2] = hzV[idx];
  3650. cmVectArrayAppendR(p->levelVa,v,vn);
  3651. }
  3652. }
  3653. void _cmFrqTrkWriteLog( cmFrqTrk* p )
  3654. {
  3655. unsigned n;
  3656. cmReal_t* vb = NULL;
  3657. if( p->logFn == NULL )
  3658. return;
  3659. if((n = _cmFrqTrkActiveChCount(p)) > 0 )
  3660. {
  3661. unsigned i,j;
  3662. // sn = count of elements in the summary sub-vector
  3663. unsigned sn = 3;
  3664. // each active channel will emit 7 values
  3665. unsigned nn = 1 + n*7 + sn;
  3666. // allocate the row vector
  3667. vb = cmMemResize(cmReal_t,vb,nn);
  3668. // row format
  3669. // [ nn idV[n] hzV[n] ... hsV[n] smV[sn] ]
  3670. // n = (nn - (1 + sn)) / 7
  3671. *vb = nn; // the first element in the vector contains the length of the row
  3672. cmReal_t* v = vb + 1;
  3673. // setup the base pointer to each sub-vector
  3674. cmReal_t* idV = v + n * 0;
  3675. cmReal_t* hzV = v + n * 1;
  3676. cmReal_t* dbV = v + n * 2;
  3677. cmReal_t* stV = v + n * 3;
  3678. cmReal_t* dsV = v + n * 4;
  3679. cmReal_t* hsV = v + n * 5;
  3680. cmReal_t* agV = v + n * 6;
  3681. cmReal_t* smV = v + n * 7; // summary information
  3682. smV[0] = p->newTrkCnt;
  3683. smV[1] = p->curTrkCnt;
  3684. smV[2] = p->deadTrkCnt;
  3685. // for each active channel
  3686. for(i=0,j=0; i<p->a.chCnt; ++i)
  3687. if( p->ch[i].activeFl )
  3688. {
  3689. assert(j < n);
  3690. // elements of each sub-vector associated with a given
  3691. // index refer to the same track record - element i therefore
  3692. // refers to active track index i.
  3693. idV[j] = p->ch[i].id;
  3694. hzV[j] = p->ch[i].hz;
  3695. dbV[j] = p->ch[i].db;
  3696. stV[j] = p->ch[i].dN;
  3697. dsV[j] = p->ch[i].db_std;
  3698. hsV[j] = p->ch[i].hz_std;
  3699. agV[j] = p->ch[i].attenGain;
  3700. ++j;
  3701. }
  3702. cmVectArrayAppendR(p->logVa, vb, nn );
  3703. }
  3704. cmMemFree(vb);
  3705. }
  3706. void _cmFrqTrkPrintChs( const cmFrqTrk* p )
  3707. {
  3708. unsigned i;
  3709. for(i=0; i<p->a.chCnt; ++i)
  3710. {
  3711. cmFrqTrkCh_t* c = p->ch + i;
  3712. printf("%i : %i tN:%i hz:%f db:%f\n",i,c->activeFl,c->tN,c->hz,c->db);
  3713. }
  3714. }
  3715. // Used to sort the channels into descending dB order.
  3716. int _cmFrqTrkChCompare( const void* p0, const void* p1 )
  3717. { return ((cmFrqTrkCh_t*)p0)->db - ((cmFrqTrkCh_t*)p1)->db; }
  3718. // Return the index of the peak associated with pkiV[i] which best matches the tracker 'c'
  3719. // or cmInvalidIdx if no valid peaks were found.
  3720. // pkiV[ pkN ] holds the indexes into dbV[] and hzV[] which are peaks.
  3721. // Some elements of pkiV[] may be set to cmInvalidIdx if the associated peak has already
  3722. // been selected by another tracker.
  3723. unsigned _cmFrqTrkFindPeak( cmFrqTrk* p, const cmFrqTrkCh_t* c, const cmReal_t* dbV, const cmReal_t* hzV, unsigned* pkiV, unsigned pkN )
  3724. {
  3725. unsigned i,pki;
  3726. cmReal_t d_max = p->a.pkThreshDb;
  3727. unsigned d_idx = cmInvalidIdx;
  3728. cmReal_t hz_min = c->hz * pow(2,-p->a.stRange/12.0);
  3729. cmReal_t hz_max = c->hz * pow(2, p->a.stRange/12.0);
  3730. // find the peak with the most energy inside the frequency range hz_min to hz_max.
  3731. for(i=0; i<pkN; ++i)
  3732. if( ((pki = pkiV[i]) != cmInvalidIdx) && hz_min <= hzV[pki] && hzV[pki] <= hz_max && dbV[pki]>d_max )
  3733. {
  3734. d_max= dbV[pki];
  3735. d_idx = i;
  3736. }
  3737. return d_idx;
  3738. }
  3739. void _cmFrqTrkScoreChs( cmFrqTrk* p )
  3740. {
  3741. unsigned i;
  3742. for(i=0; i<p->a.chCnt; ++i)
  3743. if( p->ch[i].activeFl )
  3744. {
  3745. cmFrqTrkCh_t* c = p->ch + i;
  3746. c->dbV[ c->si ] = c->db;
  3747. c->hzV[ c->si ] = c->hz;
  3748. c->si = (c->si + 1) % p->sN;
  3749. c->sn += 1;
  3750. unsigned n = cmMin(c->sn,p->sN);
  3751. c->db_mean = cmVOR_Mean(c->dbV,n);
  3752. c->db_std = sqrt(cmVOR_Variance( c->dbV,n,&c->db_mean));
  3753. c->hz_mean = cmVOR_Mean(c->hzV,n);
  3754. c->hz_std = sqrt(cmVOR_Variance( c->hzV,n,&c->hz_mean));
  3755. //c->score = c->db / ((cmMax(0.1,c->db_std) + cmMax(0.1,c->hz_std))/2);
  3756. c->score = c->db - (c->db_std * 5) - (c->hz_std/50);
  3757. //printf("%f %f %f %f %f\n",c->db,cmMin(0.1,c->db_std),c->hz,cmMin(0.1,c->hz_std),c->score);
  3758. }
  3759. }
  3760. // Generate a filter that is wider for higher frequencies than lower frequencies.
  3761. unsigned _cmFrqTrkFillMap( cmFrqTrk* p, cmReal_t* map, unsigned maxN, cmReal_t hz )
  3762. {
  3763. assert( maxN % 2 == 1 );
  3764. unsigned i;
  3765. cmReal_t maxHz = p->a.srate/2;
  3766. unsigned mapN = cmMin(maxN,ceil(hz/maxHz * maxN));
  3767. if( mapN % 2 == 0 )
  3768. mapN += 1;
  3769. mapN = cmMin(maxN,mapN);
  3770. unsigned N = floor(mapN/2);
  3771. double COEFF = 0.3;
  3772. for(i=0; i<N; ++i)
  3773. {
  3774. map[i] = pow(((double)i+1)/(N+1),COEFF);
  3775. map[mapN-(i+1)] = map[i];
  3776. }
  3777. map[N] = 1.0;
  3778. return mapN;
  3779. }
  3780. void _cmFrqTrkApplyAtten( cmFrqTrk* p, cmReal_t* aV, cmReal_t gain, cmReal_t hz )
  3781. {
  3782. int cbi = cmMin(p->a.binCnt,cmMax(0,round(hz/p->binHz)));
  3783. //cmReal_t map[] = { .25, .5, 1, .5, .25 };
  3784. //int mapN = sizeof(map)/sizeof(map[0]);
  3785. unsigned maxN = 30; // must be odd
  3786. cmReal_t map[ maxN ];
  3787. int mapN = _cmFrqTrkFillMap(p, map, maxN, hz );
  3788. int j;
  3789. int ai = cbi - mapN/2;
  3790. for(j=0; j<mapN; ++j,++ai)
  3791. if( 0 <= ai && ai < p->a.binCnt )
  3792. aV[ai] *= 1.0 - (map[j] * gain);
  3793. }
  3794. void _cmFrqTrkUpdateFilter( cmFrqTrk* p )
  3795. {
  3796. unsigned i;
  3797. cmVOR_Fill(p->aV,p->a.binCnt,1.0);
  3798. for(i=0; i<p->a.chCnt; ++i)
  3799. if( p->ch[i].activeFl )
  3800. {
  3801. cmFrqTrkCh_t* c = p->ch + i;
  3802. //
  3803. if( c->score >= p->a.attenThresh && c->state == kNoStateFrqTrkId )
  3804. {
  3805. //printf("%f\n",c->score);
  3806. c->attenPhsIdx = 0;
  3807. c->state = kDlyFrqTrkId;
  3808. }
  3809. switch( c->state )
  3810. {
  3811. case kNoStateFrqTrkId:
  3812. break;
  3813. case kDlyFrqTrkId:
  3814. c->attenPhsIdx += 1;
  3815. if( c->attenPhsIdx >= p->attenDlyPhsMax && c->dN == 0 )
  3816. c->state = kAtkFrqTrkId;
  3817. break;
  3818. case kAtkFrqTrkId:
  3819. if( c->attenPhsIdx < p->attenDlyPhsMax + p->attenPhsMax )
  3820. {
  3821. c->attenGain = cmMin(1.0,p->a.attenGain * c->attenPhsIdx / p->attenPhsMax);
  3822. _cmFrqTrkApplyAtten(p, p->aV, c->attenGain, c->hz);
  3823. }
  3824. c->attenPhsIdx += 1;
  3825. if( c->attenPhsIdx >= p->attenDlyPhsMax + p->attenPhsMax )
  3826. c->state = kSusFrqTrkId;
  3827. break;
  3828. case kSusFrqTrkId:
  3829. if( c->dN > 0 )
  3830. {
  3831. if( c->attenPhsIdx > 0 )
  3832. {
  3833. c->attenPhsIdx -= 1;
  3834. c->attenGain = cmMin(1.0,p->a.attenGain * c->attenPhsIdx / p->attenPhsMax);
  3835. }
  3836. }
  3837. _cmFrqTrkApplyAtten(p,p->aV, c->attenGain, c->hz);
  3838. if( c->dN >= p->deadN_max )
  3839. c->state = kDcyFrqTrkId;
  3840. break;
  3841. case kDcyFrqTrkId:
  3842. if( c->attenPhsIdx > 0 )
  3843. {
  3844. c->attenPhsIdx -= 1;
  3845. c->attenGain = cmMin(1.0,p->a.attenGain * c->attenPhsIdx / p->attenPhsMax);
  3846. _cmFrqTrkApplyAtten(p,p->aV, c->attenGain, c->hz);
  3847. }
  3848. if( c->attenPhsIdx == 0 )
  3849. c->activeFl = false;
  3850. break;
  3851. }
  3852. }
  3853. }
  3854. // Extend the existing trackers
  3855. void _cmFrqTrkExtendChs( cmFrqTrk* p, const cmReal_t* dbV, const cmReal_t* hzV, unsigned* pkiV, unsigned pkN )
  3856. {
  3857. unsigned i;
  3858. p->curTrkCnt = 0;
  3859. p->deadTrkCnt = 0;
  3860. // sort the channels in descending order
  3861. qsort(p->ch,p->a.chCnt,sizeof(cmFrqTrkCh_t),_cmFrqTrkChCompare);
  3862. // for each active channel
  3863. for(i=0; i<p->a.chCnt; ++i)
  3864. {
  3865. cmFrqTrkCh_t* c = p->ch + i;
  3866. if( c->activeFl )
  3867. {
  3868. unsigned pki;
  3869. // find the best peak to extend tracker 'c'.
  3870. if((pki = _cmFrqTrkFindPeak(p,c,dbV,hzV,pkiV,pkN)) == cmInvalidIdx )
  3871. {
  3872. // no valid track was found to extend tracker 'c'
  3873. c->dN += 1;
  3874. c->tN += 1;
  3875. if( c->dN >= p->deadN_max )
  3876. {
  3877. if( c->attenPhsIdx == 0 )
  3878. c->activeFl = false;
  3879. p->deadTrkCnt += 1;
  3880. }
  3881. }
  3882. else // ... update the tracker using the matching peak
  3883. {
  3884. unsigned j = pkiV[pki];
  3885. c->dN = 0;
  3886. c->db = dbV[ j ];
  3887. c->hz = hzV[ j ];
  3888. c->tN += 1;
  3889. pkiV[pki] = cmInvalidIdx; // mark the peak as unavailable.
  3890. p->curTrkCnt += 1;
  3891. }
  3892. }
  3893. }
  3894. }
  3895. // disable peaks which are within 'stRange' semitones of the frequency of active trackers.
  3896. void _cmFrqTrkDisableClosePeaks( cmFrqTrk* p, const cmReal_t* dbV, const cmReal_t* hzV, unsigned* pkiV, unsigned pkN )
  3897. {
  3898. unsigned i;
  3899. for(i=0; i<p->a.chCnt; ++i)
  3900. {
  3901. const cmFrqTrkCh_t* c = p->ch + i;
  3902. if( !c->activeFl )
  3903. continue;
  3904. cmReal_t hz_min = c->hz * pow(2,-p->a.stRange/12.0);
  3905. cmReal_t hz_max = c->hz * pow(2, p->a.stRange/12.0);
  3906. unsigned j;
  3907. // find all peaks within the frequency range hz_min to hz_max.
  3908. for(j=0; j<pkN; ++j)
  3909. if( pkiV[j] != cmInvalidIdx && hz_min <= c->hz && c->hz <= hz_max )
  3910. pkiV[j] = cmInvalidIdx;
  3911. }
  3912. }
  3913. // Return the index into pkiV[] of the maximum energy peak in dbV[]
  3914. // that is also above kAtkThreshDb.
  3915. unsigned _cmFrqTrkMaxEnergyPeakIndex( const cmFrqTrk* p, const cmReal_t* dbV, const cmReal_t* hzV, const unsigned* pkiV, unsigned pkN )
  3916. {
  3917. cmReal_t mv = p->a.pkAtkThreshDb;
  3918. unsigned mi = cmInvalidIdx;
  3919. unsigned i;
  3920. for(i=0; i<pkN; ++i)
  3921. if( pkiV[i] != cmInvalidIdx && dbV[pkiV[i]] >= mv && hzV[pkiV[i]] < p->a.pkMaxHz )
  3922. {
  3923. mi = i;
  3924. mv = dbV[pkiV[i]];
  3925. }
  3926. return mi;
  3927. }
  3928. // start new trackers
  3929. void _cmFrqTrkNewChs( cmFrqTrk* p, const cmReal_t* dbV, const cmReal_t* hzV, unsigned* pkiV, unsigned pkN )
  3930. {
  3931. p->newTrkCnt = 0;
  3932. while(1)
  3933. {
  3934. unsigned db_max_idx;
  3935. cmFrqTrkCh_t* c;
  3936. // find an inactive channel
  3937. if((c = _cmFrqTrkFindAvailCh(p)) == NULL )
  3938. break;
  3939. // find the largest peak that is above pkAtkThreshDb && less than pkAtkHz.
  3940. if((db_max_idx = _cmFrqTrkMaxEnergyPeakIndex(p,dbV,hzV,pkiV,pkN)) == cmInvalidIdx )
  3941. break;
  3942. // activate a new channel
  3943. c->activeFl = true;
  3944. c->tN = 1;
  3945. c->dN = 0;
  3946. c->hz = hzV[ pkiV[ db_max_idx ] ];
  3947. c->db = dbV[ pkiV[ db_max_idx ] ];
  3948. c->id = p->nextTrkId++;
  3949. c->si = 0;
  3950. c->sn = 0;
  3951. c->score = 0;
  3952. c->state = kNoStateFrqTrkId;
  3953. c->attenPhsIdx = cmInvalidIdx;
  3954. c->attenGain = 1.0;
  3955. // mark the peak as unavailable
  3956. pkiV[ db_max_idx ] = cmInvalidIdx;
  3957. p->newTrkCnt += 1;
  3958. }
  3959. }
  3960. void _cmFrqTrkApplyFrqBias( cmFrqTrk* p, cmReal_t* xV )
  3961. {
  3962. // convert to decibel scale (0.0 - 100.0) and then scale to (0.0 to 1.0)
  3963. unsigned i;
  3964. for(i=0; i<p->bN; ++i)
  3965. xV[i] = cmMax(0.0, (20*log10( cmMax(xV[i]/1.5,0.00001)) + 100.0)/100.0);
  3966. }
  3967. cmRC_t cmFrqTrkExec( cmFrqTrk* p, const cmReal_t* magV, const cmReal_t* phsV, const cmReal_t* hertzV )
  3968. {
  3969. cmRC_t rc = cmOkRC;
  3970. cmReal_t hzV[ p->bN ];
  3971. //cmReal_t powV[ p->bN ];
  3972. //cmReal_t yV[ p->bN];
  3973. //cmVOR_MultVVV(powV,p->bN,magV,magV);
  3974. //cmWhFiltExec(p->wf,powV,p->dbV,p->bN);
  3975. // convert magV to Decibels
  3976. //cmVOR_AmplToDbVV(p->dbV,p->bN, magV, -200.0);
  3977. // copy p->dbV to dbM[hi,:]
  3978. //cmVOR_CopyN(p->dbM + p->hi, p->bN, p->hN, p->dbV, 1 );
  3979. //cmVOR_CopyN(p->dbM + p->hi, p->bN, p->hN, whV, 1 );
  3980. if( 1 )
  3981. {
  3982. cmReal_t powV[ p->bN ];
  3983. cmVOR_MultVVV(powV,p->bN,magV,magV);
  3984. cmWhFiltExec(p->wf,powV,p->dbV,p->bN);
  3985. _cmFrqTrkApplyFrqBias(p,p->dbV);
  3986. }
  3987. else
  3988. {
  3989. // convert magV to Decibels
  3990. cmVOR_AmplToDbVV(p->dbV,p->bN, magV, -200.0);
  3991. }
  3992. // copy p->dbV to dbM[hi,:]
  3993. cmVOR_CopyN(p->dbM + p->hi, p->bN, p->hN, p->dbV, 1 );
  3994. // increment hi to next column to fill in dbM[]
  3995. p->hi = (p->hi + 1) % p->hN;
  3996. // Set dbV[] to spectral magnitude profile by taking the mean over time
  3997. // of the last hN magnitude vectors
  3998. cmVOR_MeanM2(p->dbV, p->dbM, p->hN, p->bN, 0, cmMin(p->fN+1,p->hN));
  3999. //cmVOR_MeanM(p->dbV, p->dbM, p->hN, p->bN, 0);
  4000. if( p->fN >= p->hN )
  4001. {
  4002. // set pkiV[] to the indexes of the peaks above pkThreshDb in i0[]
  4003. unsigned pkN = cmVOR_PeakIndexes(p->pkiV, p->bN, p->dbV, p->bN, p->a.pkThreshDb );
  4004. // set hzV[] to the peak frequencies assoc'd with peaks at dbV[ pkiV[] ].
  4005. _cmFrqTrkMagnToHz(p, p->dbV, p->pkiV, pkN, hzV );
  4006. // extend the existing trackers
  4007. _cmFrqTrkExtendChs(p, p->dbV, hzV, p->pkiV, pkN );
  4008. //_cmFrqTrkDisableClosePeaks(p, p->dbV, hzV, p->pkiV, pkN );
  4009. // create new trackers
  4010. _cmFrqTrkNewChs(p,p->dbV,hzV,p->pkiV,pkN);
  4011. //
  4012. _cmFrqTrkScoreChs(p);
  4013. //
  4014. _cmFrqTrkUpdateFilter(p);
  4015. /*
  4016. // write the log file
  4017. _cmFrqTrkWriteLog(p);
  4018. // write the spectrum output file
  4019. if( p->specFn != NULL )
  4020. cmVectArrayAppendR(p->specVa,p->dbV,p->bN);
  4021. // write the atten output file
  4022. if( p->attenFn != NULL )
  4023. cmVectArrayAppendR(p->attenVa,p->aV,p->bN);
  4024. // write the the level file
  4025. _cmFrqTrkWriteLevel(p,p->dbV,hzV,p->bN);
  4026. */
  4027. }
  4028. p->fN += 1;
  4029. return rc;
  4030. }
  4031. void cmFrqTrkPrint( cmFrqTrk* p )
  4032. {
  4033. printf("srate: %f\n",p->a.srate);
  4034. printf("chCnt: %i\n",p->a.chCnt);
  4035. printf("binCnt: %i (bN=%i)\n",p->a.binCnt,p->bN);
  4036. printf("hopSmpCnt: %i\n",p->a.hopSmpCnt);
  4037. printf("stRange: %f\n",p->a.stRange);
  4038. printf("wndSecs: %f (%i)\n",p->a.wndSecs,p->hN);
  4039. printf("minTrkSec: %f (%i)\n",p->a.minTrkSec,p->minTrkN);
  4040. printf("maxTrkDeadSec: %f (%i)\n",p->a.maxTrkDeadSec,p->deadN_max);
  4041. printf("pkThreshDb: %f\n",p->a.pkThreshDb);
  4042. printf("pkAtkThreshDb: %f\n",p->a.pkAtkThreshDb);
  4043. }
  4044. //------------------------------------------------------------------------------------------------------------
  4045. cmFbCtl_t* cmFbCtlAlloc( cmCtx* c, cmFbCtl_t* ap, const cmFbCtlArgs_t* a )
  4046. {
  4047. cmFbCtl_t* p = cmObjAlloc( cmFbCtl_t, c, ap );
  4048. p->sva = cmVectArrayAlloc(c,kRealVaFl);
  4049. p->uva = cmVectArrayAlloc(c,kRealVaFl);
  4050. if( a != NULL )
  4051. {
  4052. if( cmFbCtlInit( p, a ) != cmOkRC )
  4053. cmFbCtlFree(&p);
  4054. }
  4055. return p;
  4056. }
  4057. cmRC_t cmFbCtlFree( cmFbCtl_t** pp )
  4058. {
  4059. if( pp == NULL || *pp == NULL )
  4060. return cmOkRC;
  4061. cmFbCtl_t* p = *pp;
  4062. cmVectArrayWrite(p->sva, "/home/kevin/temp/frqtrk/fb_ctl_s.va");
  4063. cmVectArrayWrite(p->uva, "/home/kevin/temp/frqtrk/fb_ctl_u.va");
  4064. cmMemFree(p->bM);
  4065. cmMemFree(p->rmsV);
  4066. cmVectArrayFree(&p->sva);
  4067. cmVectArrayFree(&p->uva);
  4068. cmObjFree(pp);
  4069. return cmOkRC;
  4070. }
  4071. cmRC_t cmFbCtlInit( cmFbCtl_t* p, const cmFbCtlArgs_t* a )
  4072. {
  4073. cmRC_t rc;
  4074. if((rc = cmFbCtlFinal(p)) != cmOkRC )
  4075. return rc;
  4076. double binHz = a->srate / ((a->binCnt-1)*2);
  4077. p->a = *a;
  4078. p->frmCnt = (a->bufMs * a->srate / 1000.0) /a->hopSmpCnt;
  4079. p->binCnt = cmMin(p->a.binCnt, a->maxHz/binHz);
  4080. p->bM = cmMemResizeZ(cmReal_t, p->bM, p->binCnt*p->frmCnt);
  4081. p->rmsV = cmMemResizeZ(cmReal_t, p->rmsV, p->frmCnt);
  4082. p->sV = cmMemResizeZ(cmReal_t, p->sV, p->binCnt);
  4083. p->uV = cmMemResizeZ(cmReal_t, p->uV, p->binCnt);
  4084. printf("cmFbCtl: frmCnt:%i binCnt:%i \n",p->frmCnt,p->binCnt);
  4085. return rc;
  4086. }
  4087. cmRC_t cmFbCtlFinal(cmFbCtl_t* p )
  4088. { return cmOkRC; }
  4089. cmRC_t cmFbCtlExec( cmFbCtl_t* p, const cmReal_t* x0V )
  4090. {
  4091. unsigned i;
  4092. cmRC_t rc = cmOkRC;
  4093. cmReal_t xV[ p->binCnt ];
  4094. cmVOR_AmplToDbVV(xV, p->binCnt, x0V, -1000.0 );
  4095. cmVOR_Shift( p->rmsV, p->frmCnt, -1, 0 );
  4096. p->rmsV[0] = cmVOR_Mean(xV,p->binCnt);
  4097. cmVOR_CopyN(p->bM + p->bfi, p->binCnt, p->frmCnt, xV, 1 );
  4098. p->bfi = (p->bfi + 1) % p->frmCnt;
  4099. p->bfN = cmMin(p->bfN+1,p->frmCnt);
  4100. for(i=0; i<p->binCnt; ++i)
  4101. {
  4102. const cmReal_t* v = p->bM + i * p->frmCnt;
  4103. cmReal_t u = cmVOR_Mean(v, p->bfN );
  4104. cmReal_t s = sqrt(cmVOR_Variance(v, p->bfN,&u));
  4105. p->sV[i] = (0.0002 - s);
  4106. p->uV[i] = u;
  4107. }
  4108. cmVectArrayAppendR(p->sva,p->sV,p->binCnt);
  4109. cmVectArrayAppendR(p->uva,p->uV,p->binCnt);
  4110. return rc;
  4111. }
  4112. //=======================================================================================================================
  4113. cmExpander* cmExpanderAlloc( cmCtx* c, cmExpander* p,
  4114. double srate, unsigned procSmpCnt, double threshDb, double rlsDb,
  4115. double threshMs, double rmsMs, double atkMs, double rlsMs )
  4116. {
  4117. cmExpander* op = cmObjAlloc(cmExpander,c,p);
  4118. if( srate > 0 )
  4119. if( cmExpanderInit(op,srate, procSmpCnt, threshDb, rlsDb, threshMs, rmsMs, atkMs, rlsMs) != cmOkRC )
  4120. cmExpanderFree(&op);
  4121. return op;
  4122. }
  4123. cmRC_t cmExpanderFree( cmExpander** pp )
  4124. {
  4125. cmRC_t rc = cmOkRC;
  4126. if( pp==NULL || *pp==NULL )
  4127. return rc;
  4128. cmExpander* p = *pp;
  4129. if((rc = cmExpanderFinal(p)) != cmOkRC )
  4130. return rc;
  4131. cmMemFree(p->rmsV);
  4132. cmMemFree(p->envV);
  4133. cmObjFree(pp);
  4134. return rc;
  4135. }
  4136. cmRC_t cmExpanderInit( cmExpander* p,
  4137. double srate, unsigned procSmpCnt, double threshDb, double rlsDb,
  4138. double threshMs, double rmsMs, double atkMs, double rlsMs )
  4139. {
  4140. cmRC_t rc;
  4141. unsigned i;
  4142. if((rc = cmExpanderFinal(p)) != cmOkRC )
  4143. return rc;
  4144. unsigned atkN = cmMax(1,ceil( atkMs / (srate * 1000.0)));
  4145. unsigned rlsN = cmMax(1,ceil( rlsMs / (srate * 1000.0)));
  4146. p->rmsN = cmMax(1,ceil(rmsMs / (srate * 1000.0)));
  4147. p->rmsV = cmMemResizeZ(cmReal_t,p->rmsV,p->rmsN);
  4148. p->rmsIdx = 0;
  4149. p->envN = atkN + rlsN;
  4150. p->envV = cmMemResizeZ(cmSample_t,p->envV,p->envN);
  4151. p->envIdx = p->envN;
  4152. p->threshN = cmMax(1,ceil(threshMs / (srate * 1000.0)));
  4153. p->threshIdx = 0;
  4154. p->threshLvl = pow(10.0,(threshDb/20.0));
  4155. p->rlsLvl = pow(10.0,(rlsDb/20.0));
  4156. p->gain = 1.0;
  4157. p->atkCnt = 0;
  4158. cmSample_t G = (1.0 - p->rlsLvl);
  4159. for(i=0; i<atkN; ++i)
  4160. {
  4161. p->envV[i] = 1.0 - G*i/atkN;
  4162. }
  4163. for(i=0; i<rlsN; ++i)
  4164. {
  4165. p->envV[atkN+i] = p->rlsLvl + (G*i/rlsN);
  4166. }
  4167. printf("rmsN:%i atkN:%i rlsN:%i thr:%f %f rls:%f %f\n",p->rmsN,atkN,rlsN,threshDb,p->threshLvl,rlsDb,p->rlsLvl);
  4168. //for(i=0; i<p->envN; ++i)
  4169. // printf("%i %f\n",i,p->envV[i]);
  4170. //printf("\n");
  4171. return cmOkRC;
  4172. }
  4173. cmRC_t cmExpanderFinal( cmExpander* p )
  4174. { return cmOkRC; }
  4175. cmRC_t cmExpanderExec( cmExpander* p, cmSample_t* x, cmSample_t* y, unsigned xyN )
  4176. {
  4177. unsigned i;
  4178. // update the RMS buffer
  4179. for(i=0; i<xyN; ++i)
  4180. {
  4181. // NOTE: using abs() instead of pow(x,2)
  4182. p->rmsV[p->rmsIdx] = abs(x[i]);
  4183. if( ++p->rmsIdx >= p->rmsN )
  4184. p->rmsIdx = 0;
  4185. }
  4186. // calculate the RMS
  4187. double rms = cmVOR_Mean(p->rmsV,p->rmsN);
  4188. // update the duration that the signal has been above the threshold
  4189. if( rms > p->threshLvl )
  4190. p->threshIdx += 1;
  4191. else
  4192. p->threshIdx = 0;
  4193. // begin the atk phase?
  4194. if( p->threshIdx > p->threshN && p->envIdx >= p->envN )
  4195. {
  4196. p->envIdx = 0;
  4197. }
  4198. // update the output
  4199. if( p->envIdx >= p->envN )
  4200. {
  4201. if( y != NULL )
  4202. cmVOS_Copy(y,xyN,x);
  4203. }
  4204. else
  4205. {
  4206. if( y == NULL )
  4207. y = x;
  4208. for(i=0; i<xyN && p->envIdx<p->envN; ++i,++p->envIdx)
  4209. y[i] = p->envV[p->envIdx] * x[i];
  4210. }
  4211. return cmOkRC;
  4212. }
  4213. cmRC_t cmExpanderExecD( cmExpander* p, double* x, double* y, unsigned xyN )
  4214. {
  4215. unsigned i;
  4216. // update the RMS buffer
  4217. for(i=0; i<xyN; ++i)
  4218. {
  4219. // NOTE: using abs() instead of pow(x,2)
  4220. p->rmsV[p->rmsIdx] = x[i];
  4221. p->rmsIdx += 1;
  4222. if( p->rmsIdx >= p->rmsN )
  4223. p->rmsIdx = 0;
  4224. }
  4225. // calculate the RMS
  4226. p->rmsValue = cmVOR_Mean(p->rmsV,p->rmsN);
  4227. // update the duration that the signal has been above the threshold
  4228. if( p->rmsValue > p->threshLvl )
  4229. p->threshIdx += 1;
  4230. else
  4231. p->threshIdx = 0;
  4232. // begin the atk phase?
  4233. if( p->threshIdx > p->threshN && p->envIdx >= p->envN )
  4234. {
  4235. p->envIdx = 0;
  4236. p->atkCnt += 1;
  4237. }
  4238. /*
  4239. if( p->envIdx >= p->envN )
  4240. p->gain = 1.0;
  4241. else
  4242. {
  4243. p->gain = p->envV[p->envIdx];
  4244. p->envIdx += 1;
  4245. }
  4246. */
  4247. // update the output
  4248. if( p->envIdx >= p->envN )
  4249. {
  4250. if( y != NULL )
  4251. cmVOD_Copy(y,xyN,x);
  4252. }
  4253. else
  4254. {
  4255. if( y == NULL )
  4256. y = x;
  4257. for(i=0; i<xyN && p->envIdx<p->envN; ++i,++p->envIdx)
  4258. y[i] = p->envV[p->envIdx] * x[i];
  4259. }
  4260. return cmOkRC;
  4261. }
  4262. //=======================================================================================================================
  4263. cmExpanderBank* cmExpanderBankAlloc( cmCtx* c, cmExpanderBank* p, unsigned bandN, double srate, unsigned procSmpCnt, double threshDb, double rlsDb, double threshMs, double rmsMs, double atkMs, double rlsMs )
  4264. {
  4265. cmExpanderBank* op = cmObjAlloc(cmExpanderBank,c,p);
  4266. if( bandN > 0 )
  4267. if( cmExpanderBankInit(op,bandN,srate, procSmpCnt, threshDb, rlsDb, threshMs, rmsMs, atkMs, rlsMs) != cmOkRC )
  4268. cmExpanderBankFree(&op);
  4269. return op;
  4270. }
  4271. cmRC_t cmExpanderBankFree( cmExpanderBank** pp )
  4272. {
  4273. cmRC_t rc = cmOkRC;
  4274. if( pp==NULL || *pp==NULL )
  4275. return rc;
  4276. cmExpanderBank* p = *pp;
  4277. if((rc = cmExpanderBankFinal(p)) != cmOkRC )
  4278. return rc;
  4279. cmMemFree(p->b);
  4280. cmObjFree(pp);
  4281. return rc;
  4282. }
  4283. cmRC_t cmExpanderBankInit( cmExpanderBank* p, unsigned bandN, double srate, unsigned procSmpCnt, double threshDb, double rlsDb, double threshMs, double rmsMs, double atkMs, double rlsMs )
  4284. {
  4285. cmRC_t rc;
  4286. unsigned i;
  4287. if((rc = cmExpanderBankFinal(p)) != cmOkRC )
  4288. return rc;
  4289. p->bandN = bandN;
  4290. p->b = cmMemResizeZ(cmExpander*,p->b,p->bandN);
  4291. for(i=0; i<bandN; ++i)
  4292. p->b[i] = cmExpanderAlloc(p->obj.ctx,NULL,srate, procSmpCnt,threshDb,rlsDb,threshMs,rmsMs,atkMs,rlsMs);
  4293. return cmOkRC;
  4294. }
  4295. cmRC_t cmExpanderBankFinal( cmExpanderBank* p )
  4296. {
  4297. unsigned i;
  4298. for(i=0; i<p->bandN; ++i)
  4299. cmExpanderFree(&p->b[i]);
  4300. return cmOkRC;
  4301. }
  4302. cmRC_t cmExpanderBankExec( cmExpanderBank* p, cmSample_t* x, unsigned bandN )
  4303. {
  4304. assert( bandN <= p->bandN);
  4305. unsigned i;
  4306. for(i=0; i<bandN; ++i)
  4307. {
  4308. cmExpanderExec(p->b[i],&x[i],NULL,1);
  4309. }
  4310. return cmOkRC;
  4311. }
  4312. /*
  4313. cmRC_t cmExpanderBankExecD( cmExpanderBank* p, double* x, unsigned binN )
  4314. {
  4315. unsigned i;
  4316. p->rmsValue = 0;
  4317. for(i=0; i<p->bandN; ++i)
  4318. {
  4319. double sum = cmVOD_Sum(x,binN);
  4320. cmExpanderExecD(p->b[i],&sum,NULL,1);
  4321. //printf("%f %f\n",sum, p->b[i]->rmsValue);
  4322. p->rmsValue += p->b[i]->rmsValue;
  4323. cmVOR_MultVS(x,binN,p->b[i]->gain);
  4324. }
  4325. p->rmsValue /= p->bandN;
  4326. return cmOkRC;
  4327. }
  4328. */
  4329. cmRC_t cmExpanderBankExecD( cmExpanderBank* p, double* x, unsigned bandN )
  4330. {
  4331. unsigned i;
  4332. //unsigned n = 3;
  4333. //unsigned no2 = n/2;
  4334. double xx;
  4335. p->rmsValue = 0;
  4336. p->atkCnt = 0;
  4337. for(i=0; i<p->bandN; ++i)
  4338. {
  4339. unsigned atkCnt = p->b[i]->atkCnt;
  4340. //if( i >= no2 && i < bandN-no2 )
  4341. // xx = cmVOR_Mean(x-no2,n);
  4342. //else
  4343. xx = x[i];
  4344. cmExpanderExecD(p->b[i],&xx,NULL,1);
  4345. p->rmsValue += p->b[i]->rmsValue;
  4346. p->atkCnt += p->b[i]->atkCnt != atkCnt;
  4347. }
  4348. p->rmsValue /= p->bandN;
  4349. return cmOkRC;
  4350. }
  4351. //------------------------------------------------------------------------------------------------------------
  4352. cmSpecDist_t* cmSpecDistAlloc( cmCtx* ctx,cmSpecDist_t* ap, unsigned procSmpCnt, double srate, unsigned wndSmpCnt, unsigned hopFcmt, unsigned olaWndTypeId )
  4353. {
  4354. cmSpecDist_t* p = cmObjAlloc( cmSpecDist_t, ctx, ap );
  4355. //p->iSpecVa = cmVectArrayAlloc(ctx,kRealVaFl);
  4356. p->oSpecVa = cmVectArrayAlloc(ctx,kRealVaFl);
  4357. if( procSmpCnt != 0 )
  4358. {
  4359. if( cmSpecDistInit( p, procSmpCnt, srate, wndSmpCnt, hopFcmt, olaWndTypeId ) != cmOkRC )
  4360. cmSpecDistFree(&p);
  4361. }
  4362. return p;
  4363. }
  4364. cmRC_t cmSpecDistFree( cmSpecDist_t** pp )
  4365. {
  4366. if( pp == NULL || *pp == NULL )
  4367. return cmOkRC;
  4368. cmSpecDist_t* p = *pp;
  4369. cmSpecDistFinal(p);
  4370. //cmVectArrayFree(&p->iSpecVa);
  4371. cmVectArrayFree(&p->oSpecVa);
  4372. cmMemPtrFree(&p->hzV);
  4373. cmMemPtrFree(&p->iSpecM);
  4374. cmMemPtrFree(&p->oSpecM);
  4375. cmMemPtrFree(&p->iSpecV);
  4376. cmMemPtrFree(&p->oSpecV);
  4377. cmObjFree(pp);
  4378. return cmOkRC;
  4379. }
  4380. cmRC_t cmSpecDistInit( cmSpecDist_t* p, unsigned procSmpCnt, double srate, unsigned wndSmpCnt, unsigned hopFcmt, unsigned olaWndTypeId )
  4381. {
  4382. //cmFrqTrkArgs_t fta;
  4383. cmRC_t rc;
  4384. if((rc = cmSpecDistFinal(p)) != cmOkRC )
  4385. return rc;
  4386. unsigned flags = 0;
  4387. p->srate = srate;
  4388. p->wndSmpCnt = wndSmpCnt;
  4389. p->hopSmpCnt = (unsigned)floor(wndSmpCnt/hopFcmt);
  4390. p->procSmpCnt = procSmpCnt;
  4391. p->mode = kBasicModeSdId;
  4392. p->thresh = 60;
  4393. p->offset = 0;
  4394. p->invertFl = false;
  4395. p->uprSlope = 0.0;
  4396. p->lwrSlope = 2.0;
  4397. p->pva = cmPvAnlAlloc( p->obj.ctx, NULL, procSmpCnt, srate, wndSmpCnt, p->hopSmpCnt, flags );
  4398. p->pvs = cmPvSynAlloc( p->obj.ctx, NULL, procSmpCnt, srate, wndSmpCnt, p->hopSmpCnt, olaWndTypeId );
  4399. /*
  4400. fta.srate = srate;
  4401. fta.chCnt = 50;
  4402. fta.binCnt = p->pva->binCnt;
  4403. fta.hopSmpCnt = p->pva->hopSmpCnt;
  4404. fta.stRange = 0.25;
  4405. fta.wndSecs = 0.25;
  4406. fta.minTrkSec = 0.25;
  4407. fta.maxTrkDeadSec = 0.25;
  4408. fta.pkThreshDb = 0.1; //-110.0;
  4409. fta.pkAtkThreshDb = 0.4; //-60.0;
  4410. fta.pkMaxHz = 20000;
  4411. fta.whFiltCoeff = 0.33;
  4412. fta.attenThresh = 0.4;
  4413. fta.attenGain = 0.5;
  4414. fta.attenDlySec = 1.0;
  4415. fta.attenAtkSec = 1.0;
  4416. fta.logFn = "/home/kevin/temp/frqtrk/trk_log.va";
  4417. fta.levelFn = "/home/kevin/temp/frqtrk/level.va";
  4418. fta.specFn = "/home/kevin/temp/frqtrk/spec.va";
  4419. fta.attenFn = "/home/kevin/temp/frqtrk/atten.va";
  4420. */
  4421. //p->ft = cmFrqTrkAlloc( p->obj.ctx, NULL, &fta );
  4422. /*
  4423. cmFbCtlArgs_t fba;
  4424. fba.srate = srate;
  4425. fba.binCnt = p->pva->binCnt;
  4426. fba.hopSmpCnt = p->hopSmpCnt;
  4427. fba.bufMs = 500;
  4428. fba.maxHz = 5000;
  4429. */
  4430. //p->fbc = cmFbCtlAlloc( p->obj.ctx, NULL, &fba );
  4431. //unsigned expBandN = 1; //
  4432. unsigned expBandN = 20000.0 / (p->srate / p->pva->binCnt);
  4433. double expSrate = p->pva->hopSmpCnt / srate;
  4434. unsigned expProcSmpCnt = 1;
  4435. double expThreshDb = -80.0;
  4436. double expRlsDb = -18.0;
  4437. double expThreshMs = 250.0;
  4438. double expRmsMs = 100.0;
  4439. double expAtkMs = 25.0;
  4440. double expRlsMs = 1000.0;
  4441. p->exb = cmExpanderBankAlloc( p->obj.ctx, NULL, expBandN, expSrate, expProcSmpCnt, expThreshDb, expRlsDb, expThreshMs, expRmsMs, expAtkMs, expRlsMs );
  4442. p->spcBwHz = cmMin(srate/2,10000);
  4443. p->spcSmArg = 0.05;
  4444. p->spcMin = p->spcBwHz;
  4445. p->spcMax = 0.0;
  4446. p->spcSum = 0.0;
  4447. p->spcCnt = 0;
  4448. double binHz = srate / p->pva->wndSmpCnt;
  4449. p->spcBinCnt = (unsigned)floor(p->spcBwHz / binHz);
  4450. p->hzV = cmMemResizeZ(cmReal_t,p->hzV,p->spcBinCnt);
  4451. cmVOR_Seq( p->hzV, p->spcBinCnt, 0, 1 );
  4452. cmVOR_MultVS( p->hzV, p->spcBinCnt, binHz );
  4453. p->aeUnit = 0;
  4454. p->aeMin = 1000;
  4455. p->aeMax = -1000;
  4456. double histSecs = 0.05;
  4457. p->hN = cmMax(1,histSecs * p->srate / p->hopSmpCnt );
  4458. p->iSpecM = cmMemResizeZ(cmReal_t,p->iSpecM,p->hN*p->pva->binCnt);
  4459. p->oSpecM = cmMemResizeZ(cmReal_t,p->oSpecM,p->hN*p->pva->binCnt);
  4460. p->iSpecV = cmMemResizeZ(cmReal_t,p->iSpecV, p->pva->binCnt);
  4461. p->oSpecV = cmMemResizeZ(cmReal_t,p->oSpecV, p->pva->binCnt);
  4462. p->hi = 0;
  4463. //p->bypOut = cmMemResizeZ(cmSample_t, p->bypOut, procSmpCnt );
  4464. return rc;
  4465. }
  4466. cmRC_t cmSpecDistFinal(cmSpecDist_t* p )
  4467. {
  4468. cmRC_t rc = cmOkRC;
  4469. //cmVectArrayWrite(p->iSpecVa, "/home/kevin/temp/frqtrk/iSpec.va");
  4470. cmVectArrayWrite(p->oSpecVa, "/home/kevin/temp/expand/oSpec.va");
  4471. cmPvAnlFree(&p->pva);
  4472. cmPvSynFree(&p->pvs);
  4473. //cmFrqTrkFree(&p->ft);
  4474. //cmFbCtlFree(&p->fbc);
  4475. cmExpanderBankFree(&p->exb);
  4476. return rc;
  4477. }
  4478. void _cmSpecDistBasicMode0(cmSpecDist_t* p, cmReal_t* X1m, unsigned binCnt, cmReal_t thresh )
  4479. {
  4480. // octavez> thresh = 60;
  4481. // octave> X1m = [-62 -61 -60 -59];
  4482. // octave> -abs(abs(X1m+thresh)-(X1m+thresh)) - thresh
  4483. // octave> ans = -64 -62 -60 -60
  4484. unsigned i=0;
  4485. for(i=0; i<binCnt; ++i)
  4486. {
  4487. cmReal_t a = fabs(X1m[i]);
  4488. cmReal_t d = a - thresh;
  4489. X1m[i] = -thresh;
  4490. if( d > 0 )
  4491. X1m[i] -= 2*d;
  4492. }
  4493. }
  4494. void _cmSpecDistBasicMode(cmSpecDist_t* p, cmReal_t* X1m, unsigned binCnt, cmReal_t thresh )
  4495. {
  4496. unsigned i=0;
  4497. if( p->lwrSlope < 0.3 )
  4498. p->lwrSlope = 0.3;
  4499. for(i=0; i<binCnt; ++i)
  4500. {
  4501. cmReal_t a = fabs(X1m[i]);
  4502. cmReal_t d = a - thresh;
  4503. X1m[i] = -thresh;
  4504. if( d > 0 )
  4505. X1m[i] -= (p->lwrSlope*d);
  4506. else
  4507. X1m[i] -= (p->uprSlope*d);
  4508. }
  4509. }
  4510. cmReal_t _cmSpecDistCentMode( cmSpecDist_t* p, cmReal_t* X1m )
  4511. {
  4512. // calc the spectral centroid
  4513. double num = cmVOR_MultSumVV( p->pva->magV, p->hzV, p->spcBinCnt );
  4514. double den = cmVOR_Sum( p->pva->magV, p->spcBinCnt );
  4515. double result = 0;
  4516. if( den != 0 )
  4517. result = num/den;
  4518. // apply smoothing filter to spectral centroid
  4519. p->spc = (result * p->spcSmArg) + (p->spc * (1.0-p->spcSmArg));
  4520. // track spec. cetr. min and max
  4521. p->spcMin = cmMin(p->spcMin,p->spc);
  4522. p->spcMax = cmMax(p->spcMax,p->spc);
  4523. //-----------------------------------------------------
  4524. ++p->spcCnt;
  4525. p->spcSum += p->spc;
  4526. p->spcSqSum += p->spc * p->spc;
  4527. // use the one-pass std-dev calc. trick
  4528. //double mean = p->spcSum / p->spcCnt;
  4529. //double variance = p->spcSqSum / p->spcCnt - mean * mean;
  4530. //double std_dev = sqrt(variance);
  4531. double smin = p->spcMin;
  4532. double smax = p->spcMax;
  4533. //smin = mean - std_dev;
  4534. //smax = mean + std_dev;
  4535. //-----------------------------------------------------
  4536. // convert spec. cent. to unit range
  4537. double spcUnit = (p->spc - smin) / (smax - smin);
  4538. spcUnit = cmMin(1.0,cmMax(0.0,spcUnit));
  4539. if( p->invertFl )
  4540. spcUnit = 1.0 - spcUnit;
  4541. //if( p->spcMin==p->spc || p->spcMax==p->spc )
  4542. // printf("min:%f avg:%f sd:%f max:%f\n",p->spcMin,p->spcSum/p->spcCnt,std_dev,p->spcMax);
  4543. return spcUnit;
  4544. }
  4545. void _cmSpecDistBump( cmSpecDist_t* p, cmReal_t* x, unsigned binCnt, double thresh)
  4546. {
  4547. /*
  4548. thresh *= -1;
  4549. minDb = -100;
  4550. if db < minDb
  4551. db = minDb;
  4552. endif
  4553. if db > thresh
  4554. y = 1;
  4555. else
  4556. x = (minDb - db)/(minDb - thresh);
  4557. y = x + (x - (x.^coeff));
  4558. endif
  4559. y = minDb + abs(minDb) * y;
  4560. */
  4561. unsigned i=0;
  4562. //printf("%f %f %f\n",thresh,p->lwrSlope,x[0]);
  4563. double minDb = -100.0;
  4564. thresh = -thresh;
  4565. for(i=0; i<binCnt; ++i)
  4566. {
  4567. double y;
  4568. if( x[i] < minDb )
  4569. x[i] = minDb;
  4570. if( x[i] > thresh )
  4571. y = 1;
  4572. else
  4573. {
  4574. y = (minDb - x[i])/(minDb - thresh);
  4575. y += y - pow(y,p->lwrSlope);
  4576. }
  4577. x[i] = minDb + (-minDb) * y;
  4578. }
  4579. }
  4580. void _cmSpecDistAmpEnvMode( cmSpecDist_t* p, cmReal_t* X1m )
  4581. {
  4582. cmReal_t smCoeff = 0.1;
  4583. //
  4584. cmReal_t mx = cmVOR_Max(X1m,p->pva->binCnt,1);
  4585. p->aeSmMax = (mx * smCoeff) + (p->aeSmMax * (1.0-smCoeff));
  4586. cmReal_t a = cmVOR_Mean(X1m,p->pva->binCnt);
  4587. p->ae = (a * smCoeff) + (p->ae * (1.0-smCoeff));
  4588. p->aeMin = cmMin(p->ae,p->aeMin);
  4589. p->aeMax = cmMax(p->ae,p->aeMax);
  4590. p->aeUnit = (p->ae - p->aeMin) / (p->aeMax-p->aeMin);
  4591. p->aeUnit = cmMin(1.0,cmMax(0.0,p->aeUnit));
  4592. if( p->invertFl )
  4593. p->aeUnit = 1.0 - p->aeUnit;
  4594. //printf("%f\n",p->aeSmMax);
  4595. }
  4596. void _cmSpecDistPhaseMod( cmSpecDist_t* p, cmReal_t* phsV, unsigned binCnt )
  4597. {
  4598. unsigned i;
  4599. cmReal_t offs = sin( 0.1 * 2.0 * M_PI * (p->phaseModIndex++) / (p->srate/p->hopSmpCnt) );
  4600. //printf("offs %f %i %i %f\n",offs,p->phaseModIndex,p->hopSmpCnt,p->srate);
  4601. cmReal_t new_phs = phsV[0] + offs;
  4602. for(i=0; i<binCnt-1; ++i)
  4603. {
  4604. while( new_phs > M_PI )
  4605. new_phs -= 2.0*M_PI;
  4606. while( new_phs < -M_PI )
  4607. new_phs += 2.0*M_PI;
  4608. cmReal_t d = phsV[i+1] - phsV[i];
  4609. phsV[i] = new_phs;
  4610. new_phs += d;
  4611. }
  4612. }
  4613. cmRC_t cmSpecDistExec( cmSpecDist_t* p, const cmSample_t* sp, unsigned sn )
  4614. {
  4615. assert( sn == p->procSmpCnt );
  4616. bool recordFl = false;
  4617. // cmPvAnlExec() returns true when it calc's a new spectral output frame
  4618. if( cmPvAnlExec( p->pva, sp, sn ) )
  4619. {
  4620. cmReal_t X1m[p->pva->binCnt];
  4621. // take the mean of the the input magntitude spectrum
  4622. cmReal_t u0 = cmVOR_Mean(p->pva->magV,p->pva->binCnt);
  4623. if(recordFl)
  4624. {
  4625. // store a time windowed average of the input spectrum to p->iSpecV
  4626. cmVOR_CopyN(p->iSpecM + p->hi, p->pva->binCnt, p->hN, X1m, 1 );
  4627. cmVOR_MeanM2(p->iSpecV, p->iSpecM, p->hN, p->pva->binCnt, 0, cmMin(p->fi+1,p->hN));
  4628. }
  4629. cmVOR_AmplToDbVV(X1m, p->pva->binCnt, p->pva->magV, -1000.0 );
  4630. //cmVOR_AmplToDbVV(X1m, p->pva->binCnt, X1m, -1000.0 );
  4631. switch( p->mode )
  4632. {
  4633. case kBypassModeSdId:
  4634. break;
  4635. case kBasicModeSdId:
  4636. _cmSpecDistBasicMode(p,X1m,p->pva->binCnt,p->thresh);
  4637. break;
  4638. case kSpecCentSdId:
  4639. {
  4640. _cmSpecDistAmpEnvMode(p,X1m);
  4641. double spcUnit = _cmSpecDistCentMode(p,X1m);
  4642. double thresh = fabs(p->aeSmMax) - (spcUnit*p->offset);
  4643. _cmSpecDistBasicMode(p,X1m,p->pva->binCnt, thresh);
  4644. }
  4645. break;
  4646. case kAmpEnvSdId:
  4647. {
  4648. _cmSpecDistAmpEnvMode(p,X1m);
  4649. //double thresh = fabs(p->aeSmMax) - p->offset;
  4650. double thresh = fabs(p->aeSmMax) - (p->aeUnit*p->offset);
  4651. thresh = fabs(p->thresh) - (p->aeUnit*p->offset);
  4652. _cmSpecDistBasicMode(p,X1m,p->pva->binCnt, thresh);
  4653. }
  4654. break;
  4655. case kBumpSdId:
  4656. _cmSpecDistBump(p,X1m, p->pva->binCnt, p->offset);
  4657. _cmSpecDistBasicMode(p,X1m,p->pva->binCnt,p->thresh);
  4658. break;
  4659. case 5:
  4660. break;
  4661. default:
  4662. break;
  4663. }
  4664. cmVOR_DbToAmplVV(X1m, p->pva->binCnt, X1m );
  4665. // run and apply the tracker/supressor
  4666. //cmFrqTrkExec(p->ft, X1m, p->pva->phsV, NULL );
  4667. //cmVOR_MultVV(X1m, p->pva->binCnt,p->ft->aV );
  4668. // convert the mean input magnitude to db
  4669. cmReal_t idb = 20*log10(u0);
  4670. // get the mean output magnitude spectra
  4671. cmReal_t u1 = cmVOR_Mean(X1m,p->pva->binCnt);
  4672. if( idb > -150.0 )
  4673. {
  4674. // set the output gain such that the mean output magnitude
  4675. // will match the mean input magnitude
  4676. p->ogain = u0/u1;
  4677. }
  4678. else
  4679. {
  4680. cmReal_t a0 = 0.9;
  4681. p->ogain *= a0;
  4682. }
  4683. cmVOR_MultVS(X1m,p->pva->binCnt,cmMin(4.0,p->ogain));
  4684. //cmFbCtlExec(p->fbc,X1m);
  4685. //cmReal_t v[ p->pva->binCnt ];
  4686. //cmVOR_Copy(v,p->pva->binCnt,p->pva->phsV);
  4687. //_cmSpecDistPhaseMod(p, v, p->pva->binCnt );
  4688. if(recordFl)
  4689. {
  4690. // store a time windowed average of the output spectrum to p->iSpecV
  4691. cmVOR_CopyN(p->oSpecM + p->hi, p->pva->binCnt, p->hN, X1m, 1 );
  4692. cmVOR_MeanM2(p->oSpecV, p->oSpecM, p->hN, p->pva->binCnt, 0, cmMin(p->fi+1,p->hN));
  4693. // store iSpecV and oSpecV to iSpecVa and oSpecVa to create debugging files
  4694. //cmVectArrayAppendR(p->iSpecVa,p->iSpecV,p->pva->binCnt);
  4695. //cmVectArrayAppendR(p->oSpecVa,p->oSpecV,p->pva->binCnt);
  4696. p->hi = (p->hi + 1) % p->hN;
  4697. }
  4698. //unsigned binN = 12500.0 / (p->srate / p->pva->binCnt);
  4699. //cmExpanderBankExecD(p->exb, X1m, binN );
  4700. /*
  4701. cmExpanderBankExecD(p->exb, X1m, p->exb->bandN );
  4702. cmReal_t mean = cmVOR_Mean(X1m,p->exb->bandN);
  4703. cmReal_t arr[3] = { p->exb->rmsValue, mean, p->exb->atkCnt };
  4704. cmVectArrayAppendR(p->oSpecVa,arr,3);
  4705. */
  4706. cmPvSynExec(p->pvs, X1m, p->pva->phsV );
  4707. p->fi += 1;
  4708. }
  4709. return cmOkRC;
  4710. }
  4711. const cmSample_t* cmSpecDistOut( cmSpecDist_t* p )
  4712. {
  4713. return cmPvSynExecOut(p->pvs);
  4714. }
  4715. //------------------------------------------------------------------------------------------------------------
  4716. cmRC_t _cmBinMtxFileWriteHdr( cmBinMtxFile_t* p )
  4717. {
  4718. cmFileRC_t fileRC;
  4719. unsigned n = 3;
  4720. unsigned hdr[n];
  4721. hdr[0] = p->rowCnt;
  4722. hdr[1] = p->maxRowEleCnt;
  4723. hdr[2] = p->eleByteCnt;
  4724. if((fileRC = cmFileSeek(p->fh,kBeginFileFl,0)) != kOkFileRC )
  4725. return cmCtxRtCondition(&p->obj, fileRC, "File seek failed on matrix file:'%s'.", cmStringNullGuard(cmFileName(p->fh)));
  4726. if((fileRC = cmFileWriteUInt(p->fh,hdr,n)) != kOkFileRC )
  4727. return cmCtxRtCondition( &p->obj, fileRC, "Header write failed on matrix file:'%s'.", cmStringNullGuard(cmFileName(p->fh)) );
  4728. return cmOkRC;
  4729. }
  4730. cmBinMtxFile_t* cmBinMtxFileAlloc( cmCtx* ctx, cmBinMtxFile_t* ap, const cmChar_t* fn )
  4731. {
  4732. cmBinMtxFile_t* p = cmObjAlloc( cmBinMtxFile_t, ctx, ap );
  4733. if( fn != NULL )
  4734. if( cmBinMtxFileInit( p, fn ) != cmOkRC )
  4735. cmBinMtxFileFree(&p);
  4736. return p;
  4737. }
  4738. cmRC_t cmBinMtxFileFree( cmBinMtxFile_t** pp )
  4739. {
  4740. cmRC_t rc;
  4741. if( pp==NULL || *pp == NULL )
  4742. return cmOkRC;
  4743. cmBinMtxFile_t* p = *pp;
  4744. if((rc = cmBinMtxFileFinal(p)) == cmOkRC )
  4745. {
  4746. cmObjFree(pp);
  4747. }
  4748. return rc;
  4749. }
  4750. cmRC_t cmBinMtxFileInit( cmBinMtxFile_t* p, const cmChar_t* fn )
  4751. {
  4752. cmRC_t rc;
  4753. cmFileRC_t fileRC;
  4754. if((rc = cmBinMtxFileFinal(p)) != cmOkRC )
  4755. return rc;
  4756. // open the output file for writing
  4757. if((fileRC = cmFileOpen(&p->fh,fn,kWriteFileFl | kBinaryFileFl, p->obj.err.rpt)) != kOkFileRC )
  4758. return cmCtxRtCondition( &p->obj, fileRC, "Unable to open the matrix file:'%s'", cmStringNullGuard(fn) );
  4759. // iniitlaize the object
  4760. p->rowCnt = 0;
  4761. p->maxRowEleCnt = 0;
  4762. p->eleByteCnt = 0;
  4763. // write the blank header as place holder
  4764. if((rc = _cmBinMtxFileWriteHdr(p)) != cmOkRC )
  4765. return rc;
  4766. return rc;
  4767. }
  4768. cmRC_t cmBinMtxFileFinal( cmBinMtxFile_t* p )
  4769. {
  4770. cmRC_t rc;
  4771. cmFileRC_t fileRC;
  4772. if( p != NULL && cmFileIsValid(p->fh))
  4773. {
  4774. // re-write the file header
  4775. if((rc = _cmBinMtxFileWriteHdr(p)) != cmOkRC )
  4776. return rc;
  4777. // close the file
  4778. if((fileRC = cmFileClose(&p->fh)) != kOkFileRC )
  4779. return cmCtxRtCondition(&p->obj, fileRC, "Matrix file close failed on:'%s'",cmStringNullGuard(cmFileName(p->fh)));
  4780. }
  4781. return cmOkRC;
  4782. }
  4783. bool cmBinMtxFileIsValid( cmBinMtxFile_t* p )
  4784. { return p != NULL && cmFileIsValid(p->fh); }
  4785. cmRC_t _cmBinMtxFileWriteRow( cmBinMtxFile_t* p, const void* buf, unsigned eleCnt, unsigned eleByteCnt )
  4786. {
  4787. cmFileRC_t fileRC;
  4788. if((fileRC = cmFileWrite(p->fh,&eleCnt,sizeof(eleCnt))) != kOkFileRC )
  4789. return cmCtxRtCondition(&p->obj, fileRC, "Matrix file row at index %i element count write failed on '%s'.", p->rowCnt, cmStringNullGuard(cmFileName(p->fh)));
  4790. if((fileRC = cmFileWrite(p->fh,buf,eleCnt*eleByteCnt)) != kOkFileRC )
  4791. return cmCtxRtCondition(&p->obj, fileRC, "Matrix file row at index %i data write failed on '%s'.", p->rowCnt, cmStringNullGuard(cmFileName(p->fh)));
  4792. if( eleCnt > p->maxRowEleCnt )
  4793. p->maxRowEleCnt = eleCnt;
  4794. ++p->rowCnt;
  4795. return cmOkRC;
  4796. }
  4797. cmRC_t cmBinMtxFileExecS( cmBinMtxFile_t* p, const cmSample_t* x, unsigned xn )
  4798. {
  4799. // verify that all rows are written as cmSample_t
  4800. assert( p->eleByteCnt == 0 || p->eleByteCnt == sizeof(cmSample_t));
  4801. p->eleByteCnt = sizeof(cmSample_t);
  4802. return _cmBinMtxFileWriteRow(p,x,xn,p->eleByteCnt);
  4803. }
  4804. cmRC_t cmBinMtxFileExecR( cmBinMtxFile_t* p, const cmReal_t* x, unsigned xn )
  4805. {
  4806. // verify that all rows are written as cmReal_t
  4807. assert( p->eleByteCnt == 0 || p->eleByteCnt == sizeof(cmReal_t));
  4808. p->eleByteCnt = sizeof(cmReal_t);
  4809. return _cmBinMtxFileWriteRow(p,x,xn,p->eleByteCnt);
  4810. }
  4811. cmRC_t cmBinMtxFileWrite( const cmChar_t* fn, unsigned rowCnt, unsigned colCnt, const cmSample_t* sp, const cmReal_t* rp, cmCtx* ctx, cmRpt_t* rpt )
  4812. {
  4813. assert( sp == NULL || rp == NULL );
  4814. cmCtx* ctxp = NULL;
  4815. cmBinMtxFile_t* bp = NULL;
  4816. if( ctx == NULL )
  4817. ctx = ctxp = cmCtxAlloc(NULL,rpt,cmLHeapNullHandle,cmSymTblNullHandle);
  4818. if((bp = cmBinMtxFileAlloc(ctx,NULL,fn)) != NULL )
  4819. {
  4820. unsigned i = 0;
  4821. cmSample_t* sbp = sp == NULL ? NULL : cmMemAlloc(cmSample_t,colCnt);
  4822. cmReal_t* rbp = rp == NULL ? NULL : cmMemAlloc(cmReal_t,colCnt);
  4823. for(i=0; i<rowCnt; ++i)
  4824. {
  4825. if( sp!=NULL )
  4826. {
  4827. cmVOS_CopyN(sbp,colCnt,1,sp+i,rowCnt);
  4828. cmBinMtxFileExecS(bp,sbp,colCnt);
  4829. }
  4830. if( rp!=NULL )
  4831. {
  4832. cmVOR_CopyN(rbp,colCnt,1,rp+i,rowCnt);
  4833. cmBinMtxFileExecR(bp,rbp,colCnt);
  4834. }
  4835. }
  4836. cmMemPtrFree(&sbp);
  4837. cmMemPtrFree(&rbp);
  4838. cmBinMtxFileFree(&bp);
  4839. }
  4840. if( ctxp != NULL )
  4841. cmCtxFree(&ctxp);
  4842. return cmOkRC;
  4843. }
  4844. cmRC_t _cmBinMtxFileReadHdr( cmCtx_t* ctx, cmFileH_t h, unsigned* rowCntPtr, unsigned* colCntPtr, unsigned* eleByteCntPtr, const cmChar_t* fn )
  4845. {
  4846. cmRC_t rc = cmOkRC;
  4847. unsigned hdr[3];
  4848. if( cmFileRead(h,&hdr,sizeof(hdr)) != kOkFileRC )
  4849. {
  4850. rc = cmErrMsg(&ctx->err,cmSubSysFailRC,"Binary matrix file header read failed on '%s'.",cmStringNullGuard(fn));
  4851. goto errLabel;
  4852. }
  4853. if( rowCntPtr != NULL )
  4854. *rowCntPtr = hdr[0];
  4855. if( colCntPtr != NULL )
  4856. *colCntPtr = hdr[1];
  4857. if( eleByteCntPtr != NULL )
  4858. *eleByteCntPtr = hdr[2];
  4859. errLabel:
  4860. return rc;
  4861. }
  4862. cmRC_t cmBinMtxFileSize( cmCtx_t* ctx, const cmChar_t* fn, unsigned* rowCntPtr, unsigned* colCntPtr, unsigned* eleByteCntPtr )
  4863. {
  4864. cmFileH_t h = cmFileNullHandle;
  4865. cmRC_t rc = cmOkRC;
  4866. if(cmFileOpen(&h,fn,kReadFileFl | kBinaryFileFl, ctx->err.rpt) != kOkFileRC )
  4867. {
  4868. rc = cmErrMsg(&ctx->err,cmSubSysFailRC,"Binary matrix file:%s open failed.",cmStringNullGuard(fn));
  4869. goto errLabel;
  4870. }
  4871. rc = _cmBinMtxFileReadHdr(ctx,h,rowCntPtr,colCntPtr,eleByteCntPtr,fn);
  4872. errLabel:
  4873. cmFileClose(&h);
  4874. return rc;
  4875. }
  4876. cmRC_t cmBinMtxFileRead( cmCtx_t* ctx, const cmChar_t* fn, unsigned mRowCnt, unsigned mColCnt, unsigned mEleByteCnt, void* retBuf, unsigned* colCntV )
  4877. {
  4878. cmFileH_t h = cmFileNullHandle;
  4879. cmRC_t rc = cmOkRC;
  4880. char* rowBuf = NULL;
  4881. unsigned rowCnt,colCnt,eleByteCnt,i;
  4882. cmErr_t err;
  4883. cmErrSetup(&err,ctx->err.rpt,"Binary Matrix File Reader");
  4884. if(cmFileOpen(&h,fn,kReadFileFl | kBinaryFileFl, err.rpt) != kOkFileRC )
  4885. {
  4886. rc = cmErrMsg(&err,cmSubSysFailRC,"Binary matrix file:%s open failed.",cmStringNullGuard(fn));
  4887. goto errLabel;
  4888. }
  4889. if((rc = _cmBinMtxFileReadHdr(ctx,h,&rowCnt,&colCnt,&eleByteCnt,fn)) != cmOkRC )
  4890. goto errLabel;
  4891. if( mRowCnt != rowCnt )
  4892. rc = cmErrMsg(&err,cmArgAssertRC,"The binary matrix file row count and the return buffer row count are not the same.");
  4893. if( mColCnt != colCnt )
  4894. rc = cmErrMsg(&err,cmArgAssertRC,"The binary matrix file column count and the return buffer column count are not the same.");
  4895. if( mEleByteCnt != eleByteCnt )
  4896. rc = cmErrMsg(&err,cmArgAssertRC,"The binary matrix file element byte count and the return buffer element byte count are not the same.");
  4897. if( rc == cmOkRC )
  4898. {
  4899. rowBuf = cmMemAllocZ(char,colCnt*eleByteCnt);
  4900. for(i=0; i<rowCnt; ++i)
  4901. {
  4902. unsigned cn;
  4903. // read the row length
  4904. if( cmFileReadUInt(h,&cn,1) != kOkFileRC )
  4905. {
  4906. rc = cmErrMsg(&err,cmSubSysFailRC,"Error reading row length at row index:%i.",i);
  4907. goto errLabel;
  4908. }
  4909. if( colCntV != NULL )
  4910. colCntV[i] = cn;
  4911. // verify the actual col count does not exceed the max col count
  4912. if( cn > colCnt )
  4913. {
  4914. rc = cmErrMsg(&err,cmSubSysFailRC,"The actual column count:%i exceeds the max column count:%i.",cn,colCnt);
  4915. goto errLabel;
  4916. }
  4917. //read the row data
  4918. if( cmFileReadChar(h,rowBuf,cn*eleByteCnt) != kOkFileRC )
  4919. {
  4920. rc = cmErrMsg(&err,cmSubSysFailRC,"File read failed at row index:%i.",i);
  4921. goto errLabel;
  4922. }
  4923. char* dp = ((char*)retBuf) + i * eleByteCnt;
  4924. // the data is read in row-major order but the matrix must be
  4925. // returned on col major order - rearrange the columns here.
  4926. switch(eleByteCnt)
  4927. {
  4928. case sizeof(cmSample_t):
  4929. cmVOS_CopyN(((cmSample_t*)dp), cn, rowCnt, (cmSample_t*)rowBuf, 1 );
  4930. break;
  4931. case sizeof(cmReal_t):
  4932. cmVOR_CopyN(((cmReal_t*)dp), cn, rowCnt, (cmReal_t*)rowBuf, 1 );
  4933. break;
  4934. default:
  4935. rc = cmErrMsg(&err,cmSubSysFailRC,"Invalid element byte count:%i.",eleByteCnt);
  4936. goto errLabel;
  4937. }
  4938. }
  4939. }
  4940. errLabel:
  4941. cmMemPtrFree(&rowBuf);
  4942. cmFileClose(&h);
  4943. return rc;
  4944. }