filter

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

構文

y = filter(b,a,x)

y = filter(b,a,x,zi)

y = filter(b,a,x,zi,dim)

[y,zf] = filter(___)

説明

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

a(1) が 1 に等しくない場合、filter は a(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)=\frac{1}{windowSize}\left(x(n)+x(n-1)+...+x(n-(windowSize-1))\right).$$

ウィンドウ サイズが 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)=\frac{b(1)}{a(1)+a(2)z^{-1}}=\frac{1}{1-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 になります。

複数セクションによるデータのフィルター処理

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

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

x = randn(10000,1);

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

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

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

$$H(z)=\frac{b(1)+b(2)z^{-1}}{a(1)+a(2)z^{-1}}=\frac{2+3z^{-1}}{1+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);

y1 は x1 からフィルター処理されたデータで、y2 は x2 からフィルター処理されたデータです。フィルター処理された全体のシーケンスは、y1 と y2 の垂直連結です。

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

y = filter(b,a,x);

isequal(y,[y1;y2])

ans = logical
   1

入力引数

b — 有理伝達関数の分子係数

ベクトル

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

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

a — 有理伝達関数の分母係数

ベクトル

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

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

x — 入力データ

ベクトル | 行列 | 多次元配列

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

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

zi — フィルター遅延の初期条件

[] (既定値) | ベクトル | 行列 | 多次元配列

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

• 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
複素数のサポート: あり

dim — 演算の対象の次元

正の整数スカラー

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

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

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

  

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

  

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

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

出力引数

y — フィルター処理されたデータ

ベクトル | 行列 | 多次元配列

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

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

データ型: double | single

zf — フィルター遅延の最終状態

ベクトル | 行列 | 多次元配列

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

• 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)=\frac{b(1)+b(2)z^{-1}+\mathrm{...}+b(n_b+1)z^{-n_b}}{1+a(2)z^{-1}+\mathrm{...}+a(n_a+1)z^{-n_a}}X(z),$$

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

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

$$\begin{array}{c}a(1)y(L)=b(1)x(L)+b(2)x(L-1)+\mathrm{...}+b(n_b+1)x(L-n_b)\\ -a(2)y(L-1)-\mathrm{...}-a(n_a+1)y(L-n_a).\end{array}$$

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

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

$$\begin{array}{l}y(m)=b(1)x(m)+w_1(m-1)\\ w_1(m)=b(2)x(m)+w_2(m-1)-a(2)y(m)\\ ⋮=⋮⋮\\ w_{n-2}(m)=b(n-1)x(m)+w_{n-1}(m-1)-a(n-1)y(m)\\ w_{n-1}(m)=b(n)x(m)-a(n)y(m).\end{array}$$

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

拡張機能

tall 配列

メモリの許容量を超えるような多数の行を含む配列を計算します。

filter 関数は tall 配列をサポートしていますが、以下の使用上の注意および制限があります。

• dim > 1 の場合、2 出力の構文 [y,zf] = filter(___) はサポートされていません。

• digitalFilter (Signal Processing Toolbox) オブジェクト d を指定して y = filter(d,x) を使用する場合、filter 関数は入力信号 x で tall 配列をサポートしません。

詳細については、tall 配列を参照してください。

C/C++ コード生成

MATLAB® Coder™ を使用して C および C++ コードを生成します。

使用上の注意および制限:

• 入力引数 b と a では、次のようになります。

  • 入力ベクトルは、コード生成時に固定サイズまたは可変長のベクトルとして指定しなければなりません。ベクトルの最初の次元または 2 番目の次元を可変サイズにすることができます。他の次元はすべて固定サイズ 1 でなければなりません。

• 入力引数 x では、次のようになります。

  • 次元を指定しない場合、コード ジェネレーターは、可変サイズであるかサイズが 1 と等しくない入力配列の最初の次元に沿って演算を行います。この次元がコード生成時に可変サイズで実行時に 1 になる場合、ランタイム エラーが発生することがあります。このエラーを回避するには、次元を指定します。

• dim 引数を使用する場合、コード生成時に定数でなければなりません。

GPU コード生成

GPU Coder™ を使用して NVIDIA® GPU のための CUDA® コードを生成します。

使用上の注意および制限:

• 引数 b と a はコード生成時にベクトルでなければなりません。これは、各引数の最初の 2 つの次元のうちの少なくとも一方が固定サイズ 1 でなければならないことを意味します。

• 入力引数

  次元を指定しない場合、コード ジェネレーターは、可変サイズであるかサイズが 1 と等しくない入力配列の最初の次元に沿って演算を行います。この次元がコード生成時に可変サイズで実行時に 1 になる場合、ランタイム エラーが発生することがあります。このエラーを回避するには、次元を指定します。

• dim 引数を使用する場合、コード生成時に定数でなければなりません。

スレッドベースの環境

MATLAB® の backgroundPool を使用してバックグラウンドでコードを実行するか、Parallel Computing Toolbox™ の ThreadPool を使用してコードを高速化します。

filter 関数はスレッドベースの環境を完全にサポートしています。詳細については、スレッドベースの環境での MATLAB 関数の実行を参照してください。

GPU 配列

Parallel Computing Toolbox™ を使用してグラフィックス処理装置 (GPU) 上で実行することにより、コードを高速化します。

filter 関数は、GPU 配列を完全にサポートします。GPU 上で関数を実行するには、入力データを gpuArray (Parallel Computing Toolbox) として指定します。詳細については、GPU での MATLAB 関数の実行 (Parallel Computing Toolbox)を参照してください。

y = filter(d,x) 構文を使用して、GPU で 1 次元デジタル フィルターを適用することもできます。ここで、d は digitalFilter (Signal Processing Toolbox) オブジェクト、x は gpuArray です。digitalFilter オブジェクトの係数を gpuArray に変換する必要はありません。

分散配列

Parallel Computing Toolbox™ を使用して、クラスターの結合メモリ上で大きなアレイを分割します。

filter 関数は分散配列を完全にサポートしています。詳細については、分散配列を使用した MATLAB 関数の実行 (Parallel Computing Toolbox)を参照してください。

バージョン履歴

R2006a より前に導入

R2026a: 有限インパルス応答 (FIR) フィルターのパフォーマンスが向上

filter 関数で FIR フィルターのパフォーマンスが向上しています。FIR フィルターは、有限インパルス応答時間と、零点のみをもち極をもたない z 領域の有理伝達関数によって特徴付けられます。パフォーマンスの向上は、伝達関数に多数の係数がある場合に最も顕著になります。

たとえば、次のコードでは、伝達関数の分子に 10,001 個の係数があるローパス フィルターを使用して、長さが 1,000,000 の入力信号をフィルター処理しています。直前のリリースと比較して、このコードは約 2.3 倍速くなっています。

function ts = timingFilter
t = 1:1000000;
x = sin(5e-5*t) + sin(10*t) + 0.1*rand(size(t));

N = 10000;
freq = 0.5;
k = -(N/2):(N/2);
h_ideal = sin(freq*pi*k)./(pi*k);
h_ideal(k==0) = freq;
w = 0.54 - 0.46*cos(2*pi*(0:N)/N);
b = h_ideal.*w;

y = @() filter(b,1,x);
ts = timeit(y);
end

おおよその実行時間は以下のとおりです。

R2025b: 0.28 秒

R2026a: 0.12 秒

このコードの時間測定では、Windows^® 11、AMD EPYC™ 74F3 24 コア プロセッサ (3.19 GHz) 搭載のテスト システムで、timingFilter 関数を呼び出しました。

参考

conv | filter2

トピック

• データのフィルター処理

授業用リソース

• Digital Signal Processing: Signals and Filter Design (MathWorks Teaching Resources)