Skip to content

Instantly share code, notes, and snippets.

@hecomi
Last active August 24, 2023 09:10
Show Gist options
  • Save hecomi/7398975 to your computer and use it in GitHub Desktop.
Save hecomi/7398975 to your computer and use it in GitHub Desktop.
#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;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment