matlab IIRフィルタは、同じフィルタ/コードがC/C＃などで実装されているときに、異なる出力を出力します。

Question

今週、実際のヘッドスクラッチャーに遭遇しました。私はfilterにMathWorks社のMATLAB sourceから直接コピーして、私はC＃でIIRフィルタを実装しています自分の時間ドメインフィルタ機能（直接型II転置）：

 // direct form ii transposed 
     for (int i = 0; i < data.Length; i++) 
     { 
      Xi = data[i]; 
      Yi = b[0] * Xi + Z[0]; 

      for (int j = 1; j < order; j++) 
       Z[j - 1] = b[j] * Xi + Z[j] - a[j] * Yi; 

      Z[order - 1] = b[order] * Xi - a[order] * Yi; 

      output[i] = Yi; 
     } 
     return output; 

何奇妙なことは、私は衝動のフィルタをテストするとき、私はMatlabによって報告されたものとは若干異なる値を得る。私はMatlabからフィルタ係数も取得しています。

[b,a] = butter(3, [.0360, .1160], 'bandpass'); 
x = zeros(50,1); 
x(1) = 1.0; 
y = filter(b,a,x) 

私は私のC＃コードでcoeffientsとしてbとa内の値を使用します。ここでは、コードです。

のMatlabによって報告されるようにyのための最初のいくつかの値：

>> y(1:13) 

    ans = 

    0.0016 
    0.0084 
    0.0216 
    0.0368 
    0.0487 
    0.0537 
    0.0501 
    0.0382 
    0.0194 
    -0.0038 
    -0.0286 
    -0.0519 
    -0.0713 

これは私のC＃のポートと異なっていたので、私は直接Cファイルへfilterからコードをコピーして、同じ係数を用いて、そこに走りました。私はfilterにソースを注意深く見て、私は前の係数をマッサージの証拠を見ない

[0] 0.0016 double 
[1] 0.0086161600000000012 double 
[2] 0.022182403216000009 double 
[3] 0.038161063110481647 double 
[4] 0.051323531488129848 double 
[5] 0.05827273642334313 double 
[6] 0.057456579295617483 double 
[7] 0.048968543791003127 double 
[8] 0.034196988694833064 double 
[9] 0.015389667539999874 double 
[10] -0.0048027826346631469 double 
[11] -0.023749640330880527 double 
[12] -0.039187648694732236 double 
[13] -0.04946710058803272 double 

：出力は少し私は私のC＃実装に入ったことを、インパルス応答のバージョンオフ、まったく同じでしたフィルタ出力を計算する。 filterは、a[0]が1に等しくない場合にのみフィードフォワード係数を正規化します（この場合は最も確かにそうです）。それ以外にも、私はMatlabのようにMatlabのCコードから正確に同じフィルタ出力を得ることを期待しています。

私のフィルタが正確に正しいと確信する必要があるので、私はこの鑑別がどこから来ているのか本当に知りたいと思います。フィルター係数を確認して再確認しました。これらはC/C＃とMatlabの間では同じです。

これらの「間違った」値を取得するために使用した完全なCファイルが続きます。私は固定数の状態（この場合は6）とN状態の一般的な場合（ここでコメントアウト）として実装されたフィルタの両方を試しました。どちらも、Matlabのソースコードから来て、両方が同じで、「間違った」の出力を生成：私はフィルタ係数のすべての数字を使用している場合

# define order 6 
# define len 50 

int main(void){ 

    // filter coeffs  
    float a[7] = { 1.0000, -5.3851, 12.1978, -14.8780, 10.3077, -3.8465, 0.6041 }; 
    float b[7] = { 0.0016, 0.0, -0.0047, 0.0, 0.0047, 0, -0.0016 }; 
    float a1 = a[1], a2 = a[2], a3 = a[3], a4 = a[4], a5 = a[5], a6 = a[6]; 

    // input, output, and state arrays 
    float X[len]; 
    float Y[len]; 
    float Z[order]; 
    float Xi; 
    float Yi; 
    float z0, z1, z2, z3,z4, z5; 

    // indeces 
    int i,j; 

    // initialize arrays 
    for(i=0;i<len;i++) { 
     X[i] = 0.0; 
     Y[i] = 0.0; 
    } 
    X[0] = 1.0; 

    for(i=0;i<order;i++) 
     Z[i] = 0.0; 

    z0 = Z[0]; 
    z1 = Z[1]; 
    z2 = Z[2]; 
    z3 = Z[3]; 
    z4 = Z[4]; 
    z5 = Z[5]; 
    i = 0; 
    while (i < len) { 
     Xi = X[i]; 
     Yi = b[0] * Xi + z0; 
     z0 = b[1] * Xi + z1 - a1 * Yi; 
     z1 = b[2] * Xi + z2 - a2 * Yi; 
     z2 = b[3] * Xi + z3 - a3 * Yi; 
     z3 = b[4] * Xi + z4 - a4 * Yi; 
     z4 = b[5] * Xi + z5 - a5 * Yi; 
     z5 = b[6] * Xi  - a6 * Yi; 
     Y[i++] = Yi; 
    } 


    //// copied from matlab filter source code 
    //i=0; 
    //while (i < len) { 
//  Xi = X[i];      // Get signal 
//  Yi = b[0] * Xi + Z[0];   // Filtered value 
//  for (j = 1; j < order; j++) { // Update conditions 
//   Z[j - 1] = b[j] * Xi + Z[j] - a[j] * Yi; 
//  } 
//  Z[order - 1] = b[order] * Xi - a[order] * Yi; 
//  
//  Y[i++] = Yi;      // Write to output 
// } 


} 

Answer 1

が、私はMatlabの答えを得ます。具体的に

、

float a[7] = {1, 
       -5.3850853610906082025167052052, 
       12.1978301571107792256043467205, 
       -14.8779557262737220924009307055, 
       10.3076512098041828124905805453, 
       -3.84649525959781790618308150442, 
       0.604109699507274999774608659209}; 
    float b[7] = {0.00156701035058826889344307797813, 0, 
       -0.0047010310517648064634887994373, 0, 
       0.0047010310517648064634887994373, 0, 
       -0.00156701035058826889344307797813}; 

（フィルタ係数を使用して得られscipyのダウンロード：。私が作っていないよので、[b,a] = scipy.signal.butter(3, [.0360, .1160], 'bandpass')）代わりにこれにより

、あなたのCコードをプリントアウトすることができます：

[0] = 0.0015670103020966053009033203125 
[1] = 0.00843848474323749542236328125 
[2] = 0.02162680588662624359130859375 
[3] = 0.036844909191131591796875 
[4] = 0.048709094524383544921875 
[5] = 0.05368389189243316650390625 
[6] = 0.05014741420745849609375 
[7] = 0.0382179915904998779296875 
[8] = 0.0194064676761627197265625 
[9] = -0.003834001719951629638671875 

あなたのMatlab出力と一致します。

float/doubleとは関係ありません。このような状況では、ある言語から別の言語へ実装を移植する際には、正確な入力を保証し、相対的なエラー - ではなく、を使用して出力をコピーアンドペーストして比較するよう努力してください。

（PS。、どのように通過帯域の対称性のおかげで、注意してください、bの他のすべての要素はこれが必要なプの数を減らすために使用することができます0です！）