Add dywapitchtracker

Improvements for PtAKF (finer resolution at end)
Bugfix for WavFile

PitchTracker/Helper.cpp

@@ -0,0 +1,43 @@
+#include "Helper.h"
+#include <stdlib.h>
+#include <cmath>
+
+float* short2FloatArray(short* in, size_t len){
+	const float maxShort = 32767.0;
+
+	float* result = static_cast<float*>(malloc(len * sizeof(float)));
+	float* curOut = result;
+	for(size_t i=0; i<len; i++){
+		*curOut = *in / maxShort;
+		curOut++;
+		in++;
+	}
+	return result;
+}
+
+double* short2DoubleArray(short* in, size_t len, bool scale){
+	const double maxShort = 32767.0;
+
+	double* result = static_cast<double*>(malloc(len * sizeof(double)));
+	double* curOut = result;
+	for(size_t i=0; i<len; i++){
+		*curOut = *in;
+		if(scale)
+			*curOut /= maxShort;
+		curOut++;
+		in++;
+	}
+	return result;
+}
+
+void freeFloatArray(float* floats){
+	free(floats);
+}
+
+void freeDoubleArray(double* floats){
+	free(floats);
+}
+
+double FreqToNote(double freq){
+	return 12.0 * std::log(freq / 127.09) / std::log(2.0) + 27.5 - 16; // 16 = C2
+}

PitchTracker/Helper.h

@@ -1,47 +1,10 @@
 #pragma once
-#include "PitchWrapper.h"
-#include <stdlib.h>
-#include <cmath>
 
-float* short2FloatArray(short* in, size_t len){
-	const float maxShort = 32767.0;
-
-	float* result = static_cast<float*>(malloc(len * sizeof(float)));
-	float* curOut = result;
-	for(size_t i=0; i<len; i++){
-		*curOut = *in / maxShort;
-		curOut++;
-		in++;
-	}
-	return result;
-}
-
-double* short2DoubleArray(short* in, size_t len, bool scale = true){
-	const double maxShort = 32767.0;
-
-	double* result = static_cast<double*>(malloc(len * sizeof(double)));
-	double* curOut = result;
-	for(size_t i=0; i<len; i++){
-		*curOut = *in;
-		if(scale)
-			*curOut /= maxShort;
-		curOut++;
-		in++;
-	}
-	return result;
-}
-
-void freeFloatArray(float* floats){
-	free(floats);
-}
-
-void freeDoubleArray(double* floats){
-	free(floats);
-}
-
-double FreqToNote(double freq){
-	return 12.0 * std::log(freq / 127.09) / std::log(2.0) + 27.5 - 16; // 16 = C2
-}
+float* short2FloatArray(short* in, size_t len);
+double* short2DoubleArray(short* in, size_t len, bool scale = true);
+void freeFloatArray(float* floats);
+void freeDoubleArray(double* floats);
+double FreqToNote(double freq);
 
 /*
 namespace Native{

PitchTracker/PitchTracker.vcxproj

@@ -170,12 +170,17 @@
   </ItemDefinitionGroup>
   <ItemGroup>
     <ClCompile Include="compatibility.cpp" />
+    <ClCompile Include="dywapitchtrack\dywapitchtrack.cpp" />
+    <ClCompile Include="dywapitchtrack\ptDyWa.cpp" />
+    <ClCompile Include="Helper.cpp" />
     <ClCompile Include="performous\pitch.cc" />
     <ClCompile Include="ptAKF.cpp" />
     <ClCompile Include="PitchWrapper.cpp" />
   </ItemGroup>
   <ItemGroup>
     <ClInclude Include="compatibility.h" />
+    <ClInclude Include="dywapitchtrack\dywapitchtrack.h" />
+    <ClInclude Include="dywapitchtrack\ptDyWa.h" />
     <ClInclude Include="Helper.h" />
     <ClInclude Include="performous\libda\fft.hpp" />
     <ClInclude Include="performous\libda\sample.hpp" />

PitchTracker/PitchTracker.vcxproj.filters

@@ -22,6 +22,12 @@
     <Filter Include="Headerdateien\performous\libda">
       <UniqueIdentifier>{d71782a0-47fc-4fe1-8569-0a0ce8fe9b45}</UniqueIdentifier>
     </Filter>
+    <Filter Include="Quelldateien\dywapitchtrack">
+      <UniqueIdentifier>{e8eab624-6cc3-42db-93d9-575aaeff06f7}</UniqueIdentifier>
+    </Filter>
+    <Filter Include="Headerdateien\dywapitchtrack">
+      <UniqueIdentifier>{2a09cee1-51b7-4c24-aa9f-32bec53c3a36}</UniqueIdentifier>
+    </Filter>
   </ItemGroup>
   <ItemGroup>
     <ClCompile Include="PitchWrapper.cpp">
@@ -36,6 +42,15 @@
     <ClCompile Include="ptAKF.cpp">
       <Filter>Quelldateien</Filter>
     </ClCompile>
+    <ClCompile Include="dywapitchtrack\ptDyWa.cpp">
+      <Filter>Quelldateien\dywapitchtrack</Filter>
+    </ClCompile>
+    <ClCompile Include="dywapitchtrack\dywapitchtrack.cpp">
+      <Filter>Quelldateien\dywapitchtrack</Filter>
+    </ClCompile>
+    <ClCompile Include="Helper.cpp">
+      <Filter>Quelldateien</Filter>
+    </ClCompile>
   </ItemGroup>
   <ItemGroup>
     <ClInclude Include="compatibility.h">
@@ -56,6 +71,12 @@
     <ClInclude Include="ptAKF.h">
       <Filter>Headerdateien</Filter>
     </ClInclude>
+    <ClInclude Include="dywapitchtrack\dywapitchtrack.h">
+      <Filter>Headerdateien\dywapitchtrack</Filter>
+    </ClInclude>
+    <ClInclude Include="dywapitchtrack\ptDyWa.h">
+      <Filter>Headerdateien\dywapitchtrack</Filter>
+    </ClInclude>
   </ItemGroup>
   <ItemGroup>
     <None Include="performous\pitch.hh">

PitchTracker/PitchWrapper.cpp

@@ -58,23 +58,23 @@ bool Analyzer_OutputFloat(Analyzer* analyzer, float* data, int sampleCt, float r
     return analyzer->output(data, data + sampleCt, rate);
 }
 
-DllExport PtAKF* PtAKF_Create(unsigned step){
+PtAKF* PtAKF_Create(unsigned step){
 	return new PtAKF(step);
 }
 
-DllExport void PtAKF_Free(PtAKF* analyzer){
+void PtAKF_Free(PtAKF* analyzer){
 	if(analyzer)
 		delete analyzer;
 }
 
-DllExport int PtAKF_GetNumHalfTones(){
+int PtAKF_GetNumHalfTones(){
 	return PtAKF::GetNumHalfTones();
 }
 
-DllExport void PtAKF_InputByte(PtAKF* analyzer, char* data, int sampleCt){
+void PtAKF_InputByte(PtAKF* analyzer, char* data, int sampleCt){
 	if(sampleCt <= 0 || !analyzer)
 		return;
-	short* dataShort = static_cast<short*>(static_cast<void*>(data));
+	//short* dataShort = static_cast<short*>(static_cast<void*>(data));
 	float* dataFloat = short2FloatArray(static_cast<short*>(static_cast<void*>(data)), sampleCt);
 	analyzer->input(dataFloat, dataFloat + sampleCt);
 	freeFloatArray(dataFloat);
@@ -82,12 +82,41 @@ DllExport void PtAKF_InputByte(PtAKF* analyzer, char* data, int sampleCt){
 	//analyzer->input(dataShort, dataShort + sampleCt);
 }
 
-DllExport int PtAKF_GetNote(PtAKF* analyzer, double* maxVolume, float* weights){
+int PtAKF_GetNote(PtAKF* analyzer, double* maxVolume, float* weights){
 	if(!analyzer)
 		return -1;
 	return analyzer->GetNote(maxVolume, weights);
 }
 
+PtDyWa* PtDyWa_Create(unsigned step){
+	return new PtDyWa(step);
+}
+
+void PtDyWa_Free(PtDyWa* analyzer){
+	if(analyzer)
+		delete analyzer;
+}
+
+void PtDyWa_InputByte(PtDyWa* analyzer, char* data, int sampleCt){
+	if(sampleCt <= 0 || !analyzer)
+		return;
+	float* dataFloat = short2FloatArray(static_cast<short*>(static_cast<void*>(data)), sampleCt);
+	analyzer->input(dataFloat, dataFloat + sampleCt);
+	freeFloatArray(dataFloat);
+}
+
+void PtDyWa_Process(PtDyWa* analyzer){
+	if(!analyzer)
+		return;
+	analyzer->Process();
+}
+
+double PtDyWa_GetNote(PtDyWa* analyzer){
+	if(!analyzer)
+		return -1;
+	return analyzer->GetNote();
+}
+
 /*
 namespace Native{
 	namespace PitchTracking{

PitchTracker/PitchWrapper.h

@@ -2,6 +2,7 @@
 
 #include "performous/pitch.hh"
 #include "ptAKF.h"
+#include "dywapitchtrack/ptDyWa.h"
 
 #define DllExport extern "C" __declspec(dllexport)
 
@@ -21,6 +22,12 @@ DllExport int PtAKF_GetNumHalfTones();
 DllExport void PtAKF_InputByte(PtAKF* analyzer, char* data, int sampleCt);
 DllExport int PtAKF_GetNote(PtAKF* analyzer, double* maxVolume, float* weights);
 
+DllExport PtDyWa* PtDyWa_Create(unsigned step);
+DllExport void PtDyWa_Free(PtDyWa* analyzer);
+DllExport void PtDyWa_InputByte(PtDyWa* analyzer, char* data, int sampleCt);
+DllExport void PtDyWa_Process(PtDyWa* analyzer);
+DllExport double PtDyWa_GetNote(PtDyWa* analyzer);
+
 /*
 public ref class CTone {
 public:

PitchTracker/dywapitchtrack/DOCUMENTATION.txt

@@ -0,0 +1,29 @@
+
+dywapitchtrack pitch tracking library
+Documentation
+
+The dywapitchtrack library computes the pitch of an audio stream in real time. The pitch is the main frequency of the waveform (the 'note' being played or sung). It is expressed as a float in Hz.
+
+Unlike the human ear, pitch detection is difficult to achieve for computers. Many algorithms have been designed and experimented, but there is no 'best' algorithm. They all depend on the context and the tradeoffs acceptable in terms of speed and latency. The context includes the quality and 'cleanness' of the audio : obviously polyphonic sounds (multiple instruments playing different notes at the same time) are extremely difficult to track, percussive or noisy audio has no pitch, most real-life audio have some noisy moments, some instruments have a lot of harmonics, etc...
+
+The dywapitchtrack is based on a custom-tailored algorithm which is of very high quality: both very accurate (precision < 0.05 semitones), very low latency (< 23 ms) and very low error rate. It has been thoroughly tested on human voice.
+
+It can best be described as a dynamic wavelet algorithm (dywa):
+
+The heart of the algorithm is a very powerful wavelet algorithm, described in a paper by Eric Larson and Ross Maddox :
+"Real-Time Time-Domain Pitch Tracking Using Wavelets" http://online.physics.uiuc.edu/courses/phys406/NSF_REU_Reports/2005_reu/Real-Time_Time-Domain_Pitch_Tracking_Using_Wavelets.pdf
+
+This central algorithm has been improved by adding dynamic tracking, to reduce the common problems of frequency-halving and voiced/unvoiced errors. This dynamic tracking explains the need for a tracking structure (dywapitchtracker). The dynamic tracking assumes that the main function dywapitch_computepitch is called repeatedly, as it follows the pitch over time and makes assumptions about human voice capabilities and reallife conditions (as documented inside the code).
+
+Note : The algorithm currently assumes a 44100Hz audio sampling rate. If you use a different samplerate, you can just multiply the resulting pitch by the ratio between your samplerate and 44100.
+
+See the dywapitchtrack.h file for the library API documentation.
+
+History :
+This algorithm has been designed during a mission for a customer. As I had kept the author's rights on the source code, I eventually decided to make this code public and open source, because it is of really high quality. This algorithms has been extensibely tested, especially it has been included in an Adobe Director Xtra (plugin), asPitchDetect Xtra, that has been distributed widely for a few years.
+
+I hope that it can be used with pleasure and efficiency in your project.
+
+Thank you
+Antoine Schmitt April 2010
+

PitchTracker/dywapitchtrack/README.txt

@@ -0,0 +1,43 @@
+dywapitchtrack
+Dynamic Wavelet Algorithm Pitch Tracking library
+
+Copyright (c) 2010 Antoine Schmitt as@schmittmachine.com
+
+-- Documentation
+
+See the dywapitchtrack.h file for the library API documentation.
+
+
+-- Version history
+
+Inital release : Antoine Schmitt April 2010
+jan 2012 : Minor update : changed the link to the Eric Larson and Ross Maddox paper
+apr 2012 : Included the SpeakHere iPhone app example
+
+-- Building, Installing
+
+- the algorithm can be included directly in your project by including the .h and .c files. 
+- the iPhone/dywapitchtrack.xcodeproj builds an iPhone static library, against which one can statically link an iPhone application.
+
+
+-- MIT open source licence
+ 
+ Copyright (c) 2010 Antoine Schmitt
+  
+ Permission is hereby granted, free of charge, to any person obtaining a copy
+ of this software and associated documentation files (the "Software"), to deal
+ in the Software without restriction, including without limitation the rights
+ to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+ copies of the Software, and to permit persons to whom the Software is
+ furnished to do so, subject to the following conditions:
+ 
+ The above copyright notice and this permission notice shall be included in
+ all copies or substantial portions of the Software.
+ 
+ THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+ IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+ AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+ THE SOFTWARE.

PitchTracker/dywapitchtrack/dywapitchtrack.cpp

@@ -0,0 +1,448 @@
+/* dywapitchtrack.c
+ 
+ Dynamic Wavelet Algorithm Pitch Tracking library
+ Released under the MIT open source licence
+  
+ Copyright (c) 2010 Antoine Schmitt
+ 
+ Permission is hereby granted, free of charge, to any person obtaining a copy
+ of this software and associated documentation files (the "Software"), to deal
+ in the Software without restriction, including without limitation the rights
+ to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+ copies of the Software, and to permit persons to whom the Software is
+ furnished to do so, subject to the following conditions:
+ 
+ The above copyright notice and this permission notice shall be included in
+ all copies or substantial portions of the Software.
+ 
+ THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+ IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+ AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+ THE SOFTWARE.
+*/
+
+#include "dywapitchtrack.h"
+#include <math.h>
+#include <stdlib.h>
+#include <string.h> // for memset
+
+
+//**********************
+//       Utils
+//**********************
+
+#ifndef max
+#define max(x, y) ((x) > (y)) ? (x) : (y)
+#endif
+#ifndef min
+#define min(x, y) ((x) < (y)) ? (x) : (y)
+#endif
+
+// returns 1 if power of 2
+int _power2p(int value) {
+	if (value == 0) return 1;
+	if (value == 2) return 1;
+	if (value & 0x1) return 0;
+	return (_power2p(value >> 1));
+}
+
+// count number of bits
+int _bitcount(int value) {
+	if (value == 0) return 0;
+	if (value == 1) return 1;
+	if (value == 2) return 2;
+	return _bitcount(value >> 1) + 1;
+}
+
+// closest power of 2 above or equal
+int _ceil_power2(int value) {
+	if (_power2p(value)) return value;
+	
+	if (value == 1) return 2;
+	int j, i = _bitcount(value);
+	int res = 1;
+	for (j = 0; j < i; j++) res <<= 1;
+	return res;
+}
+
+// closest power of 2 below or equal
+int _floor_power2(int value) {
+	if (_power2p(value)) return value;
+	return _ceil_power2(value)/2;
+}
+
+// abs value
+int _iabs(int x) {
+	if (x >= 0) return x;
+	return -x;
+}
+
+// 2 power
+int _2power(int i) {
+	int res = 1, j;
+	for (j = 0; j < i; j++) res <<= 1;
+	return res;
+}
+
+//******************************
+// the Wavelet algorithm itself
+//******************************
+
+int dywapitch_neededsamplecount(int minFreq) {
+	int nbSam = 3*44100/minFreq; // 1017. for 130 Hz
+	nbSam = _ceil_power2(nbSam); // 1024
+	return nbSam;
+}
+
+typedef struct _minmax {
+	int index;
+	struct _minmax *next;
+} minmax;
+
+double _dywapitch_computeWaveletPitch(double * samples, int startsample, int samplecount) {
+	double pitchF = 0.0;
+	
+	int i, j;
+	double si, si1;
+	
+	// must be a power of 2
+	samplecount = _floor_power2(samplecount);
+	
+	double *sam = (double *)malloc(sizeof(double)*samplecount);
+	memcpy(sam, samples + startsample, sizeof(double)*samplecount);
+	int curSamNb = samplecount;
+	
+	int *distances = (int *)malloc(sizeof(int)*samplecount);
+	int *mins = (int *)malloc(sizeof(int)*samplecount);
+	int *maxs = (int *)malloc(sizeof(int)*samplecount);
+	int nbMins, nbMaxs;
+	
+	// algorithm parameters
+	int maxFLWTlevels = 6;
+	double maxF = 3000.;
+	int differenceLevelsN = 3;
+	double maximaThresholdRatio = 0.75;
+	
+	double ampltitudeThreshold;  
+	double theDC = 0.0;
+	
+	{ // compute ampltitudeThreshold and theDC
+		//first compute the DC and maxAMplitude
+		double maxValue = 0.0;
+		double minValue = 0.0;
+		for (i = 0; i < samplecount;i++) {
+			si = sam[i];
+			theDC = theDC + si;
+			if (si > maxValue) maxValue = si;
+			if (si < minValue) minValue = si;
+		}
+		theDC = theDC/samplecount;
+		maxValue = maxValue - theDC;
+		minValue = minValue - theDC;
+		double amplitudeMax = (maxValue > -minValue ? maxValue : -minValue);
+		
+		ampltitudeThreshold = amplitudeMax*maximaThresholdRatio;
+		//asLog("dywapitch theDC=%f ampltitudeThreshold=%f\n", theDC, ampltitudeThreshold);
+		
+	}
+	
+	// levels, start without downsampling..
+	int curLevel = 0;
+	double curModeDistance = -1.;
+	int delta;
+	
+	while(1) {
+		
+		// delta
+		delta = static_cast<int>(44100./(_2power(curLevel)*maxF));
+		//("dywapitch doing level=%ld delta=%ld\n", curLevel, delta);
+		
+		if (curSamNb < 2) goto cleanup;
+		
+		// compute the first maximums and minumums after zero-crossing
+		// store if greater than the min threshold
+		// and if at a greater distance than delta
+		double dv, previousDV = -1000;
+		nbMins = nbMaxs = 0;   
+		int lastMinIndex = -1000000;
+		int lastmaxIndex = -1000000;
+		int findMax = 0;
+		int findMin = 0;
+		for (i = 2; i < curSamNb; i++) {
+			si = sam[i] - theDC;
+			si1 = sam[i-1] - theDC;
+			
+			if (si1 <= 0 && si > 0) findMax = 1;
+			if (si1 >= 0 && si < 0) findMin = 1;
+			
+			// min or max ?
+			dv = si - si1;
+			
+			if (previousDV > -1000) {
+				
+				if (findMin && previousDV < 0 && dv >= 0) { 
+					// minimum
+					if (fabs(si) >= ampltitudeThreshold) {
+						if (i > lastMinIndex + delta) {
+							mins[nbMins++] = i;
+							lastMinIndex = i;
+							findMin = 0;
+							//if DEBUGG then put "min ok"&&si
+							//
+						} else {
+							//if DEBUGG then put "min too close to previous"&&(i - lastMinIndex)
+							//
+						}
+					} else {
+						// if DEBUGG then put "min "&abs(si)&" < thresh = "&ampltitudeThreshold
+						//--
+					}
+				}
+				
+				if (findMax && previousDV > 0 && dv <= 0) {
+					// maximum
+					if (fabs(si) >= ampltitudeThreshold) {
+						if (i > lastmaxIndex + delta) {
+							maxs[nbMaxs++] = i;
+							lastmaxIndex = i;
+							findMax = 0;
+						} else {
+							//if DEBUGG then put "max too close to previous"&&(i - lastmaxIndex)
+							//--
+						}
+					} else {
+						//if DEBUGG then put "max "&abs(si)&" < thresh = "&ampltitudeThreshold
+						//--
+					}
+				}
+			}
+			
+			previousDV = dv;
+		}
+		
+		if (nbMins == 0 && nbMaxs == 0) {
+			// no best distance !
+			//asLog("dywapitch no mins nor maxs, exiting\n");
+			
+			// if DEBUGG then put "no mins nor maxs, exiting"
+			goto cleanup;
+		}
+		//if DEBUGG then put count(maxs)&&"maxs &"&&count(mins)&&"mins"
+		
+		// maxs = [5, 20, 100,...]
+		// compute distances
+		int d;
+		memset(distances, 0, samplecount*sizeof(int));
+		for (i = 0 ; i < nbMins ; i++) {
+			for (j = 1; j < differenceLevelsN; j++) {
+				if (i+j < nbMins) {
+					d = _iabs(mins[i] - mins[i+j]);
+					//asLog("dywapitch i=%ld j=%ld d=%ld\n", i, j, d);
+					distances[d] = distances[d] + 1;
+				}
+			}
+		}
+		for (i = 0 ; i < nbMaxs ; i++) {
+			for (j = 1; j < differenceLevelsN; j++) {
+				if (i+j < nbMaxs) {
+					d = _iabs(maxs[i] - maxs[i+j]);
+					//asLog("dywapitch i=%ld j=%ld d=%ld\n", i, j, d);
+					distances[d] = distances[d] + 1;
+				}
+			}
+		}
+		
+		// find best summed distance
+		int bestDistance = -1;
+		int bestValue = -1;
+		for (i = 0; i< curSamNb; i++) {
+			int summed = 0;
+			for (j = -delta ; j <= delta ; j++) {
+				if (i+j >=0 && i+j < curSamNb)
+					summed += distances[i+j];
+			}
+			//asLog("dywapitch i=%ld summed=%ld bestDistance=%ld\n", i, summed, bestDistance);
+			if (summed == bestValue) {
+				if (i == 2*bestDistance)
+					bestDistance = i;
+				
+			} else if (summed > bestValue) {
+				bestValue = summed;
+				bestDistance = i;
+			}
+		}
+		//asLog("dywapitch bestDistance=%ld\n", bestDistance);
+		
+		// averaging
+		double distAvg = 0.0;
+		double nbDists = 0;
+		for (j = -delta ; j <= delta ; j++) {
+			if (bestDistance+j >=0 && bestDistance+j < samplecount) {
+				int nbDist = distances[bestDistance+j];
+				if (nbDist > 0) {
+					nbDists += nbDist;
+					distAvg += (bestDistance+j)*nbDist;
+				}
+			}
+		}
+		// this is our mode distance !
+		distAvg /= nbDists;
+		//asLog("dywapitch distAvg=%f\n", distAvg);
+		
+		// continue the levels ?
+		if (curModeDistance > -1.) {
+			double similarity = fabs(distAvg*2 - curModeDistance);
+			if (similarity <= 2*delta) {
+				//if DEBUGG then put "similarity="&similarity&&"delta="&delta&&"ok"
+ 				//asLog("dywapitch similarity=%f OK !\n", similarity);
+				// two consecutive similar mode distances : ok !
+				pitchF = 44100./(_2power(curLevel-1)*curModeDistance);
+				goto cleanup;
+			}
+			//if DEBUGG then put "similarity="&similarity&&"delta="&delta&&"not"
+		}
+		
+		// not similar, continue next level
+		curModeDistance = distAvg;
+		
+		curLevel = curLevel + 1;
+		if (curLevel >= maxFLWTlevels) {
+			// put "max levels reached, exiting"
+ 			//asLog("dywapitch max levels reached, exiting\n");
+			goto cleanup;
+		}
+		
+		// downsample
+		if (curSamNb < 2) {
+ 			//asLog("dywapitch not enough samples, exiting\n");
+			goto cleanup;
+		}
+		for (i = 0; i < curSamNb/2; i++) {
+			sam[i] = (sam[2*i] + sam[2*i + 1])/2.;
+		}
+		curSamNb /= 2;
+	}
+	
+	///
+cleanup:
+	free(distances);
+	free(mins);
+	free(maxs);
+	free(sam);
+	
+	return pitchF;
+}
+
+// ***********************************
+// the dynamic postprocess
+// ***********************************
+
+/***
+It states: 
+ - a pitch cannot change much all of a sudden (20%) (impossible humanly,
+ so if such a situation happens, consider that it is a mistake and drop it. 
+ - a pitch cannot double or be divided by 2 all of a sudden : it is an
+ algorithm side-effect : divide it or double it by 2. 
+ - a lonely voiced pitch cannot happen, nor can a sudden drop in the middle
+ of a voiced segment. Smooth the plot. 
+***/
+
+double _dywapitch_dynamicprocess(dywapitchtracker *pitchtracker, double pitch) {
+	
+	// equivalence
+	if (pitch == 0.0) pitch = -1.0;
+	
+	//
+	double estimatedPitch = -1;
+	double acceptedError = 0.2f;
+	int maxConfidence = 5;
+	
+	if (pitch != -1) {
+		// I have a pitch here
+		
+		if (pitchtracker->_prevPitch == -1) {
+			// no previous
+			estimatedPitch = pitch;
+			pitchtracker->_prevPitch = pitch;
+			pitchtracker->_pitchConfidence = 1;
+			
+		} else if (abs(pitchtracker->_prevPitch - pitch)/pitch < acceptedError) {
+			// similar : remember and increment pitch
+			pitchtracker->_prevPitch = pitch;
+			estimatedPitch = pitch;
+			pitchtracker->_pitchConfidence = min(maxConfidence, pitchtracker->_pitchConfidence + 1); // maximum 3
+			
+		} else if ((pitchtracker->_pitchConfidence >= maxConfidence-2) && abs(pitchtracker->_prevPitch - 2.*pitch)/(2.*pitch) < acceptedError) {
+			// close to half the last pitch, which is trusted
+			estimatedPitch = 2.*pitch;
+			pitchtracker->_prevPitch = estimatedPitch;
+			
+		} else if ((pitchtracker->_pitchConfidence >= maxConfidence-2) && abs(pitchtracker->_prevPitch - 0.5*pitch)/(0.5*pitch) < acceptedError) {
+			// close to twice the last pitch, which is trusted
+			estimatedPitch = 0.5*pitch;
+			pitchtracker->_prevPitch = estimatedPitch;
+			
+		} else {
+			// nothing like this : very different value
+			if (pitchtracker->_pitchConfidence >= 1) {
+				// previous trusted : keep previous
+				estimatedPitch = pitchtracker->_prevPitch;
+				pitchtracker->_pitchConfidence = max(0, pitchtracker->_pitchConfidence - 1);
+			} else {
+				// previous not trusted : take current
+				estimatedPitch = pitch;
+				pitchtracker->_prevPitch = pitch;
+				pitchtracker->_pitchConfidence = 1;
+			}
+		}
+		
+	} else {
+		// no pitch now
+		if (pitchtracker->_prevPitch != -1) {
+			// was pitch before
+			if (pitchtracker->_pitchConfidence >= 1) {
+				// continue previous
+				estimatedPitch = pitchtracker->_prevPitch;
+				pitchtracker->_pitchConfidence = max(0, pitchtracker->_pitchConfidence - 1);
+			} else {
+				pitchtracker->_prevPitch = -1;
+				estimatedPitch = -1.;
+				pitchtracker->_pitchConfidence = 0;
+			}
+		}
+	}
+	
+	// put "_pitchConfidence="&pitchtracker->_pitchConfidence
+	if (pitchtracker->_pitchConfidence >= 1) {
+		// ok
+		pitch = estimatedPitch;
+	} else {
+		pitch = -1;
+	}
+	
+	// equivalence
+	if (pitch == -1) pitch = 0.0;
+	
+	return pitch;
+}
+
+
+// ************************************
+// the API main entry points
+// ************************************
+
+void dywapitch_inittracking(dywapitchtracker *pitchtracker) {
+	pitchtracker->_prevPitch = -1.;
+	pitchtracker->_pitchConfidence = -1;
+}
+
+double dywapitch_computepitch(dywapitchtracker *pitchtracker, double * samples, int startsample, int samplecount) {
+	double raw_pitch = _dywapitch_computeWaveletPitch(samples, startsample, samplecount);
+	return _dywapitch_dynamicprocess(pitchtracker, raw_pitch);
+}
+
+
+

PitchTracker/dywapitchtrack/dywapitchtrack.h

@@ -0,0 +1,110 @@
+/* dywapitchtrack.h
+ 
+ Dynamic Wavelet Algorithm Pitch Tracking library
+ Released under the MIT open source licence
+ 
+ Copyright (c) 2010 Antoine Schmitt
+  
+ Permission is hereby granted, free of charge, to any person obtaining a copy
+ of this software and associated documentation files (the "Software"), to deal
+ in the Software without restriction, including without limitation the rights
+ to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+ copies of the Software, and to permit persons to whom the Software is
+ furnished to do so, subject to the following conditions:
+ 
+ The above copyright notice and this permission notice shall be included in
+ all copies or substantial portions of the Software.
+ 
+ THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+ IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+ AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+ THE SOFTWARE.
+*/
+
+/* Documentation
+ 
+ The dywapitchtrack library computes the pitch of an audio stream in real time.
+ 
+ The pitch is the main frequency of the waveform (the 'note' being played or sung).
+ It is expressed as a float in Hz.
+ 
+ Unlike the human ear, pitch detection is difficult to achieve for computers. Many
+ algorithms have been designed and experimented, but there is no 'best' algorithm.
+ They all depend on the context and the tradeoffs acceptable in terms of speed and
+ latency. The context includes the quality and 'cleanness' of the audio : obviously
+ polyphonic sounds (multiple instruments playing different notes at the same time)
+ are extremely difficult to track, percussive or noisy audio has no pitch, most
+ real-life audio have some noisy moments, some instruments have a lot of harmonics,
+ etc...
+ 
+ The dywapitchtrack is based on a custom-tailored algorithm which is of very high quality:
+ both very accurate (precision < 0.05 semitones), very low latency (< 23 ms) and
+ very low error rate. It has been  thoroughly tested on human voice.
+ 
+ It can best be described as a dynamic wavelet algorithm (dywa):
+ 
+ The heart of the algorithm is a very powerful wavelet algorithm, described in a paper
+ by Eric Larson and Ross Maddox : "Real-Time Time-Domain Pitch Tracking Using Wavelets"
+ http://online.physics.uiuc.edu/courses/phys498pom/NSF_REU_Reports/2005_reu/Real-Time_Time-Domain_Pitch_Tracking_Using_Wavelets.pdf
+ 
+ This central algorithm has been improved by adding dynamic tracking, to reduce the
+ common problems of frequency-halving and voiced/unvoiced errors. This dynamic tracking
+ explains the need for a tracking structure (dywapitchtracker). The dynamic tracking assumes
+ that the main function dywapitch_computepitch is called repeatedly, as it follows the pitch
+ over time and makes assumptions about human voice capabilities and reallife conditions
+ (as documented inside the code).
+ 
+ Note : The algorithm currently assumes a 44100Hz audio sampling rate. If you use a different
+ samplerate, you can just multiply the resulting pitch by the ratio between your samplerate and 44100.
+*/
+
+/* Usage
+
+ // Allocate your audio buffers and start the audio stream.
+ // Allocate a 'dywapitchtracker' structure.
+ // Start the pitch tracking by calling 'dywapitch_inittracking'.
+ dywapitchtracker pitchtracker;
+ dywapitch_inittracking(&pitchtracker);
+ 
+ // For each available audio buffer, call 'dywapitch_computepitch'
+ double thepitch = dywapitch_computepitch(&pitchtracker, samples, start, count);
+ 
+*/
+
+#ifndef dywapitchtrack__H
+#define dywapitchtrack__H
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+// structure to hold tracking data
+typedef struct _dywapitchtracker {
+	double	_prevPitch;
+	int		_pitchConfidence;
+} dywapitchtracker;
+
+// returns the number of samples needed to compute pitch for fequencies equal and above the given minFreq (in Hz)
+// useful to allocate large enough audio buffer 
+// ex : for frequencies above 130Hz, you need 1024 samples (assuming a 44100 Hz samplerate)
+int dywapitch_neededsamplecount(int minFreq);
+
+// call before computing any pitch, passing an allocated dywapitchtracker structure
+void dywapitch_inittracking(dywapitchtracker *pitchtracker);
+
+// computes the pitch. Pass the inited dywapitchtracker structure
+// samples : a pointer to the sample buffer
+// startsample : the index of the first sample to use in the sample buffer
+// samplecount : the number of samples to use to compte the pitch
+// return 0.0 if no pitch was found (sound too low, noise, etc..)
+double dywapitch_computepitch(dywapitchtracker *pitchtracker, double * samples, int startsample, int samplecount);
+
+#ifdef __cplusplus
+} // extern "C"
+#endif
+
+#endif
+

PitchTracker/dywapitchtrack/ptDyWa.cpp

@@ -0,0 +1,27 @@
+#include "ptDyWa.h"
+#include "../Helper.h"
+
+
+PtDyWa::PtDyWa(unsigned step){
+	_Step = step;
+	dywapitch_inittracking(&_State);
+}
+
+void PtDyWa::Process(){
+	double AnaylsisBuf[_SampleCt];
+	size_t size = _AnalysisBuf.size();
+	if(size > _SampleCt)
+		_AnalysisBuf.pop(size - _SampleCt);
+	if(!_AnalysisBuf.read(AnaylsisBuf, AnaylsisBuf + _SampleCt))
+		return;
+	_AnalysisBuf.pop(_Step);
+	double pitch = dywapitch_computepitch(&_State, AnaylsisBuf,0,_SampleCt);
+	if(pitch == 0.0)
+		_LastNote = -1;
+	else
+		_LastNote = FreqToNote(pitch);
+}
+
+double PtDyWa::GetNote(){
+	return _LastNote;
+}

PitchTracker/dywapitchtrack/ptDyWa.h

@@ -0,0 +1,22 @@
+#pragma once
+#include "../performous/pitch.hh"
+#include "dywapitchtrack.h"
+
+class PtDyWa{
+public:
+	PtDyWa(unsigned step);
+
+	/** Add input data to buffer. This is thread-safe (against other functions). **/
+	template <typename InIt> void input(InIt begin, InIt end) {
+		_AnalysisBuf.insert(begin, end);
+	}
+
+	void Process();
+	double GetNote();
+private:
+	dywapitchtracker _State;
+	constexpr static size_t _SampleCt = 2048;
+	RingBuffer<_SampleCt * 2> _AnalysisBuf;
+	unsigned _Step;
+	double _LastNote;
+};

PitchTracker/performous/pitch.cc

@@ -139,7 +139,7 @@ void Analyzer::calcTones() {
 		double delta = phase - m_fftLastPhase[k];
 		m_fftLastPhase[k] = static_cast<float>(phase);
 		delta -= k * phaseStep;  // subtract expected phase difference
-		delta = remainder(delta, 2.0 * M_PI);  // map delta phase into +/- M_PI interval
+		delta = remainder(delta, M_PI);  // map delta phase into +/- M_PI interval
 		delta /= phaseStep;  // calculate diff from bin center frequency
 		double freq = (k + delta) * freqPerBin;  // calculate the true frequency
 		if (freq > 1.0 && magnitude > minMagnitude) {
@@ -164,15 +164,16 @@ void Analyzer::calcTones() {
 		int bestScore = 0;
 		for (std::size_t div = 2; div <= Tone::MAXHARM && k / div > 1; ++div) {
 			double freq = peaks[k].freq / div; // Fundamental
+			double db = peaks[k].db;
 			int score = 0;
 			for (std::size_t n = 1; n < div && n < 8; ++n) {
 				Peak& p = match(peaks, k * n / div);
 				--score;
-				if (p.db < -90.0 || std::abs(p.freq / n / freq - 1.0) > .03) continue;
-				if (n == 1) score += 4; // Extra for fundamental
+				if (std::abs(p.db/db - 1.0) > 0.4 || std::abs(p.freq / n / freq - 1.0) > .04) continue;
+				if (n == 1) score += 1; // Extra for fundamental
 				score += 2;
 			}
-			if (score > bestScore) {
+			if (score >= bestScore) {
 				bestScore = score;
 				bestDiv = div;
 			}
@@ -185,12 +186,11 @@ void Analyzer::calcTones() {
 		for (std::size_t n = 1; n <= bestDiv; ++n) {
 			// Find the peak for n'th harmonic
 			Peak& p = match(peaks, k * n / bestDiv);
-			if (std::abs(p.freq / n / freq - 1.0) > .03) continue; // Does it match the fundamental freq?
-			if (p.db > t.db - 10.0) {
-				t.db = std::max(t.db, p.db);
-				++count;
-				t.freq += p.freq / n;
-			}
+			if (std::abs(p.freq / n / freq - 1.0) > .04) continue; // Does it match the fundamental freq?
+			if (std::abs(p.db/t.db - 1.0) > 0.4) continue;
+			t.db = std::max(t.db, p.db);
+			++count;
+			t.freq += p.freq / n;
 			t.harmonics[n - 1] = p.db;
 			p.clear();
 		}

PitchTracker/ptAKF.cpp

@@ -13,23 +13,37 @@
 // You should have received a copy of the GNU General Public License
 // along with Vocaluxe. If not, see <http://www.gnu.org/licenses/>.
 
+#define _USE_MATH_DEFINES
 #include "ptAKF.h"
-#include <math.h>
+#include <cmath>
 
 static constexpr double BaseToneFrequency = 65.4064; // lowest (half-)tone to analyze (C2 = 65.4064 Hz)
 static constexpr double HalftoneBase = 1.05946309436; // 2^(1/12) -> HalftoneBase^12 = 2 (one octave)
 int PtAKF::_InitCount = 0;
-float* PtAKF::_SamplesPerPeriodPerTone = NULL;
+float* restrict PtAKF::_SamplesPerPeriodPerTone = NULL;
+float* restrict PtAKF::_SamplesPerPeriodPerToneFine = NULL;
+float* restrict PtAKF::_Window = NULL;
+
+double GetFrequencyFromTone(double toneIndex){
+	return BaseToneFrequency * pow(HalftoneBase, toneIndex);
+}
 
 PtAKF::PtAKF(unsigned step){
 	_InitCount++;
 	if(_InitCount == 1){
 		_SamplesPerPeriodPerTone = new float[_MaxHalfTone + 1 + _HalfTonesAdd];
+		_SamplesPerPeriodPerToneFine = new float[(_MaxHalfTone + 1 + _HalfTonesAdd) * 2];
 		//Init Array to avoid costly calculations
 		for (int toneIndex = 0; toneIndex <= _MaxHalfTone + _HalfTonesAdd; toneIndex++)
 		{
-			double freq = BaseToneFrequency * pow(HalftoneBase, toneIndex);
-			_SamplesPerPeriodPerTone[toneIndex] = static_cast<float>(44100.0 / freq); // samples in one period
+			_SamplesPerPeriodPerTone[toneIndex] = static_cast<float>(44100.0 / GetFrequencyFromTone(toneIndex)); // samples in one period
+			_SamplesPerPeriodPerToneFine[2 * toneIndex] = static_cast<float>(44100.0 / GetFrequencyFromTone(toneIndex - 1./3.)); // samples for a bit below exact frequency
+			_SamplesPerPeriodPerToneFine[2 * toneIndex + 1] = static_cast<float>(44100.0 / GetFrequencyFromTone(toneIndex + 1./3.)); // samples for a bit above exact frequency
+		}
+		_Window = new float[_SampleCt];
+		double cosMult = 2. * M_PI / (_SampleCt - 1.); //To simplify and speed up
+		for (int i = 0; i < _SampleCt; i++) {
+			_Window[i] = static_cast<float>(0.54 - 0.46 * cos(i * cosMult));
 		}
 	}
 	_Step = step;
@@ -40,10 +54,12 @@ PtAKF::~PtAKF(){
 	if(!_InitCount){
 		delete[] _SamplesPerPeriodPerTone;
 		_SamplesPerPeriodPerTone = NULL;
+		delete[] _SamplesPerPeriodPerToneFine;
+		_SamplesPerPeriodPerToneFine = NULL;
 	}
 }
 
-int PtAKF::GetNote(double* maxVolume, float* weights){
+int PtAKF::GetNote(double* restrict maxVolume, float* restrict weights){
 	float AnaylsisBuf[_SampleCt];
 	size_t size = _AnalysisBuf.size();
 	if(size > _SampleCt)
@@ -54,9 +70,37 @@ int PtAKF::GetNote(double* maxVolume, float* weights){
 	return _GetNote(AnaylsisBuf, maxVolume, weights);
 }
 
-int PtAKF::_GetNote(float* samples, double* maxVolume, float* weights){
+void PtAKF::InitPeaks(SPeak peaks[_MaxPeaks]){
+	for(int i = 0; i < _MaxPeaks; i++){
+		peaks[i].toneIndex = -1;
+		peaks[i].weight = -9.f;
+	}
+}
+
+void PtAKF::AddPeak(SPeak peaks[], float curWeight, int toneIndex){
+	if(curWeight > peaks[0].weight){
+		int i;
+		for(i = 1; i < _MaxPeaks; i++){
+			if(curWeight < peaks[i].weight){
+				peaks[i-1].toneIndex = toneIndex;
+				peaks[i-1].weight = curWeight;
+				break;
+			}else{
+				peaks[i-1] = peaks[i];
+			}
+		}
+		if(i == _MaxPeaks){
+			peaks[i-1].toneIndex = toneIndex;
+			peaks[i-1].weight = curWeight;
+		}
+	}
+}
+
+int PtAKF::_GetNote(float* restrict samples, double* restrict maxVolume, float* restrict weights){
+	// Calculate maximum volume
+
 	float maxVolumeL = 0;
-	for(int i = 0; i < _SampleCt; i++){
+	for(int i = _SampleCt / 2; i < _SampleCt; i++){
 		float vol = abs(samples[i]);
 		if(vol > maxVolumeL)
 			maxVolumeL = vol;
@@ -64,15 +108,14 @@ int PtAKF::_GetNote(float* samples, double* maxVolume, float* weights){
 	//maxVolumeL /= 32767.;
 	*maxVolume = maxVolumeL;
 
-	//Scale the samples to have a peak of 1 to better decide between voiced/unvoiced
-	//float scaleVal = 1 / maxVolumeL;
-	//for(int i = 0; i < _SampleCt; i++){
-	//	samples[i] *= scaleVal;
-	//}
+	for(int i = 0; i < _SampleCt; i++){
+		samples[i] *= _Window[i];
+	}
 
-	float maxWeight = -1;
-	float minWeight = 1;
-	int maxTone = -1;
+	// Now analyze the samples and get peaks at the most appropriate tones
+
+	SPeak peaks[_MaxPeaks];
+	InitPeaks(peaks);
 
 	//Attention: We have a peak at lag 0 that might stretch that far, that we detect a wrong "peak" at _MaxHalfTone
 	//Because of that we filter out all tones that are past the last zero crossing from below but keep tones with decreasing weights (going towards zero crossing from above)
@@ -80,16 +123,9 @@ int PtAKF::_GetNote(float* samples, double* maxVolume, float* weights){
 	float lastWeight = 1.f;
 	for (int toneIndex = _MinHalfTone; toneIndex <= _MaxHalfTone; toneIndex++)
 	{
-		float curWeight = _AnalyzeToneFunc(samples, toneIndex);
+		float curWeight = _AnalyzeByTone(samples, toneIndex);
 
-		if (curWeight > maxWeight)
-		{
-			maxWeight = curWeight;
-			maxTone = toneIndex;
-		}
-
-		if (curWeight < minWeight)
-			minWeight = curWeight;
+		AddPeak(peaks, curWeight, toneIndex);
 
 		weights[toneIndex - _MinHalfTone] = curWeight;
 
@@ -99,54 +135,96 @@ int PtAKF::_GetNote(float* samples, double* maxVolume, float* weights){
 		lastWeight = curWeight;
 	}
 
-	if(maxTone==_MaxHalfTone){
+	// Now clear off the lag 0 peak if required
+
+	if(peaks[_MaxPeaks - 1].toneIndex ==_MaxHalfTone){
 		//We might have caught the lag 0 peak so go a bit further to check for other zero crossings (or we won't be able to detect _MaxHalfTone)
-		float lastWeight =  _FastAKF(samples, _SamplesPerPeriodPerTone[_MaxHalfTone]);
+		float lastWeight =  _AKFByTone(samples, _MaxHalfTone);
 		for(int toneIndex = _MaxHalfTone+1; toneIndex <= _MaxHalfTone + _HalfTonesAdd; toneIndex++){
-			float curWeight = _FastAKF(samples, _SamplesPerPeriodPerTone[toneIndex]);
+			float curWeight = _AKFByTone(samples, toneIndex);
 			if(lastWeight > curWeight || (lastWeight > 0.f && curWeight <= 0.f)){
 				lastValidTone = toneIndex - 1;
 				break;
 			}
 			lastWeight = curWeight;
 		}
-		if(lastValidTone < maxTone){
-			maxWeight = -1.f;
-			minWeight = 1.f;
+		if(lastValidTone < peaks[_MaxPeaks - 1].toneIndex){
+			InitPeaks(peaks);
 			for(int toneIndex = _MinHalfTone; toneIndex <= lastValidTone; toneIndex++){
-				float curWeight = weights[toneIndex-_MinHalfTone];
-				if (curWeight > maxWeight)
-				{
-					maxWeight = curWeight;
-					maxTone = toneIndex;
-				}
-
-				if (curWeight < minWeight)
-					minWeight = curWeight;
+				AddPeak(peaks, weights[toneIndex - _MinHalfTone], toneIndex);
 			}
 			//Set all invalid weights to 0
 			for(int toneIndex = lastValidTone + 1;toneIndex < _MaxHalfTone; toneIndex++){
-				weights[toneIndex - _MinHalfTone]=0.f;
+				weights[toneIndex - _MinHalfTone] = 0.f;
 			}
 		}
 	}
 
-	float energy = _FastAKF(samples, 0);
-	float maxAKF = _FastAKF(samples, _SamplesPerPeriodPerTone[maxTone]);
-	//printf("%g %g\n",energy, maxAKF);
+	// Further analyze the found peaks and find the real maximum
+	// This is necessary because we may have our real peak a bit off the exact tone frequency
+	// and a 'wrong' peak that is exactly at another tone which might become higher than the one at the 'right' tone
+
+	int maxTone = peaks[_MaxPeaks - 1].toneIndex;
+	int maxToneFine = -1;
+	float maxWeight = peaks[_MaxPeaks - 1].weight;
+	for(int i = 0; i < _MaxPeaks; i++){
+		int toneIndex = peaks[i].toneIndex;
+		float curWeightCenter = peaks[i].weight;
+		float curWeightDown = _AnalyzeBySampleCt(samples, _SamplesPerPeriodPerToneFine[toneIndex * 2]);
+		float curWeight;
+		int otherToneIndex;
+		int otherToneFineIndex;
+		if(curWeightDown > curWeightCenter){
+			otherToneIndex = toneIndex - 1;
+			otherToneFineIndex = 1;
+			curWeight = curWeightDown;
+		}else{
+			float curWeightUp = _AnalyzeBySampleCt(samples, _SamplesPerPeriodPerToneFine[toneIndex * 2 + 1]);
+			if(curWeightUp > curWeightCenter){
+				otherToneIndex = (toneIndex < _MaxHalfTone)? toneIndex + 1: -1; // If the other tone is invalid just set it to -1 which will be skipped below
+				otherToneFineIndex = 0;
+				curWeight = curWeightUp;
+			}else
+				continue;
+		}
+		weights[toneIndex - _MinHalfTone] = curWeight;
+		if(curWeight > maxWeight){
+			maxTone = toneIndex;
+			maxToneFine = 1 - otherToneFineIndex;
+			maxWeight = curWeight;
+		}
+
+		// Now check also neighbouring tone
+		if(otherToneIndex > 0){
+			float curWeightOther = _AnalyzeBySampleCt(samples, _SamplesPerPeriodPerToneFine[otherToneIndex * 2 + otherToneFineIndex]);
+			if(curWeightOther > weights[otherToneIndex - _MinHalfTone]){
+				weights[otherToneIndex - _MinHalfTone] = curWeightOther;
+				if(curWeightOther > maxWeight){
+					maxTone = otherToneIndex;
+					maxToneFine = otherToneFineIndex;
+					maxWeight = curWeightOther;
+				}
+			}
+		}
+	}
+
+	float energy = _AKFBySampleCt(samples, 0);
+	float maxAKF = (maxToneFine < 0) ? _AKFByTone(samples, maxTone):_AKFBySampleCt(samples, _SamplesPerPeriodPerToneFine[maxTone * 2 + maxToneFine]);
+	float aa = _AKFByTone(samples, maxTone);
 
 	//if(maxWeight - minWeight > 0.025){
-	if(maxAKF >= 0.3 * energy){
+	if(maxAKF >= 0.4 * energy){
 		return maxTone;
 	}else return -1;
 }
 
-float PtAKF::_AnalyzeToneFunc(float* samples, int toneIndex){
-	// Use method by Kobayashi and Shimamura (2001): Combine AKF and AMDF to a new f(z)=AKF(z)/(AMDF(z)+k) with k=1
+float PtAKF::_AnalyzeByTone(float* restrict samples, int toneIndex){
+	return _AnalyzeBySampleCt(samples, _SamplesPerPeriodPerTone[toneIndex]);
+}
 
-	float samplesPerPeriodD = _SamplesPerPeriodPerTone[toneIndex]; // samples in one period
+float PtAKF::_AnalyzeBySampleCt(float* restrict samples, float samplesPerPeriodD){
+	// Use method by Kobayashi and Shimamura (2001): Combine AKF and AMDF to a new f(z)=AKF(z)/(AMDF(z)+k) with k=1
 	int samplesPerPeriod = static_cast<int>(samplesPerPeriodD);
-	int maxSampleCt = static_cast<int>(_SamplesPerPeriodPerTone[0]);
 	float fHigh = samplesPerPeriodD - samplesPerPeriod;
 	float fLow = 1.0f - fHigh;
 
@@ -174,11 +252,13 @@ float PtAKF::_AnalyzeToneFunc(float* samples, int toneIndex){
 	//{toneIndex}: {accumDistAKF} ; {accumDistAMDF}; {result}
 }
 
-float PtAKF::_FastAKF(float* samples, float samplesPerPeriodD){
-	// Do a fast, non-scaled AKF (used for Lag 0 peak removal)
+float PtAKF::_AKFByTone(float* restrict samples, int toneIndex){
+	return _AKFBySampleCt(samples, _SamplesPerPeriodPerTone[toneIndex]);
+}
 
+float PtAKF::_AKFBySampleCt(float* restrict samples, float samplesPerPeriodD){
+	// Do a fast, non-scaled AKF (used e.g. for Lag 0 peak removal)
 	int samplesPerPeriod = static_cast<int>(samplesPerPeriodD);
-	int maxSampleCt = static_cast<int>(_SamplesPerPeriodPerTone[0]);
 	float fHigh = samplesPerPeriodD - samplesPerPeriod;
 	float fLow = 1.0f - fHigh;
 

PitchTracker/ptAKF.h

@@ -1,6 +1,23 @@
 #pragma once
 #include "performous/pitch.hh"
 
+#if __STDC__ != 1
+#    define restrict __restrict
+#else
+#    ifndef __STDC_VERSION__
+#        define restrict __restrict
+#    else
+#        if __STDC_VERSION__ < 199901L
+#            define restrict __restrict
+#        endif
+#    endif
+#endif
+
+struct SPeak{
+	int toneIndex;
+	float weight;
+};
+
 // A pitch detection that is based on the AKF and AMDF (initially from original Vocaluxe/USDx)
 class PtAKF{
 public:
@@ -12,20 +29,31 @@ public:
 		_AnalysisBuf.insert(begin, end);
 	}
 
-	int GetNote(double* maxVolume, float* weights);
+	int GetNote(double* restrict maxVolume, float* restrict weights);
 	static int GetNumHalfTones(){ return _MaxHalfTone - _MinHalfTone + 1;}
+
 private:
 	static int _InitCount;
-	static float* _SamplesPerPeriodPerTone;
+	static float* restrict _SamplesPerPeriodPerTone;
+	// Use a 3 times finer resolution for exact peak detection
+	// So this array will store for each note the number of samples 1/3 below and up
+	// (cannot use halves as it would be ambiguous)
+	static float* restrict _SamplesPerPeriodPerToneFine;
+	static float* restrict _Window;
 	static constexpr int _MinHalfTone = 0;  //C2
 	static constexpr int _MaxHalfTone = 55;//47; //B5
 	static constexpr int _HalfTonesAdd = 4; //Additonal half tones to analyze to remove the peak at lag 0
+	static constexpr int _MaxPeaks = 3;
 
 	constexpr static size_t _SampleCt = 2048;
 	RingBuffer<_SampleCt * 2> _AnalysisBuf;
 	unsigned _Step;
 
-	static int _GetNote(float* samples, double* maxVolume, float* weights);
-	static float _AnalyzeToneFunc(float* samples, int toneIndex);
-	static float _FastAKF(float* samples, float samplesPerPeriodD);
+	static int _GetNote(float* restrict samples, double* restrict maxVolume, float* restrict weights);
+	static float _AnalyzeBySampleCt(float* restrict samples, float samplesPerPeriodD);
+	static inline float _AnalyzeByTone(float* restrict samples, int toneIndex);
+	static float _AKFBySampleCt(float* restrict samples, float samplesPerPeriodD);
+	static inline float _AKFByTone(float* restrict samples, int toneIndex);
+	static inline void InitPeaks(SPeak peaks[_MaxPeaks]);
+	static inline void AddPeak(SPeak peaks[_MaxPeaks], float curWeight, int toneIndex);
 };

Vocaluxe/Lib/Sound/Record/CAnalyzer.cs

@@ -159,4 +159,69 @@ namespace Vocaluxe.Lib.Sound.Record
             GC.SuppressFinalize(this);
         }
     }
+
+    /// <summary>
+    /// Pitchtracker (Pt) that uses a dynamic wavelet (DyWa) algorithm 
+    /// Quite fast and perfect in artificial tests, but may fail in real world scenarios (e.g. missing fundamental)
+    /// </summary>
+    class CPtDyWa : IDisposable
+    {
+        #region Imports
+        [DllImport("PitchTracker.dll", CallingConvention = CallingConvention.Cdecl)]
+        private static extern IntPtr PtDyWa_Create(uint step);
+
+        [DllImport("PitchTracker.dll", CallingConvention = CallingConvention.Cdecl)]
+        private static extern void PtDyWa_Free(IntPtr analyzer);
+
+        [DllImport("PitchTracker.dll", CallingConvention = CallingConvention.Cdecl)]
+        private static extern void PtDyWa_InputByte(IntPtr analyzer, [In] byte[] data, int sampleCt);
+
+        [DllImport("PitchTracker.dll", CallingConvention = CallingConvention.Cdecl)]
+        private static extern int PtDyWa_Process(IntPtr analyzer);
+
+        [DllImport("PitchTracker.dll", CallingConvention = CallingConvention.Cdecl)]
+        private static extern double PtDyWa_GetNote(IntPtr analyzer);
+        #endregion
+
+        private IntPtr _Instance;
+
+        public CPtDyWa(uint step = 200)
+        {
+            _Instance = PtDyWa_Create(step);
+        }
+
+        ~CPtDyWa()
+        {
+            _Dispose(false);
+        }
+
+        public void Input(byte[] data)
+        {
+            PtDyWa_InputByte(_Instance, data, data.Length / 2);
+        }
+
+        public void Process()
+        {
+            PtDyWa_Process(_Instance);
+        }
+
+        public double GetNote()
+        {
+            return PtDyWa_GetNote(_Instance);
+        }
+
+        private void _Dispose(bool disposing)
+        {
+            if (_Instance == IntPtr.Zero)
+                throw new ObjectDisposedException(GetType().Name);
+            PtDyWa_Free(_Instance);
+            _Instance = IntPtr.Zero;
+        }
+
+        public void Dispose()
+        {
+            _Dispose(true);
+            GC.SuppressFinalize(this);
+        }
+    }
 }

Vocaluxe/Lib/Sound/Record/CBuffer.cs

@@ -38,7 +38,6 @@ namespace Vocaluxe.Lib.Sound.Record
         private const double _HalftoneBase = 1.05946309436; // 2^(1/12) -> HalftoneBase^12 = 2 (one octave)
 #endif
         private CPtAKF _Analyzer = new CPtAKF();
-        private CAnalyzer _Analyzer2 = new CAnalyzer();
 
         private double _MaxVolume;
 
@@ -98,12 +97,10 @@ namespace Vocaluxe.Lib.Sound.Record
             //if (assigned(fVoiceStream)) then
             //fVoiceStream.WriteData(Buffer, BufferSize);
             _Analyzer.Input(buffer);
-            _Analyzer2.Input(buffer);
         }
 
         public void AnalyzeBuffer()
         {
-            _Analyzer2.Process();
             //_MaxVolume = _Analyzer.GetPeak() / 43 + 1;
             int tone = _Analyzer.GetNote(out _MaxVolume, ToneWeigths); //(int)Math.Round(_Analyzer.FindNote()); // 
             if (tone >= 0 && _MaxVolume >= MinVolume)
@@ -136,6 +133,8 @@ namespace Vocaluxe.Lib.Sound.Record
 
         private static string _ToneToNote(int tone, bool withOctave = true)
         {
+            if (tone < 0)
+                return "inv.";
             tone += 24;
             string[] notes = {"C", "C#", "D", "D#", "E", "F", "F#", "G", "G#", "A", "A#", "B"};
             string result = notes[tone % 12];
@@ -149,7 +148,7 @@ namespace Vocaluxe.Lib.Sound.Record
             if (_PitchTestRun)
                 return;
             _PitchTestRun = true;
-            int toneFrom = 0;
+            const int toneFrom = 0;
             const int toneTo = -1; // 47; //B5
             Console.WriteLine("Testing notes " + _ToneToNote(toneFrom) + " - " + _ToneToNote(toneTo));
             byte[] data;
@@ -185,7 +184,7 @@ namespace Vocaluxe.Lib.Sound.Record
                     tests++;
                 }
             }
-            if (_TestFile("toneG3.wav", 19))
+            /*if (_TestFile("toneG3.wav", 19))
                 ok++;
             tests++;
             if (_TestFile("toneG3Miss.wav", 19))
@@ -193,9 +192,9 @@ namespace Vocaluxe.Lib.Sound.Record
             tests++;
             if (_TestFile("toneG4.wav", 31))
                 ok++;
-            tests++;
-            _TestFile("whistling3.wav", "whistling3.txt", ref tests, ref ok);
+            tests++;*/
             _TestFile("sClausVoc.wav", "sClausVoc.txt", ref tests, ref ok);
+            _TestFile("whistling3.wav", "whistling3.txt", ref tests, ref ok);
 
             _GetSineWave(_BaseToneFreq * Math.Pow(_HalftoneBase, 5), 44100, ref angle, out data);
             Stopwatch sw = new Stopwatch();
@@ -271,35 +270,52 @@ namespace Vocaluxe.Lib.Sound.Record
                 int samplesRead = 0;
                 int curTimeIndex = -1;
                 int curNote = -1;
+                int lastNote = -1; //It is ok not to detect a note change immediately so the old note is kept for a bit as ok
+                int lastNoteEndTime = 0;
+                const int lastNoteMaxTimeDiff = 1024 * 1000 / 44100; // old note is valid for 1024 more samples
                 const int maxSamplesPerBatch = 512;
-                while (wavFile.NumSamplesLeft / wavFile.NumChannels > maxSamplesPerBatch)
+                CAnalyzer analyzer2 = new CAnalyzer();
+                CPtDyWa analyzer3 = new CPtDyWa();
+                while (wavFile.NumSamplesLeft > maxSamplesPerBatch)
                 {
                     byte[] samples = wavFile.GetNextSamples16BitAsBytes(maxSamplesPerBatch, 1);
                     samplesRead += samples.Length / 2;
                     int time = samplesRead * 1000 / wavFile.SampleRate;
                     ProcessNewBuffer(samples);
                     AnalyzeBuffer();
+                    analyzer2.Input(samples);
+                    analyzer3.Input(samples);
+                    analyzer2.Process();
+                    analyzer3.Process();
                     while (curTimeIndex + 1 < tones.Count && time >= tones[curTimeIndex + 1].Time)
                     {
                         curTimeIndex++;
+                        lastNote = curNote;
+                        lastNoteEndTime = tones[curTimeIndex].Time;
                         curNote = tones[curTimeIndex].Note;
                     }
                     if (curNote < 0)
                         continue;
-                    ct++;
+                    if (lastNoteEndTime + lastNoteMaxTimeDiff < time)
+                        lastNote = curNote;
+                    int tone1 = ToneValid ? ToneAbs : -1;
+                    int tone2 = (int)Math.Round(analyzer2.FindNote(64, 1770));
+                    int tone3 = (int)Math.Round(analyzer3.GetNote());
+                    /*ct++;
                     if (!ToneValid)
                         CBase.Log.LogDebug("Note " + _ToneToNote(curNote) + " at " + time + "ms not detected");
                     else if (ToneAbs % 12 != curNote % 12)
                         CBase.Log.LogDebug("Note " + _ToneToNote(curNote) + " at " + time + "ms wrongly detected as " + _ToneToNote(ToneAbs));
                     else
-                        ok++;
-                    int note2 = (int)Math.Round(_Analyzer2.FindNote());
-                    if (note2 < 0)
-                        CBase.Log.LogDebug("Note2 " + _ToneToNote(curNote) + " at " + time + "ms not detected");
-                    else if (note2 % 12 != curNote % 12)
-                        CBase.Log.LogDebug("Note2 " + _ToneToNote(curNote) + " at " + time + "ms wrongly detected as " + _ToneToNote(note2));
-                    else
-                        ok++;
+                        ok++;*/
+                    bool ok1 = tone1 == curNote || tone1 == lastNote;
+                    bool ok2 = tone2 == curNote || tone2 == lastNote;
+                    bool ok3 = tone3 == curNote || tone3 == lastNote;
+                    if (!ok1 || !ok2 || !ok3)
+                    {
+                        CBase.Log.LogDebug("Note " + _ToneToNote(curNote) + " at " + time + "ms detected as " + _ToneToNote(tone1) + (ok1 ? "" : "(!)") + "; " + _ToneToNote(tone2) +
+                                           (ok2 ? "" : "(!)") + "; " + _ToneToNote(tone3) + (ok3 ? "" : "(!)") + "; ");
+                    }
                 }
             }
             catch (Exception e)

VocaluxeLib/Utils/CWavFile.cs

@@ -14,11 +14,12 @@ namespace VocaluxeLib.Utils
         public string FileName { get; private set; }
         //Unless this is true, all fields are to be considered as invalid
         public int DataSize { get; private set; }
-        //Number of samples (=SampesPerChannel*NumChannels)
+        //Number of samples (per channel)
         public int NumSamples
         {
-            get { return (DataSize / (BitsPerSample / 8)); }
+            get { return (DataSize / (BitsPerSample / 8 * NumChannels)); }
         }
+        //Number of samples left to be read (per channel)
         public int NumSamplesLeft { get; private set; }
         public bool IsOpen { get; private set; }
         private int _DataPos;
@@ -129,13 +130,12 @@ namespace VocaluxeLib.Utils
         {
             if (channel < 0 || channel > NumChannels || BitsPerSample != 8)
                 return null;
-            int samplesPerBlock = (channel == 0 ? NumChannels : 1);
-            if (NumSamplesLeft < numSamples * samplesPerBlock)
+            if (numSamples > NumSamplesLeft)
                 return null;
             byte[] result;
             BinaryReader reader = new BinaryReader(_File);
             if (channel == 0 || NumChannels == 1)
-                result = reader.ReadBytes(numSamples * samplesPerBlock);
+                result = reader.ReadBytes(numSamples * NumChannels);
             else
             {
                 result = new byte[numSamples];
@@ -146,7 +146,7 @@ namespace VocaluxeLib.Utils
                     _File.Seek(NumChannels - channel, SeekOrigin.Current);
                 }
             }
-            NumSamplesLeft -= numSamples * samplesPerBlock;
+            NumSamplesLeft -= numSamples;
             return result;
         }
 
@@ -160,8 +160,7 @@ namespace VocaluxeLib.Utils
         {
             if (channel < 0 || channel > NumChannels || BitsPerSample != 16)
                 return null;
-            int samplesPerBlock = (channel == 0 ? NumChannels : 1);
-            if (NumSamplesLeft < numSamples * samplesPerBlock)
+            if (numSamples > NumSamplesLeft)
                 return null;
             short[] result;
             BinaryReader reader = new BinaryReader(_File);
@@ -183,7 +182,7 @@ namespace VocaluxeLib.Utils
                         _File.Seek((NumChannels - channel) * 2, SeekOrigin.Current);
                     }
                 }
-                NumSamplesLeft -= numSamples * samplesPerBlock;
+                NumSamplesLeft -= numSamples;
             }
             catch (EndOfStreamException e)
             {