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

スクリプトを開く

この例では、簡単な固定小数点アルゴリズムの開発方法と検証方法を説明します。

アルゴリズム開発の簡単な例

この例では、簡単な固定小数点フィルターアルゴリズムの開発と検証について説明しています。以下の手順に従います。

1) 2 次フィルターアルゴリズムを実装し、倍精度浮動小数点でシミュレートします。

2) コードをインストルメント化し、出力と状態のダイナミックレンジを可視化します。

3) 変数のデータ型を変更することにより、アルゴリズムを固定小数点に変換します。アルゴリズム自体は変わりません。

4) 固定小数点と浮動小数点の結果を比較し、プロットします。

浮動小数点変数の定義

倍精度浮動小数点でアルゴリズムを開発します。2 次ローパスフィルターを使用して、入力信号の高周波数を除去します。

b = [ 0.25 0.5      0.25    ]; % Numerator coefficients
a = [ 1    0.09375  0.28125 ]; % Denominator coefficients
% Random input that has both high and low frequencies.
s = rng; rng(0,'v5uniform');
x = randn(1000,1);
rng(s); % restore RNG state
% Pre-allocate the output and state for speed.
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

% Save the Floating-Point Result
ydouble = y;

ダイナミックレンジの可視化

固定小数点に変換するには、変数の範囲を知る必要があります。アルゴリズムの複雑さによって、このタスクは簡単な場合もあれば非常に難しい場合もあります。この例では、入力値の範囲がわかっているため、適切な固定小数点データ型を選択することは簡単です。出力 (y) と状態 (z) の範囲が不明であるため、この 2 つに注目します。出力と状態のダイナミックレンジを表示するには、コードを少し変更してインストルメントします。2 つの NumericTypeScope オブジェクトを作成し、出力 (y) と状態 (z) のダイナミックレンジを同時に表示します。

浮動小数点コードのインストルメント

% Reset states
z = [0;0];

hscope1 = NumericTypeScope;
hscope2 = 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
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) となります。提案されたこのデータ型を使用すると、オーバーフローまたはアンダーフローが表示されません。

固定小数点変数の定義

変数を固定小数点に変換し、アルゴリズムを再度実行します。ログをオンにし、選択したデータ型によるオーバーフローとアンダーフローを確認します。

% Turn on logging to see overflows/underflows.
FIPREF_STATE = get(fipref);
reset(fipref)
fp = fipref;
default_loggingmode = fp.LoggingMode;
fp.LoggingMode = 'On';
% Capture the present state of and reset the global fimath to the factory
% settings.
globalFimathAtStart = fimath;
resetglobalfimath;
% Define the fixed-point types for the variables in the below format:
%   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® ファイル関数に含めることをお勧めします。データ型が不明なアルゴリズムの記述と使用の例は、filimitcycledemo.m を参照してください。

浮動小数点と固定小数点の結果の比較とプロット

浮動小数点と固定小数点の結果の振幅応答とフィルターの応答をプロットし、固定小数点への変換時にフィルターが期待どおりの動作をしているかどうかを確認します。

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 コードに合わせて設定されています。ブロックをダブルクリックして設定を確認します。

if fidemo.hasSimulinkLicense

    % 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
    fitdf2filter_demo

    % Verify that the Simulink results are the same as the MATLAB file
    isequal(y, out_sim.get('y_sim'))

end

ans =

  logical

   1

この例の前提条件

例を簡単にするために、既定の数学パラメーターである、最も近い正の整数方向への丸め、オーバーフローで飽和、全精度の積、および最大桁数の和を使用しています。これらのパラメーターはいずれも、実際のシステムに合わせて変更できます。

これらの設定は、アルゴリズム開発の出発点として選択されています。この MATLAB ファイルのコピーを保存し、これらのパラメーターを使用してみて、出力への影響を確認します。入力が異なる場合、アルゴリズムの動作はどうなるでしょうか。その他のパラメーター (丸めモードやオーバーフローモードなど) の設定方法の詳細は、fi、fimath、および numerictype のヘルプを参照してください。

close all force;
bdclose all;
% Reset the global fimath
globalfimath(globalFimathAtStart);
fipref(FIPREF_STATE);