Rによる音声のサンプリング周波数変換

書いたひと: ると
2015年2月22日

三重大学の奥村晴彦先生はRを使った統計に関する資料をウェブ上で公開なさっています。その中のtuneRパッケージによる音声処理を解説したページの中で、ハイレゾ音源を通常の音質に変換する例として「次のサンプルとの平均を取ってからサンプルを間引く」という方法を取っています。この方法は高速で多くの場合十分な音質が得られますが、音質を求める際には最適な方法とは言えず、不要な雑音が付加される恐れがあります。実はtuneRパッケージにはサンプリング周波数の変換に適した関数が用意されておらず、signalパッケージのresample関数などを使う必要があります。ここでは先生が省略された部分である「平均を取って間引く方法では雑音が出る恐れがある理由」と「Rにおいて音質を重視する場合の変換方法」を例を交えて簡単に見てみます。

平均すらしないただの間引きの場合

問題をはっきりさせるために、平均すら取らずに単に間引いた場合を見てみます。音だとわかりづらいので、画像で次のような縞模様を考えてみます。

ここで偶数番目のピクセルを奇数番目のピクセルで置き換えてみます。

右側の方になにやら模様が出ています。これは折り返し雑音エイリアシングと呼ばれる現象です(元の奥村先生のページでもこのキーワードは軽く出てきますが、読者には既知のものとして説明はされていません)

例えばサンプリング周波数が48 kHzの場合、24 kHzまでの音を表せます。この性質はサンプリング定理と呼ばれています。では、24 kHzを越える音を無理矢理48 kHzでサンプリングすると何が起きるのでしょうか。その場合、その音が消えてしまうのではなく、24 kHzを越えた分を24 kHzから引いた周波数の雑音として現れます。これはサンプリング周波数が足りない場合、次の画像のように低周波の波と高周波の波でサンプリング結果が全く同じになってしまう現象によるものです。


画像はMoxfyre氏による。Wikimedia Commons: AliasingSine.svgより。ライセンスはCC-BY-SA 3.0

そのためサンプリング周波数48 kHzで録音する場合は24 kHzを越える音を先にフィルタで消去してからサンプリングする必要があります。サンプリング周波数96 kHzのデータを48 kHzにする場合も同様に24 kHzを越える音をフィルタで消してから間引く必要があります。このようなフィルタはローパスフィルタ(LPF、低い(ロー)周波数の音だけを通す(パスする)フィルタ)と言います。

なお、tuneRパッケージに含まれるdownsample関数は単純な間引きをします。

フィルタとしての平均

実は「次の値との平均を取る」という操作も高周波数の音を弱くする効果があります。例えば1, -1, 1, -1, ...というデータで平均を取ると0, 0, 0, ...となるので高周波が消えたことになります。しかし平均を取る方法の特性はあまり理想的ではありません。周波数別の特性を見ると次のようになります(参考資料: The Scientist and Engineer's Guide to Digital Signal Processing)。

理想的なフィルタは24 kHzを越える音を全て0にして、24 kHz以下の音はそのままにします。しかし平均をとる方法だと24 kHzを越える音も大分残ってしまいますし、24 kHz以下の音も少し弱くなっています。前述の縞模様について平均を取ってから間引く処理をしてみると次のようになり、単純に間引くよりはましですがまだ模様が出ています。

一方、理想的なフィルタを近似したフィルタだと次のようになります。模様は少し出ていますが、右側を中心にかなり消えています。

理想的なフィルタは一定の周波数を越える成分を全て0にして、それ以下の成分はそのままにします。しかしそのような理想的なフィルタは実は無限にデータを積分しないといけないので実現できません。実際のアプリケーションではなるべく近いフィルタで代用します。

Rでの変換例

Rでサンプリング周波数を適切に変換するにはsignalパッケージのresample関数を使います。次のように使います。

resample(x, p, q, d)

ここでxは音のデータ(tuneRのWaveの場合は@leftや@rightで取り出したデータ)、p/qが倍率、dはフィルタのパラメータで、大きいほど精度が良くなりますが、計算は遅くなります。

実際に使ってみましょう。ハイレゾの音を扱うのは再生環境などによって上手くいかないのでここではサンプリング周波数48 kHzの音を24 kHzにしてみます。20 kHzの音を作成して、平均を取る手法とresample関数を比較してみます。

library(methods)
library(tuneR)
library(signal)
 
# 48 kHzのサンプリング周波数で、20 kHzの2秒の音を生成する。
original = sine(freq=20000,
                duration=48000*2,
                samp.rate=48000,
                bit=16,
                pcm=TRUE,
                stereo=FALSE) / 2
 
 
# 奇数番目と偶数番目の平均を取る。
len = length(original@left)
left = floor((original@left[seq(1, len, 2)]
              + original@left[seq(2, len, 2)]) / 2)
averaged = Wave(left, original@right, samp.rate=24000, bit=16)
 
# 理想的なフィルタを近似したフィルタをかけてから間引き。
left = resample(original@left, 1, 2, 16)
resampled = Wave(left, original@right, samp.rate=24000, bit=16)
 
# .wavとして保存。
writeWave(original, "original.wav")
writeWave(averaged, "averaged.wav")
writeWave(resampled, "resampled.wav")
 
# 波形(最初の100点)をプロット。
svg("original_wave.svg")
plot(extractWave(original, from=1, to=100),
     main="original",
     ylim=c(-32768, 32767))
svg("averaged_wave.svg")
plot(extractWave(averaged, from=1, to=100),
     main="averaged",
     ylim=c(-32768, 32767))
svg("resampled_wave.svg")
plot(extractWave(resampled, from=1, to=100),
     main="resampled",
     ylim=c(-32768, 32767))
 
# 周波数スペクトルをプロット。
svg("original_freq.svg")
plot(periodogram(original, normalize=TRUE, width=1024, overlap=512),
     main="original",
     ylim=c(0, 1))
svg("averaged_freq.svg")
plot(periodogram(averaged, normalize=TRUE, width=1024, overlap=512),
     main="averaged",
     ylim=c(0, 1))
svg("resampled_freq.svg")
plot(periodogram(resampled, normalize=TRUE, width=1024, overlap=512),
     main="resampled",
     ylim=c(0, 1))

波形をプロットすると次のようになります。

平均を取ったものはそれなりの大きさの音が残っているのがわかります。一方resample関数を使ったものはほとんど消えています(ほんのわずか残っています)。

周波数スペクトルをプロットするとよりはっきりします。

サンプリング周波数 24 kHzでは12 kHzまでの音しか再現できず、20 kHzの音はそれより8 kHz高いので、12 - 8 = 4 kHzの雑音として現れているのがわかります。resample関数を使うと雑音はほとんどなくなりますが、それでもまだ4 kHzの雑音がわずかに残っています。

最後に音も聞いてみましょう。

original.wav

averaged.wav

resampled.wav

元の20 kHzの音はほとんど(あるいは全く)聞こえません。しかし平均を取ったものははっきりと音が聞こえます。一方resampleによるものはほとんど聞こえなくなります。

元記事で違いが出なかった理由

ではなぜ元の記事では平均を取る方法でも違いが出なかったのでしょうか。おそらく次のような要因の組み合わせであると推測します。