固定小数点アルゴリズムの開発

ライブスクリプトを開く

この例では、簡単な固定小数点アルゴリズムの開発方法と検証方法を説明します。この例は、以下のアルゴリズム開発の手順に沿っています。

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);

NumericTypeScope オブジェクトを作成し、出力 y のダイナミックレンジを表示します。

hscope2 = NumericTypeScope;

step(hscope2,y);

% Clear the information stored in the object hscope2
reset(hscope2);

まず、変数 z (状態) について表示されている情報を解析します。ヒストグラムから、ダイナミックレンジが [(2^(1) 2^(-12)] 内にあることがわかります。既定では、スコープでは numerictype('double') データ型が使用されます。スコープの [推定されたデータ型] フィールドを調整して、データ型を numerictype(true,16,14) にするとオーバーフローがゼロになることがわかります。

次に、変数 y (出力) を確認します。ヒストグラムプロットから、ダイナミックレンジが [(2^(1) 2^(-13)] 内にあることがわかります。ここでも、データ型を 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');

Figure contains an axes object. The axes object with title Magnitude response of Floating-point and Fixed-point results, xlabel Normalized Frequency, ylabel Magnitude in dB contains 3 objects of type line. These objects represent Input, Floating point output, Fixed point output.

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');

Figure contains an axes object. The axes object with title Magnitude response of the filter, xlabel Frequency, ylabel Magnitude in dB contains an object of type line.

入力信号の高周波数はローパスフィルターで減衰されており、これは期待どおりの動作です。

誤差のプロット

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')

Figure contains an axes object. The axes object with xlabel n (samples), ylabel error contains 5 objects of type line, text.

Simulink でのアルゴリズムの実装

Simulink® と Fixed-Point Designer™ をご利用の場合は、上記のアルゴリズムに該当する fitdf2filter_demo モデルを実行できます。出力 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);