分数のおはなし

小数点以下がダラダラと続く数字を、なるべくなら分数で近似したい、という事がよくある。
嘘。ほとんどない。
でも、そういう話をしようかと思う。

唐突だけども、音声のサンプリングレートは44100Hzが一般的な値だ。
つまり、空気の振動を1秒間に44100回サンプリングして音声データに変換している。

ところで平均律音階では、隣り合う半音の周波数比はだいたい 1:1.0595 となる。
正確には 1:2^(1/12) だ。
12半音で周波数は2倍、つまり音が1オクターブ上がる。
「ラ」の音が一般的には440Hzといわれるので、3半音上の「ド」は約523.25Hzとなる。正確な値は無理数であるので、小数点以下は不規則数字の列が無限に続く。

さて、このドの音をデータにしようとする。波形は正弦波でも矩形波でも何でもいい。
サンプリングレート毎に振動の大きさを計算してもよいのだけど、その計算コストが問題になるような場合には予め計算したデータを格納しておいて、そのデータをサンプリングレート毎に呼び出すのがよいだろう。

では、ド(523.25Hz)の波を表現するのに、44100Hzでのサンプリング点はいくつあればいいだろうか。
ドの1周期は約 1/523.25 秒なので、1周期あたりの必要サンプリング回数は 1/523.25×44100 = 約54.28回となる。

当然ながらサンプリング回数は整数でなければならない。しかし、サンプリング点を54個(小数切捨)や55個(小数繰上)にすると1周期毎に音の波形が不連続な部分が現れ、その波の不連続はノイズとして聴こえる事になってしまう。

それなら、たとえば3周期分をまとめて計算して格納してはどうか。3周期では必要サンプリング回数が約162.84回となり、1周期のときより整数に近い値になる。これは、ノイズがより目立たなくなったことに相当する。
ということは、まとめる周期を適切に選べば、必要サンプリング回数がもっと整数に近づくんじゃないだろうか、という話になる。

上で目標としていることをより一般的に表現するならば、

「Nを整数として、N×(1/523.25)×44100 の値をできるだけ整数に近づけたい。」

ということになる。あるいは、

「M,Nを整数として、M/N の値が 44100/523.25 に限りなく近い組み合わせを見つけ出す。」

ということだ。さらに一般的な表現をするなら、

「ある無理数に最も近い有理数を求める。」

ということになる。

この問題の解決には「連分数(Continued Fraction)」なる概念を用いる。これは 数+分数 という式が、再帰的に第2項の分母の中に入っており、多くの場合その再帰的構造が無限に続く、というシロモノ。擬似コードでは

int[] a;
int[] b;
double continuedFraction(int idx) {
	return (double) a[idx] + (double) b[idx]/ continuedFraction(idx-1);
}

とでも表現されようか。
とくに分子(b[idx])が1であるような連分数(Simple Continued Fraction)であれば、任意の実数を連分数表現に書き直した場合のa[]を簡単に計算することが出来る。
詳しい話はwikipediaとかココ(http://mathworld.wolfram.com/ContinuedFraction.html)でも読んでもらうとして、先程のM,Nを求めるために連分数を使って何をするのかだけを簡単に述べたいと思う。

以降はAS3のコード。
まずは連分数を計算する。求めるのは先程の擬似コードのa[]に相当する。

/**
 * Calculate Simple Continued Fraction
 * @param	real 連分数で表現したい実数
 * @param	order 連分数の次数
 * @return 連分数の数列表現
 */
public static function calcSCF(real:Number, order:uint):Vector.<int>
{
	var scf :Vector.<int> = new Vector.<int>(order);
	var ir:int;
	for (var i:int = 0; i < order; i++) 
	{
		ir = int(real); //小数点以下を切り捨て
		scf[i] = ir; //連分数に格納
		real = 1 / (real - ir); //realの小数点部分の逆数が次のrealになる
	}
	return scf;
}

得られた連分数を単純な分数に戻すには以下のような計算を行う。

/**
 * Convert SCF to Fraction
 * @param	scf 数列表現の連分数
 * @return Fractionオブジェクトで表現された有理数
 */
public static function convertSCF(scf:Vector.<int>):Fraction
{
	var order:int = scf.length;
	
	// 次数がゼロだったらNaNに相当する分数を返す
	if (order == 0) return new Fraction(1, 0);
	
	var p0:int = 1;
	var q0:int = 0;
	var p1:int = scf[0];
	var q1:int = 1;
	var pTmp:int;
	var qTmp:int;
	for (var i:int = 1; i < order; i++) 
	{
		pTmp = p1 * scf[i] + p0;
		qTmp = q1 * scf[i] + q0;
		p0 = p1;
		q0 = q1;
		p1 = pTmp;
		q1 = qTmp;
	}
	// 分数クラスのコンストラクタ Fraction(分子:int, 分母:int)
	return new Fraction(p1, q1);
}

不思議なことに、この計算で得られる分数は既約分数であるので、分子と分母が必要以上に大きい値になる心配はない。


ということであらためて「ド」のサンプリングの話に戻りたい。
1周期あたりに必要なサンプリング数である 44100.0/523.2511... = 84.2807... に近い有理数はいったいいくらなのか。
計算結果は以下のようになった。

次数: 1, 連分数: [84], 有理数: ( 84 / 1 = 84 )
次数: 2, 連分数: [84,3], 有理数: ( 253 / 3 = 84.3333... )
次数: 3, 連分数: [84,3,1], 有理数: ( 337 / 4 = 84.25 )
次数: 4, 連分数: [84,3,1,1], 有理数: ( 590 / 7 = 84.2857... )
次数: 5, 連分数: [84,3,1,1,3], 有理数: ( 2107 / 25 = 84.28 )
次数: 6, 連分数: [84,3,1,1,3,1], 有理数: ( 2697 / 32 = 84.28125 )
次数: 7, 連分数: [84,3,1,1,3,1,1], 有理数: ( 4804 / 57 = 84.2807... )

つまりこの結果によれば、サンプリング数4804のデータが有れば十分に連続な（ノイズの無い）「ド」の57周期相当の波形データが得られることが分かる。


以上が連分数による任意実数の有理数近似の例であるが、他にどのような使い道があるだろう。
例えば、格子点間を結んで任意の傾きの直線を引く、などの場合が考えられる。

とりあえず終わり。