内挿、間引きによるレート変換
matlab の upfirdn() の例に、48KHzから44.1KHzへのサンプリングレートの変換がありました。
残念ながら、R言語 の signalパッケージには該当がないようです。
R言語 と matlab には、resample() があり、レート変換はこれを使うようです。
この2つのresample()は異なるもののようです。
レート変換とリサンプル
サンプリングで作られる音声データ(PCMデータ)は、サンプリング時のサンプリングレートで再生することが想定されています。
再生時に、異なるサンプリングレートで再生するなら「レート変換」する必要があります。
これは、PCMデータから振幅波形を再現し、新たなレートでサンプリングし直すことです。
この意味では、レート変換とリサンプルは同じ操作です。
サンプリング時や再生時に使われるサンプリングレートの値がどうなるかを抜きにして、サンプリングデータの数列に対する操作だけを取ると、リサンプルは、波形の形状を維持して、サンプル数を増減することです。
この操作は、「振幅波形の伸縮(sinc)」で作った fit_length() を使って行えます。
リサンプルは、再生時間を維持すると考えると、波形の伸縮になります。
リサンプルは、「レート変換」、「波形の伸縮」、「サンプル数の増減」などいろいろな見方ができる操作です。
R言語のresample()
R言語では、resample と関数名のみ入力すると、関数のスクリプトが表示されます。
使い方は、resample(x, p, q = 1, d = 5) で、48KHzから44.1KHzへの変換は、resample(x,44100,48000) のようにします。
-
p,q は p/q の値が使われるので、resample(x,44100,48000) と resample(x,147,160) も同じ結果になる。
-
dは、フィルタの長さを決めるもので、2d(前後にdずつ)のサンプルを補完に使う意図のようだ。
-
p/q < 1は、ダウンサンプリングを示す。この場合は、入力をFIRフィルタを使って、出力のナイキスト周波数以上をカットする。
このフィルタの係数は、fir1(2d+1, p/q) で求める。係数列の長さがデフォルトでは11となる。
この長さ(dの値)は、カットオフ周波数の精度に関係する。
-
アップサンプリングでもダウンサンプリングでも同じ処理でサンプル数を変換する。
入力をnサンプルとすると、出力は m = np/q。
出力のサンプル間の間隔を1とすると、入力のサンプル間の間隔は p/q。
出力の各サンプル値は、入力の時間的にもっとも近いものを取れば良い。
実際には、出力の時間を超えないもっとも近い値を使って、その±dの範囲で補完した値を使っている。
-
ダウンサンプリングの場合、FIRフィルタを通すが、最初のフィルタの係数の長さ(dに依存)分の出力は入力に依存していない。
この分も最終的に出力される。
これは、意図したことなのかどうかわからない。
Octaveのresample()
使い方は、resample( x, p, q[, h] )。
resample()は、fir1() と upfirdn() でできている。
upfirdn() は、upsample()、downsample()、filter()でできている。
-
pのアップサンプリング、qのダウンサンプリングを行う。
48000Hzから44100Hzに変換するには、p=147,q=160 とする必要がある。
-
hが指定されなければ、kaiser()を使って数列が用意される。
-
upfirdn() が呼び出される。
-
upfirdn() では、pのupsample(),filter(),qのdownsample() を行う。
-
upsample() は、ゼロを補ってサンプル数を増やす関数。
-
filter()は平準化を行う。
-
downsample()は間引きを行う。
-
octaveのresample()は、フィルタによる先頭部分の余分をカットしてある。
フィルタの係数列
octaveのresample()は、目的がアップサンプリングかダウンサンプリングかにかかわらず、upsample()とdownsample()の両方を呼び出します。
p>1であれば、qの値にかかわらず、upsample()とfilter()が機能します。
このフィルタの係数列は、upfirdn() の例(48KHzから44.1KHzへの変換)では、fir1()を使って作っています。
resample()では、fir1()を使っておらず、より長い数列が作られます。
左図は、長さが異なるので、中央を揃えて描きました。
upfirdn()の例では、以下のようにして算出しています。
N = 24*147;
h = fir1(N-1,1/160,kaiser(N,7.8562)) * 147;
フィルタの係数列の長さは 24*pで、最小が24(p=1)です。
fir1()の2番目の引数が通過帯域の上限で、ナイキスト周波数に対する割合を与えます。
48KHzから44.1KHzへの変換の場合、22050Hzが通過帯域の上限です。
1/160 は、22050/(48000 * 147 / 2)で、アップサンプリングした状態に対する割合で指定されています。
resample()の場合は、概ね以下のようになっているようです。
カットオフ周波数を元の周波数に対する割合で表して、
cutoff = 1 / (2*max(p,q))
この値の1/10をロールオフ帯域の幅に、
rolloff = cutoff / 10
減衰量を、
rejection = 60
この値から、係数列の長さを計算。
L = 2 * (rejection - 8) / 28.714*rolloff
p=147,q=160 のとき、
L = 2 * ceil(52 / (28.714* 1/3200)) + 1 = 11593
係数列は、以下のように計算されます。
L2=(L-1)/2;
t=(-L2:L2);
filt = 2 * p * cutoff * sinc(2*cutoff*t) * kaiser(L, 0.1102 * (rejection-8.7));
波形の外観は、sinc(1/q*t) で決まることになります。
R言語でoctaveのresample()をして見ると
R言語のresample()と、octaveのresample()は、異なっています。
R言語のresample()は、「resample」と打てば、スクリプトが表示されますし、短くて理解しやすいものです。
octaveの方は、操作やスクリプトが良くわかっていないこともあって難解です。
処理内容を類推して、R言語のスクリプトで表現して見ます。
upsample(),downsample()がないので、seq()でインデクスを生成して代入することで代替しました。
フィルタの出力の最初の方は、入力データと無関係な値になります。フィルタの段数分進んだところからが有効なデータになります。
R言語のresample()は、これをそのまま出力しています。このことは、終端のデータも、フィルタに残って出力に反映していないことになります。
octaveの方は、先頭をカットしてして、後端も押し出しているようです。
後端を押し出すためにフィルタの段数以上のゼロをデータの後半に加えることにします。フィルタの係数列を長くする理由はわかりません。
図は、1000Hzのサイン波のテスト結果です。
48000Hzで10msサンプリングしたデータを入力にしています。
このデータを青で書いています。このデータは、480個です。
これを44100Hzでリサンプルする場合のサンプリングポイントが赤です。同じ10msで441個のデータになります。
- # テスト1
- # 48000Hzで1000Hzのサイン波を0.01秒間サンプリングすると
- t <- seq(0,by=1/48000,length.out=480)
- x <- sin(2*pi*1000*t)
- # リサンプリング
- y <- octave_resample(x,147,160)
- # オリジナルを表示
- tx <- seq(0,by=1/48000,length.out=480)
- plot(tx,x,type="l",col="blue",xlab="秒",ylab="振幅")
- # 新しいサンプリング箇所を表示
- par(new=T)
- ty <- seq(0,by=1/44100,length.out=441)
- plot(ty,y,type="p",col="red",xlab="",ylab="")
- title(main="48K->44.1K Hz rate convert")
- #----------------------------------------------
- # テスト2
- y <- octave_resample(x,2,1)
- length(y)
- plot(y,type="l")
- windows()
- y <- octave_resample(x,1,2)
- length(y)
- plot(y,type="l")
octave_resample()
- # Octaveのresample()の動作の類推
- octave_resample <- function(x,p,q)
- {
- # フィルタ係数列の算出
- cutoff <- 1 / (2*max(p,q))
- rolloff <- cutoff / 10
- rejection <- 60
- L <- 2 * ceiling((rejection - 8) / (28.714 * rolloff)) + 1
- L2 <- (L-1)/2
- t <- (-L2:L2)
- filt <- 2 * p * cutoff * sinc(2*cutoff*t)
- filt <- filt * kaiser(L, 0.1102 * (rejection-8.7))
- filt <- ifelse(is.nan(filt),1,filt)
- # フィルタによる先頭部分の不要部分の除去のために、フィルタ係数列とデータを調整
- fp <- floor(q-(L2%%q))
- filt <- c(rep(0,fp),filt) # フィルタ係数列の長さをfpだけ伸ばす
- dp <- floor((L2+fp)/q)
- np <- 0
- m <- ceiling(length(x)*p/q)
- n <- p*(length(x)-1) + fp + L
- while( ceiling((n + np) / q) - dp < m ) {
- np <- np+1
- }
- filt <- c(filt,rep(0,np)) # フィルタ係数列の後部にゼロ追加
- x <- c(x,rep(0,dp*2)) # 押し出すためにデータ後部にゼロ追加
- # upsample()
- y <- rep(0,length(x)*p)
- y[seq(1,by=p,length.out=length(x))] <- x
- # フィルタ
- x <- filter(filt,a=1,x=y)
- # downsample()
- y <- x[seq(1,length(x),by=q)]
- y <- y[(dp+1):(dp+m)] # 先頭の無効な部分をカット
- }
|