rec.cpp

#include <fstream>
#include <cstdio>
#include <string>
#include <memory>
#include <complex>
#include <vector>
#include <cmath>
#include <sndfile.h>
#include <fftw3.h>
#include <portaudio.h>
#include <boost/format.hpp>

using std::complex;
using std::vector;
using std::pair;


// LPC 係数の次数
const int LPC_ORDER = 64;

// 解析に用いる微小区間 (sec)
const double WINDOW_DURATION = 0.04;

// FFT する
vector<double> fft(const vector<double>& input);

// LPC スペクトル包絡を得る
vector<double> lpc(const vector<double>& input, int order, double df);

// [1 : -1] に正規化
vector<double> normalize(const vector<double>& input);

// ハミング窓をかける
vector<double> hamming(const vector<double>& input);

// デジタルフィルタ
vector<double> freqz(const vector<double>& b, const vector<double>& a, double df, int N);

// フォルマント f1 / f2 を返す
pair<double, double> formant(const vector<double>& input, double df);

// フォルマントから母音を推定する
std::string vowel(double f1, double f2);

// 与えられた範囲の音の平均
double volume(const vector<double>& input);

// PortAudio のエラーチェック
void check_error(PaError code);

// PortAudio の録音のコールバック
int recordCallback(
	const void *input,
	void *output,
	unsigned long frameCount,
	const PaStreamCallbackTimeInfo* timeInfo,
	PaStreamCallbackFlags statusFlags,
	void *userData
);


int main()
{
	// 初期化
	check_error( Pa_Initialize() );

	std::shared_ptr<void> terminator(nullptr, [](void*) {
		check_error( Pa_Terminate() );
		std::cout << "terminate" << std::endl;
	});

	// ストリームを開く
	PaStream *stream;
	auto err = Pa_OpenDefaultStream(
		 &stream, 1, 0, paInt16, 44100, 1024, recordCallback, nullptr);
	check_error(err);

	std::shared_ptr<void> closer(nullptr, [&](void*) {
		check_error( Pa_CloseStream(stream) );
		std::cout << "close" << std::endl;
	});

	// 録音の開始と停止
	check_error( Pa_StartStream(stream) );
	getchar();
	check_error( Pa_StopStream(stream) );

	return 0;
}


// FFT
vector<double> fft(const vector<double>& input)
{
	const int frame = input.size();
	vector<complex<double>> in(frame), out(frame);
	auto plan = fftw_plan_dft_1d(
			frame,
			reinterpret_cast<fftw_complex*>(&in[0]),
			reinterpret_cast<fftw_complex*>(&out[0]),
			FFTW_FORWARD, FFTW_ESTIMATE
	);
	for (int i = 0; i < frame; ++i) {
		in[i] = complex<double>(input[i], 0.0);
	}
	fftw_execute(plan);

	vector<double> fft_result(frame);
	for (int i = 0; i < frame; ++i) {
		fft_result[i] = abs(out[i]);
	}

	return fft_result;
}


// LPC
vector<double> lpc(const vector<double>& input, int order, double df)
{
	const int N = input.size();

	// 自己相関関数
	const int lags_num = order + 1;
	vector<double> r(lags_num);
	for (int l = 0; l < lags_num; ++l) {
		r[l] = 0.0;
		for (int n = 0; n < N - l; ++n) {
			r[l] += input[n] * input[n + l];
		}
	}

	// Levinson-Durbin のアルゴリズムで LPC 係数を計算
	vector<double> a(order + 1, 0.0), e(order + 1, 0.0);
	a[0] = e[0] = 1.0;
	a[1] = - r[1] / r[0];
	e[1] = r[0] + r[1] * a[1];
	for (int k = 1; k < order; ++k) {
		double lambda = 0.0;
		for (int j = 0; j < k + 1; ++j) {
			lambda -= a[j] * r[k + 1 - j];
		}
		lambda /= e[k];

		vector<double> U(k + 2), V(k + 2);
		U[0] = 1.0; V[0] = 0.0;
		for (int i = 1; i < k + 1; ++i) {
			U[i] = a[i];
			V[k + 1 - i] = a[i];
		}
		U[k + 1] = 0.0; V[k + 1] = 1.0;

		for (int i = 0; i < k + 2; ++i) {
			a[i] = U[i] + lambda * V[i];
		}

		e[k + 1] = e[k] * (1.0 - lambda * lambda);
	}

	// LPC 係数から音声信号再現
	// vector<double> lpc_result(N, 0.0);
	// for (int i = 0; i < N; ++i) {
	// 	if (i < order) {
	// 		lpc_result[i] = input[i];
	// 	} else {
	// 		for (int j = 1; j < order; ++j) {
	// 			lpc_result[i] -= a[j] * input[i + 1 - j];
	// 		}
	// 	}
	// }
	// return lpc_result;

	return freqz(e, a, df, N);
}


// 正規化
vector<double> normalize(const vector<double>& input)
{
	// 最大 / 最小値
	auto max = abs( *std::max_element(input.begin(), input.end()) );
	auto min = abs( *std::min_element(input.begin(), input.end()) );
	double factor = std::max(max, min);
	vector<double> result( input.size() );
	std::transform(input.begin(), input.end(), result.begin(), [factor](double x) {
		return x / factor;
	});
	return result;
}


// ハミング窓をかける
vector<double> hamming(const vector<double>& input)
{
	const double N = input.size();
	vector<double> result(N);
	for (int i = 1; i < N - 1; ++i) {
		const double h = 0.54 - 0.46 * cos(2 * M_PI * i / (N - 1));
		result[i] = input[i] * h;
	}
	result[0] = result[N - 1] = 0;
	return result;
}


// デジタルフィルタ
vector<double> freqz(const vector<double>& b, const vector<double>& a, double df, int N)
{
	vector<double> H(N);
	for (int n = 0; n < N; ++n) {
		auto z = std::exp(complex<double>(0.0, -2.0 * M_PI * n / N));
		complex<double> numerator(0.0, 0.0), denominator(0.0, 0.0);
		for (int i = 0; i < b.size(); ++i) {
			numerator += b[b.size() - 1 - i] * pow(z, i);
		}
		for (int i = 0; i < a.size(); ++i) {
			denominator += a[a.size() - 1 - i] * pow(z, i);
		}
		H[n] = abs(numerator / denominator);
	}

	return H;
}


// フォルマント f1 / f2 を返す
pair<double, double> formant(const vector<double>& input, double df)
{
	pair<double, double> result(0.0, 0.0);
	bool is_find_first = false;
	for (int i = 1; i < input.size() - 1; ++i) {
		if (input[i] > input[i-1] && input[i] > input[i+1]) {
			if (!is_find_first) {
				result.first = df * i;
				is_find_first = true;
			} else {
				result.second = df * i;
				break;
			}
		}
	}
	return result;
}


// 母音を推定
// NOTE: 機械学習にする予定
std::string vowel(double f1, double f2)
{
	if (f1 > 600 && f1 < 1400 && f2 > 900  && f2 < 2000) return "あ";
	if (f1 > 100 && f1 < 410  && f2 > 1900 && f2 < 3500) return "い";
	if (f1 > 100 && f1 < 700  && f2 > 1100 && f2 < 2000) return "う";
	if (f1 > 400 && f1 < 800  && f2 > 1700 && f2 < 3000) return "え";
	if (f1 > 300 && f1 < 900  && f2 > 500  && f2 < 1300) return "お";
	return "-";
}


// ボリューム
double volume(const vector<double>& input)
{
	double v = 0.0;
	std::for_each(input.begin(), input.end(), [&v](double x) { v += x*x; });
	v /= input.size();
	return v;
}


// エラーチェック
void check_error(PaError code)
{
	if (code == paNoError) { return; }
	std::cout << Pa_GetErrorText(code) << std::endl;
	exit(EXIT_FAILURE);
}


// 録音のコールバック
int recordCallback(
	const void *input,
	void *output,
	unsigned long frameCount,
	const PaStreamCallbackTimeInfo* timeInfo,
	PaStreamCallbackFlags statusFlags,
	void *userData
) {
	const auto input_data = reinterpret_cast<const short*>(input);
	const double df = 44100.0 / 1024.0;

	vector<double> data(frameCount + 1);
	for (int i = 0; i < frameCount; ++i) {
		data[i] = input_data[i];
	}
	auto hamming_result = normalize( hamming(data) );
	auto lpc_result = normalize( lpc(hamming_result, LPC_ORDER, df) );
	auto fft_result = normalize( fft(hamming_result) );

	// 得られた LPC スペクトル包絡線からフォルマントを抽出
	auto formant_result = formant(lpc_result, df);
	double f1 = formant_result.first;
	double f2 = formant_result.second;

	if (volume(data) < 1e4) {
		std::cout << "-";
	} else {
		std::cout << boost::format("%1%  (f1:%2%, f2:%3%)")
			% vowel(f1, f2) % f1 % f2;
	}
	std::cout << std::endl;

	return 0;
}