変数ログを使用した固定小数点の数値的問題のデバッグ

この例では、組み込みの変換より優れた変換が必要となる場合に、固定小数点コードをデバッグするためのベスト プラクティスをいくつか示します。

必要条件

1. ローカルの作業フォルダー (例: c:\kalman_filter) を作成します。

2. ローカルの作業フォルダーで以下のファイルを作成します。

   •  関数 kalman_filter

     これはプロジェクトのメイン エントリポイント関数です。

     function [y] = kalman_filter(z,N0)
         %#codegen
         A = kalman_stm();
         
         % Measurement Matrix
         H = [1 0];
         
         % Process noise variance
         Q = 0;
         % Measurement noise variance
         R = N0 ;
         
         persistent x_est p_est
         if isempty(x_est)
             % Estimated state
             x_est = [0; 1];
             % Estimated error covariance
             p_est = N0 * eye(2, 2);
         end
             
         % Kalman algorithm
         % Predicted state and covariance
         x_prd = A * x_est;   
         p_prd = A * p_est * A' + Q;
         
         % Estimation
         S = H * p_prd' * H' + R;
         B = H * p_prd';
         klm_gain = matrix_solve(S,B)';
         
         % Estimated state and covariance
         x_est = x_prd + klm_gain * (z - H * x_prd);
         p_est = p_prd - klm_gain * H * p_prd;
         
         % Compute the estimated measurements
         y = H * x_est;    
     end

   •  関数 kalman_stm

     この関数は関数 kalman_filter から呼び出され、状態遷移マトリックスを計算します。

     function A = kalman_stm()    
         f0 = 200;
         dt = 1/1e4;
         % Kalman filter initialization
         % State transition Matrix
         A = [cos(2*pi*f0*dt), -sin(2*pi*f0*dt);
             sin(2*pi*f0*dt), cos(2*pi*f0*dt)];    
     end

   •  関数 matrix_solve

     この関数は行列の左除算をより効率的に実装しています。

     function x = matrix_solve(a,b)
         %fixed-point conversion friendly matrix solve: a * x = b
         
         % initialize x
         x = zeros(size(a,1),size(b,2));
         % compute lu decomposition of a
         [l, u] = lu_replacement(a);
         % solve x = a\b for every column of b 
         % through forward and backward substitution
         for col = 1:size(b,2)
             bcol = b(:,col);
             y = forward_substitute(l,bcol);
             x(:,col) = back_substitute(u,y);
         end
         
     end

   •  関数 lu_replacement

     この関数は関数 matrix_solve から呼び出されます。

     function [l,A]=lu_replacement(A)
         N=size(A,1);
         l = eye(N);
         for n=1:N-1
             piv = A(n,n);
             for k=n+1:N
                 mult = divide_no_zero(A(k,n),piv);
                 A(k,:) = -mult*A(n,:) + A(k,:);
                 l(k,n) = mult;
             end
         end
     end

   •  関数 forward_substitute

     この関数は関数 matrix_solve から呼び出されます。

     function y = forward_substitute(l,b)
         % forward substitution
         N = size(b,1);
         y = zeros(N,1);
         % forward substitution
         y(1) = divide_no_zero(b(1),l(1,1));
         for n = 2:N
             acc = 0;
             for k = 1:n-1
                 acc(:) = acc + y(k)*l(n,k);
             end
             y(n) = divide_no_zero((b(n)-acc),l(n,n));
         end
         
     end

   •  関数 back_substitute

     この関数は関数 matrix_solve から呼び出されます。

     function x = back_substitute(u,y)
         % backwards substitution
         N = size(u,1);
         x = zeros(N,1);
         
         % backward substitution
         x(N) = divide_no_zero(y(N),u(N,N));
         
         for n = (N-1):(-1):(1)
             acc = 0;
             for k = n:(N)
                 acc(:) = acc + x(k)*u(n,k);
             end
             x(n) = divide_no_zero((y(n) - acc),u(n,n));
         end
     end

   •  関数 divide_no_zero

     この関数は関数 lu_replacement、forward_substitute および back_substitute から呼び出されます。

     function y = divide_no_zero(num, den)
         % Divide and avoid division by zero
         if den == 0
             y = 0;
         else
             y = num/den;
         end
     end

   •  テスト ファイル kalman_filter_tb

     このスクリプトはノイズを含む正弦波を生成し、関数 kalman_filter を呼び出してノイズを含む信号をフィルターします。その後、信号をプロットして比較します。

     % KALMAN FILTER EXAMPLE TEST BENCH
     clear all
     step = ((400*pi)/1000)/10;
     TIME_STEPS = 400;
     X = 0:step:TIME_STEPS;
     rng default;
     rng(1);
     Orig_Signal = sin(X);
     Noisy_Signal = Orig_Signal + randn(size(X));
     Clean_Signal = zeros(size(X));
     for i = 1:length(X)  
     Clean_Signal(i) = kalman_filter(Noisy_Signal(i), 1);
     end
     figure
     subplot(5,1,1)
     plot(X,rand(size(X)))
     axis([1 TIME_STEPS 0 1.25]);
     title('Noise')
     
     % Plot Noisy Signal
     subplot(5,1,2)
     plot(X,Noisy_Signal)
     axis([1 TIME_STEPS -4 4]);
     title('Noisy Signal')
     
     % Plot Filtered Clean Signal
     subplot(5,1,3)
     plot(X,Clean_Signal)
     axis([1 TIME_STEPS -2 2]);
     title('Filtered Signal')
     
     % Plot Original Signal
     subplot(5,1,4)
     plot(X,Orig_Signal)
     axis([1 TIME_STEPS -2 2]);
     title('Original Signal')
     
     % Plot Error
     subplot(5,1,5)
     plot(X, (Clean_Signal - Orig_Signal))
     axis([1 TIME_STEPS -1 1]);
     title('Error')
     figure(gcf)

既定の構成を使用した固定小数点への変換

1. アプリ ギャラリーから固定小数点コンバーター アプリを開きます。

2. [選択] ページで kalman_filter.m ファイルを探して [開く] をクリックします。

3. [次へ] をクリックします。[入力の型を定義] ページで kalman_filter_tb ファイルを参照します。[入力の型の自動定義] をクリックします。

   テストファイルが実行され、ノイズを含む入力信号、フィルター処理された信号、理想的なフィルターで処理された信号、およびフィルター処理された信号と理想的なフィルターで処理された信号の間の差がプロットされます。

   

4. [次へ] をクリックします。[固定小数点に変換] ページで [解析] をクリックし、既定の設定を使用して範囲情報の収集およびデータ型の推奨を行います。

5. [変換] をクリックし、推奨されたデータ型を適用します。

6. [テスト] 矢印をクリックし、[比較プロット用の入力と出力を記録] チェック ボックスをオンにします。[テスト] をクリックします。固定小数点コンバーターはテスト ファイル kalman_filter_tb.m を実行し、生成された固定小数点のコードをテストします。浮動小数点と固定小数点のシミュレーションが実行され、入力変数と出力変数の誤差が計算されます。

   生成されたプロットによって、現在の固定小数点の実装では良好な結果が得られないことが示されます。

   

   出力変数 y の誤差は非常に大きく、282% を超えています。

   

数値的問題の発生元の特定

数値的問題の原因として考えられる任意の関数の入力変数と出力変数を、最上位関数の出力引数に記録します。

1. [ソース コード] ペインで kalman_filter をクリックし、浮動小数点コードに戻ります。

   [テスト] フェーズ中に [比較プロット用の入力と出力を記録] オプションを選択すると、最上位関数 kalman_filter の入力変数と出力変数がプロット用に自動的に記録されます。

2. 関数 kalman_filter は関数 matrix_solve を呼び出します。この関数には他の複数の関数への呼び出しが含まれます。数値的問題の発生元が関数 matrix_solve であるかどうかを調査するには、[ソース コード] ペインで kalman_filter、matrix_solve を選択します。

   [ログ データ] 列で、記録する関数の入力変数と出力変数を選択します。この例では、a、b、および x の 3 つをすべて選択します。

   

3. [テスト] をクリックします。

   生成されたプロットによって関数 matrix_solve の出力変数の大きな誤差が示されます。

   

fimath 設定の調整

1. [固定小数点に変換] ページで [設定] をクリックします。

   [fimath] で [丸め方法] を [最も近い正の整数方向] に設定します。[オーバーフロー アクション] を [飽和] に設定します。

2. [変換] をクリックして新しい設定を適用します。

3. [テスト] の隣の矢印をクリックし、[比較プロット用の入力と出力を記録] が選択されていることを確認します。関数の入力変数または出力変数のログを有効にします。[テスト] をクリックします。

   最上位関数の出力変数 y のプロットを確認します。

   

   新しい fimath の設定によって出力は向上しましたが、まだ多少の誤差が残っています。

語長設定の調整

既定の語長を調整することによって、生成される固定小数点設計の精度が向上します。

1. [設定] をクリックして既定の語長を 32 に変更します。[変換] をクリックして新しい設定を適用します。

2. [テスト] をクリックします。出力変数 y の誤差が累計されます。

   

3. 固定小数点コンバーターおよびプロット ウィンドウを閉じます。

定数関数の置換

関数 kalman_stm が状態遷移マトリックスを計算すると、結果は定数になります。この関数を固定小数点に変換する必要はありません。計算による不必要な量子化を避けるために、この関数を倍精度の定数で置き換えます。関数を定数に置き換えることにより、状態遷移マトリックスが量子化されるのは 1 回のみになります。

1. [ソース コード] ペインで関数 kalman_filter をクリックします。関数 kalman_filter を編集します。関数 kalman_stm の呼び出しを等価な double の定数で置き換えます。

   A = [0.992114701314478, -0.125333233564304;
       0.125333233564304, 0.992114701314478];

   変更内容を保存します。

2. [解析] をクリックして推奨を更新します。

3. [変換] をクリックして新しい推奨を適用します。

4. [テスト] をクリックします。関数の出力 y のプロットの誤差は 10^-6 の次数になり、この設計で許容できます。