AI エンジンへの IIR フィルターのインプリメント - パート 1a
バージョン: Vitis 2024.1
下準備
IIR フィルターの転置直接形式 II (TDF2) を使用して数値的なロバスト性を高め、双 2 次セクションのインプリメントに焦点を当て、複数のセクションをカスケードして高次の IIR フィルターを形成できるようにします。
TDF2 の IIR 双 2 次 フィルターのブロック図は次のとおりです。
TDF2 の IIR 双 2 次フィルターには、次の差分方程式があります。
最初の 8 つの入力サンプル (すなわち、x[0]、x[1]...、 x[7]) から開始し、4 つ前のステートは x[-1]= x[-2]= y[-1]= y[-2]=0 に初期化されます。この結果、対応する 8 つの出力は次のように示すことができます。
この方程式は逐次的に解くことも並列的に解くこともできます。ベクター プロセッサの SIMD 機能をフルに活用するには、このシステムの方程式を並列化する必要があります [2]。
注記: すべての出力を、現在の入力と 4 つ前の状態で表すことができます。たとえば、y[n + 0] の式を置換して、現在の入力 (x[n]~x[n + 1]) と 4 つのステート (x[n-1]、x[n-2]、y[n-1]、y[n-2]) を解きます。
wxMaxima または Mathematica ような記号ソルバーを使用すると、次のようになります。
説明
*************************************
Ky0_ym2 = -a2;
Ky0_ym1 = -a1;
Ky0_xm2 = b2*K;
Ky0_xm1 = b1*K;
Ky0_x0 = b0*K;
*************************************
Ky1_ym2 = a1*a2;
Ky1_ym1 = a1^2 - a2;
Ky1_xm2 = -(a1*b2*K);
Ky1_xm1 = (-(a1*b1) + b2)*K;
Ky1_x0 = (-(a1*b0) + b1)*K;
Ky1_x1 = b0*K;
*************************************
Ky2_ym2 = a2*(-a1^2 + a2);
Ky2_ym1 = -a1^3 + 2*a1*a2;
Ky2_xm2 = (a1^2 - a2)*b2*K;
Ky2_xm1 = (a1^2*b1 - a2*b1 - a1*b2)*K;
Ky2_x0 = (a1^2*b0 - a2*b0 - a1*b1 + b2)*K;
Ky2_x1 = (-(a1*b0) + b1)*K;
Ky2_x2 = b0*K;
*************************************
Ky3_ym2 = a1*(a1^2 - 2*a2)*a2;
Ky3_ym1 = a1^4 - 3*a1^2*a2 + a2^2;
Ky3_xm2 = -(a1*(a1^2 - 2*a2)*b2*K);
Ky3_xm1 = (-(a1^3*b1) + 2*a1*a2*b1 + a1^2*b2 - a2*b2)*K;
Ky3_x0 = -((a1^3*b0 - a1^2*b1 + a2*b1 + a1*(-2*a2*b0 + b2))*K);
Ky3_x1 = (a1^2*b0 - a2*b0 - a1*b1 + b2)*K;
Ky3_x2 = (-(a1*b0) + b1)*K;
Ky3_x3 = b0*K;
*************************************
Ky4_ym2 = -(a2*(a1^4 - 3*a1^2*a2 + a2^2));
Ky4_ym1 = -(a1*(a1^4 - 4*a1^2*a2 + 3*a2^2));
Ky4_xm2 = (a1^4 - 3*a1^2*a2 + a2^2)*b2*K;
Ky4_xm1 = (a1^4*b1 - 3*a1^2*a2*b1 + a2^2*b1 - a1^3*b2 + 2*a1*a2*b2)*K;
Ky4_x0 = (a1^4*b0 - a1^3*b1 + 2*a1*a2*b1 + a2*(a2*b0 - b2) + a1^2*(-3*a2*b0 + b2))*K;
Ky4_x1 = -((a1^3*b0 - a1^2*b1 + a2*b1 + a1*(-2*a2*b0 + b2))*K);
Ky4_x2 = (a1^2*b0 - a2*b0 - a1*b1 + b2)*K;
Ky4_x3 = (-(a1*b0) + b1)*K;
Ky4_x4 = b0*K;
*************************************
Ky5_ym2 = a1*a2*(a1^4 - 4*a1^2*a2 + 3*a2^2);
Ky5_ym1 = a1^6 - 5*a1^4*a2 + 6*a1^2*a2^2 - a2^3;
Ky5_xm2 = -(a1*(a1^4 - 4*a1^2*a2 + 3*a2^2)*b2*K);
Ky5_xm1 = (-(a1^5*b1) + 4*a1^3*a2*b1 - 3*a1*a2^2*b1 + a1^4*b2 - 3*a1^2*a2*b2 + a2^2*b2)*K;
Ky5_x0 = (-(a1^5*b0) + a1^4*b1 - 3*a1^2*a2*b1 + a2^2*b1 + a1^3*(4*a2*b0 - b2) + a1*a2*(-3*a2*b0 + 2*b2))*K;
Ky5_x1 = (a1^4*b0 - a1^3*b1 + 2*a1*a2*b1 + a2*(a2*b0 - b2) + a1^2*(-3*a2*b0 + b2))*K;
Ky5_x2 = -((a1^3*b0 - a1^2*b1 + a2*b1 + a1*(-2*a2*b0 + b2))*K);
Ky5_x3 = (a1^2*b0 - a2*b0 - a1*b1 + b2)*K;
Ky5_x4 = (-(a1*b0) + b1)*K;
Ky5_x5 = b0*K;
*************************************
Ky6_ym2 = a2*(-a1^6 + 5*a1^4*a2 - 6*a1^2*a2^2 + a2^3);
Ky6_ym1 = -a1^7 + 6*a1^5*a2 - 10*a1^3*a2^2 + 4*a1*a2^3;
Ky6_xm2 = (a1^6 - 5*a1^4*a2 + 6*a1^2*a2^2 - a2^3)*b2*K;
Ky6_xm1 = (a1^6*b1 - 5*a1^4*a2*b1 + 6*a1^2*a2^2*b1 - a2^3*b1 - a1^5*b2 + 4*a1^3*a2*b2 - 3*a1*a2^2*b2)*K;
Ky6_x0 = (a1^6*b0 - a1^5*b1 + 4*a1^3*a2*b1 - 3*a1*a2^2*b1 + 3*a1^2*a2*(2*a2*b0 - b2) + a1^4*(-5*a2*b0 + b2) + a2^2*(-(a2*b0) + b2))*K;
Ky6_x1 = (-(a1^5*b0) + a1^4*b1 - 3*a1^2*a2*b1 + a2^2*b1 + a1^3*(4*a2*b0 - b2) + a1*a2*(-3*a2*b0 + 2*b2))*K;
Ky6_x2 = (a1^4*b0 - a1^3*b1 + 2*a1*a2*b1 + a2*(a2*b0 - b2) + a1^2*(-3*a2*b0 + b2))*K;
Ky6_x3 = -((a1^3*b0 - a1^2*b1 + a2*b1 + a1*(-2*a2*b0 + b2))*K);
Ky6_x4 = (a1^2*b0 - a2*b0 - a1*b1 + b2)*K;
Ky6_x5 = (-(a1*b0) + b1)*K;
Ky6_x6 = b0*K;
*************************************
Ky7_ym2 = a1*a2*(a1^6 - 6*a1^4*a2 + 10*a1^2*a2^2 - 4*a2^3);
Ky7_ym1 = a1^8 - 7*a1^6*a2 + 15*a1^4*a2^2 - 10*a1^2*a2^3 + a2^4;
Ky7_xm2 = -(a1*(a1^6 - 6*a1^4*a2 + 10*a1^2*a2^2 - 4*a2^3)*b2*K);
Ky7_xm1 = (-(a1^7*b1) + 6*a1^5*a2*b1 - 10*a1^3*a2^2*b1 + 4*a1*a2^3*b1 + a1^6*b2 - 5*a1^4*a2*b2 + 6*a1^2*a2^2*b2 - a2^3*b2)*K;
Ky7_x0 = -((a1^7*b0 - a1^6*b1 + 5*a1^4*a2*b1 - 6*a1^2*a2^2*b1 + a2^3*b1 + 2*a1^3*a2*(5*a2*b0 - 2*b2) + a1^5*(-6*a2*b0 + b2) + a1*a2^2*(-4*a2*b0 + 3*b2))*K);
Ky7_x1 = (a1^6*b0 - a1^5*b1 + 4*a1^3*a2*b1 - 3*a1*a2^2*b1 + 3*a1^2*a2*(2*a2*b0 - b2) + a1^4*(-5*a2*b0 + b2) + a2^2*(-(a2*b0) + b2))*K;
Ky7_x2 = (-(a1^5*b0) + a1^4*b1 - 3*a1^2*a2*b1 + a2^2*b1 + a1^3*(4*a2*b0 - b2) + a1*a2*(-3*a2*b0 + 2*b2))*K;
Ky7_x3 = (a1^4*b0 - a1^3*b1 + 2*a1*a2*b1 + a2*(a2*b0 - b2) + a1^2*(-3*a2*b0 + b2))*K;
Ky7_x4 = -((a1^3*b0 - a1^2*b1 + a2*b1 + a1*(-2*a2*b0 + b2))*K);
Ky7_x5 = (a1^2*b0 - a2*b0 - a1*b1 + b2)*K;
Ky7_x6 = (-(a1*b0) + b1)*K;
Ky7_x7 = b0*K;
*************************************
注記: (4) の定数 C の行列には、8 行および 12 列が含まれます。
aie_iir_1a.jl
は、次を実行するこのチュートリアルに含まれる Julia スクリプトです。
フィルター特性を定義します。
フィルターを 2 次セクションに分解します。
Double64 データ型を使用して、各ステージの SIMD 係数を生成します。
AI エンジンのベクター プロセッサは、1 サイクルで binary32 変数に対して 8 回の乗累算演算を実行できます。8 つの 独立した IIR 双 2 次フィルターで同時に処理される 8 つの独立した入力があれば、カーネル コードは簡単です。(1) から、これらの信号は理想的な 5 サイクルの待ち時間で処理されます。
サイクル 1: yi[n] = K * b0 * xi[n]
サイクル 2: yi[n] += K * b1 * xi[n-1]
サイクル 3: yi[n] += K * b2 * xi[n-2]
サイクル 4: yi[n] += -a1 * yi[n-1]
サイクル 5: yi[n] += -a2 * yi[n-2]
実際のインプリメンテーションでは、レジスタにデータを入力し、計算を実行し、パイプラインを通り、内部ステートを更新し、データを取り出すといった作業が必要になるため、レイテンシは 5 サイクル以上かかります。
(4) を使用すると、1 つの AI エンジン コアを使って、1 つの入力信号に対して 8 つのシーケンシャル出力を計算できます。計算を視覚化する 1 つの方法は、定数 C の行列の 1 列をクロック サイクルごとに処理することです。C は 12 列なので、理想的には 8 つの出力を生成するのに 12 サイクルかかります。
カーネル コード
1 つの 2 次ステージのカーネル コード (AI エンジン API を使用) は次のとおりです。
kernel.hpp
#ifndef __KERNEL_HPP__ // include guard to prevent multiple inclusion
#define __KERNEL_HPP__
#include <adf.h> // Adaptive DataFlow header
#include <aie_api/aie.hpp> // header files for AIE API
using Vector8f = aie::vector<float, 8>; // vector of 8 floating-point elements
using Vector16f = aie::vector<float, 16>; // vector of 16 floating-point elements
using VAcc8f = aie::accum<accfloat, 8>; // accumulator with 8 floating-point elements
template<unsigned id>
void SecondOrderSection(
adf::input_buffer<float> & __restrict idata, // 8 input samples per iteration
adf::output_buffer<float> & __restrict odata, // 8 output samples per iteration
const float (&C)[96] // run-time parameter: SIMD matrix of coefficients
);
#endif // __KERNEL_HPP__
kernel.cpp
#include <aie_api/aie_adf.hpp>
// !!! Need to separate C++20 portions for kernel from C++14-processed files (e.g. graph.hpp and tb.cpp)
#include "kernel.hpp"
template<unsigned id>
void SecondOrderSection(
adf::input_buffer<float> & __restrict idata, // 8 input samples per iteration
adf::output_buffer<float> & __restrict odata, // 8 output samples per iteration
const float (&C)[96] // run-time parameter: SIMD matrix of coefficients
) {
static Vector8f state_reg = aie::zeros<float, 8>(); // clear states
// input/output iterators
auto inIter = aie::begin_vector<8>(idata);
auto outIter = aie::begin_vector<8>(odata);
Vector8f xreg_hi = *inIter++; // fetch input samples
Vector16f xreg = aie::concat(state_reg, xreg_hi); // xreg[4]: ym2; xreg[5]: ym1; xreg[6]: xm2; xreg[7]: xm1; xreg[8:15]: x0:x7
Vector8f coeff = aie::load_v<8>(&C[0]);
VAcc8f acc = aie::mul(coeff, xreg[4]); // do 1st multiplication instead of zeroing
for (auto i = 1; i < 12; i++) {
coeff = aie::load_v<8>(&C[8 * i]);
float xval = xreg[i + 4];
acc = aie::mac(acc, coeff, xval);
}
Vector8f yout = acc;
// update states
state_reg[4] = yout[6];
state_reg[5] = yout[7];
state_reg[6] = xreg_hi[6];
state_reg[7] = xreg_hi[7];
*outIter++ = yout;
} // end SecondOrderSection()
注記:
カーネルのコードは C++20 準拠のコンパイラでコンパイルされ、それ以外のコード (
graph.hpp
とシミュレーション テストベンチ) は C++14 準拠のコンパイラでコンパイルされます。カーネル ヘッダーはgraph.hpp
に含まれるため、C++20 コンストラクトを含むことはできません。テンプレート パラメーター
id
は、SecondOrderSection()
関数の複数インスタンスをインスタンシエートするために使用されます。この関数は、
graph.hpp
で定義された数のエレメントを含む入力バッファーを受け取り、出力バッファーを生成します。入力と出力には、
__restrict
キーワードを使用すると、コンパイラの最適化が実行されます (詳細は UG1079 を参照)。フィルター係数は、
C
引数を介して 1 次元配列として渡されます。フィルター ステート (
state_reg
) は、関数呼び出しの間で保持する必要があるため、static
を宣言します。(4) で示されるような通常の行列とベクターの乗算の代わりに、
for
ループの各反復は C 行列の n 番目の列を使用し、その列のすべてのエレメントと x ベクターの n 番目のエレメントを乗算、つまりベクターのスケーリング演算を実行します。
Julia スクリプトに関する注記
カーネル コードの機能をチェックするため、aie_iir_1a.jl
を使用して、1 つの 2 次セクションの係数とインパルス応答を生成します。このスクリプトは、カーネルへの入力としてユニット サンプル関数も生成します。
Julia スクリプトのユーザー設定可能なパラメーターの一部を次に示します。
# --- begin user parameters
fp = 10.0e6 # passband frequency
fs = 100.0e6 # sampling frequency
p = 2 # no. of poles
rp = 0.1 # passband ripple (dB)
rs = 80.0 # stopband attenuation (dB)
N = 256 # no. of samples for impulse response
show_plots = true # display plots?
write_cmatrix = true # write C matrix to files?
write_impulse = true # write impulse response?
# --- end user parameters
スクリプトを実行するには、Julia を開始してから、プロンプトで次のコマンドを実行します。
julia> cd("path_to_aie_iir_1a.jl")
julia> include("aie_iir_1a.jl")
注記:
path_to_aie_iir_1a.jl
をaie_iir_1a.jl
スクリプトへの実際のパスに置き換えます。パスはダブルクォーテーションで囲む必要があります。
複数パッケージを読み込むため、最初の起動は遅く感じるかもしれません。
これにより次の 5 つのプロットが生成されます。
元のフィルターの周波数特性。
元のフィルターのインパルス応答。
SOS 周波数特性。
SOS インパルス応答。
インパルス応答誤差。
次のファイルも生成されます。
C1.h
- カーネルに渡す係数の配列。input.dat
- カーネルの入力信号として使用されるユニット サンプル関数。impresponse.dat
- AI エンジンの結果と比較するために計算されたインパルス応答。
C1.h
を AI エンジン プロジェクトのsrc
ディレクトリに、input.dat
とimpresponse.dat
をdata
ディレクトリにコピーします。
適応型データフロー グラフ
適応型データフロー グラフ ファイルは次のようになります。
graph.hpp
#ifndef __GRAPH_H__ // include guard to prevent multiple inclusion
#define __GRAPH_H__
#include <adf.h> // Adaptive DataFlow header
#include "kernel.hpp"
using namespace adf;
// dataflow graph declaration
class the_graph : public graph { // inherit all properties of the adaptive dataflow graph
private:
kernel section1;
public:
input_plio in; // input port for data to enter the kernel
input_port cmtx1; // input port for SIMD matrix coefficients
output_plio out; // output port for data to leave the kernel
// constructor
the_graph() {
// associate the kernel with the function to be executed
section1 = kernel::create(SecondOrderSection<1>);
// declare data widths and files for simulation
in = input_plio::create(plio_32_bits, "data/input.dat");
out = output_plio::create(plio_32_bits, "output.dat");
const unsigned num_samples = 8;
// establish connections
connect(in.out[0], section1.in[0]);
dimensions(section1.in[0]) = {num_samples};
connect<parameter>(cmtx1, adf::async(section1.in[1]));
connect(section1.out[0], out.in[0]);
dimensions(section1.out[0]) = {num_samples};
// specify which source code file contains the kernel function
source(section1) = "kernel.cpp";
// !!! temporary value: assumes this kernel dominates the AI Engine tile !!!
runtime<ratio>(section1) = 1.0;
} // end the_graph()
}; // end class the_graph
#endif // __GRAPH_H__
注記:
section1 = kernel::create(SecondOrderSection<1>)
は、AI エンジン カーネル プログラムが、テンプレート パラメーター1
を持つテンプレート関数SecondOrderSection
を使用することをツールに伝えます。[input|output]_plio
は、デバイスのプログラマブル ロジック (PL) 部分にあるシミュレーション用の入力または出力ポートを宣言します。[input|output]_plio::create()
は、PL で使用されるデータ バスの幅と、シミュレーション中に使用される関連する入力/出力ファイルを宣言します。connect(in.out[0], section1.in[0])
は、カーネル入力がin
ポートに接続されていることをツールに伝えます。dimensions(section1.in[0])
は、カーネルを実行する前に収集する必要があるサンプル数を宣言します。connect<parameter>(cmtx1, adf::async(section1.in[1]))
は、カーネルの最初の実行に非同期ランタイム パラメーターが必要であることをツールに伝えます。その後に続く式では使用可能な最新のランタイム パラメーターが使用されます。つまり、非同期パラメーターが一度しか送信されない場合、そのパラメーターはカーネルの残りの期間、再利用されます。connect(section1.out[0], out.in[0])
は、カーネル出力がout
ポートに接続されていることをツールに伝えます。dimensions(section1.out[0])
は、カーネルが各呼び出し中に生成するサンプル数を宣言します。source(section1) = "kernel.cpp"
は、カーネルのソース コードがどこにあるかをツールに伝えます。runtime<ratio>(section1) = 1.0
は、このカーネルのみを AI エンジンに配置できることをツールに伝えます。現時点では、実際の実行時間は不明です。
テストベンチ コード
テストベンチ コードは次のようになります。
tb.cpp
#include "kernel.hpp"
#include "graph.hpp"
#include "C1.h"
using namespace std;
using namespace adf;
// specify the DFG
the_graph my_graph;
const unsigned num_pts = 256; // number of sample points in "input.dat"
const unsigned num_iterations = num_pts/8; // number of iterations to run
// main simulation program
int main() {
my_graph.init(); // load the DFG into the AI Engine array, establish connectivity, etc.
my_graph.update(my_graph.cmtx1, C1, 96);
my_graph.run(num_iterations); // run the DFG for the specified number of iterations
my_graph.end(); // housekeeping
return (0);
} // end main()
注記:
my_graph.update(my_graph.cmtx1, C1, 96)
はフィルター用の 96 個の係数をカーネルに送信します。my_graph.run(num_iterations)
は、Julia スクリプトのインパルス応答と比較できるように、カーネルを 256/8 = 32 回実行します。
プログラムのビルドおよび実行
src
とdat
にあるファイルをプロジェクトにコピーします。Top-Level File
(値) をsrc/tb.cpp
(Rx および Tx 有効) に設定します。現打開では機能を検証するだけなので、プログラムのビルドと実行には
Emulation-SW
を使用します。プログラムがビルドされ、エラーなしに実行された場合、出力は
Emulation-SW/x86simulator_output/output.dat
となるはずです。生成された
impresponse.dat
ファイルをdata
ディレクトリにコピーします。Julia を使用してカーネル出力を検証できます。
julia> using PyPlot
julia> using DelimitedFiles
julia> ref = readdlm("{specify_directory}/impresponse.dat");
julia> dut = readdlm("{specify_directory}/output.dat");
julia> err = ref - dut;
julia> plot(err);
julia> grid("on");
julia> title("Impulse Response Error");
julia> xlabel("Sampling Index");
julia> ylabel("Error");
julia> eps(Float32)
1.1920929f-7
julia> maximum(abs.(err))
1.0517072768223557e-8
また、check.jl
を修正して実行することもできます。
$ julia
julia> include("check.jl")
# type Ctrl-D to exit Julia
その結果、インパルス応答誤差の Julia プロットは次のようになります。
最大絶対誤差は binary32 の機械イプシロン (Julia では Float32
) より小さいので、カーネル コードは問題なく動作しているといえます。
完全なデザインは、data
と src
ディレクトリに含まれています。AMD Vitis™ デザインをゼロから構築することに慣れていない場合は、aie_exp/Part1 チュートリアルを参照してください。
まとめ
デザインの機能性を検証するには、Emulation-SW
を使用します。IIR フィルターの 2 次セクションに必要な係数を生成する Julia スクリプトと、それをインプリメントするカーネル コードを提供しました。AI エンジンのカーネル コードのインプリメンテーションでは、元のフィルターとの違いはごくわずかです。
パート 1b では、カスケード接続された任意の数の 2 次セクションについて、適応型データフロー グラフを作成するプロセスを示します。
参考資料
[1] Julius O. Smith III 著。「オーディオ アプリケーションを使用したデジタル フィルター入門」
サポート
GitHub 問題は、リクエストやバグの追跡に使用します。質問については、forums.xilinx.com を参照してください。
Copyright © 2020-2024 Advanced Micro Devices, Inc