C言語／C++/wave5 - Sample is best!

開発環境	Microsoft Visual C++ 2010 Express (SP1)
実行環境	Microsoft Windows XP Home Edition (SP3)
プロジェクトの種類	Win32 コンソール アプリケーション
プロジェクト名	wave
アプリケーションの種類	コンソール アプリケーション
追加のオプション	空のプロジェクト

周波数の特定に最小二乗法を用いているが、相互相関でもいいかもしれない。

• 自己相関(Autocorrelation)

波長を求めるのに適しているが、倍音の判定に工夫が必要。

• フーリエ級数展開(Fourier series expansion)

ある波長の1/nの波長のサイン波がどの程度含まれているかを調べるのに適している。

wave.cpp

// 音階判定2
 
#define _USE_MATH_DEFINES
 
#include <stdio.h>
#include <math.h>
#include <vector>
 
using namespace std;
 
typedef unsigned char BYTE;
typedef struct {
	int pos;
	int noteNum;
	double vol;
} Wave;
typedef vector<Wave> VWave;
 
// note number
#define A1	33
#define A4	69
#define A5	81
 
#define SAMPLING_RATE	11025
#define DATA_NUM	1000
#define WAVE_LEN	100.23	// 波長
 
// 関数プロトタイプ宣言
double LeastSquares(int pos, double freq, double waveLen);
double Volume(int pos, int wl);
 
// 音階
const char *scale[] = {"C","C#","D","D#","E","F","F#","G","G#","A","A#","B"};
 
// 外部変数
BYTE waveformData[DATA_NUM];
VWave vWave;
 
//==============================================================================
int main()
{
	for (int i = 0; i < DATA_NUM; i++) {
		double t = fmod(i / WAVE_LEN, 1.0);
//		double y = t < 0.5 ? t * 2 : (t - 1) * 2;
		double y = sin(2 * M_PI * t);
		waveformData[i] = BYTE(128 + 64 * y);
//		printf("%d %u\n", i, waveformData[i]);
	}
 
	for (int pos = 0; pos <= DATA_NUM - 201; ) {
		int noteNum = -1;
		double jMin = 4.0;
		int wl;
		for (int n = A1; n <= A5; n++) {
			double freq = 440.0 * pow(pow(2.0, n - A4), 1.0 / 12);
			double waveLen = SAMPLING_RATE / freq;
			double j = LeastSquares(pos, freq, waveLen);
			//printf("#%d %-2s%d %.2fHz %.2f %.6f\n",
			//	n, scale[n % 12], n / 12 - 1, freq, waveLen, j);
			if (j < jMin) {
				noteNum = n;
				jMin = j;
				wl = (int)ceil(waveLen);
			}
		}
		Wave wave;
		wave.pos	= pos;
		wave.noteNum	= noteNum;
		wave.vol	= Volume(pos, wl);
		vWave.push_back(wave);
 
		printf("%.2f Note#%d\n", WAVE_LEN, noteNum);
		pos += wl;
	}
 
	for (VWave::iterator it = vWave.begin(); it != vWave.end(); it++) {
		printf("%d #%d %.6f\n", it->pos, it->noteNum, it->vol);
	}
	return 0;
}
 
// 最小二乗法
double LeastSquares(int pos, double freq, double waveLen)
{
	double sum = 0;
	int wl = (int)ceil(waveLen);
	for (int i = 0; i < wl; i++) {
		double t = i * freq / SAMPLING_RATE;
		double y = sin(2 * M_PI * t);
		double d = (waveformData[pos + i] - 128.0) / 128 - y;
		sum += d * d;
//		printf("%d t=%.6f y=%.6f\n", i, t, y);
	}
	return sum / waveLen;
}
 
// 音量
double Volume(int pos, int wl)
{
	double sum = 0;
	for (int i = 0; i < wl; i++) {
		sum += abs((waveformData[pos + i] - 128.0) / 128);
	}
	return sum / wl;
}