3

面倒で申し訳ありませんが、数十の記事の助けを借りてコードを数回確認しましたが、それでもKFは機能しません。「機能しない」とは、KFによる見積もりが間違っていることを意味します。これは、Real、Noiseed、およびKFの推定位置の素敵なペーストです（ほんの小さなチャンク）。

私の例は、私が見つけたすべてのチュートリアルと同じです-位置と速度の状態ベクトルがあります。位置はメートル単位で、空中の垂直位置を表します。私の実際のケースはスカイダイビング（パラシュート付き）です。私のサンプル生成データでは、3000mから開始し、速度は10m/sであると想定しています。

PS：マトリックスの計算は大丈夫だと確信しています。ロジックにエラーがあるはずです。

ここでデータを生成します：

void generateData(float** inData, float** noisedData, int x, int y){
    inData[0][0]= 3000; //start position
    inData[1][0]= -10; // 10m/s velocity; minus because we assume it's falling

    noisedData[0][0]= 2998; 
    noisedData[1][0]= -10;

    for(int i=1; i<x; i++){
        inData[0][i]= inData[0][i-1] + inData[1][i-1]; 
        inData[1][i]= inData[1][i-1]; //the velocity doesn't change for simplicity's sake

        noisedData[0][i]=inData[0][i]+(rand()%6-3); //we add noise to real measurement
        noisedData[1][i]=inData[1][i]; //velocity has no noise
    }  
}

そしてこれが私の実装です（行列の初期化はWikipedia Kalmanの例に基づいています）：

int main(int argc, char** argv) {
    srand(time(NULL));

    float** inData = createMatrix(100,2); //2 rows, 100 columns
    float** noisedData = createMatrix(100,2);
    float** estData = createMatrix(100,2);

    generateData(inData, noisedData, 100, 2);

    float sampleRate=0.1; //10hz

    float** A=createMatrix(2,2);
    A[0][0]=1;
    A[0][1]=sampleRate;
    A[1][0]=0;
    A[1][1]=1;

    float** B=createMatrix(1,2);
    B[0][0]=pow(sampleRate,2)/2;
    B[1][0]=sampleRate;

    float** C=createMatrix(2,1);
    C[0][0]=1; //we measure only position
    C[0][1]=0;


    float u=1.0; //acceleration magnitude
    float accel_noise=0.2; //acceleration noise
    float measure_noise=1.5; //1.5 m standard deviation
    float R=pow(measure_noise,2); //measure covariance
    float** Q=createMatrix(2,2); //process covariance
    Q[0][0]=pow(accel_noise,2)*(pow(sampleRate,4)/4);
    Q[0][1]=pow(accel_noise,2)*(pow(sampleRate,3)/2);
    Q[1][0]=pow(accel_noise,2)*(pow(sampleRate,3)/2);
    Q[1][1]=pow(accel_noise,2)*pow(sampleRate,2);

    float** P=createMatrix(2,2); //covariance update
    P[0][0]=0;
    P[0][1]=0; 
    P[1][0]=0; 
    P[1][1]=0;

    float** P_est=createMatrix(2,2);
    P_est[0][0]=P[0][0];
    P_est[0][1]=P[0][1];
    P_est[1][0]=P[1][0];
    P_est[1][1]=P[1][1];

    float** K=createMatrix(1,2); //Kalman gain

    float** X_est=createMatrix(1,2); //our estimated state
    X_est[0][0]=3000; X_est[1][0]=10; 

    // !! KALMAN ALGORITHM START !! //
    for(int i=0; i<100; i++)
    {        
        float** temp;
        float** temp2;
        float** temp3;

        float** C_trans=matrixTranspose(C,2,1);
        temp=matrixMultiply(P_est,C_trans,2,2,1,2); //2x1
        temp2=matrixMultiply(C,P_est,2,1,2,2); //1x2
        temp3=matrixMultiply(temp2,C_trans,2,1,1,2); //1x1
        temp3[0][0]+=R;
        K[0][0]=temp[0][0]/temp3[0][0]; // 1. KALMAN GAIN
        K[1][0]=temp[1][0]/temp3[0][0];

        temp=matrixMultiply(C,X_est,2,1,1,2);
        float diff=noisedData[0][i]-temp[0][0]; //diff between meas and est

        X_est[0][0]=X_est[0][0]+(K[0][0]*diff);  // 2. ESTIMATION CORRECTION
        X_est[1][0]=X_est[1][0]+(K[1][0]*diff);

        temp=createMatrix(2,2);
        temp[0][0]=1; temp[0][1]=0; temp[1][0]=0; temp[1][1]=1;
        temp2=matrixMultiply(K,C,1,2,2,1);
        temp3=matrixSub(temp,temp2,2,2,2,2);
        P=matrixMultiply(temp3,P_est,2,2,2,2);  // 3. COVARIANCE UPDATE



        temp=matrixMultiply(A,X_est,2,2,1,2);
        X_est[0][0]=temp[0][0]+B[0][0]*u; 
        X_est[1][0]=temp[1][0]+B[1][0]*u; // 4. PREDICT NEXT STATE


        temp=matrixMultiply(A,P,2,2,2,2);
        float** A_inv=getInverse(A,2);
        temp2=matrixMultiply(temp,A_inv,2,2,2,2);
        P_est=matrixAdd(temp2,Q,2,2,2,2); // 5. PREDICT NEXT COVARIANCE


        estData[0][i]=X_est[0][0]; //just saving here for later to write out
        estData[1][i]=X_est[1][0];
    }

    for(int i=0; i<100; i++) printf("%4.2f  :  %4.2f  :  %4.2f \n", inData[0][i], noisedData[0][i], estData[0][i]); // just writing out

    return (EXIT_SUCCESS);
}
4

2 に答える 2

0

あなたは奇妙な配列のインデックス付けをたくさん行っています。

float** A=createMatrix(2,2);
A[0][0]=1;
A[0][3]=sampleRate;
A[1][0]=0;
A[1][4]=1;

配列の境界外にインデックスを作成すると、どのような結果が期待されますか?

于 2012-08-10T11:07:28.590 に答える