Main Content

filter

1 次元のデジタル フィルター

説明

y = filter(b,a,x) は、分子係数 b と分母係数 a で定義される有理伝達関数を使用して、入力データ x をフィルター処理します。

a(1)1 に等しくない場合、filtera(1) でフィルター係数を正規化します。このため a(1) は非ゼロでなければなりません。

  • x がベクトルの場合、filter はフィルター処理されたデータを x と同じサイズのベクトルとして返します。

  • x が行列の場合、filter は最初の次元に沿って機能し、各列のフィルター処理されたデータを返します。

  • x が多次元配列の場合、filter は、サイズが 1 でない最初の配列次元に沿って機能します。

y = filter(b,a,x,zi) は、フィルター遅延の初期条件 zi を使用します。zi の長さは、max(length(a),length(b))-1 と等しくなければなりません。

y = filter(b,a,x,zi,dim) は次元 dim に沿って機能します。たとえば、x が行列の場合、filter(b,a,x,zi,2) は各行のフィルター処理されたデータを返します。

[y,zf] = filter(___) は前述の任意の構文を使用し、フィルター遅延の最終状態 zf も返します。

すべて折りたたむ

移動平均フィルターは、ノイズを含むデータの平滑化に使用される一般的なメソッドです。この例では、関数 filter を使用して、データのベクトルに沿って平均を計算します。

ランダム ノイズで乱した正弦波データの 1 行 100 列の行ベクトルを作成します。

t = linspace(-pi,pi,100);
rng default  %initialize random number generator
x = sin(t) + 0.25*rand(size(t));

移動平均フィルターは、データに沿って長さ windowSize のウィンドウをスライドし、各ウィンドウに含まれるデータの平均を計算します。次の差分方程式は、ベクトル x の移動平均フィルターを定義します。

y(n)=1windowSize(x(n)+x(n-1)+...+x(n-(windowSize-1))).

ウィンドウ サイズが 5 の場合の有理伝達関数の分子係数と分母係数を計算します。

windowSize = 5; 
b = (1/windowSize)*ones(1,windowSize);
a = 1;

データの移動平均を求め、元のデータに対してプロットします。

y = filter(b,a,x);

plot(t,x)
hold on
plot(t,y)
legend('Input Data','Filtered Data')

この例では、データ行列を次の有理伝達関数でフィルター処理します。

H(z)=b(1)a(1)+a(2)z-1=11-0.2z-1

2 行 15 列の乱数入力データの行列を作成します。

rng default  %initialize random number generator
x = rand(2,15);

有理伝達関数の分子係数と分母係数を定義します。

b = 1;
a = [1 -0.2];

x の 2 番目の次元に沿って伝達関数を適用し、各行の 1 次元デジタル フィルターを返します。フィルター処理されたデータに対して元のデータの 1 行目をプロットします。

y = filter(b,a,x,[],2);

t = 0:length(x)-1;  %index vector

plot(t,x(1,:))
hold on
plot(t,y(1,:))
legend('Input Data','Filtered Data')
title('First Row')

フィルター処理されたデータに対して入力データの 2 行目をプロットします。

figure
plot(t,x(2,:))
hold on
plot(t,y(2,:))
legend('Input Data','Filtered Data')
title('Second Row')

ウィンドウ サイズを 3 として移動平均フィルターを定義します。

windowSize = 3;
b = (1/windowSize)*ones(1,windowSize);
a = 1;

1 行 6 列のデータ行ベクトルの 3 点移動平均を求めます。

x = [2 1 6 2 4 3];
y = filter(b,a,x)
y = 1×6

    0.6667    1.0000    3.0000    3.0000    4.0000    3.0000

既定では、関数 filter はフィルター遅延をゼロに初期化し、過去の入力と出力が両方ともゼロであると仮定します。この場合、y の最初の 2 つの要素はそれぞれ、x の最初の要素と最初の 2 つの要素の 3 点移動平均です。つまり、最初の要素 0.6667 は 2 の 3 点平均で、2 番目の要素 1 は 2 と 1 の 3 点平均です。

データに過去の入力と出力を追加で含めるには、初期条件をフィルター遅延として指定します。現在の入力におけるこの初期条件は、過去の入力 (および過去の出力) に同じ伝達関数を適用して取得する最終状態です。たとえば、過去の入力 [1 3] を含めるとします。フィルター遅延がない場合、過去の出力は (0+0+1)/3 および (0+1+3)/3 です。

x_past = [1 3];
y_past = filter(b,a,x_past)
y_past = 1×2

    0.3333    1.3333

ただし、これらの過去の入力の末尾はゼロであると仮定して、引き続き同じ伝達関数を適用して非ゼロの出力をさらに生成できます。これらの追加出力は (1+3+0)/3 および (3+0+0)/3 であり、過去の入力から取得した最終状態を表します。これらの最終状態を計算するには、関数 filter の 2 番目の出力引数を指定します。

[y_past,zf] = filter(b,a,x_past)
y_past = 1×2

    0.3333    1.3333

zf = 2×1

    1.3333
    1.0000

現在のデータに過去の入力を含めるには、関数 filter の 4 番目の入力引数を使用してフィルター遅延を指定します。過去のデータの最終状態を現在のデータの初期条件として使用します。

y = filter(b,a,x,zf)
y = 1×6

    2.0000    2.0000    3.0000    3.0000    4.0000    3.0000

この場合、y の最初の要素は 1、3、2 の 3 点移動平均、つまり 2 になり、y の 2 番目の要素は 3、2、1 の移動平均、つまり 2 になります。

フィルター遅延の初期条件と最終状態を使用し、分割してデータをフィルター処理します。メモリ制限に考慮が必要な場合は特に有効です。

大規模な乱数データ シーケンスを作成し、x1x2 の 2 つのセグメントに分割します。

x = randn(10000,1);

x1 = x(1:5000);
x2 = x(5001:end);

全体のシーケンス x は、x1x2 の垂直連結です。

有理伝達関数の分子係数と分母係数を定義します。

H(z)=b(1)+b(2)z-1a(1)+a(2)z-1=2+3z-11+0.2z-1.

b = [2,3];
a = [1,0.2];

サブシーケンス x1 および x2 を一度に 1 つずつフィルター処理します。最初のセグメントが終了した時点のフィルターの内部状態を保存するために、x1 のフィルター処理の最終条件を出力します。

[y1,zf] = filter(b,a,x1);

x1 のフィルター処理の最終条件を、2 番目のセグメント x2 のフィルター処理の初期条件として使用します。

y2 = filter(b,a,x2,zf);

y1x1 からフィルター処理されたデータで、y2x2 からフィルター処理されたデータです。フィルター処理された全体のシーケンスは、y1y2 の垂直連結です。

比較のために全体シーケンスを一度にフィルター処理します。

y = filter(b,a,x);

isequal(y,[y1;y2])
ans = logical
   1

入力引数

すべて折りたたむ

有理伝達関数の分子係数。ベクトルとして指定します。

データ型: double | single | int8 | int16 | int32 | int64 | uint8 | uint16 | uint32 | uint64 | logical
複素数のサポート: あり

有理伝達関数の分母係数。ベクトルとして指定します。

データ型: double | single | int8 | int16 | int32 | int64 | uint8 | uint16 | uint32 | uint64 | logical
複素数のサポート: あり

入力データ。ベクトル、行列または多次元配列として指定します。

データ型: double | single | int8 | int16 | int32 | int64 | uint8 | uint16 | uint32 | uint64 | logical
複素数のサポート: あり

フィルター遅延の初期条件。ベクトル、行列または多次元配列として指定します。

  • zi がベクトルの場合、その長さは max(length(a),length(b))-1 でなければなりません。

  • zi が行列または多次元配列の場合、最初の次元のサイズは max(length(a),length(b))-1 でなければなりません。残りの次元のサイズはそれぞれ対応する x の次元のサイズと同じでなければなりません。たとえば、3×4×5 の配列 x の 2 番目の次元 (dim = 2) に沿って filter を使用するとします。配列 zi のサイズは、[max(length(a),length(b))-1]×3×5 でなければなりません。

既定値は [] で指定され、すべてのフィルター遅延はゼロに初期化されます。

データ型: double | single | int8 | int16 | int32 | int64 | uint8 | uint16 | uint32 | uint64 | logical
複素数のサポート: あり

演算の対象の次元。正の整数のスカラーとして指定します。次元を指定しない場合、既定値はサイズが 1 より大きい最初の配列次元です。

2 次元の入力配列 x について考えます。

  • dim = 1 の場合、filter(b,a,x,zi,1)x の列に沿って動作し、各列に適用されるフィルターを返します。

    filter(b,a,x,zi,1) column-wise operation

  • dim = 2 の場合、filter(b,a,x,zi,2)x の行に沿って動作し、各行に適用されるフィルターを返します。

    filter(b,a,x,zi,2) row-wise operation

dimndims(x) より大きい場合、filter は、xdim までの追加の次元 (サイズは 1) をもつと見なします。たとえば、x が 2 行 3 列のサイズの行列であり、dim = 3 の場合、filterx のサイズが 2×3×1 であると見なし、その 3 番目の次元に沿って動作します。

データ型: double | single | int8 | int16 | int32 | int64 | uint8 | uint16 | uint32 | uint64 | logical

出力引数

すべて折りたたむ

フィルター処理されたデータ。入力データ x と同じサイズのベクトル、行列または多次元配列として返されます。

x の型が single である場合、filter はネイティブ レベルの単精度で計算し、y の型も single になります。それ以外の場合、ydouble 型として返されます。

データ型: double | single

フィルター遅延の最終状態。ベクトル、行列または多次元配列として返されます。

  • x ベクトルの場合、zf は長さ max(length(a),length(b))-1 の列ベクトルになります。

  • x が行列または多次元配列の場合、zf は長さ max(length(a),length(b))-1 の列ベクトルの配列になり、zf の列数は x の列数と等しくなります。たとえば、3×4×5 の配列 x の 2 番目の次元 (dim = 2) に沿って filter を使用するとします。配列 zf のサイズは、[max(length(a),length(b))-1]×3×5 です。

データ型: double | single

詳細

すべて折りたたむ

有理伝達関数

Z 変換領域でのベクトルの filter 操作の入出力記述が有理伝達関数です。有理伝達関数は次の形式になります。

Y(z)=b(1)+b(2)z1+...+b(nb+1)znb1+a(2)z1+...+a(na+1)znaX(z),

これは、有限インパルス応答 (FIR) フィルターと無限インパルス応答 (IIR) フィルターの両方を扱います[1]。ここで、X(z) は入力信号 x の Z 変換、Y(z) は出力信号 y の Z 変換、na はフィードバック フィルターの次数、nb はフィードフォワード フィルターの次数です。正規化により、a(1) = 1 と仮定します。

L 個の要素を持つ離散信号の場合、有理伝達関数は、次の差分方程式で表すこともできます。

a(1)y(L)=b(1)x(L)+b(2)x(L1)+...+b(nb+1)x(Lnb)a(2)y(L1)...a(na+1)y(Lna).

さらに有理伝達関数は、IIR デジタル フィルターの次の図に示すように直接形式 II 転置構成の実装で表すこともできます。この図では、na = nb = n–1 です。フィードバック フィルターとフィードフォワード フィルターの次数が異なる場合、つまり na ≠ nb の場合は、高次項を 0 として扱うことができます。たとえば、a = [1,2]b = [2,3,2,4] をもつフィルターの場合は、a = [1,2,0,0] を想定できます。

Block diagram that illustrates the direct-form II transposed implementation of an IIR digital filter with order n-1

サンプル点 m での filter の動作は、時間領域での差分方程式で与えられます。

y(m)=b(1)x(m)+w1(m1)w1(m)=b(2)x(m)+w2(m1)a(2)y(m)       =                 wn2(m)=b(n1)x(m)+wn1(m1)a(n1)y(m)wn1(m)=b(n)x(m)a(n)y(m).

既定では、関数 filter はフィルター遅延をゼロに初期化し、wk(0) = 0 となります。この初期化では、過去の入力と出力が両方ともゼロであると仮定します。現在のデータに過去の非ゼロの入力を含めるには、現在のデータの初期条件をフィルター遅延として指定します。フィルター遅延は、過去の入力 (および過去の出力) に同じ伝達関数を適用して取得する最終状態と考えることができます。filter の使用時に 4 番目の入力引数 zi を指定してフィルター遅延を設定でき、wk(0) = zi(k) となります。filter の使用時に 2 番目の出力引数 zf を指定して最終状態にアクセスすることもでき、wk(L) = zf(k) となります。

ヒント

  • FIR フィルターの係数 b を指定して関数 filter を使用するには、y = filter(b,1,x) を使用します。

  • Signal Processing Toolbox™ がある場合は、digitalFilter (Signal Processing Toolbox) オブジェクト d を指定し、y = filter(d,x) を使用して入力信号 x をフィルター処理します。周波数応答の仕様に基づいて d を生成するには、designfilt (Signal Processing Toolbox) を使用します。

  • フィルター処理関数の詳細については、デジタル フィルター処理 (Signal Processing Toolbox)を参照してください。

参照

[1] Oppenheim, Alan V., Ronald W. Schafer, and John R. Buck. Discrete-Time Signal Processing. Upper Saddle River, NJ: Prentice-Hall, 1999.

拡張機能

バージョン履歴

R2006a より前に導入