やや古い音の波形も参照。

準備(Mac)

以下ではWAVなどのオーディオファイルを扱う。MacではWAVファイルなどはおそらくiTunes.appに関連付けられているが,FinderでQuickTime Player.appまたはVLC.appに関連付けを変更しておくと便利である。ターミナルで例えばQuickTime Player.appでWAVファイルを開くには

open -a 'QuickTime Player' *.wav

と打ち込めばよい。GUIが不要ならば

afplay *.wav

のほうが便利である。Ctrl-C で止まる。ほかに名前がaf(audio file)で始まるコマンドとして afinfoafconvertafhashafida がある。いずれも --help(または -h)オプションで簡単な説明が出る。

Rでオーディオを扱う定番のパッケージはtuneRである(2013年8月のバージョン1.0で大きく仕様が変わっている)。以下の例をRで実行するにはあらかじめ

install.packages("tuneR")
library(tuneR)

と打ち込んでおく。このパッケージの play() という関数でWAVファイルなどをプレイすることができる。そのためにはあらかじめMacなら

setWavPlayer("/usr/bin/afplay")

と打ち込んでおく。

以下ではFLACというオープンなロスレス圧縮形式も必要になるので,適当な方法でインストールしておく。Windows用バイナリなどもダウンロードできる。

MacでHomebrewを使っているなら brew install flac と打ち込むだけでインストールできる。ソースからインストールするには次のようにする:

tar xvzf flac-*.tar.xz
cd flac-*
./configure
make
make install

Musopen Collection as FLAC

Internet ArchiveMusopen Collection as FLACは,パブリックドメインの音源で,しかもFLAC(オープンなロスレス圧縮)形式である。ここでは最初のGoldberg変奏曲のAria(FLACで13519283バイト)を使う。1997年のShelley Katzのピアノ演奏である。

これをまずFLACでWAVに変換する:

flac -d Bach_GoldbergVariations-JohannSebastianBach-01-GoldbergVariationsBwv.988-Aria.flac

77376140バイトのWAVファイルができる。

Rでこれを読み込む:

library(tuneR)
bach = readWave("Bach_GoldbergVariations-JohannSebastianBach-01-GoldbergVariationsBwv.988-Aria.wav")
object.size(bach)
206337672 bytes
str(bach)
Formal class 'Wave' [package "tuneR"] with 6 slots
  ..@ left     : num [1:12896016] 0 0 0 0 0 0 0 0 0 0 ...
  ..@ right    : num [1:12896016] 0 0 0 0 0 0 0 0 0 0 ...
  ..@ stereo   : logi TRUE
  ..@ samp.rate: int 44100
  ..@ bit      : int 24
  ..@ pcm      : logi TRUE

readWave() で読み出したオブジェクトは「Wave」というS4クラスであり,各スロットにアクセスするには $ ではなく @ を使う。例えば左チャンネルは bach@left または slot(bach, "left") のようにしてアクセスする。

オブジェクトのサイズを見てわかるように(あるいは leftrightnum であることからわかるように),各サンプル値は8バイトのnumeric型(double型)に格納されている。

冒頭の度数分布を調べてみる:

table(bach@left[1:100000])

-2304 -2048 -1792 -1536 -1280 -1024  -768  -512  -256     0   256   512   768 
    5   167   353   468   830  1101  1446 13575  4101 72289  1133  1124  1117 
 1024  1280  1536  1792  2048  2304  2560  2816 
  696   447   424   321   236   130    33     4 

あれあれ,すべて256の倍数だ。せっかく24ビット音だと思ったら,16ビットだ。全体を256で割ってしまおう。

bachl = as.integer(bach@left / 256)
bachr = as.integer(bach@right / 256)
rm(bach) # メモリが逼迫したら,もとのデータは消す

適当なところをプロットしてみる:

s = 44100*4    # 4秒後から
t = 44100*10   # 10秒後まで
plot  ((s:t)/44100, bachl[s:t], type="l", col="#0068b780", xlab="", ylab="")
points((s:t)/44100, bachr[s:t], type="l", col="#f3980080")

:波形がわかるように,例えば6秒〜6.01秒をプロットしてみよ。波形から音の振動数(周波数)を求めよ。

全体のヒストグラムを描いてみる。左右チャンネルはアペンドする:

a = append(bachl, bachr)
length(a)
[1] 25792032
range(a)
[1] -9322  8909
hist(a, main="", xlab="", ylab="", col="gray")

実はこのように広いヒストグラムになるのは階級の幅が広すぎるためである。階級の幅を1にして比べてみよう。もっと0に集中している。

hist(a, main="", xlab="", ylab="", col="gray", breaks=20000)

もっと詳しく見るために,プロットせず度数分布だけ取り出してみよう:

h = hist(a, breaks=20000, plot=FALSE)
str(h)  # hの構造を確認
List of 6
 $ breaks  : num [1:18232] -9322 -9321 -9320 -9319 -9318 ...
 $ counts  : int [1:18231] 2 0 0 0 0 0 0 0 0 0 ...
 $ density : num [1:18231] 7.75e-08 0.00 0.00 0.00 0.00 ...
 $ mids    : num [1:18231] -9322 -9320 -9320 -9318 -9318 ...
 $ xname   : chr "a"
 $ equidist: logi TRUE
 - attr(*, "class")= chr "histogram"
sum(h$density)
[1] 1

上では分布の範囲がほぼ2万とわかっていたので breaks=20000 としたが,一般には breaks=-32768:32767 のようにすれば16ビット整数の全範囲がカバーできる。

h$density に確率分布が入る。この分布のエントロピーを求めてみよう:

-sum(h$density * log2(h$density), na.rm=TRUE)
[1] 11.94497

約12ビット/サンプルである。

エントロピーを求めるだけなら hist() でなく table() を使ってもできる:

t = table(a)
p = t / sum(t)
-sum(p * log2(p))
[1] 11.94497

ところで,a の代わりに,その差分 diff(a) を使ってみたらどうなるか? 差分の差分(2階差分)ではどうか? …

audio というパッケージにもWAVを読み書きする関数があるが,私の環境ではうまく動作しない。

seewave パッケージに,WAVとFLACを相互変換する wav2flac() という関数があるが,単にFLACを呼び出しているだけである。

別の方法として,WAVファイルを拙作dumpwave.cでテキストに落としてから読み込む方法を紹介しておく。この方法なら,各サンプル値を4バイトのinteger型で読み込むことができる。

$ gcc dumpwave.c -o dumpwave
$ ./dumpwave Bach_GoldbergVariations-JohannSebastianBach-01-GoldbergVariationsBwv.988-Aria.wav >x.txt
[RIFF] (77376132 bytes)
[WAVEfmt ] (16 bytes)
  Data type = 1 (1 = PCM)
  Number of channels = 2 (1 = mono, 2 = stereo)
  Sampling rate = 44100Hz
  Bytes / second = 264600
  Bytes x channels = 6
  Bits / sample = 24

[data] (77376096 bytes)

できたテキストファイルは179573070バイト(約180MB)ある。これをRで読み込むには read.table() でもよいが,非常に遅いので,data.tableを使うとよい:

library("data.table")
x = fread("x.txt")
dim(x)
[1] 12896016        2
head(x)
   V1 V2
1:  0  0
2:  0  0
3:  0  0
4:  0  0
5:  0  0
6:  0  0
object.size(x)
103169360 bytes
# または
library(pryr)
object_size(x)
103 MB
str(x)
Classes ‘data.table’ and 'data.frame':	12896016 obs. of  2 variables:
 $ V1: int  0 0 0 0 0 0 0 0 0 0 ...
 $ V2: int  0 0 0 0 0 0 0 0 0 0 ...
 - attr(*, ".internal.selfref")=<externalptr>

変数名が V1V2 ではわかりにくいので,LR に直しておく:

setnames(x, c("L","R"))
# names(x) = c("L", "R") ではx全体をコピーする無駄が生じる
# data.tableではsetnames()を使うべき

SONYサンプル音源

2015年2月時点で,SONYのハイレゾ用「高音質」miniSDカードが話題になっている:

上の3番目の記事の3ページ目には次のような記述がある:

S/Nが格段に高まり、ボーカルは透明感やディティールの見晴らしが良くなる。バックバンドによる演奏は奥行きの深みが増して、楽器固有の音色がより濃く表れるようになった。音の粒立ちが明瞭で、特に高域が煌びやかでエネルギッシュだ。その違いは、何度も聴き直せばそんな気がしてくるというレベルではなく、ファーストインプレッションでクリアにわかるレベルだった。

4番目の記事によれば,「静かな環境下で聴き比べれば、大半の人が瞬時にわかるレベルの違いがある」。

真相はともかく,ここではSONYのサンプル音源ダウンロードからダウンロードしたハイレゾデータを調べてみよう。「Blue Monday FM “Bee Moved”より」,ZIPファイルで提供されているが,中身のFLACファイルのほうがサイズが若干小さい。展開してFLACでWAVに変換し,Rで読み込む:

library(tuneR)
bee = readWave("Sample_BeeMoved_96kHz24bit.wav")
str(bee)
Formal class 'Wave' [package "tuneR"] with 6 slots
  ..@ left     : num [1:3828096] -2 1 2 -8 11 -13 9 -6 7 3 ...
  ..@ right    : num [1:3828096] 0 3 -5 8 -11 16 -10 -2 -1 7 ...
  ..@ stereo   : logi TRUE
  ..@ samp.rate: int 96000
  ..@ bit      : int 24
  ..@ pcm      : logi TRUE
range(bee@left)
[1] -8292596  8292598
range(bee@right)
[1] -8292595  8292598

24ビット(±8388607)ほぼフルに使っているところがすごい。

t = table(a)
p = t / sum(t)
-sum(p * log2(p))
[1] 22.04475

これを48000サンプル/秒,16ビットに丸めてみる:

left = floor(((bee@left[seq(1,3828096,2)] + bee@left[seq(2,3828096,2)]) / 2 + 128) / 256)
right = floor(((bee@right[seq(1,3828096,2)] + bee@right[seq(2,3828096,2)]) / 2 + 128) / 256)
bee2 = Wave(left, right, samp.rate=48000, bit=16)
writeWave(bee2, "bee2.wav")  # あるいは play(bee2)

聴き比べてみたが,再生装置がハイレゾ対応していないので当然ながらまったく同じ音に聞こえた。

[追記] ると @cocoa_ruto さんが Rによる音声のサンプリング周波数変換 を書いてくださった。要は,単純に隣同士の平均をとるだけでは折り返し雑音(エイリアシング)が生じる。次のようにするのが正しい(16 は調整可能):

library(signal)
left = resample(bee@left, 1, 2, 16)
right = resample(bee@right, 1, 2, 16)
range(left)
[1] -8515382  8356196
range(right)
[1] -8352250  8552039

あれあれ,24ビットの範囲を超えてしまった。プロットしてみる:

t = 1257596
r = (t-20):(t+20)
r2 = (2*t-40):(2*t+40)
plot(r2, bee@right[r2], type="o", pch=16, xlab="", ylab="")
points(r*2, right[r+8], type="o", pch=1)

16ビット化するときには256より少し大きい値で割る必要がある:

bee3 = Wave(floor((left + 130.5) / 261), floor((right + 130.5) / 261),
            samp.rate=48000, bit=16)
writeWave(bee3, "bee3.wav")  # または play(bee3)

あるいは全体の音量を変えたくないなら pmax(pmin(x, 32767), -32768) のような関数を通して16ビットの範囲に制限するだけでもいいかもしれない。

ハイレゾの効果を評価する際には,圧縮音と聴き比べるのではなく,上記のような方法でサンプリングレートや量子化ビット数を(個別に・両方)通常の値に落としたものと聴き比べるべきであろう。また,ハイレゾとノーマルをランダムに切り替えて聴き比べ,後で正解と比較する二重盲検法を使うべきである。そういうコードはRで簡単に書ける。ハイレゾ再生装置をお持ちの若い方にぜひお願いしたい(私のような老人は20kHzなんてとても無理)。

人工音

高い音がどこまで聞こえるか確かめてみる。

library(tuneR)
s20k = sine(freq=20000, duration=96000*10, samp.rate=96000, bit=16, pcm=TRUE, stereo=TRUE)
writeWave(s20k, "s20k.wav")

のようにしてファイルに出力してから再生するか,あるいは設定さえ正しければ

play(sine(freq=10000, duration=44100*5, samp.rate=44100, bit=16, pcm=TRUE, stereo=TRUE))

でも音が出る。私の環境では,高い音は標本化周波数によって聞こえたり聞こえなかったりする。折り返し雑音みたいな音が聞こえることもある。再生装置と私の耳のどちらが効いているかわからないが,44100Hz標本化で試すと,14200〜14500Hzあたりに上限があるようだ。ちなみに,アナログテレビの時代には,雨の日は15750Hzの水平同期信号がよく聞こえた。若者だけに聞こえる17kHzあたりの「モスキート音」を使った「若者よけ」装置についてもしばしば話題になる(例:「若者よけ」のモスキート音、課題は利用ルール

[追記] この項に関しても Rによる音声のサンプリング周波数変換 を参照されたい。

周波数成分

このページの最初のほうで使った例(バッハ)の周波数成分を調べる:

library(tuneR)
bach = readWave("Bach_GoldbergVariations-JohannSebastianBach-01-GoldbergVariationsBwv.988-Aria.wav")
p = periodogram(Wave(bach@left / 256))
plot(p)

1000Hz以下に集中しているように見える。

plot(p, xlim=c(0,1000))

周波数の高い成分はほぼ0なので,対数グラフにしないとよく見えない。

plot(p, log="y", ylim=c(1e-13,1e-3))

中音域(440Hz前後)と比べて,10kHz以上は8桁の違いがある。periodogram() の縦軸はエネルギーなので,振幅にすれば4桁の違い。元の音が ±10000 に収まっているので,10kHz以上は ±1 つまり量子化ノイズなのだろう。実際,サンプリング周波数を半分に resample() してみても,区別がつかない(私の耳 and/or 再生装置の問題も多分にある。追試してください)。

その他のサンプル音源