Merge branch 'master' of http://192.168.0.5:3000/ftc/soundthermo-cpp

Measurer.cpp

@@ -67,27 +67,42 @@ float Measurer::Measure(float duration)
 float Measurer::ProcessData(QByteArray data)
 {
 	// convert byte data to sample array
-	vector<float> samples;
+	vector<float> read_samples;
 	// note that we fixed format to be little-endian, two-byte, signed
 	for (int i = 0; i < data.count(); i += 2)
 	{
 		int t = qFromLittleEndian<qint16>(data.data() + i);
-		samples.push_back(t / 32767.0);
+		read_samples.push_back(t / 32767.0);
 	}
 	if (debug)
 	{
 		ostringstream oss;
-		for (int i = 0; i < (int)samples.size(); i++)
+		for (int i = 0; i < (int)read_samples.size(); i++)
 		{
-			oss << samples[i] << " ";
+			oss << read_samples[i] << " ";
 		}
 		cerr << "Samples read: " << oss.str() << "\n";
 	}
+	// cut off some initial samples because of delays in sound subsystem 
+	// (otherwise we can hear some older sound which is not finished playing yet)
+	vector<float> samples;
+	for (int i = (int)(read_samples.size() * INITIAL_DELAY_PERCENTAGE); i < (int)read_samples.size(); i++)
+	{
+		samples.push_back(read_samples[i]);
+	}
+	// determine silence threshold (will consider signal less than SILENCE of maximum as silence)
+	float max_sample = abs(samples[0]);
+	for (int i = 1; i < (int)samples.size(); i++)
+	{
+		float q = abs(samples[i]);
+		if (q > max_sample) max_sample = q;
+	}
+	float silence_level = SILENCE_LEVEL * max_sample;
 	// cut off starting and ending silence
 	int left_pos = 0;
-	while (left_pos < (int)samples.size() && abs(samples[left_pos]) < 1e-4) left_pos++;
+	while (left_pos < (int)samples.size() && abs(samples[left_pos]) < silence_level) left_pos++;
 	int right_pos = (int)samples.size() - 1;
-	while (right_pos >= 0 && abs(samples[right_pos]) < 1e-4) right_pos--;
+	while (right_pos >= 0 && abs(samples[right_pos]) < silence_level) right_pos--;
 	int len = right_pos - left_pos + 1;
 	int pow = 1, m = 0;
 	while (2 * pow <= len)
@@ -101,24 +116,45 @@ float Measurer::ProcessData(QByteArray data)
 		d[i] = Complex(samples[left_pos + i], 0);
 	}
 	FFT(d, pow, false);
+	return ProcessSpectrum (pow, d);
+}
+
+float Measurer::ProcessSpectrum(int pow, Complex *d)
+{
+	if (debug)
+	{
+		ostringstream oss;
+		for (int i = 0; i < pow; i++)
+		{
+			float amp = (float)sqrt(d[i].real() * d[i].real() + d[i].imag() * d[i].imag());
+			oss << amp << " ";
+		}
+		cerr << "Spectrum: " << oss.str() << "\n";
+	}
+	// will search for maximum value near our frequency
+	// also will 
+	float desired = TONE_FREQUENCY * pow / SAMPLE_RATE;
+	int find_left = (int)((TONE_FREQUENCY - FREQUENCY_GAP) * pow / SAMPLE_RATE);
+	if (find_left < 0) find_left = 0;
+	int find_right = (int)((TONE_FREQUENCY + FREQUENCY_GAP) * pow / SAMPLE_RATE);
+	if (find_right > pow) find_right = pow;
 	float max = -1;
-	int argmax = 0;
-	for (int i = 0; i < pow; i++)
+	float avg = 0;
+	for (int i = find_left; i < find_right; i++)
 	{
 		float amp = (float)sqrt(d[i].real() * d[i].real() + d[i].imag() * d[i].imag());
 		if (amp > max)
 		{
 			max = amp;
-			argmax = i;
 		}
+		avg += amp;
 	}
-	float desired = TONE_FREQUENCY * pow / SAMPLE_RATE;
+	avg /= (find_right - find_left);
 	if (debug)
 	{
-		cerr << "Measure: desired = " << desired << " argmax = " << argmax << " max = " 
-		     << max << " left_pos = " << left_pos << " right_pos = " << right_pos << "\n";
+		cerr << "ProcessSpectrum: desired = " << desired << " avg = " << avg << " max = " << max << "\n";
 	}
-	if (abs(desired - argmax) < 1e-3)
+	if (abs(max / avg) - 1 > PEAK_LEVEL)
 	{
 		return max / pow;
 	}

Measurer.h

@@ -3,6 +3,8 @@
 
 #include <QAudioDeviceInfo>
 
+#include "fft.h"
+
 enum MeasureChannel
 {
 	Left,
@@ -11,14 +13,24 @@ enum MeasureChannel
 
 class Measurer
 {
-	const int TONE_FREQUENCY = 6000;
+	// frequency for generating tone
+	const int TONE_FREQUENCY = 1000;
+	// frequency window width to search for resulting tone 
+	const int FREQUENCY_GAP = 100;
+	// our resulting amplitude should differ at least this times from average
+	// to be considered as measured signal
+	const float PEAK_LEVEL = 1.0;
+	//
 	const int SAMPLE_RATE = 44100;
+	const float SILENCE_LEVEL = 0.3;
+	const float INITIAL_DELAY_PERCENTAGE = 0.3;
 
 	MeasureChannel channel;
 	QAudioDeviceInfo input, output;
 	bool debug;
 
 	float ProcessData(QByteArray data);
+	float ProcessSpectrum(int pow, Complex* d);
 public:
 	Measurer(QAudioDeviceInfo input, QAudioDeviceInfo output, MeasureChannel channel, bool debug);
 	float Measure(float duration);

README.md

@@ -37,3 +37,8 @@ really equation is like `Vout = R2 * Vin * k / (R1 + R2)`. So firstly we do meas
 channel (where we know both `R1` and `R2`), to determine `k`. After that we repeat same measurement
 using left channel. Knowing `Vin`, `Vout`, `R2`, and `k` we can determine value of thermistor. 
 After that we compute real temperature via approximation formula.
+
+# Notes
+
+For Linux users if you're running this from root user, you will need to have Pulseaudio started 
+system-wide. In OpenSuSE this is done via installing package `pulseaudio-system-wide`.

fft.cpp

@@ -0,0 +1,60 @@
+#include <complex>
+#include "fft.h"
+
+using namespace std;
+
+const long double PI = 3.1415926535897932384626433832795;
+
+int Rev(int i, int l)
+{
+	int res = 0;
+	int j;
+	for (j = 0; j < l; j++)
+	{
+		res <<= 1;
+		if (i & (1 << j)) res++;
+	}
+	return res;
+}
+
+void FFT(Complex *a, int n, bool inv)
+{
+	int i, j, k;
+	int logn = 0;
+	while ((1 << logn) < n) logn++;
+	for (i = 0; i < n; i++)
+	{
+		j = Rev(i, logn);
+		if (j > i)
+		{
+			Complex t = a[i];
+			a[i] = a[j];
+			a[j] = t;
+		}
+	}
+	Complex u, v;
+	for (i = 1; i <= logn; i++)
+	{
+		long double ang = 2 * PI / (1 << i);
+		if (inv) ang = -ang;
+		Complex w(cos(ang), sin(ang));
+		for (j = 0; j < n; j += (1 << i))
+		{
+			Complex wk(1, 0);
+			for (k = j; k < j + (1 << (i - 1)); k++)
+			{
+				u = a[k];
+				v = a[k + (1 << (i - 1))] * wk;
+				a[k] = u + v;
+				a[k + (1 << (i - 1))] = u - v;
+				wk *= w;
+			}
+			wk *= w;
+		}
+	}
+	if (inv)
+	{
+		for (i = 0; i < n; i++) a[i] /= n;
+	}
+}
+

fft.h

@@ -10,64 +10,5 @@ using namespace std;
 /* Defines Fast Fourier Transform */ 
 
 void FFT(Complex *a, int n, bool inv = false);
-
-/*
-  	IMPLEMENTATION
-*/
-
-const long double PI = 3.1415926535897932384626433832795;
-
-int Rev(int i, int l)
-{
-	int res = 0;
-	int j;
-	for (j = 0; j < l; j++)
-	{
-		res <<= 1;
-		if (i & (1 << j)) res++;
-	}
-	return res;
-}
-
-void FFT(Complex *a, int n, bool inv)
-{
-	int i, j, k;
-	int logn = 0;
-	while ((1 << logn) < n) logn++;
-	for (i = 0; i < n; i++)
-	{
-		j = Rev(i, logn);
-		if (j > i)
-		{
-			Complex t = a[i];
-			a[i] = a[j];
-			a[j] = t;
-		}
-	}
-	Complex u, v;
-	for (i = 1; i <= logn; i++)
-	{
-		long double ang = 2 * PI / (1 << i);
-		if (inv) ang = -ang;
-		Complex w(cos(ang), sin(ang));
-		for (j = 0; j < n; j += (1 << i))
-		{
-			Complex wk(1, 0);
-			for (k = j; k < j + (1 << (i - 1)); k++)
-			{
-				u = a[k];
-				v = a[k + (1 << (i - 1))] * wk;
-				a[k] = u + v;
-				a[k + (1 << (i - 1))] = u - v;
-				wk *= w;
-			}
-			wk *= w;
-		}
-	}
-	if (inv)
-	{
-		for (i = 0; i < n; i++) a[i] /= n;
-	}
-}
-
+
 #endif

main.cpp

@@ -4,6 +4,7 @@
 
 #include <iostream>
 #include <cmath>
+#include <unistd.h>
 
 #include "Measurer.h"
 
@@ -55,6 +56,7 @@ void Measure(QString input_device_name, QString output_device_name, float durati
 	auto output = GetOutputDevice(output_device_name);
 	Measurer left(input, output, MeasureChannel::Left, debug);
 	float left_result = left.Measure(duration);
+	//usleep(duration * 1000000);
 	Measurer right(input, output, MeasureChannel::Right, debug);
 	float right_result = right.Measure(duration);
 	if (debug)

soundthermo-cpp.pro

@@ -18,6 +18,7 @@ DEFINES += QT_DEPRECATED_WARNINGS
 SOURCES += \
         Generator.cpp \
         Measurer.cpp \
+        fft.cpp \
         main.cpp
 
 # Default rules for deployment.