comparison hs.cpp @ 1:6422640a802f

first upload
author Wen X <xue.wen@elec.qmul.ac.uk>
date Tue, 05 Oct 2010 10:45:57 +0100
parents
children fc19d45615d1
comparison
equal deleted inserted replaced
0:9b9f21935f24 1:6422640a802f
1 //---------------------------------------------------------------------------
2 #include <math.h>
3 #include "hs.h"
4 #include "align8.h"
5 #include "arrayalloc.h"
6 #include "procedures.h"
7 #include "SinEst.h"
8 #include "SinSyn.h"
9 #include "splines.h"
10 #include "xcomplex.h"
11 #ifdef I
12 #undef I
13 #endif
14
15 //---------------------------------------------------------------------------
16 //stiffcandid methods
17
18 /*
19 method stiffcandid::stiffcandid: creates a harmonic atom with minimal info, i.e. fundamantal frequency
20 between 0 and Nyqvist frequency, stiffness coefficient between 0 and STIFF_B_MAX.
21
22 In: Wid: DFT size (frame size, frequency resolution)
23 */
24 stiffcandid::stiffcandid(int Wid)
25 {
26 F=new double[6];
27 G=&F[3];
28 double Fm=Wid*Wid/4; //NyqvistF^2
29 F[0]=0, G[0]=0;
30 F[1]=Fm, G[1]=STIFF_B_MAX*Fm;
31 F[2]=Fm, G[2]=0;
32 N=3;
33 P=0;
34 s=0;
35 p=NULL;
36 f=NULL;
37 }//stiffcandid
38
39 /*
40 method stiffcandid::stiffcandid: creates a harmonic atom with a frequency range for the $ap-th partial,
41 without actually taking this partial as "known".
42
43 In; Wid: DFT size
44 ap: partial index
45 af, errf: centre and half-width of the frequency range of the ap-th partial, in bins
46 */
47 stiffcandid::stiffcandid(int Wid, int ap, double af, double errf)
48 {
49 F=new double[10];
50 G=&F[5];
51 double Fm=Wid*Wid/4;
52 F[0]=0, G[0]=0;
53 F[1]=Fm, G[1]=STIFF_B_MAX*Fm;
54 F[2]=Fm, G[2]=0;
55 N=3;
56 P=0;
57 s=0;
58 p=NULL;
59 f=NULL;
60 double k=ap*ap-1;
61 double g1=(af-errf)/ap, g2=(af+errf)/ap;
62 if (g1<0) g1=0;
63 g1=g1*g1, g2=g2*g2;
64 //apply constraints F+kG-g2<0 and -F-kG+g1<0
65 cutcvpoly(N, F, G, 1, k, -g2);
66 cutcvpoly(N, F, G, -1, -k, g1);
67 }//stiffcandid
68
69 /*
70 method stiffcandid::stiffcandid: creates a harmonic atom from one atom.
71
72 In; Wid: DFT size
73 ap: partial index of the atom.
74 af, errf: centre and half-width of the frequency range of the atom, in bins.
75 ds: contribution to candidate score from this atom.
76 */
77 stiffcandid::stiffcandid(int Wid, int ap, double af, double errf, double ds)
78 {
79 //the starting polygon implements the trivial constraints 0<f0<NyqvistF,
80 // 0<B<Bmax. in terms of F and G, 0<F<NF^2, 0<G<FBmax. this is a triangle
81 // with vertices (0, 0), (NF^2, Bmax*NF^2), (NF^2, 0)
82 F=new double[12];
83 G=&F[5], f=&F[10], p=(int*)&F[11];
84 double Fm=Wid*Wid/4; //NyqvistF^2
85 F[0]=0, G[0]=0;
86 F[1]=Fm, G[1]=STIFF_B_MAX*Fm;
87 F[2]=Fm, G[2]=0;
88
89 N=3;
90 P=1;
91 p[0]=ap;
92 f[0]=af;
93 s=ds;
94
95 double k=ap*ap-1;
96 double g1=(af-errf)/ap, g2=(af+errf)/ap;
97 if (g1<0) g1=0;
98 g1=g1*g1, g2=g2*g2;
99 //apply constraints F+kG-g2<0 and -F-kG+g1<0
100 cutcvpoly(N, F, G, 1, k, -g2);
101 cutcvpoly(N, F, G, -1, -k, g1);
102 }//stiffcandid
103
104 /*
105 method stiffcandid::stiffcandid: creates a harmonic atom by adding a new partial to an existing
106 harmonic atom.
107
108 In; prev: initial harmonic atom
109 ap: partial index of new atom
110 af, errf: centre and half-width of the frequency range of the new atom, in bins.
111 ds: contribution to candidate score from this atom.
112 */
113 stiffcandid::stiffcandid(stiffcandid* prev, int ap, double af, double errf, double ds)
114 {
115 N=prev->N, P=prev->P, s=prev->s;
116 int Np2=N+2, Pp1=P+1;
117 F=new double[(Np2+Pp1)*2];
118 G=&F[Np2]; f=&G[Np2]; p=(int*)&f[Pp1];
119 memcpy(F, prev->F, sizeof(double)*N);
120 memcpy(G, prev->G, sizeof(double)*N);
121 if (P>0)
122 {
123 memcpy(f, prev->f, sizeof(double)*P);
124 memcpy(p, prev->p, sizeof(int)*P);
125 }
126 //see if partial ap is already involved
127 int i=0; while (i<P && p[i]!=ap) i++;
128 //if it is not:
129 if (i>=P)
130 {
131 p[P]=ap;
132 f[P]=af;
133 s+=ds;
134 P++;
135
136 double k=ap*ap-1;
137 double g1=(af-errf)/ap, g2=(af+errf)/ap;
138 if (g1<0) g1=0;
139 g1=g1*g1, g2=g2*g2;
140 //apply constraints F+kG-g2<0 and -F-kG+g1<0
141 cutcvpoly(N, F, G, 1, k, -g2);
142 cutcvpoly(N, F, G, -1, -k, g1);
143 }
144 //otherwise, the new candid is simply a copy of prev
145 }//stiffcandid
146
147 //method stiffcandid::~stiffcandid: destructor
148 stiffcandid::~stiffcandid()
149 {
150 delete[] F;
151 }//~stiffcandid
152
153
154 //---------------------------------------------------------------------------
155 //THS methods
156
157 //method THS::THS: default constructor
158 THS::THS()
159 {
160 memset(this, 0, sizeof(THS));
161 memset(BufM, 0, sizeof(void*)*128);
162 }//THS
163
164 /*
165 method THS::THS: creates an HS with given number of partials and frames.
166
167 In: aM, aFr: number of partials and frames
168 */
169 THS::THS(int aM, int aFr)
170 {
171 memset(this, 0, sizeof(THS));
172 M=aM; Fr=aFr; Allocate2(atom, M, Fr, Partials);
173 memset(Partials[0], 0, sizeof(atom)*M*Fr);
174 memset(BufM, 0, sizeof(void*)*128);
175 }//THS
176
177 /*
178 method THS::THS: creates a duplicate HS. This does not duplicate contents of BufM.
179
180 In: HS: the HS to duplicate
181 */
182 THS::THS(THS* HS)
183 {
184 memcpy(this, HS, sizeof(THS));
185 Allocate2(atom, M, Fr, Partials);
186 for (int m=0; m<M; m++) memcpy(Partials[m], HS->Partials[m], sizeof(atom)*Fr);
187 Allocate2(double, M, st_count, startamp);
188 if (st_count) for (int m=0; m<M; m++) memcpy(startamp[m], HS->startamp[m], sizeof(double)*st_count);
189 memset(BufM, 0, sizeof(void*)*128);
190 }//THS
191
192 /*
193 method THS::THS: create a new HS which is a trunction of an existing HS.
194
195 In: HS: the exinting HS from which a sub-HS is to be created
196 start, end: truncation points. Atoms of *HS beyond these points are discarded.
197 */
198 THS::THS(THS* HS, double start, double end)
199 {
200 memset(this, 0, sizeof(THS));
201 M=HS->M; Fr=HS->Fr; isconstf=HS->isconstf;
202 int frst=0, fren=Fr;
203 while (HS->Partials[0][frst].t<start && frst<fren) frst++;
204 while (HS->Partials[0][fren-1].t>end && fren>frst) fren--;
205 Fr=fren-frst;
206 Allocate2(atom, M, Fr, Partials);
207 for (int m=0; m<M; m++)
208 {
209 memcpy(Partials[m], &HS->Partials[m][frst], sizeof(atom)*Fr);
210 for (int fr=0; fr<Fr; fr++) Partials[m][fr].t-=start;
211 }
212 if (frst>0){DeAlloc2(startamp); startamp=0;}
213 }//THS
214
215 //method THS::~THS: default descructor.
216 THS::~THS()
217 {
218 DeAlloc2(Partials);
219 DeAlloc2(startamp);
220 ClearBufM();
221 }//~THS
222
223 /*
224 method THS::ClearBufM: free buffers pointed to by members of BufM.
225 */
226 void THS::ClearBufM()
227 {
228 for (int m=0; m<128; m++) delete[] BufM[m];
229 memset(BufM, 0, sizeof(void*)*128);
230 }//ClearBufM
231
232 /*
233 method THS::EndPos: the end position of an HS.
234
235 Returns the time of the last frame of the first partial.
236 */
237 int THS::EndPos()
238 {
239 return Partials[0][Fr-1].t;
240 }//EndPos
241
242 /*
243 method THS::EndPosEx: the extended end position of an HS.
244
245 Returns the time later than the last frame of the first partial by the interval between the first and
246 second frames of that partial.
247 */
248 int THS::EndPosEx()
249 {
250 return Partials[0][Fr-1].t+Partials[0][1].t-Partials[0][0].t;
251 }//EndPosEx
252
253 /*
254 method THS::ReadAtomFromStream: reads an atom from a stream.
255
256 Returns the number of bytes read, or 0 on failure. The stream pointer is shifted by this value.
257 */
258 int THS::ReadAtomFromStream(TStream* Stream, atom* atm)
259 {
260 int StartPosition=Stream->Position, atomlen; char s[4];
261 if (Stream->Read(s, 4)!=4) goto ret0;
262 if (strcmp(s, "ATM")) goto ret0;
263 if (Stream->Read(&atomlen, sizeof(int))!=sizeof(int)) goto ret0;
264 if (Stream->Read(atm, atomlen)!=atomlen) goto ret0;
265 return atomlen+4+sizeof(int);
266 ret0:
267 Stream->Position=StartPosition;
268 return 0;
269 }//ReadAtomFromStream
270
271 /*
272 method THS::ReadFromStream: reads an HS from a stream.
273
274 Returns the number of bytes read, or 0 on failure. The stream pointer is shifted by this value.
275 */
276 int THS::ReadFromStream(TStream* Stream)
277 {
278 int StartPosition=Stream->Position, lenevt; char s[4];
279 if (Stream->Read(s, 4)!=4) goto ret0;
280 if (strcmp(s, "EVT")) goto ret0;
281 if (Stream->Read(&lenevt, sizeof(int))!=sizeof(int)) goto ret0;
282 if (!ReadHdrFromStream(Stream)) goto ret0;
283 Resize(M, Fr, false, true);
284 for (int m=0; m<M; m++) for (int fr=0; fr<Fr; fr++) if (!ReadAtomFromStream(Stream, &Partials[m][fr])) goto ret0;
285 lenevt=lenevt+4+sizeof(int);
286 if (lenevt!=Stream->Position-StartPosition) goto ret0;
287 return lenevt;
288 ret0:
289 Stream->Position=StartPosition;
290 return 0;
291 }//ReadFromStream
292
293 /*
294 method THS::ReadHdrFromStream: reads header from a stream into this HS.
295
296 Returns the number of bytes read, or 0 on failure. The stream pointer is shifted by this value.
297 */
298 int THS::ReadHdrFromStream(TStream* Stream)
299 {
300 int StartPosition=Stream->Position, hdrlen, tags; char s[4];
301 if (Stream->Read(s, 4)!=4) goto ret0;
302 if (strcmp(s, "HDR")) goto ret0;
303 if (Stream->Read(&hdrlen, sizeof(int))!=sizeof(int)) goto ret0;
304 if (Stream->Read(&Channel, sizeof(int))!=sizeof(int)) goto ret0;
305 if (Stream->Read(&M, sizeof(int))!=sizeof(int)) goto ret0;
306 if (Stream->Read(&Fr, sizeof(int))!=sizeof(int)) goto ret0;
307 if (Stream->Read(&tags, sizeof(int))!=sizeof(int)) goto ret0;
308 isconstf=tags&HS_CONSTF;
309 hdrlen=hdrlen+4+sizeof(int); //expected read count
310 if (hdrlen!=Stream->Position-StartPosition) goto ret0;
311 return hdrlen;
312 ret0:
313 Stream->Position=StartPosition;
314 return 0;
315 }//ReadHdrFromStream
316
317 /*
318 method THS::Resize: change the number of partials or frames of an HS structure. Unless forcealloc is
319 set to true, this allocates new buffers to host HS data only when the new dimensions are larger than
320 current ones.
321
322 In: aM, aFr: new number of partials and frames
323 copydata: specifies if HS data is to be copied to new buffers.
324 forcealloc: specifies if new buffers are to be allocated in all cases.
325 */
326 void THS::Resize(int aM, int aFr, bool copydata, bool forcealloc)
327 {
328 if (forcealloc || M<aM || Fr<aFr)
329 {
330 atom** Allocate2(atom, aM, aFr, NewPartials);
331 if (copydata)
332 {
333 int minM=(M<aM)?M:aM, minFr=(Fr<aFr)?Fr:aFr;
334 for (int m=0; m<minM; m++) memcpy(NewPartials[m], Partials[m], sizeof(atom)*minFr);
335 }
336 DeAlloc2(Partials);
337 Partials=NewPartials;
338 }
339 if (forcealloc || M<aM)
340 {
341 double** Allocate2(double, aM, st_count, newstartamp);
342 if (copydata)
343 {
344 for (int m=0; m<M; m++) memcpy(newstartamp[m], startamp[m], sizeof(double)*st_count);
345 }
346 DeAlloc2(startamp);
347 startamp=newstartamp;
348 }
349 M=aM, Fr=aFr;
350 }//Resize
351
352 /*
353 method THS::StartPos: the start position of an HS.
354
355 Returns the time of the first frame of the first partial.
356 */
357 int THS::StartPos()
358 {
359 return Partials[0][0].t;
360 }//StartPos
361
362 /*
363 method THS::StdOffst: standard hop size of an HS.
364 Returns the interval between the timing of the first and second frames of the first partial, which is
365 the "standard" hop size if all measurement points are uniformly positioned.
366 */
367 int THS::StdOffst()
368 {
369 return Partials[0][1].t-Partials[0][0].t;
370 }//StdOffst
371
372 /*
373 method THS::WriteAtomToStream: writes an atom to a stream. An atom is stored in a stream as an RIFF
374 chunk identified by "ATM\0".
375
376 Returns the number of bytes written. The stream pointer is shifted by this value.
377 */
378 int THS::WriteAtomToStream(TStream* Stream, atom* atm)
379 {
380 Stream->Write((void*)"ATM", 4);
381 int atomlen=sizeof(atom);
382 Stream->Write(&atomlen, sizeof(int));
383 atomlen=Stream->Write(atm, atomlen);
384 return atomlen+4+sizeof(int);
385 }//WriteAtomToStream
386
387 /*
388 method THS::WriteHdrToStream: writes header to a stream. HS header is stored in a stream as an RIFF
389 chunk, identified by "HDR\0".
390
391 Returns the number of bytes written. The stream pointer is shifted by this value.
392 */
393 int THS::WriteHdrToStream(TStream* Stream)
394 {
395 Stream->Write((void*)"HDR", 4);
396 int hdrlen=0, tags=0;
397 if (isconstf) tags=tags|HS_CONSTF;
398 Stream->Write(&hdrlen, sizeof(int));
399 Stream->Write(&Channel, sizeof(int)); hdrlen+=sizeof(int);
400 Stream->Write(&M, sizeof(int)); hdrlen+=sizeof(int);
401 Stream->Write(&Fr, sizeof(int)); hdrlen+=sizeof(int);
402 Stream->Write(&tags, sizeof(int)); hdrlen+=sizeof(int);
403 Stream->Seek(-hdrlen-sizeof(int), soFromCurrent); Stream->Write(&hdrlen, sizeof(int));
404 Stream->Seek(hdrlen, soFromCurrent);
405 return hdrlen+4+sizeof(int);
406 }//WriteHdrToStream
407
408 /*
409 method THS::WriteToStream: write an HS to a stream. An HS is stored in a stream as an RIFF chunk
410 identified by "EVT\0". Its chunk data contains an HS header chunk followed by M*Fr atom chunks, in
411 rising order of m then fr.
412
413 Returns the number of bytes written. The stream pointer is shifted by this value.
414 */
415 int THS::WriteToStream(TStream* Stream)
416 {
417 Stream->Write((void*)"EVT", 4);
418 int lenevt=0;
419 Stream->Write(&lenevt, sizeof(int));
420 lenevt+=WriteHdrToStream(Stream);
421
422 for (int m=0; m<M; m++) for (int fr=0; fr<Fr; fr++) lenevt+=WriteAtomToStream(Stream, &Partials[m][fr]);
423 Stream->Seek(-lenevt-sizeof(int), soFromCurrent);
424 Stream->Write(&lenevt, sizeof(int));
425 Stream->Seek(lenevt, soFromCurrent);
426 return lenevt+4+sizeof(int);
427 }//WriteToStream
428
429
430 //---------------------------------------------------------------------------
431 //TPolygon methods
432
433 /*
434 method TPolygon::TPolygon: create an empty polygon with a storage capacity
435
436 In: cap: number of vertices to allocate memory for.
437 */
438 TPolygon::TPolygon(int cap)
439 {
440 X=(double*)malloc8(sizeof(double)*cap*2);
441 Y=&X[cap];
442 }//TPolygon
443
444 /*
445 method TPolygon::TPolygon: create a duplicate polygon.
446
447 In: cap: number of vertices to allocate memory for.
448 R: the polygon to duplicate.
449
450 If cap is smaller than the number of vertices of R, the latter is used for memory allocation.
451 */
452 TPolygon::TPolygon(int cap, TPolygon* R)
453 {
454 if (cap<R->N) cap=R->N; X=(double*)malloc8(sizeof(double)*cap*2); Y=&X[cap]; N=R->N;
455 memcpy(X, R->X, sizeof(double)*R->N); memcpy(Y, R->Y, sizeof(double)*R->N);
456 }//TPolygon
457
458 //method TPolygon::~TPolygon: default destructor.
459 TPolygon::~TPolygon()
460 {
461 free8(X);
462 }//~TPolygon
463
464 //---------------------------------------------------------------------------
465 //TTempAtom methods
466
467 /*
468 method TTempAtom::TTempAtom: initializes an empty atom with a fundamental frequency range.
469
470 In: af, ef: centre and half width of the fundamental frequency
471 maxB: upper bound of stiffness coefficient
472 */
473 TTempAtom::TTempAtom(double af, double ef, double maxB)
474 {
475 Prev=NULL; pind=f=a=s=acce=0;
476 R=new TPolygon(4);
477 InitializeR(R, af, ef, maxB);
478 }//TTempAtom
479
480 /*
481 method TTempAtom::TTempAtom: initialize an empty atom with a frequency range for a certain partial.
482
483 In: apin: partial index
484 af, ef: centre and half width of the fundamental frequency
485 maxB: upper bound of stiffness coefficient
486 */
487 TTempAtom::TTempAtom(int apin, double af, double ef, double maxB)
488 {
489 Prev=NULL; pind=f=a=s=acce=0;
490 R=new TPolygon(4);
491 InitializeR(R, apin, af, ef, maxB);
492 }//TTempAtom
493
494 /*
495 method TTempAtom::TTempAtom: initialize an empty atom whose F-G polygon is the extension of AR
496
497 In: AR: the F-G polygon to extend
498 delf1, delf2: amount of extension, in semitones, of fundamental frequency
499 minf: minimal fundamental frequency: extension will not go below this value
500 */
501 TTempAtom::TTempAtom(TPolygon* AR, double delf1, double delf2, double minf)
502 {
503 Prev=NULL; pind=f=a=s=acce=0;
504 R=new TPolygon(AR->N+2);
505 R->N=AR->N;
506 memcpy(R->X, AR->X, sizeof(double)*AR->N); memcpy(R->Y, AR->Y, sizeof(double)*AR->N);
507 ExtendR(R, delf1, delf2, minf);
508 }//TTempAtom
509
510 /*
511 method TTempAtom::TTempAtom: initialize a atom as the only partial.
512
513 In: apind: partial index
514 af, ef: partial frequency and its error range
515 aa, as: partial amplitude and measurement scale
516 maxB: stiffness coefficient upper bound
517 */
518 TTempAtom::TTempAtom(int apind, double af, double ef, double aa, double as, double maxB)
519 {
520 Prev=NULL; pind=apind; f=af; a=aa; s=as; acce=aa*aa;
521 R=new TPolygon(4);
522 InitializeR(R, apind, af, ef, maxB);
523 }//TTempAtom
524
525 /*
526 method TTempAtom::TTempAtom: creates a duplicate atom. R can be duplicated or snatched according to
527 DupR.
528
529 In: APrev: the atom to duplicate
530 DupR: specifies if the newly created atom is to have its F-G polygon duplicated from that of APrev
531 or to snatch the F-G polygon from APrev (in the latter case APrev losts its polygon).
532 */
533 TTempAtom::TTempAtom(TTempAtom* APrev, bool DupR)
534 {
535 Prev=APrev; pind=APrev->pind; f=APrev->f; a=APrev->a; s=APrev->s; acce=APrev->acce;
536 TPolygon* PrevR=APrev->R;
537 if (DupR)
538 {//duplicate R
539 R=new TPolygon(PrevR->N);
540 R->N=PrevR->N; memcpy(R->X, PrevR->X, sizeof(double)*PrevR->N); memcpy(R->Y, PrevR->Y, sizeof(double)*PrevR->N);
541 }
542 else
543 {//snatch R from APrev
544 R=APrev->R;
545 APrev->R=0;
546 }
547 }//TTempAtom
548
549 /*
550 method TTempAtom::TTempAtom:: initializes an atom by adding a new partial to an existing group.
551
552 In: APrev: the existing atom.
553 apind: partial index
554 af, ef: partial frequency and its error range
555 aa, as: partial amplitude and measurement scale
556 updateR: specifies if the F-G polygon is updated using (af, ef).
557 */
558 TTempAtom::TTempAtom(TTempAtom* APrev, int apind, double af, double ef, double aa, double as, bool updateR)
559 {
560 Prev=APrev; pind=apind; f=af; a=aa;
561 s=as; acce=APrev->acce+aa*aa;
562 TPolygon* PrevR=APrev->R;
563 R=new TPolygon(PrevR->N+2); // create R
564 R->N=PrevR->N; //
565 memcpy(R->X, PrevR->X, sizeof(double)*PrevR->N); // copy R
566 memcpy(R->Y, PrevR->Y, sizeof(double)*PrevR->N); //
567 if (updateR) CutR(R, apind, af, ef);
568 }//TTempAtom
569
570 //method TTempAtom::~TTempAtom: default destructor
571 TTempAtom::~TTempAtom()
572 {
573 delete R;
574 }//~TTempAtom
575
576
577 //---------------------------------------------------------------------------
578 //functions
579
580 /*
581 function areaandcentroid: calculates the area and centroid of a convex polygon.
582
583 In: x[N], y[N]: x- and y-coordinates of vertices of a polygon
584 Out: A: area
585 (cx, cy): coordinate of the centroid
586
587 No return value.
588 */
589 void areaandcentroid(double& A, double& cx, double& cy, int N, double* x, double* y)
590 {
591 if (N==0) {A=0;}
592 else if (N==1) {A=0; cx=x[0], cy=y[0];}
593 else if (N==2) {A=0; cx=(x[0]+x[1])/2, cy=(y[0]+y[1])/2;}
594 A=cx=cy=0;
595 for (int i=0; i<N; i++)
596 {
597 double xi1, yi1; //x[i+1], y[i+1], index modular N
598 if (i+1==N) xi1=x[0], yi1=y[0];
599 else xi1=x[i+1], yi1=y[i+1];
600 double tmp=x[i]*yi1-xi1*y[i];
601 A+=tmp;
602 cx+=(x[i]+xi1)*tmp;
603 cy+=(y[i]+yi1)*tmp;
604 }
605 A/=2;
606 cx=cx/6/A;
607 cy=cy/6/A;
608 }//areaandcentroid
609
610 /*
611 function AtomsToPartials: sort a list of atoms as a number of partials. It is assumed that the atoms
612 are simply those from a harmonic sinusoid in arbitrary order with uniform timing and partial index
613 clearly marked, so that it is easy to sort them into a list of partials. However, the sorting process
614 does not check for harmonicity, missing or duplicate atoms, uniformity of timing, etc.
615
616 In: part[k]: a list of atoms
617 offst: interval between adjacent measurement points (hop size)
618 Out: M, Fr, partials[M][Fr]: a list of partials containing the atoms.
619
620 No return value. partials[][] is allocated anew and must be freed by caller.
621 */
622 void AtomsToPartials(int k, atom* part, int& M, int& Fr, atom**& partials, int offst)
623 {
624 if (k>0)
625 {
626 int t1, t2, i=0;
627 while (i<k && part[i].f<=0) i++;
628 if (i<k)
629 {
630 t1=part[i].t, t2=part[i].t, M=part[i].pin; i++;
631 while (i<k)
632 {
633 if (part[i].f>0)
634 {
635 if (M<part[i].pin) M=part[i].pin;
636 if (t1>part[i].t) t1=part[i].t;
637 if (t2<part[i].t) t2=part[i].t;
638 }
639 i++;
640 }
641 Fr=(t2-t1)/offst+1;
642 Allocate2(atom, M, Fr, partials); memset(partials[0], 0, sizeof(atom)*M*Fr);
643 for (int i=0; i<k; i++)
644 if (part[i].f>0)
645 {
646 int fr=(part[i].t-t1)/offst;
647 int pp=part[i].pin;
648 partials[pp-1][fr]=part[i];
649 }
650 }
651 else M=0, Fr=0;
652 }
653 else M=0, Fr=0;
654 }//AtomsToPartials
655
656 /*
657 function conta: a continuity measure between two set of amplitudes
658
659 In: a1[N], a2[N]: the amplitudes
660
661 Return the continuity measure, between 0 and 1
662 */
663 double conta(int N, double* a1, double* a2)
664 {
665 double e11=0, e22=0, e12=0;
666 for (int n=0; n<N; n++)
667 {
668 double la1=(a1[n]>0)?a1[n]:0, la2=(a2[n]>0)?a2[n]:0;
669 if (la1>1e10) la1=la2;
670 e11+=la1*la1;
671 e12+=la1*la2;
672 e22+=la2*la2;
673 }
674 return 2*e12/(e11+e22);
675 }//conta
676
677 /*
678 function cutcvpoly: severs a polygon by a straight line
679
680 In: x[N], y[N]: x- and y-coordinates of vertices of a polygon, starting from the leftmost (x[0]=
681 min x[n]) point clockwise; in case of two leftmost points, x[0] is the upper one and x[N-1] is
682 the lower one.
683 A, B, C: coefficients of a straight line Ax+By+C=0.
684 protect: specifies what to do if the whole polygon satisfy Ax+By+C>0, true=do noting,
685 false=eliminate.
686 Out: N, x[N], y[N]: the polygon severed by the line, retaining that half for which Ax+By+C<=0.
687
688 No return value.
689 */
690 void cutcvpoly(int& N, double* x, double* y, double A, double B, double C, bool protect)
691 {
692 int n=0, ns;
693 while (n<N && A*x[n]+B*y[n]+C<=0) n++;
694 if (n==N) //all points satisfies the constraint, do nothing
695 {}
696 else if (n>0) //points 0, 1, ..., n-1 satisfies the constraint, point n does not
697 {
698 ns=n;
699 n++;
700 while (n<N && A*x[n]+B*y[n]+C>0) n++;
701 //points 0, ..., ns-1 and n, ..., N-1 satisfy the constraint,
702 //points ns, ..., n-1 do not satisfy the constraint,
703 // ns>0, n>ns, however, it's possible that n-1=ns
704 int pts, pte; //indicates (by 0/1) whether or now a new point is to be added for xs/xe
705 double xs, ys, xe, ye;
706
707 //find intersection of x:y[ns-1:ns] and A:B:C
708 double x1=x[ns-1], x2=x[ns], y1=y[ns-1], y2=y[ns];
709 double z1=A*x1+B*y1+C, z2=A*x2+B*y2+C; //z1<=0, z2>0
710 if (z1==0) pts=0; //as point ns-1 IS point s
711 else
712 {
713 pts=1, xs=(x1*z2-x2*z1)/(z2-z1), ys=(y1*z2-y2*z1)/(z2-z1);
714 }
715
716 //find intersection of x:y[n-1:n] and A:B:C
717 x1=x[n-1], y1=y[n-1];
718 if (n==N) x2=x[0], y2=y[0];
719 else x2=x[n], y2=y[n];
720 z1=A*x1+B*y1+C, z2=A*x2+B*y2+C; //z1>0, z2<=0
721 if (z2==0) pte=0; //as point n is point e
722 else
723 {
724 pte=1, xe=(x1*z2-x2*z1)/(z2-z1), ye=(y1*z2-y2*z1)/(z2-z1);
725 }
726
727 //the new polygon is formed by points 0, 1, ..., ns-1, (s), (e), n, ..., N-1
728 memmove(&x[ns+pts+pte], &x[n], sizeof(double)*(N-n));
729 memmove(&y[ns+pts+pte], &y[n], sizeof(double)*(N-n));
730 if (pts) x[ns]=xs, y[ns]=ys;
731 if (pte) x[ns+pts]=xe, y[ns+pts]=ye;
732 N=ns+pts+pte+N-n;
733 }
734 else //n=0, point 0 does not satisfy the constraint
735 {
736 n++;
737 while (n<N && A*x[n]+B*y[n]+C>0) n++;
738 if (n==N) //no point satisfies the constraint
739 {
740 if (!protect) N=0;
741 }
742 else
743 {
744 //points 0, 1, ..., n-1 do not satisfy the constraint, point n does
745 ns=n;
746 n++;
747 while (n<N && A*x[n]+B*y[n]+C<=0) n++;
748 //points 0, ..., ns-1 and n, ..., N-1 do not satisfy constraint,
749 //points ns, ..., n-1 satisfy the constraint
750 // ns>0, n>ns, however ns can be equal to n-1
751 int pts, pte; //indicates (by 0/1) whether or now a new point is to be added for xs/xe
752 double xs, ys, xe, ye;
753
754 //find intersection of x:y[ns-1:ns] and A:B:C
755 double x1=x[ns-1], x2=x[ns], y1=y[ns-1], y2=y[ns];
756 double z1=A*x1+B*y1+C, z2=A*x2+B*y2+C; //z1>0, z2<=0
757 if (z2==0) pts=0; //as point ns IS point s
758 else pts=1;
759 xs=(x1*z2-x2*z1)/(z2-z1), ys=(y1*z2-y2*z1)/(z2-z1);
760
761 //find intersection of x:y[n-1:n] and A:B:C
762 x1=x[n-1], y1=y[n-1];
763 if (n==N) x2=x[0], y2=y[0];
764 else x2=x[n], y2=y[n];
765 z1=A*x1+B*y1+C, z2=A*x2+B*y2+C; //z1<=0, z2>0
766 if (z1==0) pte=0; //as point n-1 is point e
767 else pte=1;
768 xe=(x1*z2-x2*z1)/(z2-z1), ye=(y1*z2-y2*z1)/(z2-z1);
769
770 //the new polygon is formed of points (s), ns, ..., n-1, (e)
771 if (xs<=xe)
772 {
773 //the point sequence comes as (s), ns, ..., n-1, (e)
774 memmove(&x[pts], &x[ns], sizeof(double)*(n-ns));
775 memmove(&y[pts], &y[ns], sizeof(double)*(n-ns));
776 if (pts) x[0]=xs, y[0]=ys;
777 N=n-ns+pts+pte;
778 if (pte) x[N-1]=xe, y[N-1]=ye;
779 }
780 else //xe<xs
781 {
782 //the sequence comes as e, (s), ns, ..., n-2 if pte=0 (notice point e<->n-1)
783 //or e, (s), ns, ..., n-1 if pte=1
784 if (pte==0)
785 {
786 memmove(&x[pts+1], &x[ns], sizeof(double)*(n-ns-1));
787 memmove(&y[pts+1], &y[ns], sizeof(double)*(n-ns-1));
788 if (pts) x[1]=xs, y[1]=ys;
789 }
790 else
791 {
792 memmove(&x[pts+1], &x[ns], sizeof(double)*(n-ns));
793 memmove(&y[pts+1], &y[ns], sizeof(double)*(n-ns));
794 if (pts) x[1]=xs, y[1]=ys;
795 }
796 x[0]=xe, y[0]=ye;
797 N=n-ns+pts+pte;
798 }
799 }
800 }
801 }//cutcvpoly
802
803 /*
804 funciton CutR: update an F-G polygon with an new partial
805
806 In: R: F-G polygon
807 apind: partial index
808 af, ef: partial frequency and error bound
809 protect: specifies what to do if the new partial has no intersection with R, true=do noting,
810 false=eliminate R
811 Out: R: the updated F-G polygon
812
813 No return value.
814 */
815 void CutR(TPolygon* R, int apind, double af, double ef, bool protect)
816 {
817 double k=apind*apind-1;
818 double g1=(af-ef)/apind, g2=(af+ef)/apind;
819 if (g1<0) g1=0;
820 g1=g1*g1, g2=g2*g2;
821 //apply constraints F+kG-g2<0 and -F-kG+g1<0
822 cutcvpoly(R->N, R->X, R->Y, 1, k, -g2, protect); // update R
823 cutcvpoly(R->N, R->X, R->Y, -1, -k, g1, protect); //
824 }//CutR
825
826 /*
827 function DeleteByHalf: reduces a list of stiffcandid objects by half, retaining those with higher s
828 and delete the others, freeing their memory.
829
830 In: cands[pcs]: the list of stiffcandid objects
831 Out: cands[return value]: the list of retained ones
832
833 Returns the size of the new list.
834 */
835 int DeleteByHalf(stiffcandid** cands, int pcs)
836 {
837 int newpcs=pcs/2;
838 double* tmp=new double[newpcs];
839 memset(tmp, 0, sizeof(double)*newpcs);
840 for (int j=0; j<pcs; j++) InsertDec(cands[j]->s, tmp, newpcs);
841 int k=0;
842 for (int j=0; j<pcs; j++)
843 {
844 if (cands[j]->s>=tmp[newpcs-1]) cands[k++]=cands[j];
845 else delete cands[j];
846 }
847 return k;
848 }//DeleteByHalf
849
850 /*
851 function DeleteByHalf: reduces a list of TTempAtom objects by half, retaining those with higher s and
852 delete the others, freeing their memory.
853
854 In: cands[pcs]: the list of TTempAtom objects
855 Out: cands[return value]: the list of retained ones
856
857 Returns the size of the new list.
858 */
859 int DeleteByHalf(TTempAtom** cands, int pcs)
860 {
861 int newpcs=pcs/2;
862 double* tmp=new double[newpcs];
863 memset(tmp, 0, sizeof(double)*newpcs);
864 for (int j=0; j<pcs; j++) InsertDec(cands[j]->s, tmp, newpcs);
865 int k=0;
866 for (int j=0; j<pcs; j++)
867 {
868 if (cands[j]->s>=tmp[newpcs-1]) cands[k++]=cands[j];
869 else delete cands[j];
870 }
871 delete[] tmp;
872 return k;
873 }//DeleteByHalf
874
875 /*
876 function ds0: atom score 0 - energy
877
878 In: a: atom amplitude
879
880 Returns a^2, or s[0]+a^2 is s is not NULL, as atom score.
881 */
882 double ds0(double a, void* s)
883 {
884 if (s) return *(double*)s+a*a;
885 else return a*a;
886 }//ds0
887
888 /*
889 function ds2: atom score 1 - cross-correlation coefficient with another HA, e.g. from previous frame
890
891 In: a: atom amplitude
892 params: pointer to dsprams1 structure supplying other inputs.
893 Out: cross-correlation coefficient as atom score.
894
895 Returns the cross-correlation coefficient.
896 */
897 double ds1(double a, void* params)
898 {
899 dsparams1* lparams=(dsparams1*)params;
900 double hs=lparams->lastene+lparams->currentacce, hsaa=hs+a*a;
901 if (lparams->p<=lparams->lastP) return (lparams->s*hs+a*lparams->lastvfp[lparams->p-1])/hsaa;
902 else return (lparams->s*hs+a*a)/hsaa;
903 }//ds1
904
905 /*
906 function ExBStiff: finds the minimal and maximal values of stiffness coefficient B given a F-G polygon
907
908 In: F[N], G[N]: vertices of a F-G polygon
909 Out: Bmin, Bmax: minimal and maximal values of the stiffness coefficient
910
911 No reutrn value.
912 */
913 void ExBStiff(double& Bmin, double& Bmax, int N, double* F, double* G)
914 {
915 Bmax=G[0]/F[0];
916 Bmin=Bmax;
917 for (int i=1; i<N; i++)
918 {
919 double vi=G[i]/F[i];
920 if (Bmin>vi) Bmin=vi;
921 else if (Bmax<vi) Bmax=vi;
922 }
923 }//ExBStiff
924
925 /*
926 function ExFmStiff: finds the minimal and maximal frequecies of partial m given a F-G polygon
927
928 In: F[N], G[N]: vertices of a F-G polygon
929 m: partial index, fundamental=1
930 Out: Fmin, Fmax: minimal and maximal frequencies of partial m
931
932 No return value.
933 */
934 void ExFmStiff(double& Fmin, double& Fmax, int m, int N, double* F, double* G)
935 {
936 int k=m*m-1;
937 double vmax=F[0]+k*G[0];
938 double vmin=vmax;
939 for (int i=1; i<N; i++)
940 {
941 double vi=F[i]+k*G[i];
942 if (vmin>vi) vmin=vi;
943 else if (vmax<vi) vmax=vi;
944 }
945 testnn(vmin); Fmin=m*sqrt(vmin);
946 testnn(vmax); Fmax=m*sqrt(vmax);
947 }//ExFmStiff
948
949 /*
950 function ExtendR: extends a F-G polygon on the f axis (to allow a larger pitch range)
951
952 In: R: the F-G polygon
953 delp1: amount of extension to the low end, in semitones
954 delpe: amount of extension to the high end, in semitones
955 minf: minimal fundamental frequency
956 Out: R: the extended F-G polygon
957
958 No return value.
959 */
960 void ExtendR(TPolygon* R, double delp1, double delp2, double minf)
961 {
962 double mdelf1=pow(2, -delp1/12), mdelf2=pow(2, delp2/12);
963 int N=R->N;
964 double *G=R->Y, *F=R->X, slope, prevslope, f, ftmp;
965 if (N==0) {}
966 else if (N==1)
967 {
968 R->N=2;
969 slope=G[0]/F[0];
970 testnn(F[0]); f=sqrt(F[0]);
971 ftmp=f*mdelf2;
972 F[1]=ftmp*ftmp;
973 ftmp=f*mdelf1;
974 if (ftmp<minf) ftmp=minf;
975 F[0]=ftmp*ftmp;
976 G[0]=F[0]*slope;
977 G[1]=F[1]*slope;
978 }
979 else if (N==2)
980 {
981 prevslope=G[0]/F[0];
982 slope=G[1]/F[1];
983 if (fabs((slope-prevslope)/prevslope)<1e-10)
984 {
985 testnn(F[0]); ftmp=sqrt(F[0])*mdelf1; if (ftmp<0) ftmp=0; F[0]=ftmp*ftmp; G[0]=F[0]*prevslope;
986 testnn(F[1]); ftmp=sqrt(F[1])*mdelf2; F[1]=ftmp*ftmp; G[1]=F[1]*slope;
987 }
988 else
989 {
990 if (prevslope>slope)
991 {
992 R->N=4;
993 testnn(F[1]); f=sqrt(F[1]); ftmp=f*mdelf1; if (ftmp<minf) ftmp=minf; F[3]=ftmp*ftmp; G[3]=F[3]*slope;
994 ftmp=f*mdelf2; F[2]=ftmp*ftmp; G[2]=F[2]*slope;
995 testnn(F[0]); f=sqrt(F[0]); ftmp=f*mdelf1; if (ftmp<minf) ftmp=minf; F[0]=ftmp*ftmp; G[0]=F[0]*prevslope;
996 ftmp=f*mdelf2; F[1]=ftmp*ftmp; G[1]=F[1]*prevslope;
997 if (F[0]==0 && F[3]==0) R->N=3;
998 }
999 else
1000 {
1001 R->N=4;
1002 testnn(F[1]); f=sqrt(F[1]); ftmp=f*mdelf1; if (ftmp<minf) ftmp=minf; F[1]=ftmp*ftmp; G[1]=F[1]*slope;
1003 ftmp=f*mdelf2; F[2]=ftmp*ftmp; G[2]=F[2]*slope;
1004 testnn(F[0]); f=sqrt(F[0]); ftmp=f*mdelf1; if (ftmp<minf) ftmp=minf; F[0]=ftmp*ftmp; G[0]=F[0]*prevslope;
1005 ftmp=f*mdelf2; F[4]=ftmp*ftmp; G[4]=F[4]*prevslope;
1006 if (F[0]==0 && F[1]==0) {R->N=3; F[1]=F[2]; F[2]=F[3]; G[1]=G[2]; G[3]=G[3];}
1007 }
1008 }
1009 }
1010 else
1011 {
1012 //find the minimal and maximal angles of R
1013 //maximal slope is taken at Ms if Mt=1, or Ms-1 and Ms if Mt=0
1014 //minimal slope is taken at ms if mt=1, or ms-1 and ms if mt=0
1015 int Ms, ms, Mt=-1, mt=-1;
1016 double prevslope=G[0]/F[0], dslope;
1017 for (int n=1; n<N; n++)
1018 {
1019 double slope=G[n]/F[n];
1020 dslope=atan(slope)-atan(prevslope);
1021 if (Mt==-1)
1022 {
1023 if (dslope<-1e-10) Ms=n-1, Mt=1;
1024 else if (dslope<1e-10) Ms=n, Mt=0;
1025 }
1026 else if (mt==-1)
1027 {
1028 if (dslope>1e-10) ms=n-1, mt=1;
1029 else if (dslope>-1e-10) ms=n, mt=0;
1030 }
1031 prevslope=slope;
1032 }
1033 if (mt==-1)
1034 {
1035 slope=G[0]/F[0];
1036 dslope=atan(slope)-atan(prevslope);
1037 if (dslope>1e-10) ms=N-1, mt=1;
1038 else if (dslope>-1e-10) ms=0, mt=0;
1039 }
1040 R->N=N+mt+Mt;
1041 int newn=R->N-1;
1042 if (ms==0 && mt==1)
1043 {
1044 slope=G[0]/F[0];
1045 testnn(F[0]); ftmp=sqrt(F[0])*mdelf2; F[newn]=ftmp*ftmp; G[newn]=F[newn]*slope; newn--;
1046 mt=-1;
1047 }
1048 else if (ms==0)
1049 mt=-1;
1050 for (int n=N-1; n>=0; n--)
1051 {
1052 if (mt!=-1)
1053 {
1054 slope=G[n]/F[n];
1055 testnn(F[n]); f=sqrt(F[n]); ftmp=f*mdelf1; if (ftmp<minf) ftmp=minf; F[newn]=ftmp*ftmp; G[newn]=F[newn]*slope; newn--;
1056 if (ms==n && mt==1)
1057 {
1058 ftmp=f*mdelf2; F[newn]=ftmp*ftmp; G[newn]=F[newn]*slope; newn--;
1059 mt=-1;
1060 }
1061 else if (ms==n) //mt=0
1062 mt=-1;
1063 }
1064 else if (Mt!=-1)
1065 {
1066 slope=G[n]/F[n];
1067 testnn(F[n]); f=sqrt(F[n]); ftmp=f*mdelf2; F[newn]=ftmp*ftmp; G[newn]=F[newn]*slope; newn--;
1068 if (Ms==n && Mt==1)
1069 {
1070 ftmp=f*mdelf1; if (ftmp<minf) ftmp=minf; F[newn]=ftmp*ftmp; G[newn]=F[newn]*slope; newn--;
1071 Mt=-1;
1072 }
1073 else if (Ms==n)
1074 Mt=-1;
1075 }
1076 else
1077 {
1078 slope=G[n]/F[n];
1079 testnn(F[n]); f=sqrt(F[n]); ftmp=f*mdelf1; if (ftmp<minf) ftmp=minf; F[newn]=ftmp*ftmp; G[newn]=F[newn]*slope; newn--;
1080 }
1081 }
1082 //merge multiple vertices at the origin
1083 N=R->N;
1084 while (N>0 && F[N-1]==0) N--;
1085 R->N=N;
1086 int n=1;
1087 while (n<N && F[n]==0) n++;
1088 if (n!=1)
1089 {
1090 R->N=N-(n-1);
1091 memcpy(&F[1], &F[n], sizeof(double)*(N-n));
1092 memcpy(&G[1], &G[n], sizeof(double)*(N-n));
1093 }
1094 }
1095 }//ExtendR
1096
1097 /*
1098 function FGFrom2: solves the following equation system given m[2], f[2], norm[2]. This is interpreted
1099 as searching from (F0, G0) down direction (dF, dG) until the first equation is satisfied, where
1100 k[*]=m[*]^2-1.
1101
1102 m[0](F+k[0]G)^0.5-f[0] m[1](F+k[1]G)^0.5-f[1]
1103 ------------------------ = ------------------------ >0
1104 norm[0] norm[1]
1105
1106 F=F0+lmd*dF
1107 G=G0+lmd*dG
1108
1109 In: (F0, G0): search origin
1110 (dF, dG): search direction
1111 m[2], f[2], norm[2]: coefficients in the first equation
1112 Out: F[return value], G[return value]: solutions.
1113
1114 Returns the number of solutions for (F, G).
1115 */
1116 int FGFrom2(double* F, double* G, double F0, double dF, double G0, double dG, int* m, double* f, double* norm)
1117 {
1118 double m1=m[0]/norm[0], m2=m[1]/norm[1],
1119 f1=f[0]/norm[0], f2=f[1]/norm[1],
1120 k1=m[0]*m[0]-1, k2=m[1]*m[1]-1;
1121 double A=m1*m1-m2*m2*(dF+k2*dG)/(dF+k1*dG),
1122 B=2*m1*(f2-f1),
1123 C=(f2-f1)*(f2-f1)-(k1-k2)*(F0*dG-G0*dF)/(dF+k1*dG)*m2*m2;
1124 if (A==0)
1125 {
1126 double x=-C/B;
1127 if (x<0) return 0;
1128 else if (m1*x-f1<0) return 0;
1129 else
1130 {
1131 double lmd=(x*x-(F0+k1*G0))/(dF+k1*dG);
1132 F[0]=F0+lmd*dF;
1133 G[0]=G0+lmd*dG;
1134 return 1;
1135 }
1136 }
1137 else
1138 {
1139 double delta=B*B-4*A*C;
1140 if (delta<0) return 0;
1141 else if (delta==0)
1142 {
1143 double x=-B/2/A;
1144 if (x<0) return 0;
1145 else if (m1*x-f1<0) return 0;
1146 else
1147 {
1148 double lmd=(x*x-(F0+k1*G0))/(dF+k1*dG);
1149 F[0]=F0+lmd*dF;
1150 G[0]=G0+lmd*dG;
1151 return 1;
1152 }
1153 }
1154 else
1155 {
1156 int roots=0;
1157 double x=(-B+sqrt(delta))/2/A;
1158 if (x<0) {}
1159 else if (m1*x-f1<0) {}
1160 else
1161 {
1162 double lmd=(x*x-(F0+k1*G0))/(dF+k1*dG);
1163 F[0]=F0+lmd*dF;
1164 G[0]=G0+lmd*dG;
1165 roots++;
1166 }
1167 x=(-B-sqrt(delta))/2/A;
1168 if (x<0) {}
1169 else if (m1*x-f1<0) {}
1170 else
1171 {
1172 double lmd=(x*x-(F0+k1*G0))/(dF+k1*dG);
1173 F[roots]=F0+lmd*dF;
1174 G[roots]=G0+lmd*dG;
1175 roots++;
1176 }
1177 return roots;
1178 }
1179 }
1180 }//FGFrom2
1181
1182 /*
1183 function FGFrom2: solves the following equation system given m[2], f[2], k[2]. This is interpreted as
1184 searching from (F0, G0) down direction (dF, dG) until the first equation is satisfied. Functionally
1185 this is the same as the version using norm[2], with m[] anf f[] normalized by norm[] beforehand so
1186 that norm[] is no longer needed. However, k[2] cannot be computed from normalized m[2] so that it must
1187 be specified explicitly.
1188
1189 m[0](F+k[0]G)^0.5-f[0] = m[1](F+k[1]G)^0.5-f[1] > 0
1190
1191 F=F0+lmd*dF
1192 G=G0+lmd*dG
1193
1194 In: (F0, G0): search origin
1195 (dF, dG): search direction
1196 m[2], f[2], k[2]: coefficients in the first equation
1197 Out: F[return value], G[return value]: solutions.
1198
1199 Returns the number of solutions for (F, G).
1200 */
1201 int FGFrom2(double* F, double* G, double* lmd, double F0, double dF, double G0, double dG, double* m, double* f, int* k)
1202 {
1203 double m1=m[0], m2=m[1],
1204 f1=f[0], f2=f[1],
1205 k1=k[0], k2=k[1];
1206 double A=m1*m1-m2*m2*(dF+k2*dG)/(dF+k1*dG),
1207 B=2*m1*(f2-f1),
1208 C=(f2-f1)*(f2-f1)-(k1-k2)*(F0*dG-G0*dF)/(dF+k1*dG)*m2*m2;
1209 //x is obtained by solving Axx+Bx+C=0
1210 if (fabs(A)<1e-6) //A==0
1211 {
1212 double x=-C/B;
1213 if (x<0) return 0;
1214 else if (m1*x-f1<0) return 0;
1215 else
1216 {
1217 lmd[0]=(x*x-(F0+k1*G0))/(dF+k1*dG);
1218 F[0]=F0+lmd[0]*dF;
1219 G[0]=G0+lmd[0]*dG;
1220 return 1;
1221 }
1222 }
1223 else
1224 {
1225 int roots=0;
1226 double delta=B*B-4*A*C;
1227 if (delta<0) return 0;
1228 else if (delta==0)
1229 {
1230 double x=-B/2/A;
1231 if (x<0) return 0;
1232 else if (m1*x-f1<0) return 0;
1233 else if ((m1*x-f1+f2)*m2<0) {}
1234 else
1235 {
1236 lmd[0]=(x*x-(F0+k1*G0))/(dF+k1*dG);
1237 F[0]=F0+lmd[0]*dF;
1238 G[0]=G0+lmd[0]*dG;
1239 roots++;
1240 }
1241 return roots;
1242 }
1243 else
1244 {
1245 int roots=0;
1246 double x=(-B+sqrt(delta))/2/A;
1247 if (x<0) {}
1248 else if (m1*x-f1<0) {}
1249 else if ((m1*x-f1+f2)*m2<0) {}
1250 else
1251 {
1252 lmd[0]=(x*x-(F0+k1*G0))/(dF+k1*dG);
1253 F[0]=F0+lmd[0]*dF;
1254 G[0]=G0+lmd[0]*dG;
1255 roots++;
1256 }
1257 x=(-B-sqrt(delta))/2/A;
1258 if (x<0) {}
1259 else if (m1*x-f1<0) {}
1260 else if ((m1*x-f1+f2)*m2<0) {}
1261 else
1262 {
1263 lmd[roots]=(x*x-(F0+k1*G0))/(dF+k1*dG);
1264 F[roots]=F0+lmd[roots]*dF;
1265 G[roots]=G0+lmd[roots]*dG;
1266 roots++;
1267 }
1268 return roots;
1269 }
1270 }
1271 }//FGFrom2
1272
1273 /*
1274 function FGFrom3: solves the following equation system given m[3], f[3], norm[3].
1275
1276 m[0](F+k[0]G)^0.5-f[0] m[1](F+k[1]G)^0.5-f[1] m[2](F+k[2]G)^0.5-f[2]
1277 ------------------------ = ------------------------ = ------------------------ > 0
1278 norm[0] norm[1] norm[2]
1279
1280 In: m[3], f[3], norm[3]: coefficients in the above
1281 Out: F[return value], G[return value]: solutions.
1282
1283 Returns the number of solutions for (F, G).
1284 */
1285 int FGFrom3(double* F, double* G, int* m, double* f, double* norm)
1286 {
1287 double m1=m[0]/norm[0], m2=m[1]/norm[1], m3=m[2]/norm[2],
1288 f1=f[0]/norm[0], f2=f[1]/norm[1], f3=f[2]/norm[2],
1289 k1=m[0]*m[0]-1, k2=m[1]*m[1]-1, k3=m[2]*m[2]-1;
1290 double h21=k2-k1, h31=k3-k1; //h21 and h31
1291 double h12=m1*m1-m2*m2, h13=m1*m1-m3*m3; //h12' and h13'
1292 double a=m2*m2*h13*h21-m3*m3*h12*h31;
1293 if (a==0)
1294 {
1295 double x=h12*(f3-f1)*(f3-f1)-h13*(f2-f1)*(f2-f1),
1296 b=h13*(f2-f1)-h12*(f3-f1);
1297 x=x/(2*m1*b);
1298 if (x<0) return 0; //x must the square root of something
1299 else if (m1*x-f1<0) return 0; //discarded because we are solving e1=e2=e3>0
1300 else
1301 {
1302 if (h21!=0)
1303 {
1304 G[0]=(h12*x*x+2*m1*(f2-f1)*x+(f2-f1)*(f2-f1))/(m2*m2*h21);
1305 }
1306 else
1307 {
1308 G[0]=(h13*x*x+2*m1*(f3-f1)*x+(f3-f1)*(f3-f1))/(m3*m3*h31);
1309 }
1310 F[0]=x*x-k1*G[0];
1311 return 1;
1312 }
1313 }
1314 else
1315 {
1316 double b=(h13*(f2-f1)*(f2-f1)-h12*(f3-f1)*(f3-f1))/a;
1317 a=2*m1*(h13*(f2-f1)-h12*(f3-f1))/a;
1318 double A=h12,
1319 B=2*m1*(f2-f1)-m2*m2*h21*a,
1320 C=(f2-f1)*(f2-f1)-m2*m2*h21*b;
1321 double delta=B*B-4*A*C;
1322 if (delta<0) return 0;
1323 else if (delta==0)
1324 {
1325 double x=-B/2/A;
1326 if (x<0) return 0; //x must the square root of something
1327 else if (m1*x-f1<0) return 0; //discarded because we are solving e1=e2=e3>0
1328 else
1329 {
1330 G[0]=a*x+b;
1331 F[0]=x*x-k1*G[0];
1332 return 1;
1333 }
1334 }
1335 else
1336 {
1337 int roots=0;
1338 double x=(-B+sqrt(delta))/2/A;
1339 if (x<0) {}
1340 else if (m1*x-f1<0) {}
1341 else
1342 {
1343 G[0]=a*x+b;
1344 F[0]=x*x-k1*G[0];
1345 roots++;
1346 }
1347 x=(-B-sqrt(delta))/2/A;
1348 if (x<0) {}
1349 else if (m1*x-f1<0) {}
1350 else
1351 {
1352 G[roots]=a*x+b;
1353 F[roots]=x*x-k1*G[roots];
1354 roots++;
1355 }
1356 return roots;
1357 }
1358 }
1359 }//FGFrom3
1360
1361 /*
1362 function FGFrom3: solves the following equation system given m[3], f[3], k[3]. Functionally this is
1363 the same as the version using norm[3], with m[] anf f[] normalized by norm[] beforehand so that norm[]
1364 is no longer needed. However, k[3] cannot be computed from normalized m[3] so that it must be
1365 specified explicitly.
1366
1367 m[0](F+k[0]G)^0.5-f[0] = m[1](F+k[1]G)^0.5-f[1] = m[2](F+k[2]G)^0.5-f[2] >0
1368
1369 In: m[3], f[3], k[3]: coefficients in the above
1370 Out: F[return value], G[return value]: solutions.
1371
1372 Returns the number of solutions for (F, G).
1373 */
1374 int FGFrom3(double* F, double* G, double* e, double* m, double* f, int* k)
1375 {
1376 double m1=m[0], m2=m[1], m3=m[2],
1377 f1=f[0], f2=f[1], f3=f[2],
1378 k1=k[0], k2=k[1], k3=k[2];
1379 double h21=k2-k1, h31=k3-k1; //h21 and h31
1380 double h12=m1*m1-m2*m2, h13=m1*m1-m3*m3; //h12' and h13'
1381 double a=m2*m2*h13*h21-m3*m3*h12*h31;
1382 if (a==0)
1383 {
1384 //a==0 implies two the e's are conjugates
1385 double x=h12*(f3-f1)*(f3-f1)-h13*(f2-f1)*(f2-f1),
1386 b=h13*(f2-f1)-h12*(f3-f1);
1387 x=x/(2*m1*b);
1388 if (x<0) return 0; //x must the square root of something
1389 else if (m1*x-f1<-1e-6) return 0; //discarded because we are solving e1=e2=e3>0
1390 else if ((m1*x-f1+f2)*m2<0) return 0;
1391 else if ((m1*x-f1+f3)*m3<0) return 0;
1392 else
1393 {
1394 if (h21!=0)
1395 {
1396 G[0]=(h12*x*x+2*m1*(f2-f1)*x+(f2-f1)*(f2-f1))/(m2*m2*h21);
1397 }
1398 else
1399 {
1400 G[0]=(h13*x*x+2*m1*(f3-f1)*x+(f3-f1)*(f3-f1))/(m3*m3*h31);
1401 }
1402 F[0]=x*x-k1*G[0];
1403 double tmp=F[0]+G[0]*k[0]; testnn(tmp);
1404 e[0]=m1*sqrt(tmp)-f[0];
1405 return 1;
1406 }
1407 }
1408 else
1409 {
1410 double b=(h13*(f2-f1)*(f2-f1)-h12*(f3-f1)*(f3-f1))/a;
1411 a=2*m1*(h13*(f2-f1)-h12*(f3-f1))/a;
1412 double A=h12,
1413 B=2*m1*(f2-f1)-m2*m2*h21*a,
1414 C=(f2-f1)*(f2-f1)-m2*m2*h21*b;
1415 double delta=B*B-4*A*C;
1416 if (delta<0) return 0;
1417 else if (delta==0)
1418 {
1419 int roots=0;
1420 double x=-B/2/A;
1421 if (x<0) return 0; //x must the square root of something
1422 else if (m1*x-f1<0) return 0; //discarded because we are solving e1=e2=e3>0
1423 else if ((m1*x-f1+f2)*m2<0) return 0;
1424 else if ((m1*x-f1+f3)*m3<0) return 0;
1425 else
1426 {
1427 G[0]=a*x+b;
1428 F[0]=x*x-k1*G[0];
1429 double tmp=F[0]+G[0]*k[0]; testnn(tmp);
1430 e[0]=m1*sqrt(tmp)-f[0];
1431 roots++;
1432 }
1433 return roots;
1434 }
1435 else
1436 {
1437 int roots=0;
1438 double x=(-B+sqrt(delta))/2/A;
1439 if (x<0) {}
1440 else if (m1*x-f1<0) {}
1441 else if ((m1*x-f1+f2)*m2<0) {}
1442 else if ((m1*x-f1+f3)*m3<0) {}
1443 else
1444 {
1445 G[0]=a*x+b;
1446 F[0]=x*x-k1*G[0];
1447 double tmp=F[0]+G[0]*k1; testnn(tmp);
1448 e[0]=m1*sqrt(tmp)-f1;
1449 roots++;
1450 }
1451 x=(-B-sqrt(delta))/2/A;
1452 if (x<0) {}
1453 else if (m1*x-f1<0) {}
1454 else if ((m1*x-f1+f2)*m2<0) {}
1455 else if ((m1*x-f1+f3)*m3<0) {}
1456 else
1457 {
1458 G[roots]=a*x+b;
1459 F[roots]=x*x-k1*G[roots];
1460 double tmp=F[roots]+G[roots]*k1; testnn(tmp);
1461 e[roots]=m1*sqrt(tmp)-f1;
1462 roots++;
1463 }
1464 return roots;
1465 }
1466 }
1467 }//FGFrom3
1468
1469 /*
1470 function FindNote: harmonic sinusoid tracking from a starting point in time-frequency plane forward
1471 and backward.
1472
1473 In: _t, _f: start time and frequency
1474 frst, fren: tracking range, in frames
1475 Spec: spectrogram
1476 wid, offst: atom scale and hop size, must be consistent with Spec
1477 settings: note match settings
1478 brake: tracking termination threshold
1479 Out: M, Fr, partials[M][Fr]: HS partials
1480
1481 Returns 0 if tracking is done by FindNoteF(), 1 if tracking is done by FindNoteConst()
1482 */
1483 int FindNote(int _t, double _f, int& M, int& Fr, atom**& partials, int frst, int fren, int wid, int offst, TQuickSpectrogram* Spec, NMSettings settings)
1484 {
1485 double brake=0.02;
1486 if (settings.delp==0) return FindNoteConst(_t, _f, M, Fr, partials, frst, fren, wid, offst, Spec, settings, brake*5);
1487
1488 atom* part=new atom[(fren-frst)*settings.maxp];
1489 double fpp[1024]; double *vfpp=&fpp[256], *pfpp=&fpp[512]; atomtype *ptype=(atomtype*)&fpp[768]; NMResults results={fpp, vfpp, pfpp, ptype};
1490
1491 int trst=floor((_t-wid/2)*1.0/offst+0.5);
1492 if (trst<0) trst=0;
1493
1494 TPolygon* R=new TPolygon(1024); R->N=0;
1495
1496 cmplx<QSPEC_FORMAT>* spec=Spec->Spec(trst);
1497 double f0=_f*wid;
1498 cdouble *x=new cdouble[wid/2+1]; for (int i=0; i<=wid/2; i++) x[i]=spec[i];
1499
1500 settings.forcepin0=true; settings.pcount=0; settings.pin0asanchor=true;
1501 double B, starts=NoteMatchStiff3(R, f0, B, 1, &x, wid, 0, &settings, &results, 0, 0, ds0);
1502
1503 int k=0;
1504 if (starts>0)
1505 {
1506 int startp=NMResultToAtoms(settings.maxp, part, trst*offst+wid/2, wid, results);
1507 settings.pin0asanchor=false;
1508 double* startvfp=new double[startp]; memcpy(startvfp, vfpp, sizeof(double)*startp);
1509 k=startp;
1510 k+=FindNoteF(&part[k], starts, R, startp, startvfp, trst, fren, wid, offst, Spec, settings, brake);
1511 k+=FindNoteF(&part[k], starts, R, startp, startvfp, trst, frst-1, wid, offst, Spec, settings, brake);
1512 delete[] startvfp;
1513 }
1514 delete R; delete[] x;
1515
1516 AtomsToPartials(k, part, M, Fr, partials, offst);
1517 delete[] part;
1518
1519 // ReEst1(M, frst, fren, partials, WV->Data16[channel], wid, offst);
1520
1521 return 0;
1522 }//Findnote*/
1523
1524 /*
1525 function FindNoteConst: constant-pitch harmonic sinusoid tracking
1526
1527 In: _t, _f: start time and frequency
1528 frst, fren: tracking range, in frames
1529 Spec: spectrogram
1530 wid, offst: atom scale and hop size, must be consistent with Spec
1531 settings: note match settings
1532 brake: tracking termination threshold
1533 Out: M, Fr, partials[M][Fr]: HS partials
1534
1535 Returns 1.
1536 */
1537 int FindNoteConst(int _t, double _f, int& M, int& Fr, atom**& partials, int frst, int fren, int wid, int offst, TQuickSpectrogram* Spec, NMSettings settings, double brake)
1538 {
1539 atom* part=new atom[(fren-frst)*settings.maxp];
1540 double fpp[1024]; double *vfpp=&fpp[256], *pfpp=&fpp[512]; atomtype *ptype=(atomtype*)&fpp[768]; NMResults results={fpp, vfpp, pfpp, ptype};
1541
1542 //trst: the frame index to start tracking
1543 int trst=floor((_t-wid/2)*1.0/offst+0.5), maxp=settings.maxp;
1544 if (trst<0) trst=0;
1545
1546 //find a note candidate for the starting frame at _t
1547 TPolygon* R=new TPolygon(1024); R->N=0;
1548 cmplx<QSPEC_FORMAT>* spec=Spec->Spec(trst);
1549 double f0=_f*wid;
1550 cdouble *x=new cdouble[wid/2+1]; for (int i=0; i<=wid/2; i++) x[i]=spec[i];
1551
1552 settings.forcepin0=true; settings.pcount=0; settings.pin0asanchor=true;
1553 double B, *starts=new double[maxp];
1554 NoteMatchStiff3(R, f0, B, 1, &x, wid, 0, &settings, &results, 0, 0);
1555
1556 //read the energy vector of this candidate HA to starts[]. starts[] will contain the highest single-frame
1557 //energy of each partial during the tracking.
1558 int P=0; while(P<maxp && f0*(P+1)*sqrt(1+B*((P+1)*(P+1)-1))<wid/2) P++;
1559 for (int m=0; m<P; m++) starts[m]=results.vfp[m]*results.vfp[m];
1560 int atomcount=0;
1561
1562 cdouble* ipw=new cdouble[maxp];
1563
1564 //find the ends of tracking by constant-pitch extension of the starting HA until
1565 //at a frame at least half of the partials fall below starts[]*brake.
1566
1567 //first find end of the note by forward extension from the frame at _t
1568 int fr1=trst;
1569 while (fr1<fren)
1570 {
1571 spec=Spec->Spec(fr1);
1572 for (int i=0; i<=wid/2; i++) x[i]=spec[i];
1573 double ls=0;
1574 for (int m=0; m<P; m++)
1575 {
1576 double fm=f0*(m+1)*sqrt(1+B*((m+1)*(m+1)-1));
1577 int K1=floor(fm-settings.hB), K2=ceil(fm+settings.hB);
1578 if (K1<0) K1=0; if (K2>=wid/2) K2=wid/2-1;
1579 ipw[m]=IPWindowC(fm, x, wid, settings.M, settings.c, settings.iH2, K1, K2);
1580 ls+=~ipw[m];
1581 if (starts[m]<~ipw[m]) starts[m]=~ipw[m];
1582 }
1583 double sumstart=0, sumstop=0;
1584 for (int m=0; m<P; m++)
1585 {
1586 if (~ipw[m]<starts[m]*brake) sumstop+=1;// sqrt(starts[m]);
1587 sumstart+=1; //sqrt(starts[m]);
1588 }
1589 double lstop=sumstop/sumstart;
1590 if (lstop>0.5) break;
1591
1592 fr1++;
1593 }
1594 //note find the start of note by backward extension from the frame at _t
1595 int fr0=trst-1;
1596 while(fr0>frst)
1597 {
1598 spec=Spec->Spec(fr0);
1599 for (int i=0; i<=wid/2; i++) x[i]=spec[i];
1600 double ls=0;
1601 for (int m=0; m<P; m++)
1602 {
1603 double fm=f0*(m+1)*sqrt(1+B*((m+1)*(m+1)-1));
1604 int K1=floor(fm-settings.hB), K2=ceil(fm+settings.hB);
1605 if (K1<0) K1=0; if (K2>=wid/2) K2=wid/2-1;
1606 ipw[m]=IPWindowC(fm, x, wid, settings.M, settings.c, settings.iH2, K1, K2);
1607 ls+=~ipw[m];
1608 if (starts[m]<~ipw[m]) starts[m]=~ipw[m];
1609 }
1610
1611 double sumstart=0, sumstop=0;
1612 for (int m=0; m<P; m++)
1613 {
1614 if (~ipw[m]<starts[m]*brake) sumstop+=1; //sqrt(starts[m]);
1615 sumstart+=1; //sqrt(starts[m]);
1616 }
1617 double lstop=sumstop/sumstart;
1618 if (lstop>0.5) {fr0++; break;}
1619
1620 fr0--;
1621 }
1622
1623 //now fr0 and fr1 are the first and last (excl.) frame indices
1624 Fr=fr1-fr0;
1625 cdouble* ipfr=new cdouble[Fr];
1626 double *as=new double[Fr*2], *phs=&as[Fr];
1627 cdouble** Allocate2(cdouble, Fr, wid/2+1, xx);
1628 for (int fr=0; fr<Fr; fr++)
1629 {
1630 spec=Spec->Spec(fr0+fr);
1631 for (int i=0; i<=wid/2; i++) xx[fr][i]=spec[i];
1632 }
1633
1634 //reestimate partial frequencies, amplitudes and phase angles using all frames. In case of frequency estimation
1635 //failure, the original frequency is used and all atoms of that partial are typed "atInfered".
1636 for (int m=0; m<P; m++)
1637 {
1638 double fm=f0*(m+1)*sqrt(1+B*((m+1)*(m+1)-1)), dfm=(settings.delm<1)?settings.delm:1;
1639 double errf=LSESinusoidMP(fm, fm-dfm, fm+dfm, xx, Fr, wid, 3, settings.M, settings.c, settings.iH2, as, phs, 1e-6);
1640 // LSESinusoidMPC(fm, fm-1, fm+1, xx, Fr, wid, offst, 3, settings.M, settings.c, settings.iH2, as, phs, 1e-6);
1641 atomtype type=(errf<0)?atInfered:atPeak;
1642 for (int fr=0; fr<Fr; fr++)
1643 {
1644 double tfr=(fr0+fr)*offst+wid/2;
1645 atom tmpatom={tfr, wid, fm/wid, as[fr], phs[fr], m+1, type, 0};
1646 part[atomcount++]=tmpatom;
1647 }
1648 }
1649
1650 //Tag the atom at initial input (_t, _f) as anchor.
1651 AtomsToPartials(atomcount, part, M, Fr, partials, offst);
1652 partials[settings.pin0-1][trst-fr0].type=atAnchor;
1653 delete[] ipfr;
1654 delete[] ipw;
1655 delete[] part;
1656 delete R; delete[] x; delete[] as; DeAlloc2(xx);
1657 delete[] starts;
1658 return 1;
1659 }//FindNoteConst
1660
1661 /*
1662 function FindNoteF: forward harmonic sinusoid tracking starting from a given harmonic atom until an
1663 endpoint is detected or a search boundary is reached.
1664
1665 In: starts: harmonic atom score of the start HA
1666 startvfp[startp]: amplitudes of partials of start HA
1667 R: F-G polygon of the start HA
1668 frst, frst: frame index of the start and end frames of tracking.
1669 Spec: spectrogram
1670 wid, offst: atom scale and hop size, must be consistent with Spec
1671 settings: note match settings
1672 brake: tracking termination threshold.
1673 Out: part[return value]: list of atoms tracked.
1674 starts: maximal HA score of tracked frames. Its product with brake is used as the tracking threshold.
1675
1676 Returns the number of atoms found.
1677 */
1678 int FindNoteF(atom* part, double& starts, TPolygon* R, int startp, double* startvfp, int frst, int fren, int wid, int offst, TQuickSpectrogram* Spec, NMSettings settings, double brake)
1679 {
1680 settings.forcepin0=false, settings.pin0=0;
1681
1682 int k=0, maxp=settings.maxp;
1683
1684 TPolygon* RR=new TPolygon(1024, R);
1685 int lastp=startp;
1686 double *lastvfp=new double[settings.maxp]; memcpy(lastvfp, startvfp, sizeof(double)*startp);
1687
1688 cdouble *x=new cdouble[wid/2+1];
1689 double *fpp=new double[maxp*4], *vfpp=&fpp[maxp], *pfpp=&fpp[maxp*2]; atomtype* ptype=(atomtype*)&fpp[maxp*3];
1690 NMResults results={fpp, vfpp, pfpp, ptype};
1691 int step=(fren>frst)?1:(-1);
1692 for (int h=frst+step; (h-fren)*step<0; h+=step)
1693 {
1694 cmplx<QSPEC_FORMAT>* spec=Spec->Spec(h);
1695 for (int i=0; i<=wid/2; i++) x[i]=spec[i];
1696 double f0=0, B=0;
1697 double tmp=NoteMatchStiff3(RR, f0, B, 1, &x, wid, 0, &settings, &results, lastp, lastvfp);
1698 if (starts<tmp) starts=tmp;
1699 if (tmp<starts*brake) break; //end condition: power drops by ??dB
1700
1701 lastp=NMResultToAtoms(settings.maxp, &part[k], h*offst+wid/2, wid, results); k+=lastp;
1702 memcpy(lastvfp, vfpp, sizeof(double)*lastp);
1703 }
1704 delete RR; delete[] lastvfp; delete[] x; delete[] fpp;
1705 return k;
1706 }//FindNoteF
1707
1708 /*
1709 function FindNoteFB: forward-backward harmonic sinusid tracking.
1710
1711 In: frst, fren: index of start and end frames
1712 Rst, Ren: F-G polygons of start and end harmonic atoms
1713 vfpst[M], vfpen[M]: amplitude vectors of start and end harmonic atoms
1714 Spec: spectrogram
1715 wid, offst: atom scale and hop size, must be consistent with Spec
1716 settings: note match settings
1717 Out: partials[0:fren-frst][M], hosting tracked HS between frst and fren (inc).
1718
1719 Returns 0 if successful, 1 if failure in pitch tracking stage (no feasible HA candidate at some frame).
1720 On start partials[0:fren-frst][M] shall be prepared to receive fren-frst+1 harmonic atoms
1721 */
1722 int FindNoteFB(int frst, TPolygon* Rst, double* vfpst, int fren, TPolygon* Ren, double* vfpen, int M, atom** partials, int wid, int offst, TQuickSpectrogram* Spec, NMSettings settings)
1723 { //*
1724 if (frst>fren) return -1;
1725 int result=0;
1726 if (settings.maxp>M) settings.maxp=M;
1727 else M=settings.maxp;
1728 int Fr=fren-frst+1;
1729 double B, fpp[1024], delm=settings.delm, delp=settings.delp, iH2=settings.iH2;
1730 double *vfpp=&fpp[256], *pfpp=&fpp[512]; atomtype *ptype=(atomtype*)&fpp[768];
1731 int maxpitch=50;
1732
1733 NMResults results={fpp, vfpp, pfpp, ptype};
1734
1735 cmplx<QSPEC_FORMAT>* spec;
1736
1737 cdouble *x=new cdouble[wid/2+1];
1738 TPolygon **Ra=new TPolygon*[maxpitch], **Rb=new TPolygon*[maxpitch], ***R; memset(Ra, 0, sizeof(TPolygon*)*maxpitch), memset(Rb, 0, sizeof(TPolygon*)*maxpitch);
1739 double *sca=new double[maxpitch], *scb=new double[maxpitch], **sc; sca[0]=scb[0]=0;
1740 double** Allocate2(double, maxpitch, M, vfpa); memcpy(vfpa[0], vfpst, sizeof(double)*M);
1741 double** Allocate2(double, maxpitch, M, vfpb); memcpy(vfpb[0], vfpen, sizeof(double)*M);
1742 double*** vfp;
1743
1744 double** Allocate2(double, Fr, wid/2, fps);
1745 double** Allocate2(double, Fr, wid/2, vps);
1746 int* pc=new int[Fr];
1747
1748 Ra[0]=new TPolygon(M*2+4, Rst);
1749 Rb[0]=new TPolygon(M*2+4, Ren);
1750 int pitchcounta=1, pitchcountb=1, pitchcount;
1751
1752 //pitch tracking
1753 int** Allocate2(int, Fr, maxpitch, prev);
1754 int** Allocate2(int, Fr, maxpitch, pitches);
1755 int* pitchres=new int[Fr];
1756 int fra=frst, frb=fren, fr;
1757 while (fra<frb-1)
1758 {
1759 if (pitchcounta<pitchcountb){fra++; fr=fra; pitchcount=pitchcounta; vfp=&vfpa; sc=&sca; R=&Ra;}
1760 else {frb--; fr=frb; pitchcount=pitchcountb; vfp=&vfpb; sc=&scb, R=&Rb;}
1761 int fr_frst=fr-frst;
1762
1763 int absfr=(partials[0][fr].t-wid/2)/offst;
1764 spec=Spec->Spec(absfr); for (int i=0; i<=wid/2; i++) x[i]=spec[i];
1765 pc[fr_frst]=QuickPeaks(fps[fr_frst], vps[fr_frst], wid, x, settings.M, settings.c, iH2, 0.0005);
1766 NoteMatchStiff3FB(pitchcount, *R, *vfp, *sc, pitches[fr_frst], prev[fr_frst], pc[fr_frst], fps[fr_frst], vps[fr_frst], x, wid, maxpitch, &settings);
1767 if (fr==fra) pitchcounta=pitchcount;
1768 else pitchcountb=pitchcount;
1769 if (pitchcount<=0)
1770 {result=1; goto cleanup;}
1771 }
1772 {
1773 //now fra=frb-1
1774 int fra_frst=fra-frst, frb_frst=frb-frst, maxpitcha=0, maxpitchb=0;
1775 double maxs=0;
1776 for (int i=0; i<pitchcounta; i++)
1777 {
1778 double f0a, f0b;
1779 if (fra==frst) f0a=partials[0][frst].f*wid;
1780 else
1781 {
1782 int peak0a=((__int16*)&pitches[fra_frst][i])[0], pin0a=((__int16*)&pitches[fra_frst][i])[1];
1783 f0a=fps[fra_frst][peak0a]/pin0a;
1784 }
1785 for (int j=0; j<pitchcountb; j++)
1786 {
1787 if (frb==fren) f0b=partials[0][fren].f*wid;
1788 else
1789 {
1790 int peak0b=((__int16*)&pitches[frb_frst][j])[0], pin0b=((__int16*)&pitches[frb_frst][j])[1];
1791 f0b=fps[frb_frst][peak0b]/pin0b;
1792 }
1793 double pow2delp=pow(2, delp/12);
1794 if (f0b<f0a*pow2delp+delm && f0b>f0a/pow2delp-delm)
1795 {
1796 double s=sca[i]+scb[j]+conta(M, vfpa[i], vfpb[j]);
1797 if (maxs<s) maxs=s, maxpitcha=i, maxpitchb=j;
1798 }
1799 }
1800 }
1801 //get pitch tracking result
1802 pitchres[fra_frst]=pitches[fra_frst][maxpitcha];
1803 for (int fr=fra; fr>frst; fr--)
1804 {
1805 maxpitcha=prev[fr-frst][maxpitcha];
1806 pitchres[fr-frst-1]=pitches[fr-frst-1][maxpitcha];
1807 }
1808 pitchres[frb_frst]=pitches[frb_frst][maxpitchb];
1809 for (int fr=frb; fr<fren; fr++)
1810 {
1811 maxpitchb=prev[fr-frst][maxpitchb];
1812 pitchres[fr-frst+1]=pitches[fr-frst+1][maxpitchb];
1813 }
1814 }
1815 {
1816 double f0s[1024]; memset(f0s, 0, sizeof(double)*1024);
1817 for (int i=frst+1; i<fren; i++)
1818 {
1819 int pitch=pitchres[i-frst];
1820 int peak0=((__int16*)&pitch)[0], pin0=((__int16*)&pitch)[1];
1821 f0s[i]=fps[i-frst][peak0]/pin0;
1822 }
1823 f0s[frst]=sqrt(Rst->X[0]), f0s[fren]=sqrt(Ren->X[0]);
1824
1825 //partial tracking
1826 int lastp=0; while (lastp<M && vfpst[lastp]>0) lastp++;
1827 double* lastvfp=new double[M]; memcpy(lastvfp, vfpst, sizeof(double)*M);
1828 delete Ra[0]; Ra[0]=new TPolygon(M*2+4, Rst);
1829 for (int h=frst+1; h<fren; h++)
1830 {
1831 int absfr=(partials[0][h].t-wid/2)/offst;
1832 int h_frst=h-frst, pitch=pitchres[h_frst];
1833 int peak0=((__int16*)&pitch)[0], pin0=((__int16*)&pitch)[1];
1834 spec=Spec->Spec(absfr);
1835 double f0=fpp[0];//, B=Form11->BEdit->Text.ToDouble();
1836 //double deltaf=f0*(pow(2, settings.delp/12)-1);
1837 for (int i=0; i<=wid/2; i++) x[i]=spec[i];
1838 memset(fpp, 0, sizeof(double)*1024);
1839 //settings.epf0=deltaf,
1840 settings.forcepin0=true; settings.pin0=pin0; f0=fps[h_frst][peak0]; settings.pin0asanchor=false;
1841 if (NoteMatchStiff3(Ra[0], f0, B, pc[h_frst], fps[h_frst], vps[h_frst], 1, &x, wid, 0, &settings, &results, lastp, lastvfp)>0)
1842 {
1843 lastp=NMResultToPartials(M, h, partials, partials[0][h].t, wid, results);
1844 memcpy(lastvfp, vfpp, sizeof(double)*lastp);
1845 }
1846 }
1847 delete[] lastvfp;
1848 }
1849 cleanup:
1850 delete[] x;
1851 for (int i=0; i<pitchcounta; i++) delete Ra[i]; delete[] Ra;
1852 for (int i=0; i<pitchcountb; i++) delete Rb[i]; delete[] Rb;
1853 delete[] sca; delete[] scb;
1854 DeAlloc2(vfpa); DeAlloc2(vfpb);
1855 DeAlloc2(fps); DeAlloc2(vps);
1856 DeAlloc2(pitches); DeAlloc2(prev);
1857 delete[] pitchres; delete[] pc;
1858 return result; //*/
1859 }//FindNoteFB
1860
1861 /*
1862 function GetResultsSingleFr: used by NoteMatchStiff3 to obtain final note match results after harmonic
1863 grouping of peaks with rough frequency estimates, including a screening of found peaks based on shape
1864 and consistency with other peak frequencies, reestimating sinusoid parameters using LSE, estimating
1865 atoms without peaks by inferring their frequencies, and translating internal peak type to atom type.
1866
1867 In: R0: F-G polygon of primitive (=rough estimates) partial frequencies
1868 x[N]: spectrum
1869 ps[P], fs[P]; partial indices and frequencies, may miss partials
1870 rsrs[P]: peak shape factors, used for evaluating peak validity
1871 psb[P]: peak type flags, related to atom::ptype.
1872 settings: note match settings
1873 numsam: number of partials to evaluate (=number of atoms to return)
1874 Out: results: note match results as a NMResult structure
1875 f0, B: fundamental frequency and stiffness coefficient
1876 Rret: F-G polygon of reestimated set of partial frequencies
1877
1878 Return the total energy of the harmonic atom.
1879 */
1880 double GetResultsSingleFr(double& f0, double& B, TPolygon* Rret, TPolygon* R0, int P, int* ps, double* fs, double* rsrs, cdouble* x, int N, int* psb, int numsam, NMSettings* settings, NMResults* results)
1881 {
1882 Rret->N=0;
1883 double delm=settings->delm, *c=settings->c, iH2=settings->iH2, epf=settings->epf, maxB=settings->maxB, hB=settings->hB;
1884 int M=settings->M, maxp=settings->maxp;
1885 double *fp=results->fp; memset(fp, 0, sizeof(double)*maxp);
1886
1887 double result=0;
1888 if (P>0)
1889 {
1890 double *vfp=results->vfp, *pfp=results->pfp;
1891 atomtype *ptype=results->ptype; //memset(ptype, 0, sizeof(int)*maxp);
1892
1893 double *f1=new double[(maxp+1)*4], *f2=&f1[maxp+1], *ft=&f2[maxp+1], *fdep=&ft[maxp+1], tmpa, cF, cG;
1894 areaandcentroid(tmpa, cF, cG, R0->N, R0->X, R0->Y);
1895 for (int p=1; p<=maxp; p++) ExFmStiff(f1[p], f2[p], p, R0->N, R0->X, R0->Y), ft[p]=p*sqrt(cF+(p*p-1)*cG);
1896 //sort peaks by rsr and departure from model so that most reliable ones are found at the start
1897 double* values=new double[P]; int* indices=new int[P];
1898 for (int i=0; i<P; i++)
1899 {
1900 values[i]=1e200;
1901 int lp=ps[i]; double lf=fs[i];
1902 if (lf>=f1[lp] && lf<=f2[lp]) fdep[lp]=0;
1903 else if (lf<f1[lp]) fdep[lp]=f1[lp]-lf;
1904 else if (lf>f2[lp]) fdep[lp]=lf-f2[lp];
1905 double tmpv=(fdep[lp]>0.5)?(fdep[lp]-0.5)*2:0;
1906 tmpv+=rsrs?rsrs[i]:0;
1907 tmpv*=pow(lp, 0.2);
1908 InsertInc(tmpv, i, values, indices, i+1);
1909 }
1910
1911 for (int i=0; i<P; i++)
1912 {
1913 int ind=indices[i];
1914 int lp=ps[ind]; //partial index
1915 //Get LSE estimation of atoms
1916 double lf=fs[ind], f1, f2, la, lph;//, lrsr=rsrs?rsrs[ind]:0
1917 if (lp==0 || lf==0) continue;
1918 ExFmStiff(f1, f2, lp, R0->N, R0->X, R0->Y);
1919 LSESinusoid(lf, f1-delm, f2+delm, x, N, 3, M, c, iH2, la, lph, epf);
1920 //lrsr=PeakShapeC(lf, 1, N, &x, 6, M, c, iH2);
1921 if (Rret->N>0) ExFmStiff(f1, f2, lp, Rret->N, Rret->X, Rret->Y);
1922 if (Rret->N==0 || (lf>=f1 && lf<=f2))
1923 {
1924 fp[lp-1]=lf;
1925 vfp[lp-1]=la;
1926 pfp[lp-1]=lph;
1927 if (psb[lp]==1 || (psb[lp]==2 && settings->pin0asanchor)) ptype[lp-1]=atAnchor; //0 for anchor points
1928 else ptype[lp-1]=atPeak; //1 for local maximum
1929 //update R using found partails with amplitude>1
1930 if (Rret->N==0) InitializeR(Rret, lp, lf, delm, maxB);
1931 else if (la>1) CutR(Rret, lp, lf, delm, true);
1932 }
1933 }
1934 //*
1935 for (int p=1; p<=maxp; p++)
1936 {
1937 if (fp[p-1]>0)
1938 {
1939 double f1, f2;
1940 ExFmStiff(f1, f2, p, Rret->N, Rret->X, Rret->Y);
1941 if (fp[p-1]<f1-0.3 || fp[p-1]>f2+0.3)
1942 fp[p-1]=0;
1943 }
1944 }// */
1945
1946
1947 //estimate f0 and B
1948 double norm[1024]; for (int i=0; i<1024; i++) norm[i]=1;
1949 areaandcentroid(tmpa, cF, cG, Rret->N, Rret->X, Rret->Y); //minimaxstiff(cF, cG, P, ps, fs, norm, R->N, R->X, R->Y);
1950 testnn(cF); f0=sqrt(cF); B=cG/cF;
1951
1952 //Get LSE estimates for unfound partials
1953 for (int i=0; i<numsam; i++)
1954 {
1955 if (fp[i]==0) //no peak is found for this partial in lcand
1956 {
1957 int m=i+1;
1958 double tmp=cF+(m*m-1)*cG; testnn(tmp);
1959 double lf=m*sqrt(tmp);
1960 if (lf<N/2.1)
1961 {
1962 fp[i]=lf;
1963 ptype[i]=atInfered; //2 for non peak
1964 cdouble r=IPWindowC(lf, x, N, M, c, iH2, lf-hB, lf+hB);
1965 vfp[i]=sqrt(r.x*r.x+r.y*r.y);
1966 pfp[i]=Atan2(r.y, r.x);//*/
1967 }
1968 }
1969 }
1970 if (f0>0) for (int i=0; i<numsam; i++) if (fp[i]>0) result+=vfp[i]*vfp[i];
1971
1972 delete[] f1; delete[] values; delete[] indices;
1973 }
1974 return result;
1975 }//GetResultsSingleFr
1976
1977 /*
1978 function GetResultsMultiFr: the constant-pitch multi-frame version of GetResultsSingleFr.
1979
1980 In: x[Fr][N]: spectrogram
1981 ps[P], fs[P]; partial indices and frequencies, may miss partials
1982 psb[P]: peak type flags, related to atom::ptype.
1983 settings: note match settings
1984 forceinputlocalfr: specifies if partial settings->pin0 is taken for granted ("pinned")
1985 numsam: number of partials to evaluate (=number of atoms to return)
1986 Out: results[Fr]: note match results as Fr NMResult structures
1987 f0, B: fundamental frequency and stiffness coefficient
1988 Rret: F-G polygon of reestimated set of partial frequencies
1989
1990 Return the total energy of the harmonic sinusoid.
1991 */
1992 double GetResultsMultiFr(double& f0, double& B, TPolygon* Rret, TPolygon* R0, int P, int* ps, double* fs, double* rsrs, int Fr, cdouble** x, int N, int offst, int* psb, int numsam, NMSettings* settings, NMResults* results, int forceinputlocalfr)
1993 {
1994 Rret->N=0;
1995 double delm=settings->delm, *c=settings->c, iH2=settings->iH2, epf=settings->epf, maxB=settings->maxB, hB=settings->hB;// *pinf=settings->pinf;
1996
1997 int M=settings->M, maxp=settings->maxp, pincount=settings->pcount, *pin=settings->pin, *pinfr=settings->pinfr;
1998 // double *fp=results->fp, *vfp=results->vfp, *pfp=results->pfp;
1999 // int *ptype=results->ptype;
2000 for (int fr=0; fr<Fr; fr++) memset(results[fr].fp, 0, sizeof(double)*maxp);
2001 int* pclass=new int[maxp+1]; memset(pclass, 0, sizeof(int)*(maxp+1));
2002 for (int i=0; i<pincount; i++) if (pinfr[i]>=0) pclass[pin[i]]=1;
2003 if (forceinputlocalfr>=0) pclass[settings->pin0]=1;
2004 double result=0;
2005 if (P>0)
2006 {
2007 double tmpa, cF, cG;
2008 //found atoms
2009 double *as=new double[Fr*2], *phs=&as[Fr];
2010 for (int i=P-1; i>=0; i--)
2011 {
2012 int lp=ps[i];
2013 double lf=fs[i];
2014 if (lp==0 || lf==0) continue;
2015
2016 if (!pclass[lp])
2017 {
2018 //Get LSE estimation of atoms
2019 double startlf=lf;
2020 LSESinusoidMP(lf, lf-1, lf+1, x, Fr, N, 3, M, c, iH2, as, phs, epf);
2021 //LSESinusoidMPC(lf, lf-1, lf+1, x, Fr, N, offst, 3, M, c, iH2, as, phs, epf);
2022 if (fabs(lf-startlf)>0.6)
2023 {
2024 lf=startlf;
2025 for (int fr=0; fr<Fr; fr++)
2026 {
2027 cdouble r=IPWindowC(lf, x[fr], N, M, c, iH2, lf-settings->hB, lf+settings->hB);
2028 as[fr]=abs(r);
2029 phs[fr]=arg(r);
2030 }
2031 }
2032 }
2033 else //frequencies of local anchor atoms are not re-estimated
2034 {
2035 for (int fr=0; fr<Fr; fr++)
2036 {
2037 cdouble r=IPWindowC(lf, x[fr], N, M, c, iH2, lf-settings->hB, lf+settings->hB);
2038 as[fr]=abs(r);
2039 phs[fr]=arg(r);
2040 }
2041 }
2042
2043 //LSESinusoid(lf, f1-delm, f2+delm, x, N, 3, M, c, iH2, la, lph, epf);
2044 //fp[lp-1]=lf; vfp[lp-1]=la; pfp[lp-1]=lph;
2045 for (int fr=0; fr<Fr; fr++) results[fr].fp[lp-1]=lf, results[fr].vfp[lp-1]=as[fr], results[fr].pfp[lp-1]=phs[fr];
2046 if (psb[lp]==1 || (psb[lp]==2 && settings->pin0asanchor)) for (int fr=0; fr<Fr; fr++) results[fr].ptype[lp-1]=atAnchor; //0 for anchor points
2047 else for (int fr=0; fr<Fr; fr++) results[fr].ptype[lp-1]=atPeak; //1 for local maximum
2048 //update R using found partails with amplitude>1
2049 if (Rret->N==0)
2050 {
2051 //temporary treatment: +0.5 should be +rsr or something similar
2052 InitializeR(Rret, lp, lf, delm+0.5, maxB);
2053 areaandcentroid(tmpa, cF, cG, Rret->N, Rret->X, Rret->Y); //minimaxstiff(cF, cG, P, ps, fs, norm, R->N, R->X, R->Y);
2054 }
2055 else //if (la>1)
2056 {
2057 CutR(Rret, lp, lf, delm+0.5, true);
2058 areaandcentroid(tmpa, cF, cG, Rret->N, Rret->X, Rret->Y); //minimaxstiff(cF, cG, P, ps, fs, norm, R->N, R->X, R->Y);
2059 }
2060 }
2061 //estimate f0 and B
2062 double norm[1024]; for (int i=0; i<1024; i++) norm[i]=1;
2063 areaandcentroid(tmpa, cF, cG, Rret->N, Rret->X, Rret->Y); //minimaxstiff(cF, cG, P, ps, fs, norm, R->N, R->X, R->Y);
2064 testnn(cF); f0=sqrt(cF); B=cG/cF;
2065
2066 //Get LSE estimates for unfound partials
2067 for (int i=0; i<numsam; i++)
2068 {
2069 if (results[0].fp[i]==0) //no peak is found for this partial in lcand
2070 {
2071 int m=i+1;
2072 double tmp=cF+(m*m-1)*cG; testnn(tmp);
2073 double lf=m*sqrt(tmp);
2074 if (lf<N/2.1)
2075 {
2076 for (int fr=0; fr<Fr; fr++)
2077 {
2078 results[fr].fp[i]=lf, results[fr].ptype[i]=atInfered; //2 for non peak
2079 int k1=ceil(lf-hB); if (k1<0) k1=0;
2080 int k2=floor(lf+hB); if (k2>=N/2) k2=N/2-1;
2081 cdouble r=IPWindowC(lf, x[fr], N, M, c, iH2, k1, k2);
2082 results[fr].vfp[i]=abs(r);
2083 results[fr].pfp[i]=arg(r);
2084 }
2085 }
2086 }
2087 }
2088 delete[] as;
2089 if (f0>0) for (int fr=0; fr<Fr; fr++) for (int i=0; i<numsam; i++) if (results[fr].fp[i]>0) result+=results[fr].vfp[i]*results[fr].vfp[i];
2090 result/=Fr;
2091 }
2092 delete[] pclass;
2093 return result;
2094 }//GetResultsMultiFr
2095
2096 /*
2097 function InitailizeR: initializes a F-G polygon with a fundamental frequency range and stiffness coefficient bound
2098
2099 In: af, ef: centre and half width of the fundamental frequency range
2100 maxB: maximal value of stiffness coefficient (the minimal is set to 0)
2101 Out: R: the initialized F-G polygon.
2102
2103 No reutrn value.
2104 */
2105 void InitializeR(TPolygon* R, double af, double ef, double maxB)
2106 {
2107 R->N=4;
2108 double *X=R->X, *Y=R->Y;
2109 double g1=af-ef, g2=af+ef;
2110 if (g1<0) g1=0;
2111 g1=g1*g1, g2=g2*g2;
2112 X[0]=X[3]=g1, X[1]=X[2]=g2;
2113 Y[0]=X[0]*maxB, Y[1]=X[1]*maxB, Y[2]=Y[3]=0;
2114 }//InitializeR
2115
2116 /*
2117 function InitialzeR: initializes a F-G polygon with a frequency range for a given partial and
2118 stiffness coefficient bound
2119
2120 In: apind: partial index
2121 af, ef: centre and half width of the frequency range of the apind-th partial
2122 maxB; maximal value of stiffness coefficient (the minimal is set to 0)
2123 Out: R: the initialized F-G polygon.
2124
2125 No return value.
2126 */
2127 void InitializeR(TPolygon* R, int apind, double af, double ef, double maxB)
2128 {
2129 R->N=4;
2130 double *X=R->X, *Y=R->Y;
2131 double k=apind*apind-1;
2132 double g1=(af-ef)/apind, g2=(af+ef)/apind;
2133 if (g1<0) g1=0;
2134 g1=g1*g1, g2=g2*g2;
2135 double kb1=1+k*maxB;
2136 X[0]=g1/kb1, X[1]=g2/kb1, X[2]=g2, X[3]=g1;
2137 Y[0]=X[0]*maxB, Y[1]=X[1]*maxB, Y[2]=Y[3]=0;
2138 }//InitializeR
2139
2140 /*
2141 function maximalminimum: finds the point within a polygon that maximizes its minimal distance to the
2142 sides.
2143
2144 In: sx[N], sy[N]: x- and y-coordinates of vertices of a polygon
2145 Out: (x, y): point within the polygon with maximal minimal distance to the sides
2146
2147 Returns the maximial minimal distance. A circle centred at (x, y) with the return value as the radius
2148 is the maximum inscribed circle of the polygon.
2149 */
2150 double maximalminimum(double& x, double& y, int N, double* sx, double* sy)
2151 {
2152 //at the beginning let (x, y) be (sx[0], sy[0]), then the mininum distance d is 0.
2153 double dm=0;
2154 x=sx[0], y=sy[0];
2155 double *A=new double[N*3];
2156 double *B=&A[N], *C=&A[N*2];
2157 //calcualte equations of all sides A[k]x+B[k]y+C[k]=0, with signs adjusted so that for
2158 // any (x, y) within the polygon, A[k]x+B[k]y+C[k]>0. A[k]^2+B[k]^2=1.
2159 for (int n=0; n<N; n++)
2160 {
2161 double x1=sx[n], y1=sy[n], x2, y2, AA, BB, CC;
2162 if (n+1!=N) x2=sx[n+1], y2=sy[n+1];
2163 else x2=sx[0], y2=sy[0];
2164 double dx=x1-x2, dy=y1-y2;
2165 double ds=sqrt(dx*dx+dy*dy);
2166 AA=dy/ds, BB=-dx/ds;
2167 CC=-AA*x1-BB*y1;
2168 //adjust signs
2169 if (n+2<N) x2=sx[n+2], y2=sy[n+2];
2170 else x2=sx[n+2-N], y2=sy[n+2-N];
2171 if (AA*x2+BB*y2+CC<0) A[n]=-AA, B[n]=-BB, C[n]=-CC;
2172 else A[n]=AA, B[n]=BB, C[n]=CC;
2173 }
2174
2175 //during the whole process (x, y) is always equal-distance to at least two sides,
2176 // namely l1--(l1+1) and l2--(l2+1). there equations are A1x+B1y+C1=0 and A2x+B2y+C2=0,
2177 // where A^2+B^2=1.
2178 int l1=0, l2=N-1;
2179
2180 double a, b;
2181 b=A[l1]-A[l2], a=B[l2]-B[l1];
2182 double ds=sqrt(a*a+b*b);
2183 a/=ds, b/=ds;
2184 //the line (x+at, y+bt) passes (x, y), and points on this line are equal-distance to l1 and l2.
2185 //along this line at point (x, y), we have dx=ads, dy=bds
2186 //now find the signs so that dm increases as ds>0
2187 double ddmds=A[l1]*a+B[l1]*b;
2188 if (ddmds<0) a=-a, b=-b, ddmds=-ddmds;
2189
2190 while (true)
2191 {
2192 //now the vector starting from (x, y) pointing in (a, b) is equi-distance-to-l1-and-l2 and
2193 // dm-increasing. actually at s from (x,y), d=dm+ddmds*s.
2194
2195 //it is now guaranteed that the distance of (x, y) to l1 (=l2) is smaller than to any other
2196 // sides. along direction (A, B) the distance of (x, y) to l1 (=l2) is increasing and
2197 // the distance to at least one other sides is increasing, so that at some value of s the
2198 // distances equal. we find the smallest s>0 where this happens.
2199 int l3=-1;
2200 double s=-1;
2201 for (int n=0; n<N; n++)
2202 {
2203 if (n==l1 || n==l2) continue;
2204 //distance of (x,y) to the side
2205 double ldm=A[n]*x+B[n]*y+C[n]; //ldm>dm
2206 double dldmds=A[n]*a+B[n]*b;
2207 //so ld=ldm+lddmds*s, the equality is dm+ddmds*s=ldm+lddmds*s
2208 if (ddmds-dldmds>0)
2209 {
2210 ds=(ldm-dm)/(ddmds-dldmds);
2211 if (l3==-1) l3=n, s=ds;
2212 else if (ds<s) l3=n, s=ds;
2213 }
2214 }
2215 if (ddmds==0) s/=2;
2216 x=x+a*s, y=y+b*s;
2217 dm=A[l1]*x+B[l1]*y+C[l1];
2218 // Form1->Canvas->Ellipse(x-dm, y-dm, x+dm, y+dm);
2219 //(x, y) is equal-distance to l1, l2 and l3
2220 //try use l3 to substitute l2
2221 b=A[l1]-A[l3], a=B[l3]-B[l1];
2222 ds=sqrt(a*a+b*b);
2223 a/=ds, b/=ds;
2224 ddmds=A[l1]*a+B[l1]*b;
2225 if (ddmds<0) a=-a, b=-b, ddmds=-ddmds;
2226 if (ddmds==0 || A[l2]*a+B[l2]*b>0)
2227 {
2228 l2=l3;
2229 }
2230 else //l1<-l3 fails, then try use l3 to substitute l2
2231 {
2232 b=A[l3]-A[l2], a=B[l2]-B[l3];
2233 ds=sqrt(a*a+b*b);
2234 a/=ds, b/=ds;
2235 ddmds=A[l3]*a+B[l3]*b;
2236 if (ddmds<0) a=-a, b=-b, ddmds=-ddmds;
2237 if (ddmds==0 || A[l1]*a+B[l1]*b>0)
2238 {
2239 l1=l3;
2240 }
2241 else break;
2242 }
2243 }
2244
2245 delete[] A;
2246 return dm;
2247 }//maximalminimum
2248
2249 /*
2250 function minimaxstiff: finds the point in polygon (N; F, G) that minimizes the maximal value of the
2251 function
2252
2253 | m[l]sqrt(F+(m[l]^2-1)*G)-f[l] |
2254 e_l = | ----------------------------- | regarding l=0, ..., L-1
2255 | norm[l] |
2256
2257 In: _m[L], _f[L], norm[L]: coefficients in the above
2258 F[N], G[N]: vertices of a F-G polygon
2259 Out: (F1, G1): the mini-maximum.
2260
2261 Returns the minimized maximum value.
2262
2263 Further reading: Wen X, "Estimating model parameters", Ch.3.2.6 from "Harmonic sinusoid modeling of
2264 tonal music events", PhD thesis, University of London, 2007.
2265 */
2266 double minimaxstiff(double& F1, double& G1, int L, int* _m, double* _f, double* norm, int N, double* F, double* G)
2267 {
2268 if (L==0 || N<=2)
2269 {
2270 F1=F[0], G1=G[0];
2271 return 0;
2272 }
2273 //normalizing
2274 double* m=(double*)malloc(sizeof(double)*L*6);//new double[L*6];
2275 double* f=&m[L*2];
2276 int* k=(int*)&m[L*4];
2277 for (int l=0; l<L; l++)
2278 {
2279 k[2*l]=_m[l]*_m[l]-1;
2280 k[2*l+1]=k[2*l];
2281 m[2*l]=_m[l]/norm[l];
2282 m[2*l+1]=-m[2*l];
2283 f[2*l]=_f[l]/norm[l];
2284 f[2*l+1]=-f[2*l];
2285 }
2286 //From this point on the L distance functions with absolute value is replace by 2L distance functions
2287 L*=2;
2288 double* vmnl=new double[N*2];
2289 int* mnl=(int*)&vmnl[N];
2290 start:
2291 //Initialize (F0, G0) to be the polygon vertex that has the minimal max_l(e_l)
2292 // maxn: the vertex index
2293 // maxl: the l that maximizes e_l at that vertex
2294 // maxsg: the sign of e_l before taking the abs. value
2295 int nc=-1, nd;
2296 int l1=-1, l2=0, l3=0;
2297 double vmax=0;
2298
2299 for (int n=0; n<N; n++)
2300 {
2301 int lmax=-1;
2302 double lvmax=0;
2303 for (int l=0; l<L; l++)
2304 {
2305 double tmp=F[n]+k[l]*G[n]; testnn(tmp);
2306 double e=m[l]*sqrt(tmp)-f[l];
2307 if (e>lvmax) lvmax=e, lmax=l;
2308 }
2309 mnl[n]=lmax, vmnl[n]=lvmax;
2310 if (n==0)
2311 {
2312 vmax=lvmax, nc=n, l1=lmax;
2313 }
2314 else
2315 {
2316 if (lvmax<vmax) vmax=lvmax, nc=n, l1=lmax;
2317 }
2318 }
2319 double F0=F[nc], G0=G[nc];
2320
2321
2322 // start searching the the minimal maximum from (F0, G0)
2323 //
2324 // Each searching step starts from (F0, G0), ends at (F1, G1)
2325 //
2326 // Starting conditions of one step:
2327 //
2328 // (F0, G0) can be
2329 // (1)inside polygon R;
2330 // (2)on one side of R (nc:(nc+1) gives the side)
2331 // (3)a vertex of R (nc being the vertex index)
2332 //
2333 // The maximum at (F0, G0) can be
2334 // (1)vmax=e1=e2=e3>...; (l1, l2, l3)
2335 // (2)vmax=e1=e2>...; (l1, l2)
2336 // (3)vmax=e1>.... (l1)
2337 //
2338 // More complication arise if we have more than 3 equal maxima, i.e.
2339 // vmax=e1=e2=e3=e4>=.... CURRENTLY WE ASSUME THIS NEVER HAPPENS.
2340 //
2341 // These are also the ending conditions of one step.
2342 //
2343 // There are types of basic searching steps, i.e.
2344 //
2345 // (1) e1=e2 search: starting with e1=e2, search down the decreasing direction
2346 // of e1=e2 until at point (F1, G1) there is another l3 so that e1(F1, G1)
2347 // =e2(F1, G1)=e3(F1, G1), or until at point (F1, G1) the search is about
2348 // to go out of R.
2349 // Arguments: l1, l2, (F0, G0, vmax)
2350 // (2) e1!=e2 search: starting with e1=e2 and (F0, G0) being on one side, search
2351 // down the side in the decreasing direction of both e1 and e1 until at point
2352 // (F1, G1) there is a l3 so that e3(F1, G1)=max(e1, e2)(F1, G1), or
2353 // at a the search reaches a vertex of R.
2354 // Arguments: l1, l2, dF, dG, (F0, G0, vmax)
2355 // (3) e1 free search: starting with e1 being the only maximum, search down the decreasing
2356 // direction of e1 until at point (F1, G1) there is another l2 so that
2357 // e1(F1, G1)=e2(F1, G1), or until at point (F1, G1) the search is about
2358 // to go out of R.
2359 // Arguments: l1, dF, dG, (F0, G0, vmax)
2360 // (4) e1 side search: starting with e1 being the only maximum, search down the
2361 // a side of R in the decreasing direction of e1 until at point (F1, G1) there
2362 // is another l2 so that e1=e2, or until point (F1, G1) the search reaches
2363 // a vertex.
2364 // Arguments: l1, destimation vertex nd, (F0, G0, vmax)
2365 //
2366 // At the beginning of each searching step we check the conditions to see which
2367 // basic step to perform, and provide the starting conditions. At the very start of
2368 // the search, (F0, G0) is a vertex of R, e_l1 being the maximum at this point.
2369 //
2370
2371 int condpos=3; //inside, on side, vertex
2372 int condmax=3; //triple max, duo max, solo max
2373 int searchmode;
2374
2375 bool minimax=false; //set to true when the minimal maximum is met
2376
2377 int iter=L*2;
2378 while (!minimax && iter>0)
2379 {
2380 iter--;
2381 double m1=m[l1], m2=m[l2], m3=m[l3], k1=k[l1], k2=k[l2], k3=k[l3];
2382 double tmp, tmp1, tmp2, tmp3;
2383
2384 switch (condmax)
2385 {
2386 case 1:
2387 tmp=F0+k3*G0; testnn(tmp);
2388 tmp3=sqrt(tmp);
2389 case 2:
2390 tmp=F0+k2*G0; testnn(tmp)
2391 tmp2=sqrt(tmp);
2392 case 3:
2393 tmp=F0+k1*G0; testnn(tmp);
2394 tmp1=sqrt(tmp);
2395 }
2396 double dF, dG;
2397 int n0=(nc==0)?(N-1):(nc-1), n1=(nc==N-1)?0:(nc+1);
2398 double x0, y0, x1, y1;
2399 if (n0>=0) x0=F[nc]-F[n0], y0=G[nc]-G[n0];
2400 if (nc>=0) x1=F[n1]-F[nc], y1=G[n1]-G[nc];
2401
2402 if (condpos==1) //(F0, G0) being inside polygon
2403 {
2404 if (condmax==1) //e1=e2=e3 being the maximum
2405 {
2406 //vmax holds the maximum
2407 //l1, l2, l3
2408
2409 //now choose a searching direction, either e1=e2 or e1=e3 or e2=e3.
2410 // choose e1=e2 if (e3-e1) decreases along the decreasing direction of e1=e2
2411 // choose e1=e3 if (e2-e1) decreases along the decreasing direction of e1=e3
2412 // choose e2=e3 if (e1-e2) decreases along the decreasing directino of e2=e3
2413 // if no condition is satisfied, then (F0, G0) is the minimal maximum.
2414 //calculate the decreasing direction of e1=e3 as (dF, dG)
2415 double den=m1*m3*(k1-k3)/2;
2416 dF=-(k1*m1*tmp3-k3*m3*tmp1)/den,
2417 dG=(m1*tmp3-m3*tmp1)/den;
2418 //the negative gradient of e2-e1 is calculated as (gF, gG)
2419 double gF=-(m2/tmp2-m1/tmp1)/2,
2420 gG=-(k2*m2/tmp2-k1*m1/tmp1)/2;
2421 if (dF*gF+dG*gG>0) //so that e2-e1 decreases in the decreasing direction of e1=e3
2422 {
2423 l2=l3;
2424 searchmode=1;
2425 }
2426 else
2427 {
2428 //calculate the decreasing direction of e2=e3 as (dF, dG)
2429 den=m2*m3*(k2-k3)/2;
2430 dF=-(k2*m2*tmp3-k3*m3*tmp2)/den,
2431 dG=(m2*tmp3-m3*tmp2)/den;
2432 //calculate the negative gradient of e1-e2 as (gF, gG)
2433 gF=-gF, gG=-gG;
2434 if (dF*gF+dG*gG>0) //so that e1-e2 decreases in the decreasing direction of e2=e3
2435 {
2436 l1=l3;
2437 searchmode=1;
2438 }
2439 else
2440 {
2441 //calculate the decreasing direction of e1=e2 as (dF, dG)
2442 den=m1*m2*(k1-k2)/2;
2443 dF=-(k1*m1*tmp2-k2*m2*tmp1)/den,
2444 dG=(m1*tmp2-m2*tmp1)/den;
2445 //calculate the negative gradient of (e3-e1) as (gF, gG)
2446 gF=-(m3/tmp3-m1/tmp1)/2,
2447 gG=-(k3*m3/tmp3-k1*m1/tmp1)/2;
2448 if (dF*gF+dG*gG>0) //so that e3-e1 decreases in the decreasing direction of e1=e2
2449 {
2450 searchmode=1;
2451 }
2452 else
2453 {
2454 F1=F0, G1=G0;
2455 searchmode=0; //no search
2456 minimax=true; //quit loop
2457 }
2458 }
2459 }
2460 } //
2461 else if (condmax==2) //e1=e2 being the maximum
2462 {
2463 //vmax holds the maximum
2464 //l1, l2
2465 searchmode=1;
2466 }
2467 else if (condmax==3) //e1 being the maximum
2468 {
2469 //the negative gradient of e1
2470 dF=-0.5*m1/tmp1;
2471 dG=k1*dF;
2472 searchmode=3;
2473 }
2474 }
2475 else if (condpos==2) //(F0, G0) being on side nc:(nc+1)
2476 {
2477 //the vector nc->(nc+1) as (x1, y1)
2478 if (condmax==1) //e1=e2=e3 being the maximum
2479 {
2480 //This case rarely happens.
2481
2482 //First see if a e1=e2 search is possible
2483 //calculate the decreasing direction of e1=e3 as (dF, dG)
2484 double den=m1*m3*(k1-k3)/2;
2485 double dF=-(k1*m1*tmp3-k3*m3*tmp1)/den, dG=(m1*tmp3-m3*tmp1)/den;
2486 //the negative gradient of e2-e1 is calculated as (gF, gG)
2487 double gF=-(m2/tmp2-m1/tmp1)/2, gG=-(k2*m2/tmp2-k1*m1/tmp1)/2;
2488 if (dF*gF+dG*gG>0 && x1*dG-y1*dF<0) //so that e2-e1 decreases in the decreasing direction of e1=e3
2489 { //~~~~~~~~~~~~~and this direction points inward
2490 l2=l3;
2491 searchmode=1;
2492 }
2493 else
2494 {
2495 //calculate the decreasing direction of e2=e3 as (dF, dG)
2496 den=m2*m3*(k2-k3)/2;
2497 dF=-(k2*m2*tmp3-k3*m3*tmp2)/den,
2498 dG=(m2*tmp3-m3*tmp2)/den;
2499 //calculate the negative gradient of e1-e2 as (gF, gG)
2500 gF=-gF, gG=-gG;
2501 if (dF*gF+dG*gG>0 && x1*dG-y1*dF<0) //so that e1-e2 decreases in the decreasing direction of e2=e3
2502 {
2503 l1=l3;
2504 searchmode=1;
2505 }
2506 else
2507 {
2508 //calculate the decreasing direction of e1=e2 as (dF, dG)
2509 den=m1*m2*(k1-k2)/2;
2510 dF=-(k1*m1*tmp2-k2*m2*tmp1)/den,
2511 dG=(m1*tmp2-m2*tmp1)/den;
2512 //calculate the negative gradient of (e3-e1) as (gF, gG)
2513 gF=-(m3/tmp3-m1/tmp1)/2,
2514 gG=-(k3*m3/tmp3-k1*m1/tmp1)/2;
2515 if (dF*gF+dG*gG>0 && x1*dG-y1*dF<0) //so that e3-e1 decreases in the decreasing direction of e1=e2
2516 {
2517 searchmode=1;
2518 }
2519 else
2520 {
2521 //see the possibility of a e1!=e2 search
2522 //calcualte the dot product of the gradients and (x1, y1)
2523 double d1=m1/2/tmp1*(x1+k1*y1),
2524 d2=m2/2/tmp2*(x1+k2*y1),
2525 d3=m3/2/tmp3*(x1+k3*y1);
2526 //we can prove that if there is a direction pointing inward R in which
2527 // e1, e2, e2 decrease, and another direction pointing outside R in
2528 // which e1, e2, e3 decrease, then on one direction along the side
2529 // all the three decrease. (Even more, this direction must be inside
2530 // the <180 angle formed by the two directions.)
2531 //
2532 // On the contrary, if there is a direction
2533 // in which all the three decrease, with two equal, it has to point
2534 // outward for the program to get here. Then if along neither direction
2535 // of side R can the three all descend, then there doesn't exist any
2536 // direction inward R in which the three descend. In that case we
2537 // have a minimal maximum at (F0, G0).
2538 if (d1*d2<=0 || d1*d3<=0 || d2*d3<=0) //so that the three don't decrease in the same direction
2539 {
2540 F1=F0, G1=G0;
2541 searchmode=0; //no search
2542 minimax=true; //quit loop
2543 }
2544 else
2545 {
2546 if (d1>0) //so that d2>0, d3>0, all three decrease in the direction -(x1, y1) towards nc
2547 {
2548 if (d1>d2 && d1>d3) //e1 decreases fastest
2549 {
2550 l1=l3; //keep e2, e3
2551 }
2552 else if (d2>=d1 && d2>d3) //e2 decreases fastest
2553 {
2554 l2=l3; //keep e1, e3
2555 }
2556 else //d3>=d1 && d3>=d2, e3 decreases fastest
2557 {
2558 //keep e1, e2
2559 }
2560 nd=nc;
2561 }
2562 else //d1<0, d2<0, d3<0, all three decrease in the direction (x1, y1)
2563 {
2564 if (d1<d2 && d1<d3) //e1 decreases fastest
2565 {
2566 l1=l3;
2567 }
2568 else if (d2<=d1 && d2<d3) //e2 decreases fastest
2569 {
2570 l2=l3;
2571 }
2572 else //d3<=d1 && d3<=d2
2573 {
2574 //keep e1, e2
2575 }
2576 nd=n1;
2577 }
2578 searchmode=2;
2579 }
2580 }
2581 }
2582 }
2583 }
2584 else if (condmax==2) //e1=e2 being the maximum
2585 {
2586 //first see if the decreasing direction of e1=e2 points to the inside
2587 //calculate the decreasing direction of e1=e2 as (dF, dG)
2588 double den=m1*m2*(k1-k2)/2;
2589 dF=-(k1*m1*tmp2-k2*m2*tmp1)/den,
2590 dG=(m1*tmp2-m2*tmp1)/den;
2591
2592 if (x1*dG-y1*dF<0) //so that (dF, dG) points inward R
2593 {
2594 searchmode=1;
2595 }
2596 else
2597 {
2598 //calcualte the dot product of the gradients and (x1, y1)
2599 double d1=m1/2/tmp1*(x1+k1*y1),
2600 d2=m2/2/tmp2*(x1+k2*y1);
2601 if (d1*d2<=0) //so that along the side e1, e2 descend in opposite directions
2602 {
2603 F1=F0, G1=G0;
2604 searchmode=0;
2605 minimax=true;
2606 }
2607 else
2608 {
2609 if (d1>0) //so that both decrease in direction -(x1, y1)
2610 nd=nc;
2611 else
2612 nd=n1;
2613 searchmode=2;
2614 }
2615 }
2616 }
2617 else if (condmax==3) //e1 being the maximum
2618 {
2619 //calculate the negative gradient of e1 as (dF, dG)
2620 dF=-0.5*m1/tmp1;
2621 dG=k1*dF;
2622
2623 if (x1*dG-y1*dF<0) //so the gradient points inward R
2624 searchmode=3;
2625 else
2626 {
2627 //calculate the dot product of the gradient and (x1, y1)
2628 double d1=m1/2/tmp1*(x1+k1*y1);
2629 if (d1>0) //so that e1 decreases in direction -(x1, y1)
2630 {
2631 nd=nc;
2632 searchmode=4;
2633 }
2634 else if (d1<0)
2635 {
2636 nd=n1;
2637 searchmode=4;
2638 }
2639 else //so that e1 does not change along side nc:(nc+1)
2640 {
2641 F1=F0, G1=G0;
2642 searchmode=0;
2643 minimax=true;
2644 }
2645 }
2646 }
2647 }
2648 else //condpos==3, (F0, G0) being vertex nc
2649 {
2650 //the vector nc->(nc+1) as (x1, y1)
2651 //the vector (nc-1)->nc as (x0, y0)
2652
2653 if (condmax==1) //e1=e2=e3 being the maximum
2654 {
2655 //This case rarely happens.
2656
2657 //First see if a e1=e2 search is possible
2658 //calculate the decreasing direction of e1=e3 as (dF, dG)
2659 double den=m1*m3*(k1-k3)/2;
2660 double dF=-(k1*m1*tmp3-k3*m3*tmp1)/den, dG=(m1*tmp3-m3*tmp1)/den;
2661 //the negative gradient of e2-e1 is calculated as (gF, gG)
2662 double gF=-(m2/tmp2-m1/tmp1)/2, gG=-(k2*m2/tmp2-k1*m1/tmp1)/2;
2663 if (dF*gF+dG*gG>0 && x1*dG-y1*dF<0 && x0*dG-y0*dF<0) //so that e2-e1 decreases in the decreasing direction of e1=e3
2664 { //~~~~~~~~~~~~~and this direction points inward
2665 l2=l3;
2666 searchmode=1;
2667 }
2668 else
2669 {
2670 //calculate the decreasing direction of e2=e3 as (dF, dG)
2671 den=m2*m3*(k2-k3)/2;
2672 dF=-(k2*m2*tmp3-k3*m3*tmp2)/den,
2673 dG=(m2*tmp3-m3*tmp2)/den;
2674 //calculate the negative gradient of e1-e2 as (gF, gG)
2675 gF=-gF, gG=-gG;
2676 if (dF*gF+dG*gG>0 && x1*dG-y1*dF<0 && x0*dG-y0*dF<0) //so that e1-e2 decreases in the decreasing direction of e2=e3
2677 {
2678 l1=l3;
2679 searchmode=1;
2680 }
2681 else
2682 {
2683 //calculate the decreasing direction of e1=e2 as (dF, dG)
2684 den=m1*m2*(k1-k2)/2;
2685 dF=-(k1*m1*tmp2-k2*m2*tmp1)/den,
2686 dG=(m1*tmp2-m2*tmp1)/den;
2687 //calculate the negative gradient of (e3-e1) as (gF, gG)
2688 gF=-(m3/tmp3-m1/tmp1)/2,
2689 gG=-(k3*m3/tmp3-k1*m1/tmp1)/2;
2690 if (dF*gF+dG*gG>0 && x1*dG-y1*dF<0 && x0*dG-y0*dF<0) //so that e3-e1 decreases in the decreasing direction of e1=e2
2691 {
2692 searchmode=1;
2693 }
2694 else
2695 {
2696 //see the possibility of a e1!=e2 search
2697 //calcualte the dot product of the gradients and (x1, y1)
2698 double d1=m1/2/tmp1*(x1+k1*y1),
2699 d2=m2/2/tmp2*(x1+k2*y1),
2700 d3=m3/2/tmp3*(x1+k3*y1);
2701 //we can also prove that if there is a direction pointing inward R in which
2702 // e1, e2, e2 decrease, and another direction pointing outside R in
2703 // which e1, e2, e3 decrease, all the three decrease either on nc->(nc+1)
2704 // direction along side nc:(nc+1), or on nc->(nc-1) direction along side
2705 // nc:(nc-1).
2706 //
2707 // On the contrary, if there is a direction
2708 // in which all the three decrease, with two equal, it has to point
2709 // outward for the program to get here. Then if along neither direction
2710 // of side R can the three all descend, then there doesn't exist any
2711 // direction inward R in which the three descend. In that case we
2712 // have a minimal maximum at (F0, G0).
2713 if (d1<0 && d2<0 && d3<0) //so that all the three decrease in the direction (x1, y1)
2714 {
2715 if (d1<d2 && d1<d3) //e1 decreases fastest
2716 {
2717 l1=l3;
2718 }
2719 else if (d2<=d1 && d2<d3) //e2 decreases fastest
2720 {
2721 l2=l3;
2722 }
2723 else //d3<=d1 && d3<=d2
2724 {
2725 //keep e1, e2
2726 }
2727 nd=n1;
2728 searchmode=2;
2729 }
2730 else
2731 {
2732 d1=m1/2/tmp1*(x0+k1*y0),
2733 d2=m2/2/tmp2*(x0+k2*y0),
2734 d3=m3/2/tmp3*(x0+k3*y0);
2735 if (d1>0 && d2>0 && d3>0) //so that all the three decrease in the direction -(x0, y0)
2736 {
2737 if (d1>d2 && d1>d3) //e1 decreases fastest
2738 {
2739 l1=l3; //keep e2, e3
2740 }
2741 else if (d2>=d1 && d2>d3) //e2 decreases fastest
2742 {
2743 l2=l3; //keep e1, e3
2744 }
2745 else //d3>=d1 && d3>=d2, e3 decreases fastest
2746 {
2747 //keep e1, e2
2748 }
2749 nd=n0;
2750 searchmode=2;
2751 }
2752 else
2753 {
2754 F1=F0, G1=G0;
2755 searchmode=0;
2756 minimax=true;
2757 }
2758 }
2759 }
2760 }
2761 }
2762 }
2763 else if (condmax==2) //e1=e2 being the maximum
2764 {
2765 //first see if the decreasing direction of e1=e2 points to the inside
2766 //calculate the decreasing direction of e1=e2 as (dF, dG)
2767 double den=m1*m2*(k1-k2)/2;
2768 dF=-(k1*m1*tmp2-k2*m2*tmp1)/den,
2769 dG=(m1*tmp2-m2*tmp1)/den;
2770
2771 if (x1*dG-y1*dF<0 && x0*dG-y1*dF<0) //so that (dF, dG) points inward R
2772 {
2773 searchmode=1;
2774 }
2775 else
2776 {
2777 //calcualte the dot product of the gradients and (x1, y1)
2778 double d1=m1/2/tmp1*(x1+k1*y1),
2779 d2=m2/2/tmp2*(x1+k2*y1);
2780 if (d1<0 && d2<0) //so that along side nc:(nc+1) e1, e2 descend in direction (x1, y1)
2781 {
2782 nd=n1;
2783 searchmode=2;
2784 }
2785 else
2786 {
2787 d1=m1/2/tmp1*(x0+k1*y0),
2788 d2=m2/2/tmp2*(x0+k2*y0);
2789 if (d1>0 && d2>0) //so that slong the side (nc-1):nc e1, e2 decend in direction -(x0, y0)
2790 {
2791 nd=n0;
2792 searchmode=2;
2793 }
2794 else
2795 {
2796 F1=F0, G1=G0;
2797 searchmode=0;
2798 minimax=true;
2799 }
2800 }
2801 }
2802 }
2803 else //condmax==3, e1 being the solo maximum
2804 {
2805 //calculate the negative gradient of e1 as (dF, dG)
2806 dF=-0.5*m1/tmp1;
2807 dG=k1*dF;
2808
2809 if (x1*dG-y1*dF<0 && x0*dG-y0*dF<0) //so the gradient points inward R
2810 searchmode=3;
2811 else
2812 {
2813 //calculate the dot product of the gradient and (x1, y1)
2814 double d1=m1/2/tmp1*(x1+k1*y1),
2815 d0=-m1/2/tmp1*(x0+k1*y0);
2816 if (d1<d0) //so that e1 decreases in direction (x1, y1) faster than in -(x0, y0)
2817 {
2818 if (d1<0)
2819 {
2820 nd=n1;
2821 searchmode=4;
2822 }
2823 else
2824 {
2825 F1=F0, G1=G0;
2826 searchmode=0;
2827 minimax=true;
2828 }
2829 }
2830 else //so that e1 decreases in direction -(x0, y0) faster than in (x1, y1)
2831 {
2832 if (d0<0)
2833 {
2834 nd=n0;
2835 searchmode=4;
2836 }
2837 else
2838 {
2839 F1=F0, G1=G0;
2840 searchmode=0;
2841 minimax=true;
2842 }
2843 }
2844 }
2845 }
2846 }
2847 // the conditioning stage ends here
2848 // the searching step starts here
2849 if (searchmode==0)
2850 {}
2851 else if (searchmode==1)
2852 {
2853 //Search mode 1: e1=e2 search
2854 // starting with e1=e2, search down the decreasing direction of
2855 // e1=e2 until at point (F1, G1) there is another l3 so that e1(F1, G1)
2856 // =e2(F1, G1)=e3(F1, G1), or until at point (F1, G1) the search is about
2857 // to go out of R.
2858 //Arguments: l1, l2, (F0, G0, vmax)
2859 double tmp=F0+k[l1]*G0; testnn(tmp);
2860 double e1=m[l1]*sqrt(tmp)-f[l1];
2861 tmp=F0+k[l2]*G0; testnn(tmp);
2862 double e2=m[l2]*sqrt(tmp)-f[l2];
2863 if (fabs(e1-e2)>1e-4)
2864 {
2865 // int kk=1;
2866 goto start;
2867 }
2868
2869 //find intersections of e1=e2 with other ek's
2870 l3=-1;
2871 loopl3:
2872 double e=-1;
2873 for (int l=0; l<L; l++)
2874 {
2875 if (l==l1 || l==l2) continue;
2876 double lF1[2], lG1[2], lE[2];
2877 double lm[3]={m[l1], m[l2], m[l]};
2878 double lf[3]={f[l1], f[l2], f[l]};
2879 int lk[3]={k[l1], k[l2], k[l]};
2880 int r=FGFrom3(lF1, lG1, lE, lm, lf, lk);
2881 for (int i=0; i<r; i++)
2882 {
2883 if (lE[i]>e && lE[i]<vmax) l3=l, e=lE[i], F1=lF1[i], G1=lG1[i];
2884 }
2885 }
2886 //l3 shall be >=0 at this point
2887 //find intersections of e1=e2 with polygon sides
2888 if (l3<0)
2889 {
2890 ///int kkk=1;
2891 goto loopl3;
2892 }
2893 nc=-1;
2894 double F2, G2;
2895 for (int n=0; n<N; n++)
2896 {
2897 int nn=(n==N-1)?0:(n+1);
2898 double xn1=F[nn]-F[n], yn1=G[nn]-G[n], xn2=F1-F[n], yn2=G1-G[n];
2899 if (xn1*yn2-xn2*yn1>=0) //so that (F1, G1) is on or out of side n:(n+1)
2900 {
2901 nc=n;
2902 break;
2903 }
2904 }
2905
2906 if (nc<0) //(F1, G1) being inside R
2907 {
2908 vmax=e;
2909 condpos=1;
2910 condmax=1;
2911 //l3 shall be >=0 at this point, so l1, l2, l3 are all ok
2912 }
2913 else //(F1, G1) being outside nc:(nc+1) //(F2, G2) being on side nc:(nc+1)
2914 {
2915 //find where e1=e2 crosses a side of R between (F0, G0) and (F1, G1)
2916 double lF2[2], lG2[2], lmd[2];
2917 double lm[2]={m[l1], m[l2]};
2918 double lf[2]={f[l1], f[l2]};
2919 int lk[2]={k[l1], k[l2]};
2920
2921 int ncc=-1;
2922 while (ncc<0)
2923 {
2924 n1=(nc==N-1)?0:(nc+1);
2925 double xn1=F[n1]-F[nc], yn1=G[n1]-G[nc];
2926
2927 int r=FGFrom2(lF2, lG2, lmd, F[nc], xn1, G[nc], yn1, lm, lf, lk);
2928
2929 bool once=false;
2930 for (int i=0; i<r; i++)
2931 {
2932 double tmp=lF2[i]+k[l1]*lG2[i]; testnn(tmp);
2933 double le=m[l1]*sqrt(tmp)-f[l1];
2934 if (le<=vmax) //this shall be satisfied once
2935 {
2936 once=true;
2937 if (lmd[i]>=0 && lmd[i]<=1) //so that curve e1=e2 does cross the boundry within it
2938 ncc=nc, F2=lF2[i], G2=lG2[i], e=le;
2939 else if (lmd[i]>1) //so that the crossing point is out of vertex n1
2940 nc=n1;
2941 else //lmd[i]<0, so that the crossing point is out of vertex nc-1
2942 nc=(nc==0)?(N-1):(nc-1);
2943 }
2944 }
2945 if (!once)
2946 {
2947 //int kk=0;
2948 goto start;
2949 }
2950 }
2951 nc=ncc;
2952
2953 vmax=e;
2954 condpos=2;
2955 if (F1==F2 && G1==G2)
2956 condmax=1;
2957 else
2958 condmax=2;
2959 F1=F2, G1=G2;
2960 }
2961 }
2962 else if (searchmode==2)
2963 {
2964 // Search mode 2 (e1!=e2 search): starting with e1=e2, and a linear equation
2965 // constraint, search down the decreasing direction until at point (F1, G1)
2966 // there is a l3 so that e3(F1, G1)=max(e1, e2)(F1, G1), or at point (F1, G1)
2967 // the search is about to go out of R.
2968 // Arguments: l1, l2, dF, dG, (F0, G0, vmax)
2969
2970 //The searching direction (dF, dG) is always along some polygon side
2971 // If conpos==2, it is along nc:(nc+1)
2972 // If conpos==3, it is along nc:(nc+1) or (nc-1):nc
2973
2974 //
2975
2976 //first check whether e1>e2 or e2>e1 down (dF, dG)
2977 dF=F[nd]-F0, dG=G[nd]-G0;
2978 //
2979 //the negative gradient of e2-e1 is calculated as (gF, gG)
2980 double gF=-(m2/tmp2-m1/tmp1)/2,
2981 gG=-(k2*m2/tmp2-k1*m1/tmp1)/2;
2982 if (gF*dF+gG*dG>0) //e2-e1 decrease down direction (dF, dG), so that e1>e2
2983 {}
2984 else //e1<e2 down (dF, dG)
2985 {
2986 int ll=l2; l2=l1; l1=ll;
2987 m1=m[l1], m2=m[l2], k1=k[l1], k2=k[l2];
2988 tmp=F0+k1*G0; testnn(tmp); tmp1=sqrt(tmp);
2989 tmp=F0+k2*G0; testnn(tmp); tmp2=sqrt(tmp);
2990 }
2991 //now e1>e2 down (dF, dG) from (F0, G0)
2992 tmp=F[nd]+k1*G[nd]; testnn(tmp);
2993 double e1=m1*sqrt(tmp)-f[l1];
2994 tmp=F[nd]+k2*G[nd]; testnn(tmp);
2995 double e2=m2*sqrt(tmp)-f[l2], FEx, GEx, eEx;
2996 int maxlmd=1;
2997 if (e1>=e2)
2998 {
2999 // then e1 is larger than e2 for the whole searching interval
3000 FEx=F[nd], GEx=G[nd];
3001 eEx=e1;
3002 }
3003 else // the they swap somewhere in between, now find it and make it maxlmd
3004 {
3005 double lF1[2], lG1[2], lmd[2];
3006 double lm[2]={m[l1], m[l2]};
3007 double lf[2]={f[l1], f[l2]};
3008 int lk[2]={k[l1], k[l2]};
3009 int r=FGFrom2(lF1, lG1, lmd, F0, dF, G0, dG, lm, lf, lk);
3010 //process the results
3011 if (r==2) //must be so
3012 {
3013 if (lmd[0]<lmd[1]) lmd[0]=lmd[1], FEx=lF1[1], GEx=lG1[1];
3014 else FEx=lF1[0], GEx=lG1[0];
3015 }
3016 if (lmd[0]>0 && lmd[0]<1) //must be so
3017 maxlmd=lmd[0];
3018 tmp=FEx+k1*GEx; testnn(tmp);
3019 eEx=m1*sqrt(tmp)-f[l1];
3020 }
3021
3022 l3=-1;
3023 // int l4;
3024 double e, llmd=maxlmd;
3025 int *tpl=new int[L], ctpl=0;
3026 for (int l=0; l<L; l++)
3027 {
3028 if (l==l1 || l==l2) continue;
3029 tmp=FEx+k[l]*GEx; testnn(tmp);
3030 double le=m[l]*sqrt(tmp)-f[l];
3031 if (le<eEx)
3032 {
3033 tpl[ctpl]=l;
3034 ctpl++;
3035 continue;
3036 }
3037 //solve for the equation system e1(F1, G1)=el(F1, G1), with (F1, G1) on nc:n1
3038 double lF1[2], lG1[2], lmd[2];
3039 double lm[2]={m[l1], m[l]};
3040 double lf[2]={f[l1], f[l]};
3041 int lk[2]={k[l1], k[l]};
3042 int r=FGFrom2(lF1, lG1, lmd, F0, dF, G0, dG, lm, lf, lk);
3043 //process the results
3044 for (int i=0; i<r; i++)
3045 {
3046 if (lmd[i]<0 || lmd[i]>maxlmd) continue; //the solution is in the wrong direction
3047 if (lmd[i]<llmd)
3048 l3=l, F1=lF1[i], G1=lG1[i], llmd=lmd[i];
3049 }
3050 }
3051 if (ctpl==L-2) //so that at (FEx, GEx) e1 is maximal, equivalent to "l3<0"
3052 {}
3053 else
3054 {
3055 //at this point F1 and G1 must have valid values already
3056 tmp=F1+k[l1]*G1; testnn(tmp);
3057 e=m[l1]*sqrt(tmp)-f[l1];
3058 //test if e1 is maximal at (F1, G1)
3059 bool lmax=true;
3060 for (int tp=0; tp<ctpl; tp++)
3061 {
3062 tmp=F1+k[tpl[tp]]*G1; testnn(tmp);
3063 double le=m[tpl[tp]]*sqrt(tmp)-f[tpl[tp]];
3064 if (le>e) {lmax=false; break;}
3065 }
3066 if (!lmax)
3067 {
3068 for (int tp=0; tp<ctpl; tp++)
3069 {
3070 double lF1[2], lG1[2], lmd[2];
3071 double lm[2]={m[l1], m[tpl[tp]]};
3072 double lf[2]={f[l1], f[tpl[tp]]};
3073 int lk[2]={k[l1], k[tpl[tp]]};
3074 int r=FGFrom2(lF1, lG1, lmd, F0, dF, G0, dG, lm, lf, lk);
3075 //process the results
3076 for (int i=0; i<r; i++)
3077 {
3078 if (lmd[i]<0 || lmd[i]>maxlmd) continue; //the solution is in the wrong direction
3079 if (lmd[i]<=llmd)
3080 l3=tpl[tp], F1=lF1[i], G1=lG1[i], llmd=lmd[i];
3081 }
3082 }
3083 tmp=F1+k[l1]*G1; testnn(tmp);
3084 e=m[l1]*sqrt(tmp)-f[l1];
3085 }
3086 }
3087 delete[] tpl;
3088 if (l3>=0) //the solution is found in the right direction
3089 {
3090 vmax=e;
3091 if (llmd==1) //(F1, G1) is a vertex of R
3092 {
3093 nc=nd;
3094 condpos=3;
3095 condmax=2;
3096 l2=l3;
3097 }
3098 else
3099 {
3100 if (condpos==3 && nd==n0) nc=n0;
3101 condpos=2;
3102 condmax=2;
3103 l2=l3;
3104 }
3105 }
3106 else if (maxlmd<1) //so that e1=e2 remains maximal at maxlmd
3107 {
3108 if (condpos==3 && nd==n0) nc=n0;
3109 condpos=2;
3110 condmax=2;
3111 //l1, l2 remain unchanged
3112 F1=FEx, G1=GEx;
3113 vmax=eEx;
3114 }
3115 else //so that e1 remains maximal at vertex nd
3116 {
3117 condpos=3;
3118 condmax=3;
3119 nc=nd;
3120 F1=FEx, G1=GEx; //this shall equal to F0+dF, G0+dG
3121 vmax=eEx;
3122 }
3123 }
3124 else if (searchmode==3)
3125 {
3126 //Search mode 3 (e1 search): starting with e1 being the only maximum, search down
3127 // the decreasing direction of e1 until at point (F1, G1) there is another l2 so
3128 // that e1(F1, G1)=e2(F1, G1), or until at point (F1, G1) the search is about
3129 // to go out of R.
3130 //
3131 // Arguments: l1, dF, dG, (F0, G0, vmax)
3132 //
3133
3134 //first calculate the value of lmd from (F0, G0) to get out of R
3135 double maxlmd=-1, FEx, GEx;
3136 double fdggdf=-F0*dG+G0*dF;
3137 nd=-1;
3138 for (int n=0; n<N; n++)
3139 {
3140 if (condpos==2 && n==nc) continue;
3141 if ((condpos==3 && n==nc) || n==n0) continue;
3142 int nn=(n==N-1)?0:n+1;
3143 double xn1=F[n], xn2=F[nn], yn1=G[n], yn2=G[nn];
3144 double test1=xn1*dG-yn1*dF+fdggdf,
3145 test2=xn2*dG-yn2*dF+fdggdf;
3146 if (test1*test2<=0)
3147 {
3148 //the two following are equivalent. we use the larger denominator.
3149 FEx=(xn1*test2-xn2*test1)/(test2-test1);
3150 GEx=(yn1*test2-yn2*test1)/(test2-test1);
3151 if (fabs(dF)>fabs(dG)) maxlmd=(FEx-F0)/dF;
3152 else maxlmd=(GEx-G0)/dG;
3153 nd=n;
3154 }
3155 }
3156
3157 //maxlmd must be >0 at this point.
3158 l2=-1;
3159 tmp=FEx+k[l1]*GEx; testnn(tmp);
3160 double e, eEx=m[l1]*sqrt(tmp)-f[l1], llmd=maxlmd;
3161 int *tpl=new int[L], ctpl=0;
3162 for (int l=0; l<L; l++)
3163 {
3164 if (l==l1) continue;
3165
3166 //if el does not catch up with e1 at (FEx, GEx), then there is no hope this
3167 // happens in between (F0, G0) and (FEx, GEx)
3168 tmp=FEx+k[l]*GEx; testnn(tmp);
3169 double le=m[l]*sqrt(tmp)-f[l];
3170 if (le<eEx)
3171 {
3172 tpl[ctpl]=l;
3173 ctpl++;
3174 continue;
3175 }
3176 //now the existence of (F1, G1) is guaranteed
3177 double lF1[2], lG1[2], lmd[2];
3178 double lm[2]={m[l1], m[l]};
3179 double lf[2]={f[l1], f[l]};
3180 int lk[2]={k[l1], k[l]};
3181 int r=FGFrom2(lF1, lG1, lmd, F0, dF, G0, dG, lm, lf, lk);
3182 for (int i=0; i<r; i++)
3183 {
3184 if (lmd[i]<0 || lmd[i]>maxlmd) continue;
3185 if (lmd[i]<=llmd) llmd=lmd[i], l2=l, F1=lF1[i], G1=lG1[i];
3186 }
3187 }
3188 if (ctpl==L-1) //e1 is maximal at eEx, equivalent to "l2<0"
3189 {}//l2 must equal -1 at this point
3190 else
3191 {
3192 tmp=F1+k[l1]*G1; testnn(tmp);
3193 e=m[l1]*sqrt(tmp)-f[l1];
3194 //test if e1 is maximal at (F1, G1)
3195 //e1 is already maximal of those l's not in tpl[ctpl]
3196 bool lmax=true;
3197 for (int tp=0; tp<ctpl; tp++)
3198 {
3199 tmp=F1+k[tpl[tp]]*G1; testnn(tmp);
3200 double le=m[tpl[tp]]*sqrt(tmp)-f[tpl[tp]];
3201 if (le>e) {lmax=false; break;}
3202 }
3203 if (!lmax)
3204 {
3205 for (int tp=0; tp<ctpl; tp++)
3206 {
3207 double lF1[2], lG1[2], lmd[2];
3208 double lm[2]={m[l1], m[tpl[tp]]};
3209 double lf[2]={f[l1], f[tpl[tp]]};
3210 int lk[2]={k[l1], k[tpl[tp]]};
3211 int r=FGFrom2(lF1, lG1, lmd, F0, dF, G0, dG, lm, lf, lk);
3212 for (int i=0; i<r; i++)
3213 {
3214 if (lmd[i]<0 || lmd[i]>maxlmd) continue;
3215 if (llmd>=lmd[i]) llmd=lmd[i], l2=tpl[tp], F1=lF1[i], G1=lG1[i];
3216 }
3217 }
3218 tmp=F1+k[l1]*G1; testnn(tmp);
3219 e=m[l1]*sqrt(tmp)-f[l1];
3220 }
3221 }
3222 delete[] tpl;
3223
3224 if (l2>=0) //so that e1 meets some other ek during the search
3225 {
3226 vmax=e; //this has to equal m[l2]*sqrt(F1+k[l2]*G1)-f[l2]
3227 if (vmax==eEx)
3228 {
3229 condpos=2;
3230 condmax=2;
3231 nc=nd;
3232 }
3233 else
3234 {
3235 condpos=1;
3236 condmax=2;
3237 }
3238 }
3239 else //l2==-1
3240 {
3241 vmax=eEx;
3242 F1=FEx, G1=GEx;
3243 condpos=2;
3244 condmax=3;
3245 nc=nd;
3246 }
3247 }
3248 else //searchmode==4
3249 {
3250 //Search mode 4: a special case of search mode 3 in which the boundry constraint
3251 // is simply 0<=lmd<=1.
3252 dF=F[nd]-F0, dG=G[nd]-G0;
3253 l2=-1;
3254 int *tpl=new int[L], ctpl=0;
3255 tmp=F[nd]+k[l1]*G[nd]; testnn(tmp);
3256 double e, eEx=m[l1]*sqrt(tmp)-f[l1], llmd=1;
3257 for (int l=0; l<L; l++)
3258 {
3259 if (l==l1) continue;
3260 //if el does not catch up with e1 at (F[nd], G[nd]), then there is no hope this
3261 // happens in between (F0, G0) and vertex nd.
3262 tmp=F[nd]+k[l]*G[nd]; testnn(tmp);
3263 double le=m[l]*sqrt(tmp)-f[l];
3264 if (le<eEx)
3265 {
3266 tpl[ctpl]=l;
3267 ctpl++;
3268 continue;
3269 }
3270 //now the existence of (F1, G1) is guaranteed
3271 double lF1[2], lG1[2], lmd[2];
3272 double lm[2]={m[l1], m[l]};
3273 double lf[2]={f[l1], f[l]};
3274 int lk[2]={k[l1], k[l]};
3275 loop1:
3276 int r=FGFrom2(lF1, lG1, lmd, F0, dF, G0, dG, lm, lf, lk);
3277 for (int i=0; i<r; i++)
3278 {
3279 if (lmd[i]<0 || lmd[i]>1) continue;
3280 if (llmd>=lmd[i])
3281 {
3282 tmp=lF1[i]+k[l1]*lG1[i]; testnn(tmp);
3283 double e1=m[l1]*sqrt(tmp)-f[l1];
3284 tmp=lF1[i]+k[l]*lG1[i]; testnn(tmp);
3285 double e2=m[l]*sqrt(tmp)-f[l];
3286 if (fabs(e1-e2)>1e-4)
3287 {
3288 //int kk=0;
3289 goto loop1;
3290 }
3291 llmd=lmd[i], l2=l, F1=lF1[i], G1=lG1[i];
3292 }
3293 }
3294 }
3295 if (ctpl==N-1) //so e1 is maximal at (F[nd], G[nd])
3296 {}
3297 else
3298 {
3299 tmp=F1+k[l1]*G1; testnn(tmp);
3300 e=m[l1]*sqrt(tmp)-f[l1];
3301 //test if e1 is maximal at (F1, G1)
3302 bool lmax=true;
3303 for (int tp=0; tp<ctpl; tp++)
3304 {
3305 tmp=F1+k[tpl[tp]]*G1; testnn(tmp);
3306 double le=m[tpl[tp]]*sqrt(tmp)-f[tpl[tp]];
3307 if (le>e) {lmax=false; break;}
3308 }
3309 if (!lmax)
3310 {
3311 for (int tp=0; tp<ctpl; tp++)
3312 {
3313 double lF1[2], lG1[2], lmd[2];
3314 double lm[2]={m[l1], m[tpl[tp]]};
3315 double lf[2]={f[l1], f[tpl[tp]]};
3316 int lk[2]={k[l1], k[tpl[tp]]};
3317 int r=FGFrom2(lF1, lG1, lmd, F0, dF, G0, dG, lm, lf, lk);
3318 for (int i=0; i<r; i++)
3319 {
3320 if (lmd[i]<0 || lmd[i]>1) continue;
3321 if (llmd>=lmd[i]) llmd=lmd[i], l2=tpl[tp], F1=lF1[i], G1=lG1[i];
3322 }
3323 }
3324 // e=m[l1]*sqrt(F1+k[l1]*G1)-f[l1];
3325 }
3326 }
3327 tmp=F1+k[l1]*G1; testnn(tmp);
3328 e=m[l1]*sqrt(tmp)-f[l1];
3329 delete[] tpl;
3330 if (l2>=0)
3331 {
3332 vmax=e;
3333 if (vmax==eEx)
3334 {
3335 condpos=3;
3336 condmax=2;
3337 nc=nd;
3338 }
3339 else
3340 {
3341 condpos=2;
3342 condmax=2;
3343 //(F1, G1) lies on side n0:nc (nd=n0) or nc:n1 (nd=nc or nd=n1)
3344 if (nd==n0) nc=n0;
3345 //else nc=nc;
3346 }
3347 }
3348 else //l2==-1,
3349 {
3350 vmax=eEx;
3351 condpos=3;
3352 condmax=3;
3353 F1=F[nd], G1=G[nd];
3354 nc=nd; //
3355 }
3356 }
3357 F0=F1, G0=G1;
3358
3359 if (condmax==1 && ((l1^l2)==1 || (l1^l3)==1 || (l2^l3)==1))
3360 minimax=true;
3361 else if (condmax==2 && ((l1^l2)==1))
3362 minimax=true;
3363 else if (vmax<1e-6)
3364 minimax=true;
3365 }
3366
3367 free(m);//delete[] m;
3368 delete[] vmnl;
3369
3370 return vmax;
3371 }//minimaxstiff*/
3372
3373 /*
3374 function NMResultToAtoms: converts the note-match result hosted in a NMResults structure, i.e.
3375 parameters of a harmonic atom, to a list of atoms. The process retrieves atom parameters from low
3376 partials until a required number of atoms have been retrieved or a non-positive frequency is
3377 encountered (non-positive frequencies are returned by note-match to indicate high partials for which
3378 no informative estimates can be obtained).
3379
3380 In: results: the NMResults structure hosting the harmonic atom
3381 M: number of atoms to request from &results
3382 t: time of this harmonic atom
3383 wid: scale of this harmonic atom with which the note-match was performed
3384 Out: HP[return value]: list of atoms retrieved from $result.
3385
3386 Returns the number of atoms retrieved.
3387 */
3388 int NMResultToAtoms(int M, atom* HP, int t, int wid, NMResults results)
3389 {
3390 double *fpp=results.fp, *vfpp=results.vfp, *pfpp=results.pfp;
3391 atomtype* ptype=results.ptype;
3392 for (int i=0; i<M; i++)
3393 {
3394 if (fpp[i]<=0) {memset(&HP[i], 0, sizeof(atom)*(M-i)); return i;}
3395 HP[i].t=t, HP[i].s=wid, HP[i].f=fpp[i]/wid, HP[i].a=vfpp[i];
3396 HP[i].p=pfpp[i]; HP[i].type=ptype[i]; HP[i].pin=i+1;
3397 }
3398 return M;
3399 }//NMResultToAtoms
3400
3401 /*
3402 function NMResultToPartials: reads atoms of a harmonic atom in a NMResult structure into HS partials.
3403
3404 In: results: the NMResults structure hosting the harmonic atom
3405 M: number of atoms to request from &results
3406 fr: frame index of this harmonic atom
3407 t: time of this harmonic atom
3408 wid: scale of this harmonic atom with which the note-match was performed
3409 Out: Partials[M][]: HS partials whose fr-th frame is updated to the harmonic partial read from $results
3410
3411 Returns the number of partials retrieved.
3412 */
3413 int NMResultToPartials(int M, int fr, atom** Partials, int t, int wid, NMResults results)
3414 {
3415 double *fpp=results.fp, *vfpp=results.vfp, *pfpp=results.pfp;
3416 atomtype* ptype=results.ptype;
3417 for (int i=0; i<M; i++)
3418 {
3419 atom* part=&Partials[i][fr];
3420 if (fpp[i]<=0) {for (int j=i; j<M; j++) memset(&Partials[j][fr], 0, sizeof(atom)); return i;}
3421 part->t=t, part->s=wid, part->f=fpp[i]/wid, part->a=vfpp[i];
3422 part->p=pfpp[i]; part->type=ptype[i]; part->pin=i+1;
3423 }
3424 return M;
3425 }//NMResultToPartials
3426
3427 /*
3428 function NoteMatchStiff3: finds harmonic atom from spectrum if Fr=1, or constant-pitch harmonic
3429 sinusoid from spectrogram if Fr>1.
3430
3431 In: x[Fr][N/2+1]: spectrogram
3432 fps[pc], vps[pc]: primitive (rough) peak frequencies and amplitudes
3433 N, offst: atom scale and hop size
3434 R: initial F-G polygon constraint, optional
3435 settings: note match settings
3436 computes: pointer to a function that computes HA score, must be ds0 or ds1
3437 lastvfp[lastp]: amplitude of the previous harmonic atom
3438 forceinputlocalfr: specifies if partial settings->pin0 is taken for granted ("pinned")
3439 Out: results: note match results
3440 f0, B: fundamental frequency and stiffness coefficient
3441 R: F-G polygon of harmonic atom
3442
3443 Returns the total energy of HA or constant-pitch HS.
3444 */
3445 double NoteMatchStiff3(TPolygon* R, double &f0, double& B, int pc, double* fps, double* vps,
3446 int Fr, cdouble** x, int N, int offst, NMSettings* settings, NMResults* results, int lastp,
3447 double* lastvfp, double (*computes)(double a, void* params), int forceinputlocalfr)
3448 {
3449 double result=0;
3450
3451 double maxB=settings->maxB, minf0=settings->minf0, maxf0=settings->maxf0,
3452 delm=settings->delm, delp=settings->delp, *c=settings->c, iH2=settings->iH2,
3453 *pinf=settings->pinf, hB=settings->hB, epf0=settings->epf0;// epf=settings->epf,
3454 int maxp=settings->maxp, *pin=settings->pin, pin0=settings->pin0, *pinfr=settings->pinfr,
3455 pcount=settings->pcount, M=settings->M;
3456
3457 // int *ptype=results->ptype;
3458
3459 //calculate the energy of the last frame (hp)
3460 double lastene=0; for (int i=0; i<lastp; i++) lastene+=lastvfp[i]*lastvfp[i];
3461 //fill frequency and amplitude buffers with 0
3462 //now determine numsam (the maximal number of partials)
3463 bool inputpartial=(f0>0 && settings->pin0>0);
3464 if (!inputpartial)
3465 {
3466 if (pcount>0) f0=pinf[0], pin0=pin[0];
3467 else f0=sqrt(R->X[0]), pin0=1;
3468 }
3469 int numsam=N*pin0/2.1/f0;
3470 if (numsam>maxp) numsam=maxp;
3471 //allocate buffer for candidate atoms
3472 TTempAtom*** Allocate2(TTempAtom*, numsam+1, MAX_CAND, cands);
3473 //pcs[p] records the number of candidates for partial index p
3474 //psb[p] = 1: anchor, 2: user input, 0: normal
3475 int *psb=new int[(numsam+1)*2], *pcs=&psb[numsam+1];
3476 memset(psb, 0, sizeof(int)*(numsam+1)*2);
3477 //if R is not specified, initialize a trivial candidate with maxB
3478 if (R->N==0) cands[0][0]=new TTempAtom(1, (minf0+maxf0)*0.5, (maxf0-minf0)*0.5, maxB);
3479 //if R is, initialize a trivial candidate by extending R by delp
3480 else cands[0][0]=new TTempAtom(R, delp, delp, minf0);
3481
3482 int pm=0, pM=0;
3483 int ind=1; pcs[0]=1;
3484 //anchor partials: highest priority atoms, must have spectral peaks
3485 for (int i=0; i<pcount; i++)
3486 {
3487 //skip out-of-scope atoms
3488 if (pin[i]>=numsam) continue;
3489 //skip user input atom
3490 if (inputpartial && pin[i]==pin0) continue;
3491
3492 void* dsparams;
3493 if (computes==ds0) dsparams=&cands[ind-1][0]->s;
3494 else
3495 {
3496 dsparams1 params={cands[ind-1][0]->s, lastene, cands[ind-1][0]->acce, pin[i], lastp, lastvfp};
3497 dsparams=&params;
3498 }
3499
3500 if (pinfr[i]>0) //local anchor
3501 {
3502 double ls=0;
3503 for (int fr=0; fr<Fr; fr++)
3504 {
3505 cdouble r=IPWindowC(pinf[i], x[fr], N, M, c, iH2, pinf[i]-hB, pinf[i]+hB);
3506 ls+=~r;
3507 }
3508 ls=sqrt(ls/Fr);
3509 cands[ind][0]=new TTempAtom(cands[ind-1][0], pin[i], pinf[i], delm, ls, computes(ls, dsparams));
3510 }
3511 else
3512 {
3513 //find the nearest peak to pinf[i]
3514 while (pm<pc-1 && fps[pm]<pinf[i]) pm++; while (pm>0 && fps[pm]>=pinf[i]) pm--;
3515 if (pm<pc-1 && pinf[i]-fps[pm]>fps[pm+1]-pinf[i]) pm++;
3516 cands[ind][0]=new TTempAtom(cands[ind-1][0], pin[i], fps[pm], delm, vps[pm], computes(vps[pm], dsparams));
3517 }
3518 delete cands[ind-1][0]->R; cands[ind-1][0]->R=0; psb[pin[i]]=1; pcs[ind]=1; ind++;
3519 }
3520 //harmonic atom grouping fails if anchors conflict
3521 if (cands[ind-1][0]->R->N==0) {f0=0; goto cleanup;}
3522
3523 //user input partial: must have spectral peak
3524 if (inputpartial)
3525 {
3526 //harmonic atom grouping fails if user partial is out of scope
3527 if (pin0>numsam){f0=0; goto cleanup;}
3528
3529 TTempAtom* lcand=cands[ind-1][0];
3530 //feasible frequency range for partial pin0
3531 double f1, f2; ExFmStiff(f1, f2, pin0, lcand->R->N, lcand->R->X, lcand->R->Y);
3532 f1-=delm, f2+=delm;
3533 if (f1<0) f1=0; if (f1>N/2.1) f1=N/2.1; if (f2>N/2.1) f2=N/2.1;
3534 //forcepin0 being true means this partial have to be the peak nearest to f0
3535 bool inputfound=false;
3536 if (forceinputlocalfr>=0)
3537 {
3538 int k1=floor(f0-epf0-1), k2=ceil(f0+epf0+1), k12=k2-k1+1, pk=-1;
3539
3540 cdouble *lx=x[forceinputlocalfr], *lr=new cdouble[k12];
3541 for (int k=0; k<k12; k++) lr[k]=IPWindowC(k1+k, lx, N, M, c, iH2, k1+k-hB, k1+k+hB);
3542 for (int k=1; k<k12-1; k++)
3543 if (~lr[k]>~lr[k-1] && ~lr[k]>~lr[k+1])
3544 {
3545 if (pk==-1 || fabs(k1+pk-f0)>fabs(k1+k-f0)) pk=k;
3546 }
3547 if (pk>0)
3548 {
3549 double lf=k1+pk, la, lph, oldlf=lf;
3550 LSESinusoid(lf, lf-1, lf+1, lx, N, hB, M, c, iH2, la, lph, settings->epf);
3551 if (fabs(lf-oldlf)>0.6)
3552 {
3553 lf=oldlf;
3554 cdouble r=IPWindowC(lf, lx, N, M, c, iH2, lf-hB, lf+hB);
3555 la=abs(r); lph=arg(r);
3556 }
3557 if (lf>f1 && lf<f2)
3558 {
3559 void* dsparams;
3560 if (computes==ds0) dsparams=&lcand->s;
3561 else
3562 {
3563 dsparams1 params={lcand->s, lastene, lcand->acce, pin0, lastp, lastvfp};
3564 dsparams=&params;
3565 }
3566 cands[ind][0]=new TTempAtom(lcand, pin0, lf, delm, la, computes(la, dsparams));
3567 pcs[ind]=1;
3568 inputfound=true;
3569 }
3570 }
3571 }
3572 if (!inputfound && settings->forcepin0)
3573 { //find the nearest peak to f0
3574 while (pm<pc-1 && fps[pm]<f0) pm++; while (pm>0 && fps[pm]>=f0) pm--;
3575 if (pm<pc-1 && f0-fps[pm]>fps[pm+1]-f0) pm++;
3576 if (fps[pm]>f1 && fps[pm]<f2)
3577 {
3578 void* dsparams;
3579 if (computes==ds0) dsparams=&lcand->s;
3580 else
3581 {
3582 dsparams1 params={lcand->s, lastene, lcand->acce, pin0, lastp, lastvfp};
3583 dsparams=&params;
3584 }
3585 cands[ind][0]=new TTempAtom(lcand, pin0, fps[pm], delm, vps[pm], computes(vps[pm], dsparams));
3586 pcs[ind]=1;
3587 }
3588 }
3589 else if (!inputfound)
3590 {
3591 double epf0=settings->epf0;
3592 //lpcs records number of candidates found for partial pin0
3593 int lpcs=0;
3594 //lcoate peaks with the feasible frequency range, specified by f0 and epf0
3595 while (pm>0 && fps[pm]>=f0-epf0) pm--; while (pm<pc-1 && fps[pm]<f0-epf0) pm++;
3596 while (pM<pc-1 && fps[pM]<=f0+epf0) pM++; while (pM>0 && fps[pM]>f0+epf0) pM--;
3597 //add peaks as candidates
3598 for (int j=pm; j<=pM; j++) if (fps[j]>f1 && fps[j]<f2)
3599 {
3600 void* dsparams;
3601 if (computes==ds0) dsparams=&lcand->s;
3602 else
3603 {
3604 dsparams1 params={lcand->s, lastene, lcand->acce, pin0, lastp, lastvfp};
3605 dsparams=&params;
3606 }
3607 cands[ind][lpcs]=new TTempAtom(lcand, pin0, fps[j], delm, vps[j], computes(vps[j], dsparams));
3608 lpcs++;
3609 }
3610 pcs[ind]=lpcs;
3611 }
3612 //harmonic atom grouping fails if user partial is not found
3613 if (pcs[ind]==0){f0=0; goto cleanup;}
3614 else {delete lcand->R; lcand->R=0;}
3615 psb[pin0]=2; ind++;
3616 }
3617 //normal partials (non-anchor, non-user)
3618 //p: partial index
3619 {
3620 int p=0;
3621 while (ind<=numsam)
3622 {
3623 p++;
3624 //skip anchor or user input partials
3625 if (psb[p]) continue;
3626 //loop over all candidates
3627 for (int i=0; i<pcs[ind-1]; i++)
3628 {
3629 //the ith candidate of last partial
3630 TTempAtom* lcand=cands[ind-1][i];
3631 //lpcs records number of candidates found for partial p
3632 int lpcs=0;
3633 //the feasible frequency range for partial p
3634 double f1, f2; ExFmStiff(f1, f2, p, lcand->R->N, lcand->R->X, lcand->R->Y); //calculate the frequency range to search for peak
3635 f1-=delm, f2+=delm; //allow a frequency error for the new partial
3636 if (f1<0) f1=0; if (f1>N/2.1) f1=N/2.1; if (f2>N/2.1) f2=N/2.1;
3637 //locate peaks within this range
3638 while (pm>0 && fps[pm]>=f1) pm--; while (pm<pc-1 && fps[pm]<f1) pm++;
3639 while (pM<pc-1 && fps[pM]<=f2) pM++; while (pM>0 && fps[pM]>f2) pM--; //an index range for peaks
3640 //add peaks as candidates
3641 for (int j=pm; j<=pM; j++)
3642 {
3643 if (fps[j]>f1 && fps[j]<f2) //this peak frequency falls in the frequency range
3644 {
3645 void* dsparams;
3646 if (computes==ds0) dsparams=&lcand->s;
3647 else
3648 {
3649 dsparams1 params={lcand->s, lastene, lcand->acce, pin0, lastp, lastvfp};
3650 dsparams=&params;
3651 }
3652 cands[ind][pcs[ind]+lpcs]=new TTempAtom(lcand, p, fps[j], delm, vps[j], computes(vps[j], dsparams));
3653 lpcs++;
3654 if (pcs[ind]+lpcs>=MAX_CAND) {int newpcs=DeleteByHalf(cands[ind], pcs[ind]); memcpy(&cands[ind][newpcs], &cands[ind][pcs[ind]], sizeof(TTempAtom)*lpcs); pcs[ind]=newpcs;}
3655 }
3656 }
3657 //if there is no peak found for partial p
3658 if (lpcs==0)
3659 {
3660 cands[ind][pcs[ind]+lpcs]=lcand; lpcs++;
3661 cands[ind-1][i]=0;
3662 if (pcs[ind]+lpcs>=MAX_CAND){int newpcs=DeleteByHalf(cands[ind], pcs[ind]); memcpy(&cands[ind][newpcs], &cands[ind][pcs[ind]], sizeof(TTempAtom*)*lpcs); pcs[ind]=newpcs;}
3663 }
3664 else{delete lcand->R; lcand->R=0;}
3665
3666 pcs[ind]+=lpcs;
3667 }
3668 ind++;
3669 }
3670
3671 //select the candidate with maximal s
3672 int maxs=0; double vmax=cands[ind-1][0]->s;
3673 for (int i=1; i<pcs[ind-1]; i++) if (vmax<cands[ind-1][i]->s) maxs=i, vmax=cands[ind-1][i]->s;
3674 TTempAtom* lcand=cands[ind-1][maxs];
3675
3676 //get results
3677 //read frequencies of found partials
3678 int P=0; int* ps=new int[maxp]; double* fs=new double[maxp]; //these are used for evaluating f0 and B
3679 while (lcand){if (lcand->pind) {ps[P]=lcand->pind; fs[P]=lcand->f; P++;} lcand=lcand->Prev;}
3680 //read R of candidate with maximal s as R0
3681 TPolygon* R0=cands[ind-1][maxs]->R;
3682
3683 if (Fr==1) result=GetResultsSingleFr(f0, B, R, R0, P, ps, fs, 0, x[0], N, psb, numsam, settings, results);
3684 else result=GetResultsMultiFr(f0, B, R, R0, P, ps, fs, 0, Fr, x, N, offst, psb, numsam, settings, results, forceinputlocalfr);
3685
3686 delete[] ps; delete[] fs;
3687 }
3688 cleanup:
3689 for (int i=0; i<=numsam; i++)
3690 for (int j=0; j<pcs[i]; j++)
3691 delete cands[i][j];
3692 DeAlloc2(cands);
3693 delete[] psb;
3694
3695 return result;
3696 }//NoteMatchStiff3
3697
3698 /*
3699 function NoteMatchStiff3: finds harmonic atom from spectrum if Fr=1, or constant-pitch harmonic
3700 sinusoid from spectrogram if Fr>1. This version uses residue-sinusoid ratio ("rsr") that measures how
3701 "good" a spectral peak is in terms of its shape. peaks with large rsr[] is likely to be contaminated
3702 and their frequency estimates are regarded unreliable.
3703
3704 In: x[Fr][N/2+1]: spectrogram
3705 N, offst: atom scale and hop size
3706 fps[pc], vps[pc], rsr[pc]: primitive (rough) peak frequencies, amplitudes and shape factor
3707 R: initial F-G polygon constraint, optional
3708 settings: note match settings
3709 computes: pointer to a function that computes HA score, must be ds0 or ds1
3710 lastvfp[lastp]: amplitude of the previous harmonic atom
3711 forceinputlocalfr: specifies if partial settings->pin0 is taken for granted ("pinned")
3712 Out: results: note match results
3713 f0, B: fundamental frequency and stiffness coefficient
3714 R: F-G polygon of harmonic atom
3715
3716 Returns the total energy of HA or constant-pitch HS.
3717 */
3718 double NoteMatchStiff3(TPolygon* R, double &f0, double& B, int pc, double* fps, double* vps,
3719 double* rsr, int Fr, cdouble** x, int N, int offst, NMSettings* settings, NMResults* results,
3720 int lastp, double* lastvfp, double (*computes)(double a, void* params), int forceinputlocalfr)
3721 {
3722 double result=0;
3723
3724 double maxB=settings->maxB, minf0=settings->minf0, maxf0=settings->maxf0,
3725 delm=settings->delm, delp=settings->delp, *c=settings->c, iH2=settings->iH2,
3726 *pinf=settings->pinf, hB=settings->hB, epf0=settings->epf0;// epf=settings->epf,
3727 int maxp=settings->maxp, *pin=settings->pin, pin0=settings->pin0, *pinfr=settings->pinfr,
3728 pcount=settings->pcount, M=settings->M;
3729
3730 // int *ptype=results->ptype;
3731
3732 //calculate the energy of the last frame (hp)
3733 double lastene=0; for (int i=0; i<lastp; i++) lastene+=lastvfp[i]*lastvfp[i];
3734 //fill frequency and amplitude buffers with 0
3735 //now determine numsam (the maximal number of partials)
3736 bool inputpartial=(f0>0 && settings->pin0>0);
3737 if (!inputpartial)
3738 {
3739 if (pcount>0) f0=pinf[0], pin0=pin[0];
3740 else f0=sqrt(R->X[0]), pin0=1;
3741 }
3742 int numsam=N*pin0/2.1/f0;
3743 if (numsam>maxp) numsam=maxp;
3744 //allocate buffer for candidate atoms
3745 TTempAtom*** Allocate2(TTempAtom*, numsam+1, MAX_CAND, cands);
3746 //pcs[p] records the number of candidates for partial index p
3747 //psb[p] = 1: anchor, 2: user input, 0: normal
3748 int *psb=new int[(numsam+1)*2], *pcs=&psb[numsam+1];
3749 memset(psb, 0, sizeof(int)*(numsam+1)*2);
3750 //if R is not specified, initialize a trivial candidate with maxB
3751 if (R->N==0) cands[0][0]=new TTempAtom(1, (minf0+maxf0)*0.5, (maxf0-minf0)*0.5, maxB);
3752 //if R is, initialize a trivial candidate by extending R by delp
3753 else cands[0][0]=new TTempAtom(R, delp, delp, minf0);
3754
3755 int pm=0, pM=0;
3756 int ind=1; pcs[0]=1;
3757 //anchor partials: highest priority atoms, must have spectral peaks
3758 for (int i=0; i<pcount; i++)
3759 {
3760 //skip out-of-scope atoms
3761 if (pin[i]>=numsam) continue;
3762 //skip user input atom
3763 if (inputpartial && pin[i]==pin0) continue;
3764
3765 void* dsparams;
3766 if (computes==ds0) dsparams=&cands[ind-1][0]->s;
3767 else
3768 {
3769 dsparams1 params={cands[ind-1][0]->s, lastene, cands[ind-1][0]->acce, pin[i], lastp, lastvfp};
3770 dsparams=&params;
3771 }
3772
3773 if (pinfr[i]>0) //local anchor
3774 {
3775 double ls=0, lrsr=PeakShapeC(pinf[i], 1, N, &x[pinfr[i]], hB*2, M, c, iH2);
3776 for (int fr=0; fr<Fr; fr++)
3777 {
3778 cdouble r=IPWindowC(pinf[i], x[fr], N, M, c, iH2, pinf[i]-hB, pinf[i]+hB);
3779 ls+=~r;
3780 }
3781 ls=sqrt(ls/Fr);
3782 cands[ind][0]=new TTempAtom(cands[ind-1][0], pin[i], pinf[i], delm+lrsr, ls, computes(ls, dsparams)); cands[ind][0]->rsr=lrsr;
3783 }
3784 else
3785 {
3786 //find the nearest peak to pinf[i]
3787 while (pm<pc-1 && fps[pm]<pinf[i]) pm++; while (pm>0 && fps[pm]>=pinf[i]) pm--;
3788 if (pm<pc-1 && pinf[i]-fps[pm]>fps[pm+1]-pinf[i]) pm++;
3789 cands[ind][0]=new TTempAtom(cands[ind-1][0], pin[i], fps[pm], delm+rsr[pm], vps[pm], computes(vps[pm], dsparams)); cands[ind][0]->rsr=rsr[pm];
3790 }
3791 delete cands[ind-1][0]->R; cands[ind-1][0]->R=0; psb[pin[i]]=1; pcs[ind]=1; ind++;
3792 }
3793 //harmonic atom grouping fails if anchors conflict
3794 if (cands[ind-1][0]->R->N==0) {f0=0; goto cleanup;}
3795
3796 //user input partial: must have spectral peak
3797 if (inputpartial)
3798 {
3799 //harmonic atom grouping fails if user partial is out of scope
3800 if (pin0>numsam){f0=0; goto cleanup;}
3801
3802 TTempAtom* lcand=cands[ind-1][0];
3803 //feasible frequency range for partial pin0
3804 double f1, f2; ExFmStiff(f1, f2, pin0, lcand->R->N, lcand->R->X, lcand->R->Y);
3805 f1-=delm, f2+=delm;
3806 if (f1<0) f1=0; if (f1>N/2.1) f1=N/2.1; if (f2>N/2.1) f2=N/2.1;
3807 bool inputfound=false;
3808 if (forceinputlocalfr>=0)
3809 {
3810 int k1=floor(f0-epf0-1), k2=ceil(f0+epf0+1), k12=k2-k1+1, pk=-1;
3811
3812 cdouble *lx=x[forceinputlocalfr], *lr=new cdouble[k12];
3813 for (int k=0; k<k12; k++) lr[k]=IPWindowC(k1+k, lx, N, M, c, iH2, k1+k-hB, k1+k+hB);
3814 for (int k=1; k<k12-1; k++)
3815 if (~lr[k]>~lr[k-1] && ~lr[k]>~lr[k+1])
3816 {
3817 if (pk==-1 || fabs(k1+pk-f0)>fabs(k1+k-f0)) pk=k;
3818 }
3819 if (pk>0)
3820 {
3821 double lf=k1+pk, la, lph, oldlf=lf;
3822 LSESinusoid(lf, lf-1, lf+1, lx, N, hB, M, c, iH2, la, lph, settings->epf);
3823 if (fabs(lf-oldlf)>0.6) //by which we judge that the LS result is biased by interference
3824 {
3825 lf=oldlf;
3826 cdouble r=IPWindowC(lf, lx, N, M, c, iH2, lf-hB, lf+hB);
3827 la=abs(r); lph=arg(r);
3828 }
3829 if (lf>f1 && lf<f2)
3830 {
3831 void* dsparams;
3832 if (computes==ds0) dsparams=&lcand->s;
3833 else
3834 {
3835 dsparams1 params={lcand->s, lastene, lcand->acce, pin0, lastp, lastvfp};
3836 dsparams=&params;
3837 }
3838 double lrsr=PeakShapeC(lf, 1, N, &x[forceinputlocalfr], hB*2, M, c, iH2);
3839 cands[ind][0]=new TTempAtom(lcand, pin0, lf, delm+lrsr, la, computes(la, dsparams)); cands[ind][0]->rsr=lrsr;
3840 pcs[ind]=1;
3841 inputfound=true;
3842 }
3843 }
3844 }
3845 //forcepin0 being true means this partial have to be the peak nearest to f0
3846 if (!inputfound && settings->forcepin0)
3847 { //find the nearest peak to f0
3848 while (pm<pc-1 && fps[pm]<f0) pm++; while (pm>0 && fps[pm]>=f0) pm--;
3849 if (pm<pc-1 && f0-fps[pm]>fps[pm+1]-f0) pm++;
3850 if (fps[pm]>f1 && fps[pm]<f2)
3851 {
3852 void* dsparams;
3853 if (computes==ds0) dsparams=&lcand->s;
3854 else
3855 {
3856 dsparams1 params={lcand->s, lastene, lcand->acce, pin0, lastp, lastvfp};
3857 dsparams=&params;
3858 }
3859 cands[ind][0]=new TTempAtom(lcand, pin0, fps[pm], delm+rsr[pm], vps[pm], computes(vps[pm], dsparams)); cands[ind][0]->rsr=rsr[pm];
3860 pcs[ind]=1;
3861 }
3862 }
3863 else if (!inputfound)
3864 {
3865 double epf0=settings->epf0;
3866 //lpcs records number of candidates found for partial pin0
3867 int lpcs=0;
3868 //lcoate peaks with the feasible frequency range, specified by f0 and epf0
3869 while (pm>0 && fps[pm]>=f0-epf0) pm--; while (pm<pc-1 && fps[pm]<f0-epf0) pm++;
3870 while (pM<pc-1 && fps[pM]<=f0+epf0) pM++; while (pM>0 && fps[pM]>f0+epf0) pM--;
3871 //add peaks as candidates
3872 for (int j=pm; j<=pM; j++) if (fps[j]>f1 && fps[j]<f2)
3873 {
3874 void* dsparams;
3875 if (computes==ds0) dsparams=&lcand->s;
3876 else
3877 {
3878 dsparams1 params={lcand->s, lastene, lcand->acce, pin0, lastp, lastvfp};
3879 dsparams=&params;
3880 }
3881 cands[ind][lpcs]=new TTempAtom(lcand, pin0, fps[j], delm+rsr[j], vps[j], computes(vps[j], dsparams)); cands[ind][lpcs]->rsr=rsr[j];
3882 lpcs++;
3883 }
3884 pcs[ind]=lpcs;
3885 }
3886 //harmonic atom grouping fails if user partial is not found
3887 if (pcs[ind]==0){f0=0; goto cleanup;}
3888 else {delete lcand->R; lcand->R=0;}
3889 psb[pin0]=2; ind++;
3890 }
3891 //normal partials (non-anchor, non-user)
3892 //p: partial index
3893 {
3894 int p=0;
3895 while (ind<=numsam)
3896 {
3897 p++;
3898 //skip anchor or user input partials
3899 if (psb[p]) continue;
3900 //loop over all candidates
3901 for (int i=0; i<pcs[ind-1]; i++)
3902 {
3903 int tightpks=0;
3904 //the ith candidate of last partial
3905 TTempAtom* lcand=cands[ind-1][i];
3906 //lpcs records number of candidates found for partial p
3907 int lpcs=0;
3908 //the feasible frequency range for partial p
3909 double tightf1, tightf2; ExFmStiff(tightf1, tightf2, p, lcand->R->N, lcand->R->X, lcand->R->Y); //calculate the frequency range to search for peak
3910 double tightf=(tightf1+tightf2)*0.5, f1=tightf1-delm, f2=tightf2+delm; //allow a frequency error for the new partial
3911 if (f1<0) f1=0; if (f1>N/2.1) f1=N/2.1; if (f2>N/2.1) f2=N/2.1;
3912 //locate peaks within this range
3913 while (pm>0 && fps[pm]>=f1) pm--; while (pm<pc-1 && fps[pm]<f1) pm++;
3914 while (pM<pc-1 && fps[pM]<=f2) pM++; while (pM>0 && fps[pM]>f2) pM--; //an index range for peaks
3915 //add peaks as candidates
3916 for (int j=pm; j<=pM; j++)
3917 {
3918 if (fps[j]>f1 && fps[j]<f2) //this peak frequency falls in the frequency range
3919 {
3920 if ((fps[j]>tightf1 && fps[j]<tightf2) || fabs(fps[j]-tightf)<0.5) tightpks++; //indicates there are "tight" peaks found for this partial
3921 void* dsparams;
3922 if (computes==ds0) dsparams=&lcand->s;
3923 else
3924 {
3925 dsparams1 params={lcand->s, lastene, lcand->acce, pin0, lastp, lastvfp};
3926 dsparams=&params;
3927 }
3928 cands[ind][pcs[ind]+lpcs]=new TTempAtom(lcand, p, fps[j], delm+rsr[j], vps[j], computes(vps[j], dsparams)); cands[ind][pcs[ind]+lpcs]->rsr=rsr[j];
3929 lpcs++;
3930 if (pcs[ind]+lpcs>=MAX_CAND) {int newpcs=DeleteByHalf(cands[ind], pcs[ind]); memcpy(&cands[ind][newpcs], &cands[ind][pcs[ind]], sizeof(TTempAtom*)*lpcs); pcs[ind]=newpcs;}
3931 }
3932 }
3933 //if there is no peak found for partial p.
3934 if (!lpcs)
3935 {
3936 //use the lcand directly
3937 //BUT UPDATE accumulative energy (acce) and s using induced partial
3938 if (tightf<N/2-hB)
3939 {
3940 double la=0; for (int fr=0; fr<Fr; fr++) la+=~IPWindowC(tightf, x[fr], N, M, c, iH2, tightf-hB, tightf+hB); la=sqrt(la/Fr);
3941 void* dsparams;
3942 if (computes==ds0) dsparams=&lcand->s;
3943 else
3944 {
3945 dsparams1 params={lcand->s, lastene, lcand->acce, pin0, lastp, lastvfp};
3946 dsparams=&params;
3947 }
3948 double ls=computes(la, dsparams);
3949
3950 lcand->acce+=la*la;
3951 lcand->s=ls;
3952 }
3953 cands[ind][pcs[ind]+lpcs]=lcand;
3954 lpcs++;
3955 cands[ind-1][i]=0;
3956 if (pcs[ind]+lpcs>=MAX_CAND){int newpcs=DeleteByHalf(cands[ind], pcs[ind]); memcpy(&cands[ind][newpcs], &cands[ind][pcs[ind]], sizeof(TTempAtom*)*lpcs); pcs[ind]=newpcs;}
3957 }
3958 //if there are peaks but no tight peak found for partial p
3959 else if (!tightpks)
3960 {
3961 if (tightf<N/2-hB)
3962 {
3963 double la=0; for (int fr=0; fr<Fr; fr++) la+=~IPWindowC(tightf, x[fr], N, M, c, iH2, tightf-hB, tightf+hB); la=sqrt(la/Fr);
3964 void* dsparams;
3965 if (computes==ds0) dsparams=&lcand->s;
3966 else
3967 {
3968 dsparams1 params={lcand->s, lastene, lcand->acce, pin0, lastp, lastvfp};
3969 dsparams=&params;
3970 }
3971 double ls=computes(la, dsparams);
3972 TTempAtom* lnewcand=new TTempAtom(lcand, true);
3973 lnewcand->f=0; lnewcand->acce+=la*la; lnewcand->s=ls;
3974 cands[ind][pcs[ind]+lpcs]=lnewcand;
3975 lpcs++;
3976 if (pcs[ind]+lpcs>=MAX_CAND)
3977 {int newpcs=DeleteByHalf(cands[ind], pcs[ind]); memcpy(&cands[ind][newpcs], &cands[ind][pcs[ind]], sizeof(TTempAtom*)*lpcs); pcs[ind]=newpcs;}
3978 }
3979 delete lcand->R; lcand->R=0;
3980 }
3981 else{delete lcand->R; lcand->R=0;}
3982
3983 pcs[ind]+=lpcs;
3984 }
3985 ind++;
3986 }
3987
3988 //select the candidate with maximal s
3989 int maxs=0; double vmax=cands[ind-1][0]->s;
3990 for (int i=1; i<pcs[ind-1]; i++) if (vmax<cands[ind-1][i]->s) maxs=i, vmax=cands[ind-1][i]->s;
3991 TTempAtom* lcand=cands[ind-1][maxs];
3992
3993 //get results
3994 //read frequencies of found partials
3995 int P=0; int* ps=new int[maxp]; double* fs=new double[maxp*2], *rsrs=&fs[maxp]; //these are used for evaluating f0 and B
3996 while (lcand){if (lcand->pind && lcand->f) {ps[P]=lcand->pind; fs[P]=lcand->f; rsrs[P]=lcand->rsr; P++;} lcand=lcand->Prev;}
3997 //read R of candidate with maximal s as R0
3998 TPolygon* R0=cands[ind-1][maxs]->R;
3999
4000 if (Fr==1) result=GetResultsSingleFr(f0, B, R, R0, P, ps, fs, rsrs, x[0], N, psb, numsam, settings, results);
4001 else result=GetResultsMultiFr(f0, B, R, R0, P, ps, fs, rsrs, Fr, x, N, offst, psb, numsam, settings, results, forceinputlocalfr);
4002
4003 delete[] ps; delete[] fs;
4004 }
4005 cleanup:
4006 for (int i=0; i<=numsam; i++)
4007 for (int j=0; j<pcs[i]; j++)
4008 delete cands[i][j];
4009 DeAlloc2(cands);
4010 delete[] psb;
4011
4012 return result;
4013 }//NoteMatchStiff3
4014
4015 /*
4016 function NoteMatchStiff3: wrapper function of the above that does peak picking and HA grouping (or
4017 constant-pitch HS finding) in a single call.
4018
4019 In: x[Fr][N/2+1]: spectrogram
4020 N, offst: atom scale and hop size
4021 R: initial F-G polygon constraint, optional
4022 startfr, validfrrange: centre and half width, in frames, of the interval used for peak picking if Fr>1
4023 settings: note match settings
4024 computes: pointer to a function that computes HA score, must be ds0 or ds1
4025 lastvfp[lastp]: amplitude of the previous harmonic atom
4026 forceinputlocalfr: specifies if partial settings->pin0 is taken for granted ("pinned")
4027 Out: results: note match results
4028 f0, B: fundamental frequency and stiffness coefficient
4029 R: F-G polygon of harmonic atom
4030
4031 Returns the total energy of HA or constant-pitch HS.
4032 */
4033 double NoteMatchStiff3(TPolygon* R, double &f0, double& B, int Fr, cdouble** x, int N, int offst, NMSettings* settings,
4034 NMResults* results, int lastp, double* lastvfp, double (*computes)(double a, void* params),
4035 bool forceinputlocalfr, int startfr, int validfrrange)
4036 {
4037 //find spectral peaks
4038 double *fps=(double*)malloc8(sizeof(double)*N*2*Fr), *vps=&fps[N/2*Fr], *rsr=&fps[N*Fr];
4039 int pc;
4040 if (Fr>1) pc=QuickPeaks(fps, vps, Fr, N, x, startfr, validfrrange, settings->M, settings->c, settings->iH2, 0.0005, -1, -1, settings->hB*2, rsr);
4041 else pc=QuickPeaks(fps, vps, N, x[0], settings->M, settings->c, settings->iH2, 0.0005, -1, -1, settings->hB*2, rsr);
4042 //if no peaks are present, return 0, indicating empty hp
4043 if (pc==0){free8(fps); return 0;}
4044
4045 double result=NoteMatchStiff3(R, f0, B, pc, fps, vps, rsr, Fr, x, N, offst, settings, results, lastp, lastvfp, computes, forceinputlocalfr?startfr:-1);
4046 free8(fps);
4047 return result;
4048 }//NoteMatchStiff3
4049
4050 /*
4051 function NoteMatchStiff3: finds harmonic atom from spectrum given the pitch as a (partial index,
4052 frequency) pair. This is used internally by NoteMatch3FB().
4053
4054 In: x[N/2+1]: spectrum
4055 fps[pc], vps[pc]: primitive (rough) peak frequencies and amplitudes
4056 R: initial F-G polygon constraint, optional
4057 pin0: partial index in the pitch specifier pair
4058 peak0: index to frequency in fps[] in the pitch specifier pair
4059 settings: note match settings
4060 Out: vfp[]: partial amplitudes
4061 R: F-G polygon of harmonic atom
4062
4063 Returns the total energy of HA.
4064 */
4065 double NoteMatchStiff3(TPolygon* R, int peak0, int pin0, cdouble* x, int pc, double* fps, double* vps, int N,
4066 NMSettings* settings, double* vfp, int** pitchind, int newpc)
4067 {
4068 double maxB=settings->maxB, minf0=settings->minf0, maxf0=settings->maxf0,
4069 delm=settings->delm, delp=settings->delp, *c=settings->c, iH2=settings->iH2,
4070 hB=settings->hB;
4071 int maxp=settings->maxp, M=settings->M; // pcount=settings->pcount;
4072
4073 double f0=fps[peak0];
4074
4075 //determine numsam
4076 int numsam=N*pin0/2.1/f0; if (numsam>maxp) numsam=maxp;
4077 if (pin0>=numsam) return 0;
4078 //pcs[p] contains the number of candidates for partial index p
4079 TTempAtom*** Allocate2(TTempAtom*, numsam+1, MAX_CAND, cands);
4080 int *pcs=new int[numsam+1]; memset(pcs, 0, sizeof(int)*(numsam+1));
4081 if (R->N==0) cands[0][0]=new TTempAtom(1, (minf0+maxf0)*0.5, (maxf0-minf0)*0.5, maxB); //initialization: trivial candidate with maxB
4082 else cands[0][0]=new TTempAtom(R, delp, delp, minf0); //initialization: trivial candidate by expanding R
4083
4084 cands[1][0]=new TTempAtom(cands[0][0], pin0, fps[peak0], delm, vps[peak0], vps[peak0]); cands[1][0]->tag[0]=peak0;
4085 delete cands[0][0]->R; cands[0][0]->R=0;
4086
4087 int pm=0, pM=0, ind=2, p=0; pcs[0]=pcs[1]=1;
4088
4089 while (ind<=numsam)
4090 {
4091 p++;
4092 if (p==pin0) continue;
4093 for (int i=0; i<pcs[ind-1]; i++)
4094 {
4095 TTempAtom* lcand=cands[ind-1][i]; //the the ith candidate of last partial
4096 int lpcs=0;
4097 {
4098 double f1, f2; ExFmStiff(f1, f2, p, lcand->R->N, lcand->R->X, lcand->R->Y); //calculate the frequency range to search for peak
4099 f1-=delm, f2+=delm; //allow a frequency error for the new partial
4100 if (f1<0) f1=0; if (f1>N/2.1) f1=N/2.1; if (f2>N/2.1) f2=N/2.1;
4101 while (pm>0 && fps[pm]>=f1) pm--; while (pm<pc-1 && fps[pm]<f1) pm++;
4102 while (pM<pc-1 && fps[pM]<=f2) pM++; while (pM>0 && fps[pM]>f2) pM--; //an index range for peaks
4103 for (int j=pm; j<=pM; j++)
4104 {
4105 if (fps[j]>f1 && fps[j]<f2 && (p>pin0 || pitchind[j][p]<0))
4106 {
4107 //this peak frequency falls in the frequency range
4108 cands[ind][pcs[ind]+lpcs]=new TTempAtom(lcand, p, fps[j], delm, vps[j], lcand->s+vps[j]); cands[ind][pcs[ind]+lpcs]->tag[0]=j; lpcs++;
4109 if (pcs[ind]+lpcs>=MAX_CAND){int newpcs=DeleteByHalf(cands[ind], pcs[ind]); memcpy(&cands[ind][newpcs], &cands[ind][pcs[ind]], sizeof(TTempAtom)*lpcs); pcs[ind]=newpcs;}
4110 }
4111 }
4112 if (lpcs==0) //there is no peak found for the partial i, the last level
4113 {
4114 cands[ind][pcs[ind]+lpcs]=lcand; lpcs++; //new TTempAtom(lcand, false);
4115 cands[ind-1][i]=0;
4116 if (pcs[ind]+lpcs>=MAX_CAND){int newpcs=DeleteByHalf(cands[ind], pcs[ind]); memcpy(&cands[ind][newpcs], &cands[ind][pcs[ind]], sizeof(TTempAtom*)*lpcs); pcs[ind]=newpcs;}
4117 }
4118 else
4119 {
4120 delete lcand->R; lcand->R=0;
4121 }
4122 }
4123 pcs[ind]+=lpcs;
4124 }
4125 ind++;
4126 }
4127
4128 //select the candidate with maximal s
4129 int maxs=0; double vmax=cands[ind-1][0]->s;
4130 for (int i=1; i<pcs[ind-1]; i++) if (vmax<cands[ind-1][i]->s) maxs=i, vmax=cands[ind-1][i]->s;
4131 TTempAtom* lcand=cands[ind-1][maxs];
4132
4133 //get fp, vfp, pfp, ptype
4134 memset(vfp, 0, sizeof(double)*maxp);
4135
4136 int P=0; int *ps=new int[maxp], *peaks=new int[maxp]; double* fs=new double[maxp]; //these are used for evaluating f0 and B
4137
4138 while (lcand){if (lcand->pind) {ps[P]=lcand->pind; fs[P]=lcand->f; peaks[P]=lcand->tag[0]; P++;} lcand=lcand->Prev;}
4139 R->N=0;
4140
4141 if (P>0)
4142 {
4143 for (int i=P-1; i>=0; i--)
4144 {
4145 int lp=ps[i];
4146 double lf=fs[i];
4147 vfp[lp-1]=vps[peaks[i]];
4148 if (R->N==0) InitializeR(R, lp, lf, delm, maxB);
4149 else if (vfp[lp-1]>1) CutR(R, lp, lf, delm, true);
4150 if (pitchind[peaks[i]][lp]>=0) {R->N=0; goto cleanup;}
4151 else pitchind[peaks[i]][lp]=newpc;
4152 }
4153
4154 //estimate f0 and B
4155 double tmpa, cF, cG;
4156 // double norm[1024]; for (int i=0; i<1024; i++) norm[i]=1;
4157 areaandcentroid(tmpa, cF, cG, R->N, R->X, R->Y);
4158 testnn(cF); f0=sqrt(cF); //B=cG/cF;
4159
4160
4161 for (int i=0; i<numsam; i++)
4162 {
4163 if (vfp[i]==0) //no peak is found for this partial in lcand
4164 {
4165 int m=i+1;
4166 double tmp=cF+(m*m-1)*cG; testnn(tmp);
4167 double lf=m*sqrt(tmp);
4168 if (lf<N/2.1)
4169 {
4170 cdouble r=IPWindowC(lf, x, N, M, c, iH2, lf-hB, lf+hB);
4171 vfp[i]=-sqrt(r.x*r.x+r.y*r.y);
4172 // vfp[i]=-IPWindowC(lf, xr, xi, N, M, c, iH2, lf-hB, lf+hB, 0);
4173 }
4174 }
4175 }
4176 }
4177
4178 cleanup:
4179 delete[] ps; delete[] fs; delete[] peaks;
4180 for (int i=0; i<=numsam; i++) for (int j=0; j<pcs[i]; j++) delete cands[i][j];
4181 DeAlloc2(cands); delete[] pcs;
4182
4183 double result=0;
4184 if (f0>0) for (int i=0; i<numsam; i++) result+=vfp[i]*vfp[i];
4185 return result;
4186 }//NoteMatchStiff3*/
4187
4188 /*
4189 function NoteMatchStiff3FB: does one dynamic programming step in forward-background pitch tracking of
4190 harmonic sinusoids. This is used internally by FindNoteFB().
4191
4192 In: R[pitchcount]: initial F-G polygons associated with pitch candidates at last frame
4193 vfp[pitchcount][maxp]: amplitude vectors associated with pitch candidates at last frame
4194 sc[pitchcount]: accumulated scores associated with pitch candidates at last frame
4195 fps[pc], vps[pc]: primitive (rough) peak frequencies and amplitudes at this frame
4196 maxpitch: maximal number of pitch candidates
4197 x[wid/2+1]: spectrum
4198 settings: note match settings
4199 Out: pitchcount: number of pitch candidiate at this frame
4200 newpitches[pitchcount]: pitch candidates at this frame, whose lower word is peak index, higher word is partial index
4201 prev[pitchcount]: indices to predecessors of the pitch candidates at this frame
4202 R[pitchcount]: F-G polygons associated with pitch candidates at this frame
4203 vfp[pitchcount][maxp]: amplitude vectors associated with pitch candidates at this frame
4204 sc[pitchcount]: accumulated scores associated with pitch candidates at this frame
4205
4206 No return value.
4207 */
4208 void NoteMatchStiff3FB(int& pitchcount, TPolygon**& R, double**& vfp, double*& sc, int* newpitches, int* prev, int pc, double* fps, double* vps, cdouble* x, int wid, int maxpitch, NMSettings* settings)
4209 {
4210 int maxpin0=6; //this specifies the maximal value that a pitch candidate may have as partial index
4211 double minf0=settings->minf0, delm=settings->delm, delp=settings->delp;
4212 int maxp=settings->maxp;
4213
4214 //pc is the number of peaks at this frame, maxp is the maximal number of partials in the model
4215 //pitchind is initialized to a sequence of -1
4216 int** Allocate2(int, pc, maxp+1, pitchind); memset(pitchind[0], 0xff, sizeof(int)*pc*(maxp+1));
4217
4218 //extend F-G polygons to allow pitch variation
4219 for (int i=0; i<pitchcount; i++) ExtendR(R[i], delp, delp, minf0);
4220
4221 double *f1=new double[pitchcount], *f2=new double[pitchcount];
4222 int pm=0, newpc=0;
4223 TPolygon** newR=new TPolygon*[maxpitch]; memset(newR, 0, sizeof(TPolygon*)*maxpitch);
4224 double* newsc=new double[maxpitch]; memset(newsc, 0, sizeof(double)*maxpitch);
4225 double** Allocate2(double, maxpitch, maxp, newvfp);
4226
4227 for (int pin0=1; pin0<=maxpin0; pin0++) //pin0: the partial index of a candidate pitch
4228 {
4229 //find out the range [pm, pM) of indices into fps[], so that for a * in this range (pin0, fps[*]) may fall
4230 //in the feasible pitch range succeeding one of the candidate pitches of the previous frame
4231 for (int i=0; i<pitchcount; i++) {ExFmStiff(f1[i], f2[i], pin0, R[i]->N, R[i]->X, R[i]->Y); f1[i]-=delm, f2[i]+=delm;}
4232 double f1a=f1[0], f2a=f2[0]; for (int i=1; i<pitchcount; i++) {if (f1a>f1[i]) f1a=f1[i]; if (f2a<f2[i]) f2a=f2[i];}
4233 while (pm<pc-1 && fps[pm]<f1a) pm++;
4234 int pM=pm; while (pM<pc-1 && fps[pM]<f2a) pM++;
4235
4236 for (int p=pm; p<pM; p++) //loop through all peaks in this range
4237 {
4238 if (pitchind[p][pin0]>=0) continue; //skip this peak if it is already ...
4239 int max; double maxs; TPolygon* lnewR=0;
4240
4241 for (int i=0; i<pitchcount; i++) //loop through candidate pitches of the previous frame
4242 {
4243 if (fps[p]>f1[i] && fps[p]<f2[i]) //so that this peak is a feasible successor of the i-th candidate pitch of the previous frame
4244 {
4245 if (pitchind[p][pin0]<0) //if this peak has not been registered as a candidate pitch with pin0 being the partial index
4246 {
4247 lnewR=new TPolygon(maxp*2+4, R[i]); //create a F-G polygon for (this peak, pin0) pair
4248 if (newpc==maxpitch) //maximal number of candidate pitches is reached
4249 {
4250 //delete the candidate with the lowest score without checking if this lowest score is below the score
4251 //of the current pitch candidate (which is not yet computed). of course there is a risk that the new
4252 //score is even lower so that $maxpitch buffered scores may not be the highest $maxpitch scores, but
4253 //the (maxpitch-1) highest of the $maxpitch buffered scores should be the highest (maxpitch01) scores.
4254 int minsc=0; for (int j=0; j<maxpitch; j++) if (newsc[j]<newsc[minsc]) minsc=j;
4255 delete newR[minsc];
4256 if (minsc!=newpc-1) {newR[minsc]=newR[newpc-1]; memcpy(newvfp[minsc], newvfp[newpc-1], sizeof(double)*maxp); newsc[minsc]=newsc[newpc-1]; prev[minsc]=prev[newpc-1]; newpitches[minsc]=newpitches[newpc-1];}
4257 }
4258 else newpc++;
4259 //try to find harmonic atom with this candidate pitch
4260 if (NoteMatchStiff3(lnewR, p, pin0, x, pc, fps, vps, wid, settings, newvfp[newpc-1], pitchind, newpc-1)>0)
4261 {//if successful, buffer its F-G polygon and save its score
4262 newR[newpc-1]=lnewR;
4263 max=i; maxs=sc[i]+conta(maxp, newvfp[newpc-1], vfp[i]);
4264 }
4265 else
4266 {//if not, discard this candidate pitch
4267 delete lnewR; lnewR=0; newpc--;
4268 }
4269 }
4270 else //if this peak has already been registered as a candidate pitch with pin0 being the partial index
4271 {
4272 //compute it score as successor to the i-th candidate of the previous frame
4273 double ls=sc[i]+conta(maxp, newvfp[newpc-1], vfp[i]);
4274 //if the score is higher than mark the i-th candidate of previous frame as its predecessor
4275 if (ls>maxs) maxs=ls, max=i;
4276 }
4277 }
4278 }
4279 if (lnewR) //i.e. a HA is found for this pitch candidate
4280 {
4281 ((__int16*)&newpitches[newpc-1])[0]=p; ((__int16*)&newpitches[newpc-1])[1]=pin0;
4282 newsc[newpc-1]=maxs; //take note of its score
4283 prev[newpc-1]=max; //take note of its predecessor
4284 }
4285 }
4286 }
4287 DeAlloc2(pitchind);
4288 delete[] f1; delete[] f2;
4289
4290 for (int i=0; i<pitchcount; i++) delete R[i]; delete[] R; R=newR;
4291 for (int i=0; i<newpc; i++) for (int j=0; j<maxp; j++) if (newvfp[i][j]<0) newvfp[i][j]=-newvfp[i][j];
4292 DeAlloc2(vfp); vfp=newvfp;
4293 delete[] sc; sc=newsc;
4294 pitchcount=newpc;
4295 }//NoteMatchStiff3FB
4296
4297 /*
4298 function PeakShapeC: residue-sinusoid-ratio for a given (hypothesis) sinusoid frequency
4299
4300 In: x[Fr][N/2+1]: spectrogram
4301 M, c[], iH2: cosine-family window specifiers
4302 f: reference frequency, in bins
4303 B: spectral truncation width
4304
4305 Returns the residue-sinusoid-ratio.
4306 */
4307 double PeakShapeC(double f, int Fr, int N, cdouble** x, int B, int M, double* c, double iH2)
4308 {
4309 cdouble* w=new cdouble[B];
4310 int fst=floor(f+1-B/2.0);
4311 if (fst<0) fst=0;
4312 if (fst+B>N/2) fst=N/2-B;
4313 Window(w, f, N, M, c, fst, fst+B-1);
4314 cdouble xx=0, ww=Inner(B, w, w);
4315 double xw2=0;
4316 for (int fr=0; fr<Fr; fr++)
4317 xw2+=~Inner(B, &x[fr][fst], w), xx+=Inner(B, &x[fr][fst], &x[fr][fst]);
4318 delete[] w;
4319 if (xx.x==0) return 1;
4320 return 1-xw2/(xx.x*ww.x);
4321 }//PeakShapeC
4322 //version using cmplx<float> as spectrogram data type
4323 double PeakShapeC(double f, int Fr, int N, cfloat** x, int B, int M, double* c, double iH2)
4324 {
4325 cdouble* w=new cdouble[B];
4326 int fst=floor(f+1-B/2.0);
4327 if (fst<0) fst=0;
4328 if (fst+B>N/2) fst=N/2-B;
4329 Window(w, f, N, M, c, fst, fst+B-1);
4330 cdouble xx=0, ww=Inner(B, w, w);
4331 double xw2=0;
4332 for (int fr=0; fr<Fr; fr++)
4333 xw2+=~Inner(B, &x[fr][fst], w), xx+=Inner(B, &x[fr][fst], &x[fr][fst]);
4334 delete[] w;
4335 if (xx.x==0) return 1;
4336 return 1-xw2/(xx.x*ww.x);
4337 }//PeakShapeC
4338
4339 /*
4340 function QuickPeaks: finds rough peaks in the spectrum (peak picking)
4341
4342 In: x[N/2+1]: spectrum
4343 M, c[], iH2: cosine-family window function specifiers
4344 mina: minimal amplitude to spot a spectral peak
4345 [binst, binen): frequency range, in bins, to look for peaks
4346 B: spectral truncation width
4347 Out; f[return value], a[return value]: frequencies (in bins) and amplitudes of found peaks
4348 rsr[return value]: residue-sinusoid-ratio of found peaks, optional
4349
4350 Returns the number of peaks found. f[] and a[] must be allocated enough space before calling.
4351 */
4352 int QuickPeaks(double* f, double* a, int N, cdouble* x, int M, double* c, double iH2, double mina, int binst, int binen, int B, double* rsr)
4353 {
4354 double hB=B*0.5;
4355 if (binst<2) binst=2;
4356 if (binst<hB) binst=hB;
4357 if (binen<0) binen=N/2-1;
4358 if (binen<0 || binen+hB>N/2-1) binen=N/2-1-hB;
4359 double a0=~x[binst-1], a1=~x[binst], a2;
4360 double minA=mina*mina*2/iH2;
4361 int p=0, n=binst;
4362 cdouble* w=new cdouble[B];
4363 while (n<binen)
4364 {
4365 a2=~x[n+1];
4366 if (a1>0 && a1>=a0 && a1>=a2)
4367 {
4368 if (a1>minA)
4369 {
4370 double A0=sqrt(a0), A1=sqrt(a1), A2=sqrt(a2);
4371 f[p]=n+(A0-A2)/2/(A0+A2-2*A1);
4372 int fst=floor(f[p]+1-hB);
4373 Window(w, f[p], N, M, c, fst, fst+B-1);
4374 double xw2=~Inner(B, &x[fst], w);
4375 a[p]=sqrt(xw2)*iH2;
4376 if (rsr)
4377 {
4378 cdouble xx=Inner(B, &x[fst], &x[fst]);
4379 if (xx.x==0) rsr[p]=1;
4380 else rsr[p]=1-xw2/xx.x*iH2;
4381 }
4382 p++;
4383 }
4384 }
4385 a0=a1;
4386 a1=a2;
4387 n++;
4388 }
4389 delete[] w;
4390 return p;
4391 }//QuickPeaks
4392
4393 /*
4394 function QuickPeaks: finds rough peaks in the spectrogram (peak picking) for constant-frequency
4395 sinusoids
4396
4397 In: x[Fr][N/2+1]: spectrogram
4398 fr0, r0: centre and half width of interval (in frames) to use for peak picking
4399 M, c[], iH2: cosine-family window function specifiers
4400 mina: minimal amplitude to spot a spectral peak
4401 [binst, binen): frequency range, in bins, to look for peaks
4402 B: spectral truncation width
4403 Out; f[return value], a[return value]: frequencies (in bins) and summary amplitudes of found peaks
4404 rsr[return value]: residue-sinusoid-ratio of found peaks, optional
4405
4406 Returns the number of peaks found. f[] and a[] must be allocated enough space before calling.
4407 */
4408 int QuickPeaks(double* f, double* a, int Fr, int N, cdouble** x, int fr0, int r0, int M, double* c, double iH2, double mina, int binst, int binen, int B, double* rsr)
4409 {
4410 double hB=B*0.5;
4411 int hN=N/2;
4412 if (binst<hB) binst=hB;
4413 if (binen<0 || binen>hN-hB) binen=hN-hB;
4414 double minA=mina*mina*2/iH2;
4415
4416 double *a0s=new double[hN*3*r0], *a1s=&a0s[hN*r0], *a2s=&a0s[hN*2*r0];
4417 int *idx=new int[hN*r0], *frc=new int[hN*r0];
4418 int** Allocate2(int, (hN*r0), (r0*2+1), frs);
4419
4420 int pc=0;
4421
4422 int lr=0;
4423 while (lr<=r0)
4424 {
4425 int fr[2]; //indices to frames to process for this lr
4426 if (lr==0) fr[0]=fr0, fr[1]=-1;
4427 else
4428 {
4429 fr[0]=fr0-lr, fr[1]=fr0+lr;
4430 if (fr[1]>=Fr) fr[1]=-1;
4431 if (fr[0]<0) {fr[0]=fr[1]; fr[1]=-1;}
4432 }
4433 for (int i=0; i<2; i++)
4434 {
4435 if (fr[i]>=0)
4436 {
4437 cdouble* xfr=x[fr[i]];
4438 for (int k=binst; k<binen; k++)
4439 {
4440 double a0=~xfr[k-1], a1=~xfr[k], a2=~xfr[k+1];
4441 if (a1>=a0 && a1>=a2)
4442 {
4443 if (a1>minA)
4444 {
4445 double A1=sqrt(a1), A2=sqrt(a2), A0=sqrt(a0);
4446 double lf=k+(A0-A2)/2/(A0+A2-2*A1);
4447 if (lr==0) //starting frame
4448 {
4449 f[pc]=lf;
4450 a0s[pc]=a0, a1s[pc]=a1, a2s[pc]=a2;
4451 idx[pc]=pc; frs[pc][0]=fr[i]; frc[pc]=1;
4452 pc++;
4453 }
4454 else
4455 {
4456 int closep, indexp=InsertInc(lf, f, pc, false); //find the closest peak ever found in previous (i.e. closer to fr0) frames
4457 if (indexp==-1) indexp=pc, closep=pc-1;
4458 else if (indexp-1>=0 && fabs(lf-f[indexp-1])<fabs(lf-f[indexp])) closep=indexp-1;
4459 else closep=indexp;
4460
4461 if (fabs(lf-f[closep])<0.2) //i.e. lf is very close to previous peak at fs[closep]
4462 {
4463 int idxp=idx[closep];
4464 a0s[idxp]+=a0, a1s[idxp]+=a1, a2s[idxp]+=a2;
4465 frs[idxp][frc[idxp]]=fr[i]; frc[idxp]++;
4466 double A0=sqrt(a0s[idxp]), A1=sqrt(a1s[idxp]), A2=sqrt(a2s[idxp]);
4467 f[closep]=k+(A0-A2)/2/(A0+A2-2*A1);
4468 }
4469 else
4470 {
4471 memmove(&f[indexp+1], &f[indexp], sizeof(double)*(pc-indexp)); f[indexp]=lf;
4472 memmove(&idx[indexp+1], &idx[indexp], sizeof(int)*(pc-indexp)); idx[indexp]=pc;
4473 a0s[pc]=a0, a1s[pc]=a1, a2s[pc]=a2, frs[pc][0]=fr[i], frc[pc]=1;
4474 pc++;
4475 }
4476 }
4477 }
4478 }
4479 }
4480 }
4481 }
4482 lr++;
4483 }
4484
4485 cdouble* w=new cdouble[B];
4486 int* frused=new int[Fr];
4487 for (int p=0; p<pc; p++)
4488 {
4489 int fst=floor(f[p]+1-hB);
4490 memset(frused, 0, sizeof(int)*Fr);
4491 int idxp=idx[p];
4492
4493 Window(w, f[p], N, M, c, fst, fst+B-1);
4494 double xw2=0, xw2r=0, xxr=0;
4495
4496 for (int ifr=0; ifr<frc[idxp]; ifr++)
4497 {
4498 int fr=frs[idxp][ifr];
4499 frused[fr]=1;
4500 xw2r+=~Inner(B, &x[fr][fst], w);
4501 xxr+=Inner(B, &x[fr][fst], &x[fr][fst]).x;
4502 }
4503 xw2=xw2r;
4504 for (int fr=0; fr<Fr; fr++)
4505 if (!frused[fr]) xw2+=~Inner(B, &x[fr][fst], w);
4506 a[p]=sqrt(xw2/Fr)*iH2;
4507 if (xxr==0) rsr[p]=1;
4508 else rsr[p]=1-xw2r/xxr*iH2;
4509 }
4510 delete[] w;
4511 delete[] frused;
4512 delete[] a0s;
4513 delete[] idx;
4514 delete[] frc;
4515 DeAlloc2(frs);
4516 return pc;
4517 }//QuickPeaks
4518
4519 /*
4520 function ReEstHS1: reestimates a harmonic sinusoid by one multiplicative reestimation using phasor
4521 multiplier
4522
4523 In: partials[M][Fr]: HS partials
4524 [frst, fren): frame (measurement point) range on which to do reestimation
4525 Data16[]: waveform data
4526 Out: partials[M][Fr]: updated HS partials
4527
4528 No return value.
4529 */
4530 void ReEstHS1(int M, int Fr, int frst, int fren, atom** partials, __int16* Data16)
4531 {
4532 double* fs=new double[Fr*3], *as=&fs[Fr], *phs=&fs[Fr*2];
4533 int wid=partials[0][0].s, offst=partials[0][1].t-partials[0][0].t;
4534 int ldatastart=partials[0][0].t-wid/2;
4535 __int16* ldata16=&Data16[ldatastart];
4536 for (int m=0; m<M; m++)
4537 {
4538 for (int fr=0; fr<Fr; fr++) {fs[fr]=partials[m][fr].f; if (fs[fr]<=0){delete[] fs; return;} as[fr]=partials[m][fr].a, phs[fr]=partials[m][fr].p;}
4539 MultiplicativeUpdateF(fs, as, phs, ldata16, Fr, frst, fren, wid, offst);
4540 for (int fr=0; fr<Fr; fr++) partials[m][fr].f=fs[fr], partials[m][fr].a=as[fr], partials[m][fr].p=phs[fr];
4541 }
4542 delete[] fs;
4543 }//ReEstHS1
4544
4545 /*
4546 function ReEstHS1: wrapper function.
4547
4548 In: HS: a harmonic sinusoid
4549 Data16: its waveform data
4550 Out: HS: updated harmonic sinusoid
4551
4552 No return value.
4553 */
4554 void ReEstHS1(THS* HS, __int16* Data16)
4555 {
4556 ReEstHS1(HS->M, HS->Fr, 0, HS->Fr, HS->Partials, Data16);
4557 }//ReEstHS1
4558
4559 /*
4560 function SortCandid: inserts a candid object into a listed of candid objects sorted first by f, then
4561 (for identical f's) by df.
4562
4563 In: cand: the candid object to insert to the list
4564 cands[newN]: the sorted list of candid objects
4565 Out: cands[newN+1]: the sorted list after the insertion
4566
4567 Returns the index of $cand in the new list.
4568 */
4569 int SortCandid(candid cand, candid* cands, int newN)
4570 {
4571 int lnN=newN-1;
4572 while (lnN>=0 && cands[lnN].f>cand.f) lnN--;
4573 while (lnN>=0 && cands[lnN].f==cand.f && cands[lnN].df>cand.df) lnN--;
4574 lnN++; //now insert cantmp as cands[lnN]
4575 memmove(&cands[lnN+1], &cands[lnN], sizeof(candid)*(newN-lnN));
4576 cands[lnN]=cand;
4577 return lnN;
4578 }//SortCandid
4579
4580 /*
4581 function SynthesisHS: synthesizes a harmonic sinusoid without aligning the phases
4582
4583 In: partials[M][Fr]: HS partials
4584 terminatetag: external termination flag. Function SynthesisHS() polls *terminatetag and exits with
4585 0 when it is set.
4586 Out: [dst, den): time interval synthesized
4587 xrec[den-dst]: resynthesized harmonic sinusoid
4588
4589 Returns pointer to xrec on normal finish, or 0 on external termination by setting $terminatetag. In
4590 the first case xrec is created anew with malloc8() and must be freed by caller using free8().
4591 */
4592 double* SynthesisHS(int M, int Fr, atom** partials, int& dst, int& den, bool* terminatetag)
4593 {
4594 int wid=partials[0][0].s, hwid=wid/2;
4595 double *as=(double*)malloc8(sizeof(double)*Fr*13);
4596 double *fs=&as[Fr], *f3=&as[Fr*2], *f2=&as[Fr*3], *f1=&as[Fr*4], *f0=&as[Fr*5],
4597 *a3=&as[Fr*6], *a2=&as[Fr*7], *a1=&as[Fr*8], *a0=&as[Fr*9], *xs=&as[Fr*11];
4598 int* ixs=(int*)&as[Fr*12];
4599
4600 dst=partials[0][0].t-hwid, den=partials[0][Fr-1].t+hwid;
4601 double* xrec=(double*)malloc8(sizeof(double)*(den-dst));
4602 memset(xrec, 0, sizeof(double)*(den-dst));
4603
4604 for (int m=0; m<M; m++)
4605 {
4606 atom* part=partials[m]; bool fzero=false;
4607 for (int fr=0; fr<Fr; fr++)
4608 {
4609 if (part[fr].f<=0){fzero=true; break;}
4610 ixs[fr]=part[fr].t;
4611 xs[fr]=part[fr].t;
4612 as[fr]=part[fr].a*2;
4613 fs[fr]=part[fr].f;
4614 if (part[fr].type==atMuted) as[fr]=0;
4615 }
4616 if (fzero) break;
4617 if (terminatetag && *terminatetag) {free8(xrec); free8(as); return 0;}
4618
4619 CubicSpline(Fr-1, f3, f2, f1, f0, xs, fs, 1, 1);
4620 CubicSpline(Fr-1, a3, a2, a1, a0, xs, as, 1, 1);
4621 double ph=0, ph0=0;
4622 for (int fr=0; fr<Fr-1; fr++)
4623 {
4624 part[fr].p=ph;
4625 ALIGN8(Sinusoid(&xrec[ixs[fr]-dst], 0, ixs[fr+1]-ixs[fr], a3[fr], a2[fr], a1[fr], a0[fr], f3[fr], f2[fr], f1[fr], f0[fr], ph, true);)
4626 if (terminatetag && *terminatetag) {free8(xrec); free8(as); return 0;}
4627 }
4628 part[Fr-1].p=ph;
4629 ALIGN8(Sinusoid(&xrec[ixs[Fr-2]-dst], ixs[Fr-1]-ixs[Fr-2], den-ixs[Fr-2], a3[Fr-2], a2[Fr-2], a1[Fr-2], a0[Fr-2], f3[Fr-2], f2[Fr-2], f1[Fr-2], f0[Fr-2], ph, true);
4630 Sinusoid(&xrec[ixs[0]-dst], dst-ixs[0], 0, a3[0], a2[0], a1[0], a0[0], f3[0], f2[0], f1[0], f0[0], ph0, true);)
4631 }
4632 free8(as);
4633 return xrec;
4634 }//SynthesisHS
4635
4636 /*
4637 function SynthesisHS: synthesizes a perfectly harmonic sinusoid without aligning the phases.
4638 Frequencies of partials above the fundamental are not used in this synthesis process.
4639
4640 In: partials[M][Fr]: HS partials
4641 terminatetag: external termination flag. This function polls *terminatetag and exits with 0 when
4642 it is set.
4643 Out: [dst, den) time interval synthesized
4644 xrec[den-dst]: resynthesized harmonic sinusoid
4645
4646 Returns pointer to xrec on normal finish, or 0 on external termination by setting $terminatetag. In
4647 the first case xrec is created anew with malloc8() and must be freed by caller using free8().
4648 */
4649 double* SynthesisHS2(int M, int Fr, atom** partials, int& dst, int& den, bool* terminatetag)
4650 {
4651 int wid=partials[0][0].s, hwid=wid/2;
4652 double *as=(double*)malloc8(sizeof(double)*Fr*7); //*as=new double[Fr*8],
4653 double *xs=&as[Fr], *f3=&as[Fr*2], *f2=&as[Fr*3], *f1=&as[Fr*4], *f0=&as[Fr*5];
4654 int *ixs=(int*)&as[Fr*6];
4655 double** a0=new double*[M*4], **a1=&a0[M], **a2=&a0[M*2], **a3=&a0[M*3];
4656 a0[0]=(double*)malloc8(sizeof(double)*M*Fr*4); //a0[0]=new double[M*Fr*4];
4657 for (int m=0; m<M; m++) a0[m]=&a0[0][m*Fr*4], a1[m]=&a0[m][Fr], a2[m]=&a0[m][Fr*2], a3[m]=&a0[m][Fr*3];
4658
4659 dst=partials[0][0].t-hwid, den=partials[0][Fr-1].t+hwid;
4660 double* xrec=(double*)malloc8(sizeof(double)*(den-dst));
4661 memset(xrec, 0, sizeof(double)*(den-dst));
4662 atom* part=partials[0];
4663
4664 for (int fr=0; fr<Fr; fr++)
4665 {
4666 xs[fr]=ixs[fr]=part[fr].t;
4667 as[fr]=part[fr].f;
4668 }
4669 CubicSpline(Fr-1, f3, f2, f1, f0, xs, as, 1, 1);
4670
4671 for (int m=0; m<M; m++)
4672 {
4673 part=partials[m];
4674 for (int fr=0; fr<Fr; fr++)
4675 {
4676 if (part[fr].f<=0){M=m; break;}
4677 as[fr]=part[fr].a*2;
4678 }
4679 if (terminatetag && *terminatetag) {free8(xrec); free8(as); free8(a0[0]); delete[] a0; return 0;}
4680 if (m>=M) break;
4681
4682 CubicSpline(Fr-1, a3[m], a2[m], a1[m], a0[m], xs, as, 1, 1);
4683 }
4684
4685 double *ph=(double*)malloc8(sizeof(double)*M*2), *ph0=&ph[M]; memset(ph, 0, sizeof(double)*M*2);
4686 for (int fr=0; fr<Fr-1; fr++)
4687 {
4688 // double *la0=new double[M*4], *la1=&la0[M], *la2=&la0[M*2], *la3=&la0[M*3]; for (int m=0; m<M; m++) la0[m]=a0[m][fr], la1[m]=a1[m][fr], la2[m]=a2[m][fr], la3[m]=a3[m][fr], part[fr].p=ph[m]; Sinusoids(M, &xrec[ixs[fr]-dst], 0, ixs[fr+1]-ixs[fr], la3, la2, la1, la0, f3[fr], f2[fr], f1[fr], f0[fr], ph, true); delete[] la0;
4689 for (int m=0; m<M; m++) ALIGN8(Sinusoid(&xrec[ixs[fr]-dst], 0, ixs[fr+1]-ixs[fr], a3[m][fr], a2[m][fr], a1[m][fr], a0[m][fr], (m+1)*f3[fr], (m+1)*f2[fr], (m+1)*f1[fr], (m+1)*f0[fr], ph[m], true);)
4690 }
4691 for (int m=0; m<M; m++)
4692 {
4693 part[Fr-1].p=ph[m];
4694 ALIGN8(Sinusoid(&xrec[ixs[Fr-2]-dst], ixs[Fr-1]-ixs[Fr-2], den-ixs[Fr-2], a3[m][Fr-2], a2[m][Fr-2], a1[m][Fr-2], a0[m][Fr-2], (m+1)*f3[Fr-2], (m+1)*f2[Fr-2], (m+1)*f1[Fr-2], (m+1)*f0[Fr-2], ph[m], true);
4695 Sinusoid(&xrec[ixs[0]-dst], dst-ixs[0], 0, a3[m][0], a2[m][0], a1[m][0], a0[m][0], (m+1)*f3[0], (m+1)*f2[0], (m+1)*f1[0], (m+1)*f0[0], ph0[m], true);)
4696 if (terminatetag && *terminatetag) {free8(xrec); free8(as); free8(a0[0]); delete[] a0; free8(ph); return 0;}
4697 }
4698 free8(as); free8(ph); free8(a0[0]); delete[] a0;
4699 return xrec;
4700 }//SynthesisHS2
4701
4702 /*
4703 function SynthesisHSp: synthesizes a harmonic sinusoid with phase alignment
4704
4705 In: partials[pm][pfr]: HS partials
4706 startamp[pm][st_count]: onset amplifiers, optional
4707 st_start, st_offst: start of and interval between onset amplifying points, optional
4708 Out: [dst, den): time interval synthesized
4709 xrec[den-dst]: resynthesized harmonic sinusoid
4710
4711 Returns pointer to xrec. xrec is created anew with malloc8() and must be freed by caller with free8().
4712 */
4713 double* SynthesisHSp(int pm, int pfr, atom** partials, int& dst, int& den, double** startamp, int st_start, int st_offst, int st_count)
4714 {
4715 int wid=partials[0][0].s, hwid=wid/2, offst=partials[0][1].t-partials[0][0].t;
4716 double *a1=new double[pfr*13];
4717 double *f1=&a1[pfr], *fa=&a1[pfr*2], *fb=&a1[pfr*3], *fc=&a1[pfr*4], *fd=&a1[pfr*5],
4718 *aa=&a1[pfr*6], *ab=&a1[pfr*7], *ac=&a1[pfr*8], *ad=&a1[pfr*9], *p1=&a1[pfr*10], *xs=&a1[pfr*11];
4719 int* ixs=(int*)&a1[pfr*12];
4720
4721 dst=partials[0][0].t-hwid, den=partials[0][pfr-1].t+hwid;
4722 double *xrec=(double*)malloc8(sizeof(double)*(den-dst)*2), *xrecm=&xrec[den-dst];
4723 memset(xrec, 0, sizeof(double)*(den-dst));
4724
4725 for (int p=0; p<pm; p++)
4726 {
4727 atom* part=partials[p];
4728 bool fzero=false;
4729 for (int fr=0; fr<pfr; fr++)
4730 {
4731 if (part[fr].f<=0)
4732 {
4733 fzero=true;
4734 break;
4735 }
4736 ixs[fr]=part[fr].t;
4737 xs[fr]=part[fr].t;
4738 a1[fr]=part[fr].a*2;
4739 f1[fr]=part[fr].f;
4740 p1[fr]=part[fr].p;
4741 if (part[fr].type==atMuted) a1[fr]=0;
4742 }
4743 if (fzero) break;
4744
4745 CubicSpline(pfr-1, fa, fb, fc, fd, xs, f1, 1, 1);
4746 CubicSpline(pfr-1, aa, ab, ac, ad, xs, a1, 1, 1);
4747
4748 for (int fr=0; fr<pfr-1; fr++) Sinusoid(&xrecm[ixs[fr]-dst], 0, offst, aa[fr], ab[fr], ac[fr], ad[fr], fa[fr], fb[fr], fc[fr], fd[fr], p1[fr], p1[fr+1], false);
4749 // Sinusoid(&xrecm[ixs[0]-dst], -hwid, 0, aa[0], ab[0], ac[0], ad[0], fa[0], fb[0], fc[0], fd[0], p1[0], p1[1], false);
4750 // Sinusoid(&xrecm[ixs[pfr-2]-dst], offst, offst+hwid, aa[pfr-2], ab[pfr-2], ac[pfr-2], ad[pfr-2], fa[pfr-2], fb[pfr-2], fc[pfr-2], fd[pfr-2], p1[pfr-2], p1[pfr-1], false);
4751 Sinusoid(&xrecm[ixs[0]-dst], -hwid, 0, 0, 0, 0, ad[0], fa[0], fb[0], fc[0], fd[0], p1[0], p1[1], false);
4752 Sinusoid(&xrecm[ixs[pfr-2]-dst], offst, offst+hwid, 0, 0, 0, ad[pfr-2], fa[pfr-2], fb[pfr-2], fc[pfr-2], fd[pfr-2], p1[pfr-2], p1[pfr-1], false);
4753
4754 if (st_count)
4755 {
4756 double* amp=startamp[p];
4757 for (int c=0; c<st_count; c++)
4758 {
4759 double a1=amp[c], a2=(c+1<st_count)?amp[c+1]:1, da=(a2-a1)/st_offst;
4760 int lst=ixs[0]-dst+st_start+c*st_offst;
4761 double *lxrecm=&xrecm[lst];
4762 if (lst>0) for (int i=0; i<st_offst; i++) lxrecm[i]*=(a1+da*i);
4763 else for (int i=-lst; i<st_offst; i++) lxrecm[i]*=(a1+da*i);
4764 }
4765 for (int i=0; i<ixs[0]-dst+st_start; i++) xrecm[i]*=amp[0];
4766 }
4767 else
4768 {
4769 /*
4770 for (int i=0; i<=hwid; i++)
4771 {
4772 double tmp=0.5+0.5*cos(M_PI*i/hwid);
4773 xrecm[ixs[0]-dst-i]*=tmp;
4774 xrecm[ixs[pfr-2]-dst+offst+i]*=tmp;
4775 } */
4776 }
4777
4778 for (int n=0; n<den-dst; n++) xrec[n]+=xrecm[n];
4779 }
4780
4781 delete[] a1;
4782 return xrec;
4783 }//SynthesisHSp
4784
4785 /*
4786 function SynthesisHSp: wrapper function.
4787
4788 In: HS: a harmonic sinusoid.
4789 Out: [dst, den): time interval synthesized
4790 xrec[den-dst]: resynthesized harmonic sinusoid
4791
4792 Returns pointer to xrec, which is created anew with malloc8() and must be freed by caller using free8().
4793 */
4794 double* SynthesisHSp(THS* HS, int& dst, int& den)
4795 {
4796 return SynthesisHSp(HS->M, HS->Fr, HS->Partials, dst, den, HS->startamp, HS->st_start, HS->st_offst, HS->st_count);
4797 }//SynthesisHSp
4798