固定小数点アルゴリズムの開発
この例では、簡単な固定小数点アルゴリズムの開発方法と検証方法を説明します。この例は、以下のアルゴリズム開発の手順に沿っています。
1) 2 次フィルター アルゴリズムを実装し、倍精度浮動小数点でシミュレートします。
2) コードをインストルメント化し、出力と状態のダイナミック レンジを可視化します。
3) 変数のデータ型を変更することにより、アルゴリズムを固定小数点に変換します。アルゴリズム自体は変わりません。
4) 固定小数点と浮動小数点の結果を比較し、プロットします。
倍精度浮動小数点変数の定義
倍精度浮動小数点でアルゴリズムを開発します。この例で使用するアルゴリズムは、入力信号の高周波数を除去する 2 次ローパス フィルターです。
分子係数 a
と分母係数 b
を定義します。
b = [ 0.25 0.5 0.25 ]; a = [ 1 0.09375 0.28125 ];
高周波数と低周波数の両方をもつ乱数入力を生成します。
s = rng; rng(0,'v5uniform'); x = randn(1000,1); rng(s); % restore |rng| state
高速化するために、出力 y
と状態 z
を事前に割り当てます。
y = zeros(size(x)); z = [0;0];
データ型に依存しないアルゴリズムの実装
これは、標準の差分方程式を実装する 2 次フィルターです。
|y(n) = b(1)*x(n) + b(2)*x(n-1) + b(3)*x(n-2) - a(2)*y(n-1) - a(3)*y(n-2)|
for k=1:length(x) y(k) = b(1)*x(k) + z(1); z(1) = (b(2)*x(k) + z(2)) - a(2)*y(k); z(2) = b(3)*x(k) - a(3)*y(k); end
浮動小数点の結果を保存します。
ydouble = y;
浮動小数点コードのインストルメントによるダイナミック レンジの可視化
固定小数点に変換するには、変数の範囲を知る必要があります。アルゴリズムの複雑さによって、このタスクは簡単な場合もあれば非常に難しい場合もあります。この例では、入力値の範囲がわかっているため、適切な固定小数点データ型を選択することは簡単です。出力 (y
) と状態 (z
) の範囲が不明であるため、この 2 つに注目します。出力と状態のダイナミック レンジを表示するには、コードを少し変更してインストルメントします。
NumericTypeScope
オブジェクトを作成し、状態 (z
) のダイナミック レンジを表示します。
z = [0;0]; % Reset states hscope1 = NumericTypeScope; for k=1:length(x) y(k) = b(1)*x(k) + z(1); z(1) = (b(2)*x(k) + z(2)) - a(2)*y(k); z(2) = b(3)*x(k) - a(3)*y(k); % Process the data and update the visual step(hscope1,z); end % clear the information stored in the object hscope1 reset(hscope1); % Create a |NumericTypeScope| object and view the dynamic range of the % output (|y|) hscope2 = NumericTypeScope; step(hscope2,y);
まず、変数 z
(状態) について表示されている情報を解析します。ヒストグラムから、ダイナミック レンジが (2^(1) 2^(-12)]
内にあることがわかります。
既定の設定では、スコープは 16 ビットの語長を使用し、オーバーフローを許容しません。これにより、オーバーフローを回避するために少なくとも 2 整数ビットを必要とするため、データ型は numerictype(true,16,14)
となります。[入力データ] および [結果のタイプ] パネルには、統計データについての詳細情報が表示されます。[入力データ] パネルでは、データに正の値と負の値の両方が含まれており、符号付きの数量であること、そしてこれらが提案された numerictype
に反映されていることがわかります。また、データの最大値は 1.51 であり、これは提案されたデータ型で表現できることもわかります。
次に、変数 y
(出力) を確認します。ヒストグラム プロットから、ダイナミック レンジが (2^(1) 2^(-13)]
内にあることがわかります。
既定の設定では、スコープは 16 ビットの語長を使用し、オーバーフローを許容しません。これにより、オーバーフローを回避するために少なくとも 2 整数ビットを必要とするため、データ型は numerictype(true,16,14)
となります。提案されたこのデータ型を使用すると、オーバーフローまたはアンダーフローが表示されません。
固定小数点変数の定義
固定小数点データ型を使用するように変数を変換し、アルゴリズムを再度実行します。
ログを有効にし、選択したデータ型によるオーバーフローとアンダーフローを確認します。
FIPREF_STATE = get(fipref);
reset(fipref)
fp = fipref;
default_loggingmode = fp.LoggingMode;
fp.LoggingMode = 'On';
グローバル fimath の現在の状態を取得し、出荷時の設定にリセットします。
globalFimathAtStart = fimath; resetglobalfimath;
変数の固定小数点型を fi(Data, Signed, WordLength, FractionLength)
という形式で定義します。
b = fi(b, 1, 8, 6); a = fi(a, 1, 8, 6); x = fi(x, 1, 16, 13); y = fi(zeros(size(x)), 1, 16, 13); z = fi([0;0], 1, 16, 14);
データ型に依存しない 1 つのアルゴリズムの実装
for k=1:length(x) y(k) = b(1)*x(k) + z(1); z(1) = (b(2)*x(k) + z(2)) - a(2)*y(k); z(2) = b(3)*x(k) - a(3)*y(k); end % Reset the logging mode. fp.LoggingMode = default_loggingmode;
この例では、アルゴリズム コードをインライン化してわかりやすくするために、固定小数点変数を浮動小数点と同じ名前で再定義しています。ただし、このアルゴリズム コードを、浮動小数点変数と固定小数点変数のどちらでも呼び出せる MATLAB® ファイル関数に含めることをお勧めします。
浮動小数点と固定小数点の結果の比較とプロット
浮動小数点と固定小数点の結果の振幅応答とフィルターの応答をプロットし、固定小数点への変換時にフィルターが期待どおりの動作をしているかどうかを確認します。
n = length(x); f = linspace(0,0.5,n/2); x_response = 20*log10(abs(fft(double(x)))); ydouble_response = 20*log10(abs(fft(ydouble))); y_response = 20*log10(abs(fft(double(y)))); plot(f,x_response(1:n/2),'c-',... f,ydouble_response(1:n/2),'bo-',... f,y_response(1:n/2),'gs-'); ylabel('Magnitude in dB'); xlabel('Normalized Frequency'); legend('Input','Floating point output','Fixed point output','Location','Best'); title('Magnitude response of Floating-point and Fixed-point results');
h = fft(double(b),n)./fft(double(a),n); h = h(1:end/2); clf hax = axes; plot(hax,f,20*log10(abs(h))); set(hax,'YLim',[-40 0]); title('Magnitude response of the filter'); ylabel('Magnitude in dB') xlabel('Frequency');
入力信号の高周波数はローパス フィルターで減衰されており、これは期待どおりの動作です。
誤差のプロット
clf n = (0:length(y)-1)'; e = double(lsb(y)); plot(n,double(y)-ydouble,'.-r', ... [n(1) n(end)],[e/2 e/2],'c', ... [n(1) n(end)],[-e/2 -e/2],'c') text(n(end),e/2,'+1/2 LSB','HorizontalAlignment','right','VerticalAlignment','bottom') text(n(end),-e/2,'-1/2 LSB','HorizontalAlignment','right','VerticalAlignment','top') xlabel('n (samples)'); ylabel('error')
Simulink® でのアルゴリズムの実装
Simulink® と Fixed-Point Designer™ をご利用の場合は、上記のアルゴリズムに該当する、次のモデルを実行できます。出力 y_sim
は、上記で MATLAB コードで計算された変数 y
に等しい固定小数点変数です。
MATLAB コード内と同様に、ブロック内の固定小数点パラメーターは実際のシステムに合うように変更できます。これらは、上記の例では MATLAB コードに合わせて設定されています。ブロックをダブルクリックして設定を確認します。
% Set up the From Workspace variable x_sim.time = n; x_sim.signals.values = x; x_sim.signals.dimensions = 1; % Run the simulation out_sim = sim('fitdf2filter_demo', 'SaveOutput', 'on', ... 'SrcWorkspace', 'current'); % Open the model open_system('fitdf2filter_demo') % Verify that the Simulink results are the same as the MATLAB file isequal(y, out_sim.get('y_sim'))
ans = logical 1
この例の前提条件
例を簡単にするために、既定の数学パラメーターである、最も近い正の整数方向への丸め、オーバーフローで飽和、全精度の積、および最大桁数の和を使用しています。これらのパラメーターはいずれも、実際のシステムに合わせて変更できます。
これらの設定は、アルゴリズム開発の出発点として選択されています。この MATLAB ファイルのコピーを保存し、これらのパラメーターを使用してみて、出力への影響を確認します。入力が異なる場合、アルゴリズムの動作はどうなるでしょうか。その他のパラメーター (丸めモードやオーバーフロー モードなど) の設定方法の詳細は、fi、fimath、および numerictype のヘルプを参照してください。
% Reset the global fimath
globalfimath(globalFimathAtStart);
fipref(FIPREF_STATE);