Mercurial > hg > x
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=¶ms; | |
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=¶ms; | |
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=¶ms; | |
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=¶ms; | |
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=¶ms; | |
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=¶ms; | |
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=¶ms; | |
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=¶ms; | |
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=¶ms; | |
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=¶ms; | |
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=¶ms; | |
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=¶ms; | |
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 |