デジタルフィルタの構成に直接型と転置形がありますが、juliaにて転置形の速度改善効果を実験してみました。
Down sample フィルタには高速処理が必要
近年のソフトウェア無線(SDR)では、ADコンバータの高速化が進んだことにより、IQデータをオーバーサンプリングした信号をダウンサンプルする構成が一般的になっている。なぜわざわざ必要以上に高いサンプリング周波数でサンプルするかというと、ADコンバータの量子ビット数が少なくても、フィルタで通過帯域を制限したのちダウンサンプルすると、量子化雑音が抑圧され、有効ビット数を増加させることができるようになるからである。
ダウンサンプルを行うにはそのまえにアンチエイリアスフィルタをかける必要があるのだが、フィルタの計算負荷はサンプリング周波数に比例して増加する。そこで計算負荷の小さなフィルタが必要となる。
FIRの実装
設計が簡単で直線位相特性をもつことから、まずはFIRの実装を検討する
直接型構成
（タップ長ー１）個のデータを蓄積し、入力の都度、データとタップの積和を計算し出力するとともに、次の計算に備えて蓄積データをシフトする。
$$ y_i = \sum_{k=0}^{N-1} b_{N-k} x_{i-k} $$
function direct_step(x,b,z)
# assert(size(b)[1] = size(z)[1]+1)
N=size(b)[1]
y=0
@inbounds for i=1:N-1
y=muladd(z[i],b[i],y)
end
y=muladd(x,b[N],y)
return y
end
function direct_filter!(out,x,b,z)
Nb=size(b)[1]
Nx=size(x)[1]
@inbounds for i=1:Nx
out[i]=direct_step(x[i],b,z)
@inbounds for j=1:Nb-2
z[j]=z[j+1]
end
z[Nb-1]=x[i]
end
end
転置形構成
入力毎に、i個あとの出力信号に対するタップをかけて蓄積し、タップ数分入力ののちに出力の値の計算が完了する。
$$ y_i = b_1 x_{i-(N-1)} + b_2 x_{i-(N-2)} + ... + b_N x_{i} $$
$$ y_i = ( ...( (b_1 x_{i-(N-1)}) + b_2 x_{i-(N-2)}) + ... b_{N-1} x_{i-1}) + b_N x_{i} $$
$$ y_{i+1} = ( ...( (b_1 x_{i-(N-2)}) + b_2 x_{i-(N-3)}) + ... b_{N-1} x_{i}) + b_N x_{i+1} $$
$ x_i $に対する積和を先に計算する。
function transpose_step(x,b,z)
# assert(size(b)[1] = size(z)[1]+1)
N=size(b)[1]
y=muladd(x,b[N],z[1])
@inbounds for i=1:N-2
z[i]=muladd(x,b[N-i],z[i+1])
end
z[N-1]=muladd(x,b[1],0)
return y
end
function transpose_filter!(out,x,b,z)
Nx=size(x)[1]
@inbounds for i=1:Nx
out[i]=transpose_step(x[i],b,z)
end
end
FIFフィルタの計算時間の比較
１万サンプルのAWGN信号に63タップのハーブバンドフィルタをかける時間で比較しました。
HBF=vec([-10 0 38 0 -102 0 232 0 -467 0 862 0 -1489 0 2440 0 -3833 0 5831 0 -8679 0 12803 0 -19086 0 29814 0 -53421 0 166138 262144 166138 0 -53421 0 29814 0 -19086 0 12803 0 -8679 0 5831 0 -3833 0 2440 0 -1489 0 862 0 -467 0 232 0 -102 0 38 0 -10])/524286;
Ntap=size(HBF)[1]
Nsmp=10000
xin=randn(Nsmp)
z=zeros(Ntap-1)
zt=zeros(Ntap-1)
out=zeros(Nsmp)
outt=zeros(Nsmp)
@time direct_filter!(out,xin,HBF,z)
@time transpose_filter!(outt,xin,HBF,zt)
（Juliaは1回目の実行時にコンパイルを行うため）２回目以降の実行時間で比較すると、ほぼ１０倍高速でした。
0.001164 seconds
0.000120 seconds
直接形構成では、１サンプル処理の都度データシフトする処理を行っている。この部分は、入力配列の参照点を変えることで省略することで、速度改善が期待できる。しかし、実際実装して速度を計測してみるとかえって遅くなった。
おそらく、バッファ処理の最初と最後に例外処理を行う必要があり、その処理の切り分けなどで最適化がうまくいかないのかもしれない。
転置形構成では、データの積和と保存処理がスムーズであり、高速化に貢献しているのではないか。
移動平均フィルタ
FIRの係数をすべて１とするとFIRの積和計算を加減算だけでフィルタ処理を行うことができ、計算負荷を小さくすることができる。この処理は信号の移動平均を計算することと同じであるため移動平均フィルタと呼ばれ、周波数特性はsinc関数となる。sinc関数は周期的にヌルとなる性質があり、ダウンサンプルで発生するalias信号を抑圧するのに都合がよい。ただし切れはよくないため、移動平均フィルタを直列接続することが一般的である。このように移動平均フィルタを直列接続したフィルタをCICと呼ぶ。
移動平均フィルタを転置構成で実装する。
function moving_average!(out, nR, x, si, col)# カラム指定することで二次元データをカラム毎に計算する
silen = length(si) # == nR-1
@inbounds for i=1:size(x, 1) #境界チェックをさぼる
xi = x[i,col]/nR
val = xi + si[1]
for j=1:(silen-1)
si[j] = xi + si[j+1]
end
si[silen] = xi
out[i,col] = val
end
out
end
計算速度は通常のFIRフィルタより約２０倍はやい。
移動平均フィルタの高速化
このアルゴリズムでは１サンプルあたり、N回の加算が発生する。移動平均は入力を加算し、N個前の値を減算することでも計算できる。この場合は１サンプル当たり加算１回減算１回の計算ですむ。
function moving_average!(out, nR, x, si, col)
silen = length(si) # == nR-1
nSmp=size(x)[1]
z = sum(si) +x[1,col]
out[1]=z/nR
@inbounds for i=2:nR
z += x[i,col] - si[i-1]
out[i]=z/nR
end
@inbounds for i=nR+1:size(x, 1) #境界チェックをさぼる
z += x[i,col] - x[i-nR,col]
out[i] = z /nR
end
@inbounds for i=1:silen
si[i]=x[end-nR+1+i]
end
out
end
区間長32ではさらに10倍高速になりました。