Main Content

このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。

castzeros を使用した浮動小数点型および固定小数点型の FIR フィルター アルゴリズムの実装

この例では、固定小数点型の指定をアルゴリズム コードから分離することによって、有限インパルス応答 (FIR) フィルターを固定小数点に変換する方法を説明します。

データ型の指定をアルゴリズム コードから分離すると、以下のことができます。

  • 異なるデータ型でアルゴリズム コードを再利用

  • データ型指定を切り離してアルゴリズムを整理し、データ型ごとにステートメントを切り替え

  • アルゴリズム コードの読みやすさを向上

  • 固定小数点と浮動小数点を切り替えてベースラインを比較

  • アルゴリズム コードを変更しないで固定小数点設定のバリエーションを切り替え

元の FIR フィルター アルゴリズム

この例では、有限インパルス応答 (FIR) フィルター用の MATLAB® コードを固定小数点に変換します。

フィルター係数が b、入力が x と仮定すると、FIR フィルターの n 番目の出力 y(n) の式は次のようになります。

y(n) = b(1)*x(n) + b(2)*x(n-1) + ... + b(end)*x(n-length(b)+1)

線形バッファーの実装

FIR フィルターを記述する方法はいくつかあります。そのうちの 1 つは、次の関数のように線形バッファーを使用する方法です。ここで、b は行ベクトル、zb と同じ長さの列ベクトルです。

function [y,z] = fir_filt_linear_buff(b,x,z)
    y = zeros(size(x));
    for n=1:length(x)
        z = [x(n); z(1:end-1)];
        y(n) = b * z;
    end
end

ベクトル z は、x の現在のサンプルと以前のサンプルで構成されます。

z = [x(n); x(n-1); .. ; x(n-length(b)+1)]

ベクトル z は状態ベクトルと呼ばれます。入力と出力の両方として使用されます。入力の場合は、フィルターの初期状態になります。出力の場合は、フィルターの最終状態になります。関数の外で状態ベクトルを定義することにより、リアルタイム システムでフィルターを使用して継続的にデータをストリーミングしたり、同じプロジェクト内でフィルター アルゴリズムを再利用したりできます。状態ベクトルは、Z 変換に関連する遅延 z-1 と同義であることから、よく z という呼び名で表現されます。Z 変換という名前は、この分野に大きな影響を与えた Lotfi Zadeh の名前から付けられたものです [1]。

線形バッファーの実装は、便利な MATLAB 行列構文に基づいており、読みやすさとわかりやすさを特長としています。ただし、この実装では、入力のサンプルごとに状態バッファーのフル コピーが導入されます。

リング バッファーの実装

FIR フィルターを効率的に実装するには、要素が z(p) = x(n) であるリング バッファー z に状態を格納します。ここで、n = 1,2,3,... の場合に p = mod(n-1,length(b))+1 です。

たとえば、length(b) = 3 として、pz をゼロに初期化します。

p = 0, z = [0 0 0]

最初のサンプルから開始し、循環しながら状態バッファー z を埋めます。

n = 1, p = 1, z(1) = x(1), z = [x(1)  0    0  ]
y(1) = b(1)*z(1) + b(2)*z(3) + b(3)*z(2)
n = 2, p = 2, z(2) = x(2), z = [x(1) x(2)  0  ]
y(2) = b(1)*z(2) + b(2)*z(1) + b(3)*z(3)
n = 3, p = 3, z(3) = x(3), z = [x(1) x(2) x(3)]
y(3) = b(1)*z(3) + b(2)*z(2) + b(3)*z(1)
n = 4, p = 1, z(1) = x(4), z = [x(4) x(2) x(3)]
y(4) = b(1)*z(1) + b(2)*z(3) + b(3)*z(2)
n = 5, p = 2, z(2) = x(5), z = [x(4) x(5) x(3)]
y(5) = b(1)*z(2) + b(2)*z(1) + b(3)*z(3)
n = 6, p = 3, z(3) = x(6), z = [x(4) x(5) x(6)]
y(6) = b(1)*z(3) + b(2)*z(2) + b(3)*z(1)
...

次の MATLAB 関数のようにリング バッファーを使用して、FIR フィルターを実装できます。

function [y,z,p] = fir_filt_circ_buff_original(b,x,z,p)
    y = zeros(size(x));
    nx = length(x);
    nb = length(b);
    for n=1:nx
        p=p+1; if p>nb, p=1; end
        z(p) = x(n);
        acc = 0;
        k = p;
        for j=1:nb
            acc = acc + b(j)*z(k);
            k=k-1; if k<1, k=nb; end
        end
        y(n) = acc;
    end
end

テスト ファイルの作成

固定小数点に変換する前に浮動小数点アルゴリズムが期待どおりに機能することを検証するテスト ファイルを作成します。これと同じテスト ファイルを使用して、固定小数点データ型の推奨を得て、変換後に固定小数点の結果を浮動小数点のベースラインと比較することができます。

テスト ベクトルは、システムによって想定されている値範囲の全体をカバーするような現実的な入力を表す必要があります。インパルス、正弦波の和およびチャープ信号を現実的な入力として使用すると、これらに対して線形理論を使用して出力が正しいことを検証できます。最大出力を生成する信号は、システムがオーバーフローしないことを検証するうえで役立ちます。

設定

グローバルな固定小数点数学設定と固定小数点基本設定の現在の状態を取得してリセットするには、次のコードを実行します。

resetglobalfimath;
FIPREF_STATE = get(fipref);
resetfipref;

フィルター係数

Signal Processing Toolbox™ の関数 fir1 を使用して計算された、以下のローパス フィルター係数を使用します。

b = fir1(11,0.25);
b = [-0.004465461051254
    -0.004324228005260
    +0.012676739550326
    +0.074351188907780
    +0.172173206073645
    +0.249588554524763
    +0.249588554524763
    +0.172173206073645
    +0.074351188907780
    +0.012676739550326
    -0.004324228005260
    -0.004465461051254]';

時間ベクトル

テスト信号を作成するには、次の時間ベクトルを使用します。

nx = 256;
t = linspace(0,10*pi,nx)';

インパルス入力

インパルス入力に対する FIR フィルターの応答はフィルター係数そのものになります。

x_impulse = zeros(nx,1); x_impulse(1) = 1;

最大出力を生成する信号

フィルターの最大出力が発生するのは、入力の符号がフィルターのインパルス応答の符号と一致する場合です。

x_max_output = sign(fliplr(b))';
x_max_output = repmat(x_max_output,ceil(nx/length(b)),1);
x_max_output = x_max_output(1:nx);

出力の最大振幅はそのインパルス応答の 1 ノルムです。ここで、norm(b,1) = sum(abs(b)) です。

maximum_output_magnitude = norm(b,1) %#ok<*NOPTS>
maximum_output_magnitude = 1.0352

正弦波の和

正弦波の和はフィルターの代表的な入力です。プロットでフィルター処理による高周波数の遮断を容易に確認できます。

f0 = 0.1; f1 = 2;
x_sines = sin(2*pi*t*f0) + 0.1*sin(2*pi*t*f1);

チャープ

チャープによって、低周波数を通過させて高周波数を減衰させるローパス フィルター動作が視覚的にわかりやすくなります。

f_chirp = 1/16;                  % Target frequency
x_chirp = sin(pi*f_chirp*t.^2);  % Linear chirp

titles = {'Impulse', 'Max output', 'Sum of sines', 'Chirp'};
x = [x_impulse, x_max_output, x_sines, x_chirp];

元の関数の呼び出し

固定小数点への変換を開始する前に、テスト ファイルの入力を使用して元の関数を呼び出し、後続の出力と比較するベースラインを確立します。

y0 = zeros(size(x));
for i=1:size(x,2)
    % Initialize the states for each column of input
    p = 0;
    z = zeros(size(b));
    y0(:,i) = fir_filt_circ_buff_original(b,x(:,i),z,p);
end

ベースライン出力

fir_filt_circ_buff_plot(1,titles,t,x,y0)

Figure contains 4 axes objects. Axes object 1 with title Impulse, xlabel t contains 2 objects of type line. These objects represent Input, Baseline Output. Axes object 2 with title Max output, xlabel t contains 2 objects of type line. These objects represent Input, Baseline Output. Axes object 3 with title Sum of sines, xlabel t contains 2 objects of type line. These objects represent Input, Baseline Output. Axes object 4 with title Chirp, xlabel t contains 2 objects of type line. These objects represent Input, Baseline Output.

インストルメンテーションとコード生成の準備

アルゴリズムが MATLAB で機能するようになった後で最初に行う手順は、インストルメンテーションの準備をすることであり、これにはコード生成が必要となります。変換の前に関数 coder.screener を使用してコードを解析し、サポートされていない関数や言語機能を特定できます。

エントリポイント関数

インストルメンテーションとコード生成を行う場合は、固定小数点に変換する関数を呼び出すエントリポイント関数があると便利です。FIR フィルターの入力をさまざまなデータ型にキャストし、さまざまなバリエーションのフィルターに呼び出しを追加して比較できます。エントリポイント関数を使用すると、固定小数点と浮動小数点の両方のバリアントのフィルターを実行するだけでなく、固定小数点のさまざまなバリアントを実行することもできます。これにより、コードをより高速に反復して、最適な固定小数点設計を得ることができます。

function y = fir_filt_circ_buff_original_entry_point(b,x,reset)
    if nargin<3, reset = true; end
    % Define the circular buffer z and buffer position index p.
    % They are declared persistent so the filter can be called in a streaming
    % loop, each section picking up where the last section left off.
    persistent z p
    if isempty(z) || reset
        p = 0;
        z = zeros(size(b));
    end
    [y,z,p] = fir_filt_circ_buff_original(b,x,z,p);
end

テスト ファイル

テスト ファイルは、コンパイルされたエントリポイント関数を呼び出します。

function y = fir_filt_circ_buff_test(b,x)
     y = zeros(size(x));
     for i=1:size(x,2)
         reset = true;
         y(:,i) = fir_filt_circ_buff_original_entry_point_mex(b,x(:,i),reset);
     end
end

元の関数のビルド

buildInstrumentedMex を使用して元のエントリポイント関数をコンパイルします。これによりコードがログ作成用にインストルメント化され、シミュレーションから最小値と最大値を収集して、データ型の推奨を得ることができます。

reset = true;
buildInstrumentedMex fir_filt_circ_buff_original_entry_point -args {b, x(:,1), reset}

元の関数の実行

最大値と最小値のログを作成するアルゴリズムを使用して、テスト ファイルの入力を実行します。

y1 = fir_filt_circ_buff_test(b,x);

型の表示

テスト ファイルの実行時にログが作成された最小値と最大値およびすべての変数のデータ型を表示するには、showInstrumentationResults を使用します。出力変数 y とアキュムレータ変数 acc についてログが作成された最大値を調べると、以前に計算した理論上の最大出力値に達していることがわかります。

showInstrumentationResults fir_filt_circ_buff_original_entry_point_mex

インストルメント化された [コード生成レポート] でこれらの結果を確認する方法は以下です。

  • 関数 fir_filt_circ_buff_original を選択します。

  • [変数] タブを選択します。

元の関数の検証

関数を変更するたびに、関数の結果がベースラインと引き続き一致することを検証します。

fir_filt_circ_buff_plot2(2,titles,t,x,y0,y1)

Figure contains 8 axes objects. Axes object 1 with title Impulse, xlabel t contains 2 objects of type line. Axes object 2 with title Impulse error, xlabel t contains an object of type line. Axes object 3 with title Max output, xlabel t contains 2 objects of type line. Axes object 4 with title Max output error, xlabel t contains an object of type line. Axes object 5 with title Sum of sines, xlabel t contains 2 objects of type line. Axes object 6 with title Sum of sines error, xlabel t contains an object of type line. Axes object 7 with title Chirp, xlabel t contains 2 objects of type line. Axes object 8 with title Chirp error, xlabel t contains an object of type line.

型テーブルを使用するように関数を変換

アルゴリズムからデータ型を分離する方法は以下です。

  1. データ型定義のテーブルを作成します。

  2. そのテーブルのデータ型を使用するようアルゴリズム コードを変更します。

この例では、さまざまなファイルを作成することによって反復的な手順について説明します。実際には、同じファイルに反復的な変更を加えることができます。

元の型テーブル

変数のプロトタイプが元の型に設定された構造体を使用して、型テーブルを作成します。ベースライン型を使用して、初回の変換が正しく行われたことを検証すると共に、プログラムによって浮動小数点と固定小数点の間で関数を切り替えます。インデックス変数 jknnbnx は MATLAB Coder™ によって自動的に整数に変換されるため、この table で型を指定する必要はありません。

プロトタイプ値を空白 ([ ]) として指定します。これはこの値ではなくデータ型を使用するためです。

function T = fir_filt_circ_buff_original_types()
    T.acc=double([]);
    T.b=double([]);
    T.p=double([]);
    T.x=double([]);
    T.y=double([]);
    T.z=double([]);
end

型認識フィルター関数

関数 castzeros および型テーブルを使用して、フィルター関数とエントリポイント関数で型認識できるようにします。

添字を使用した代入 acc(:)=...p(:)=1 および k(:)=nb を使用して、代入時にデータ型を維持します。添字を使用した代入とデータ型の維持の詳細については、fi オブジェクトのキャストを参照してください。

関数呼び出し y = zeros(size(x),'like',T.y) は、変数 T.y のプロパティを使用して、x と同じサイズのゼロの配列を作成します。T.y は初期状態では関数 fir_filt_circ_buff_original_types で定義された double ですが、この例では後で固定小数点型として再定義します。

関数呼び出し acc = cast(0,'like',T.acc) は、変数 T.acc と同じプロパティを使用して値 0 をキャストします。T.acc は初期状態では関数 fir_filt_circ_buff_original_types で定義された double ですが、この例では後で固定小数点型として再定義します。

function [y,z,p] = fir_filt_circ_buff_typed(b,x,z,p,T)
    y = zeros(size(x),'like',T.y);
    nx = length(x);
    nb = length(b);
    for n=1:nx
        p(:)=p+1; if p>nb, p(:)=1; end
        z(p) = x(n);
        acc = cast(0,'like',T.acc);
        k = p;
        for j=1:nb
            acc(:) = acc + b(j)*z(k);
            k(:)=k-1; if k<1, k(:)=nb; end
        end
        y(n) = acc;
    end
end

型認識エントリポイント関数

関数呼び出し p1 = cast(0,'like',T1.p) は、変数 T1.p と同じプロパティを使用して値 0 をキャストします。T1.p は初期状態では関数 fir_filt_circ_buff_original_types で定義された double ですが、この例では後で整数型として再定義します。

関数呼び出し z1 = zeros(size(b),'like',T1.z) は、変数 T1.z のプロパティを使用して、b と同じサイズのゼロの配列を作成します。T1.z は初期状態では関数 fir_filt_circ_buff_original_types で定義された double ですが、この例では後で固定小数点型として再定義します。

function y1 = fir_filt_circ_buff_typed_entry_point(b,x,reset)
   if nargin<3, reset = true; end
   %
   % Baseline types
   %
   T1 = fir_filt_circ_buff_original_types();
   % Each call to the filter needs to maintain its own states.
   persistent z1 p1
   if isempty(z1) || reset
       p1 = cast(0,'like',T1.p);
       z1 = zeros(size(b),'like',T1.z);
   end
   b1 = cast(b,'like',T1.b);
   x1 = cast(x,'like',T1.x);
   [y1,z1,p1] = fir_filt_circ_buff_typed(b1,x1,z1,p1,T1);
end

変更後の関数の検証

関数を変更するたびに、関数の結果がベースラインと引き続き一致することを検証します。型テーブルで元の型を使用したため、出力は同一になるはずです。これにより、アルゴリズムから型を正しく分離するための変換が行われたことが検証されます。

buildInstrumentedMex fir_filt_circ_buff_typed_entry_point -args {b, x(:,1), reset}

y1 = fir_filt_circ_buff_typed_test(b,x);
fir_filt_circ_buff_plot2(3,titles,t,x,y0,y1)

Figure contains 8 axes objects. Axes object 1 with title Impulse, xlabel t contains 2 objects of type line. Axes object 2 with title Impulse error, xlabel t contains an object of type line. Axes object 3 with title Max output, xlabel t contains 2 objects of type line. Axes object 4 with title Max output error, xlabel t contains an object of type line. Axes object 5 with title Sum of sines, xlabel t contains 2 objects of type line. Axes object 6 with title Sum of sines error, xlabel t contains an object of type line. Axes object 7 with title Chirp, xlabel t contains 2 objects of type line. Axes object 8 with title Chirp error, xlabel t contains an object of type line.

シミュレーション最小値/最大値ログからのデータ型の推奨

既定の符号付き固定小数点型と 16 ビットの語長を仮定し、関数 showInstrumentationResults を使用して固定小数点の小数部の長さを推奨します。

showInstrumentationResults fir_filt_circ_buff_original_entry_point_mex ...
    -defaultDT numerictype(1,16) -proposeFL

インストルメント化された [コード生成レポート] で、関数 fir_filt_circ_buff_original[変数] タブを選択して以下の結果を表示します。

固定小数点型テーブルの作成

[コード生成レポート] から推奨された型を参考にして固定小数点型を選択し、変数のプロトタイプをもつ構造体を使用して、固定小数点型テーブルを作成します。

アルゴリズムに関する知識に基づいて、推奨された型を改善します。たとえば、変数 acc をアキュムレータとして使用する場合は、32 ビットを指定します。[コード生成レポート] を見ると、オーバーフローを防ぐために少なくとも 2 整数ビットが acc に必要なことがわかるため、小数部の長さを 30 に設定します。

変数 p はインデックスとして使用されるため、組み込み関数を使って 16 ビット整数にすることができます。

プロトタイプ値を空白 ([ ]) として指定します。これはこの値ではなくデータ型を使用するためです。

function T = fir_filt_circ_buff_fixed_point_types()
    T.acc=fi([],true,32,30);
    T.b=fi([],true,16,17);
    T.p=int16([]);
    T.x=fi([],true,16,14);
    T.y=fi([],true,16,14);
    T.z=fi([],true,16,14);
end

固定小数点をエントリポイント関数に追加

固定小数点型テーブルへの呼び出しをエントリポイント関数に追加します。

T2 = fir_filt_circ_buff_fixed_point_types();
persistent z2 p2
if isempty(z2) || reset
    p2 = cast(0,'like',T2.p);
    z2 = zeros(size(b),'like',T2.z);
end
b2 = cast(b,'like',T2.b);
x2 = cast(x,'like',T2.x);
[y2,z2,p2] = fir_filt_circ_buff_typed(b2,x2,z2,p2,T2);

固定小数点データ型を使用したアルゴリズムのビルドと実行

buildInstrumentedMex fir_filt_circ_buff_typed_entry_point -args {b, x(:,1), reset}

[y1,y2] = fir_filt_circ_buff_typed_test(b,x);

showInstrumentationResults fir_filt_circ_buff_typed_entry_point_mex

インストルメント化された [コード生成レポート] でこれらの結果を確認する方法は以下です。

  • エントリポイント関数 fir_filt_circ_buff_typed_entry_point を選択します。

  • 以下のコード行で fir_filt_circ_buff_typed を選択します。

[y2,z2,p2] = fir_filt_circ_buff_typed(b2,x2,z2,p2,T2);
  • [変数] タブを選択します。

16 ビット語長の完全精度演算

結果がベースラインの許容誤差内に収まっていることを検証します。

fir_filt_circ_buff_plot2(4,titles,t,x,y1,y2);

Figure contains 8 axes objects. Axes object 1 with title Impulse, xlabel t contains 2 objects of type line. Axes object 2 with title Impulse error, xlabel t contains an object of type line. Axes object 3 with title Max output, xlabel t contains 2 objects of type line. Axes object 4 with title Max output error, xlabel t contains an object of type line. Axes object 5 with title Sum of sines, xlabel t contains 2 objects of type line. Axes object 6 with title Sum of sines error, xlabel t contains an object of type line. Axes object 7 with title Chirp, xlabel t contains 2 objects of type line. Axes object 8 with title Chirp error, xlabel t contains an object of type line.

アルゴリズムが固定小数点 MATLAB コードに変換されました。C コードに変換する必要がある場合は、次の節に進みます。

C コードの生成

この節では、前の節で使用した固定小数点 MATLAB コードから効率的な C コードを生成する方法について説明します。

必要な製品

C コードの生成には MATLAB Coder が必要であり、この例で使用するハードウェア実行設定には Embedded Coder® が必要となります。

C コードの効率性を最大限に高めるよう調整されたアルゴリズム

出力変数 y はゼロに初期化された後、完全に上書きされてから使用されます。したがって、y にすべてゼロを代入することは不要です。関数 coder.nullcopy を使用すると、実際に値を代入しないで変数を宣言できます。これにより、この場合はコードが効率化します。ただし、変数の要素を割り当てる前にその要素にアクセスすると、初期化されていないメモリにアクセスすることになり、メモリの内容は予測できないため、coder.nullcopy は慎重に使用する必要があります。

経験則として、coder.nullcopy の使用が適しているのは、アルゴリズムの他の部分と比べて初期化に相当な時間がかかる場合です。使用すべきかどうかが定かでない場合は、使用しないことが確実に安全な選択です。

function [y,z,p] = fir_filt_circ_buff_typed_codegen(b,x,z,p,T)
    % Use coder.nullcopy only when you are certain that every value of
    % the variable is overwritten before it is used.
    y = coder.nullcopy(zeros(size(x),'like',T.y));
    nx = length(x);
    nb = length(b);
    for n=1:nx
        p(:)=p+1; if p>nb, p(:)=1; end
        z(p) = x(n);
        acc = cast(0,'like',T.acc);
        k = p;
        for j=1:nb
            acc(:) = acc + b(j)*z(k);
            k(:)=k-1; if k<1, k(:)=nb; end
        end
        y(n) = acc;
    end
end

ネイティブ C コード型

C のネイティブ アクションと一致するように固定小数点演算プロパティを設定できます。これにより、効率性を最大限に高めた C コードを生成できますが、この例に示すとおり、オーバーフローによって問題が生じ、不正確な結果が生成される可能性があります (この問題の修正方法は次の節で示します)。ただし、この問題は必ず生じるわけではないため、クリーンな C コードを生成できるかどうかまず試してみる価値はあります。

C の既定アクションである負方向の丸めとオーバーフローのラップを使用するように、固定小数点演算プロパティを設定します。

ネイティブ C 32 ビット整数型と一致し、数学演算の最下位ビット (LSB) を維持するように、積と和の固定小数点演算プロパティを設定します。

これらの設定を固定小数点型テーブルに追加します。

function T = fir_filt_circ_buff_dsp_types()
    F = fimath('RoundingMethod','Floor',...
               'OverflowAction','Wrap',...
               'ProductMode','KeepLSB',...
               'ProductWordLength',32,...
               'SumMode','KeepLSB',...
               'SumWordLength',32);
    T.acc=fi([],true,32,30,F);
    T.p=int16([]);
    T.b=fi([],true,16,17,F);
    T.x=fi([],true,16,14,F);
    T.y=fi([],true,16,14,F);
    T.z=fi([],true,16,14,F);
end

ネイティブ C コード型のテスト

型テーブルへの呼び出しをエントリポイント関数に追加して、テスト ファイルを実行します。

[y1,y2,y3] = fir_filt_circ_buff_typed_test(b,x); %#ok<*ASGLU>

プロットの 2 行目を見てわかるとおり、最大出力誤差が入力のサイズの 2 倍になっており、本来は正であるはずの値が負にオーバーフローしていることが示されています。また、他の出力はオーバーフローしていないこともわかります。テスト ファイルでは代表的な入力だけでなく、値範囲の全体をカバーすることが重要なのはこのためです。

fir_filt_circ_buff_plot2(5,titles,t,x,y1,y3);

Figure contains 8 axes objects. Axes object 1 with title Impulse, xlabel t contains 2 objects of type line. Axes object 2 with title Impulse error, xlabel t contains an object of type line. Axes object 3 with title Max output, xlabel t contains 2 objects of type line. Axes object 4 with title Max output error, xlabel t contains an object of type line. Axes object 5 with title Sum of sines, xlabel t contains 2 objects of type line. Axes object 6 with title Sum of sines error, xlabel t contains an object of type line. Axes object 7 with title Chirp, xlabel t contains 2 objects of type line. Axes object 8 with title Chirp error, xlabel t contains an object of type line.

スケーリングされた double 型を使用したオーバーフローの検索

スケーリングされた倍精度変数は、倍精度浮動小数点のデータを格納するため、演算を全範囲で実行します。また、固定小数点設定も保持するため、計算が固定小数点型の範囲から外れた場合に報告することが可能です。

スケーリングされた double にデータ型を変更して、スケーリングされた double 型テーブルにこれらの設定を追加します。

function T = fir_filt_circ_buff_scaled_double_types()
    F = fimath('RoundingMethod','Floor',...
               'OverflowAction','Wrap',...
               'ProductMode','KeepLSB',...
               'ProductWordLength',32,...
               'SumMode','KeepLSB',...
               'SumWordLength',32);
    DT = 'ScaledDouble';
    T.acc=fi([],true,32,30,F,'DataType',DT);
    T.p=int16([]);
    T.b=fi([],true,16,17,F,'DataType',DT);
    T.x=fi([],true,16,14,F,'DataType',DT);
    T.y=fi([],true,16,14,F,'DataType',DT);
    T.z=fi([],true,16,14,F,'DataType',DT);
end

スケーリングされた double 型 table への呼び出しをエントリポイント関数に追加して、テスト ファイルを実行します。

[y1,y2,y3,y4] = fir_filt_circ_buff_typed_test(b,x); %#ok<*NASGU>

スケーリングされた double 型を使用してインストルメンテーションの結果を表示します。

showInstrumentationResults fir_filt_circ_buff_typed_entry_point_mex

インストルメント化された [コード生成レポート] でこれらの結果を確認する方法は以下です。

  • エントリポイント関数 fir_filt_circ_buff_typed_entry_point を選択します。

  • 以下のコード行で fir_filt_circ_buff_typed_codegen を選択します。

[y4,z4,p4] = fir_filt_circ_buff_typed_codegen(b4,x4,z4,p4,T4);
  • [変数] タブを選択します。

  • テーブルの変数を確認します。いずれの変数もオーバーフローしていないため、このオーバーフローは演算の結果として発生したことがわかります。

  • レポートで演算子 (+-*=) の上にマウス ポインターを置きます。

  • インストルメント化された [コード生成レポート] で、次の MATLAB コード行の + の上にマウス ポインターを置きます。

acc(:) = acc + b(j)*z(k);

和がオーバーフローしたことがレポートに次のように示されます。

和がオーバーフローした理由は、bnumerictype(true,16,17) で、znumerictype(true,16,14) であるため、b(j)*z(k) の完全精度の積によって numerictype(true,32,31) が出力されるためです。和の型は "最下位ビットを保持" (KeepLSB) に設定されているため、和は numerictype(true,32,31) になります。ただし、シミュレーションされた最小値と最大値である -1.0045+1.035 をそれぞれ格納するには、2 整数ビットが必要となります。

オーバーフローが発生しないようにするための調整

b(j)*z(k)numerictype(true,32,30) になるように、b の小数部の長さを 17 ではなく 16 に設定します。これにより、和の KeepLSB ルールに従って、この和も numerictype(true,32,30) になります。

その他の設定はすべてそのまま保持して、以下を設定します。

T.b=fi([],true,16,16,F);

これにより、この MATLAB コード行の和がオーバーフローしなくなります。

acc(:) = acc + b(j)*z(k);

新しい設定を使用してテスト ファイルを実行し、結果をプロットします。

[y1,y2,y3,y4,y5] = fir_filt_circ_buff_typed_test(b,x);

オーバーフローが発生していないことがわかります。ただし、C の固有の丸め処理である負方向の丸めを使用したことによってバイアスが生じ、誤差が大きくなっていることがプロットに示されます。このバイアスが許容可能な場合は、これ以上何も行わなくても、非常にクリーンな C コードが生成されているはずです。

fir_filt_circ_buff_plot2(6,titles,t,x,y1,y5);

Figure contains 8 axes objects. Axes object 1 with title Impulse, xlabel t contains 2 objects of type line. Axes object 2 with title Impulse error, xlabel t contains an object of type line. Axes object 3 with title Max output, xlabel t contains 2 objects of type line. Axes object 4 with title Max output error, xlabel t contains an object of type line. Axes object 5 with title Sum of sines, xlabel t contains 2 objects of type line. Axes object 6 with title Sum of sines error, xlabel t contains an object of type line. Axes object 7 with title Chirp, xlabel t contains 2 objects of type line. Axes object 8 with title Chirp error, xlabel t contains an object of type line.

バイアスの除去

生じたバイアスがアプリケーションの許容範囲を超えている場合は、丸め手法を 'Nearest' に変更してバイアスを除去します。最も近い正の整数方向に丸めると、多少複雑な C コードが生成されますが、バイアスを除去して誤差を小さくするためにはこのようにしなければならない場合があります。

最も近い正の整数方向への丸めを行い、係数の小数部の長さを調整した固定小数点型テーブルは、最終的に以下のようになります。

function T = fir_filt_circ_buff_dsp_nearest_types()
    F = fimath('RoundingMethod','Nearest',...
               'OverflowAction','Wrap',...
               'ProductMode','KeepLSB',...
               'ProductWordLength',32,...
               'SumMode','KeepLSB',...
               'SumWordLength',32);
    T.acc=fi([],true,32,30,F);
    T.p=int16([]);
    T.b=fi([],true,16,16,F);
    T.x=fi([],true,16,14,F);
    T.y=fi([],true,16,14,F);
    T.z=fi([],true,16,14,F);
end

エントリポイント関数からこの型テーブルを呼び出して実行し、出力をプロットします。

[y1,y2,y3,y4,y5,y6] = fir_filt_circ_buff_typed_test(b,x);
fir_filt_circ_buff_plot2(7,titles,t,x,y1,y6);

Figure contains 8 axes objects. Axes object 1 with title Impulse, xlabel t contains 2 objects of type line. Axes object 2 with title Impulse error, xlabel t contains an object of type line. Axes object 3 with title Max output, xlabel t contains 2 objects of type line. Axes object 4 with title Max output error, xlabel t contains an object of type line. Axes object 5 with title Sum of sines, xlabel t contains 2 objects of type line. Axes object 6 with title Sum of sines error, xlabel t contains an object of type line. Axes object 7 with title Chirp, xlabel t contains 2 objects of type line. Axes object 8 with title Chirp error, xlabel t contains an object of type line.

コード生成コマンドの実行

C コードを生成するには、次のビルド関数を実行します。ベスト プラクティスとして、エントリポイント関数やテスト ファイルがなくてもコア アルゴリズムの C コードを生成できるようにビルド関数を作成しておくと、より大きなプロジェクトにもコア アルゴリズムの C コードを含められるようになります。

function fir_filt_circ_buff_build_function()
    %
    % Declare input arguments
    %
    T = fir_filt_circ_buff_dsp_nearest_types();
    b = zeros(1,12,'like',T.b);
    x = zeros(256,1,'like',T.x);
    z = zeros(size(b),'like',T.z);
    p = cast(0,'like',T.p);
    %
    % Code generation configuration
    %
    h = coder.config('lib');
    h.PurelyIntegerCode = true;
    h.SaturateOnIntegerOverflow = false;
    h.SupportNonFinite = false;
    h.HardwareImplementation.ProdBitPerShort = 8;
    h.HardwareImplementation.ProdBitPerInt = 16;
    h.HardwareImplementation.ProdBitPerLong = 32;
    %
    % Generate C-code
    %
    codegen fir_filt_circ_buff_typed_codegen -args {b,x,z,p,T} -config h -launchreport
end

生成した C コード

MATLAB Coder はこれらの設定を使用して以下の C コードを生成します。

void fir_filt_circ_buff_typed_codegen(const int16_T b[12], const int16_T x[256],
  int16_T z[12], int16_T *p, int16_T y[256])
{
  int16_T n;
  int32_T acc;
  int16_T k;
  int16_T j;
  for (n = 0; n < 256; n++) {
    (*p)++;
    if (*p > 12) {
      *p = 1;
    }
    z[*p - 1] = x[n];
    acc = 0L;
    k = *p;
    for (j = 0; j < 12; j++) {
      acc += (int32_T)b[j] * z[k - 1];
      k--;
      if (k < 1) {
        k = 12;
      }
    }
    y[n] = (int16_T)((acc >> 16) + ((acc & 32768L) != 0L));
  }
}

以下のコードを実行してグローバル状態を復元します。

fipref(FIPREF_STATE);
clearInstrumentationResults fir_filt_circ_buff_original_entry_point_mex
clearInstrumentationResults fir_filt_circ_buff_typed_entry_point_mex
clear fir_filt_circ_buff_original_entry_point_mex
clear fir_filt_circ_buff_typed_entry_point_mex

参考文献:

[1] J. R. Ragazzini and L. A. Zadeh. "The analysis of sampled-data systems". In:Transactions of the American Institute of Electrical Engineers 71 (II) (1952), pp. 225–234.