周波数領域での線形フィルターの設計

このトピックでは、周波数領域でフィルター処理を実行する関数について説明します。空間領域でのフィルターの設計方法の詳細については、空間領域でのイメージのフィルター処理とはを参照してください。

2 次元有限インパルス応答 (FIR) フィルター

Image Processing Toolbox™ ソフトウェアは、1 つのクラスの線形フィルターをサポートしています。それは、 2 次元有限インパルス応答 (FIR) フィルターです。FIR フィルターは、単点またはインパルスに対して有限の範囲をもっています。すべての Image Processing Toolbox フィルター設計関数が FIR フィルターを返します。

FIR フィルターは、MATLAB^® 環境で最適なイメージ処理を行ういくつかの特徴をもっています。

• FIR フィルターは、係数の行列として表現することが簡単です。

• 2 次元の FIR フィルターは、1 次元の FIR フィルターを自然に拡張したものです。

• FIR フィルター設計には、よく知られた信頼性の高い方法が複数存在します。

• FIR フィルターは実装が簡単です。

• FIR フィルターは、歪みを避けるために、線形位相で設計できます。

別のクラスのフィルター、無限インパルス応答 (IIR) フィルターは、それほどイメージ処理アプリケーションには適していません。FIR フィルターの設計や実現が簡単で、安定性があるのに対し、IIR フィルターは、この部分に欠陥があります。そのためこのツールボックスでは、IIR フィルターをサポートしていません。

メモ

この節で説明しているほとんどの設計法は、Signal Processing Toolbox™ の関数で作成される 1 次元のフィルターまたはウィンドウから 2 次元フィルターを作成する方法を使用しています。このツールボックスは必須ではありませんが、Signal Processing Toolbox ソフトウェアがない場合はフィルターの設計が困難になることがあります。

1 次元フィルターの周波数変換を使用した 2 次元フィルターの作成

この例では、関数 ftrans2 を使用して 1 次元 FIR フィルターを 2 次元 FIR フィルターに変換する方法を説明します。特定の特性をもつ 1 次元フィルターは、対応する 2 次元フィルターよりも簡単に設計できるため、この関数は役に立ちます。周波数変換法では、特に遷移帯域幅とリップル特性など、1 次元フィルターのほとんどの特性が保存されます。1 次元周波数応答の形状は、2 次元応答でも明らかです。

この関数は、周波数変換を定義する一連の要素である "変換行列" を使用します。この関数の既定の変換行列は、ほぼ円対称なフィルターを作成します。自分で変換行列を定義することで、別の対称を得ることが可能です (詳細については、Jae S. Lim,『Two-Dimensional Signal and Image Processing』, 1990 を参照してください)。

Signal Processing Toolbox の関数 firpm を使用して 1 次元 FIR フィルターを作成します。

b = firpm(10,[0 0.4 0.6 1],[1 1 0 0])

b =

  Columns 1 through 9

    0.0537   -0.0000   -0.0916   -0.0001    0.3131    0.4999    0.3131   -0.0001   -0.0916

  Columns 10 through 11

   -0.0000    0.0537

1 次元フィルターを 2 次元フィルターに変換します。

h = ftrans2(b);

h =

  Columns 1 through 9

    0.0001    0.0005    0.0024    0.0063    0.0110    0.0132    0.0110    0.0063    0.0024
    0.0005    0.0031    0.0068    0.0042   -0.0074   -0.0147   -0.0074    0.0042    0.0068
    0.0024    0.0068   -0.0001   -0.0191   -0.0251   -0.0213   -0.0251   -0.0191   -0.0001
    0.0063    0.0042   -0.0191   -0.0172    0.0128    0.0259    0.0128   -0.0172   -0.0191
    0.0110   -0.0074   -0.0251    0.0128    0.0924    0.1457    0.0924    0.0128   -0.0251
    0.0132   -0.0147   -0.0213    0.0259    0.1457    0.2021    0.1457    0.0259   -0.0213
    0.0110   -0.0074   -0.0251    0.0128    0.0924    0.1457    0.0924    0.0128   -0.0251
    0.0063    0.0042   -0.0191   -0.0172    0.0128    0.0259    0.0128   -0.0172   -0.0191
    0.0024    0.0068   -0.0001   -0.0191   -0.0251   -0.0213   -0.0251   -0.0191   -0.0001
    0.0005    0.0031    0.0068    0.0042   -0.0074   -0.0147   -0.0074    0.0042    0.0068
    0.0001    0.0005    0.0024    0.0063    0.0110    0.0132    0.0110    0.0063    0.0024

  Columns 10 through 11

    0.0005    0.0001
    0.0031    0.0005
    0.0068    0.0024
    0.0042    0.0063
   -0.0074    0.0110
   -0.0147    0.0132
   -0.0074    0.0110
    0.0042    0.0063
    0.0068    0.0024
    0.0031    0.0005
    0.0005    0.0001

フィルターの周波数応答を表示します。

[H,w] = freqz(b,1,64,"whole");
colormap(jet(64))
plot(w/pi-1,fftshift(abs(H)))
figure, freqz2(h,[32 32])

1 次元の周波数応答

対応する 2 次元の周波数応答

周波数サンプリング法を使用したフィルターの作成

この例では、周波数サンプリング法を使用して、必要な周波数応答をベースに 2 次元フィルターを作成する方法を説明します。

周波数サンプリング法では、周波数応答の形状を設定する点を行列として与え、これらの点を周波数応答が通るようにフィルターを作成します。周波数サンプリングでは、与えられた点の間での周波数応答の動作に制約がありません。通常、応答はこの部分でリップルになります。(リップルは定数値周りの振動量です。実際のフィルターの周波数応答は、理想的なフィルターの周波数応答がフラットになる部分で、リップルをもつことがよくあります)。

ターゲットとする、11 行 11 列フィルターの 2 次元周波数応答を定義し、表示します。

Hd = zeros(11,11); 
Hd(4:8,4:8) = 1;
[f1,f2] = freqspace(11,"meshgrid");
mesh(f1,f2,Hd)
colormap(jet(64))
title("Desired Frequency Response")

関数 fsamp2 を使用して、ターゲットの周波数応答に基づいてフィルターを作成します。fsamp2 は、入力行列 Hd の点を通過する周波数応答をもつフィルター h を返します。

h = fsamp2(Hd);

フィルターの周波数応答をプロットします。

freqz2(h,[32 32])
title("Actual Frequency Response")

実際の周波数応答のリップルを理想の周波数応答と比較して見てください。これらのリップルは、周波数サンプリング設計法の基本的な問題です。これらのリップルは、必要な応答内の鋭い遷移部で必ず生じます。

フィルターを大きくすればするほど、リップルの空間的な広がりは小さくなります。しかしフィルターを大きくすると、リップルの高さは低くなりますが、フィルター処理の計算時間は多くなります。必要な周波数応答に対してよりスムーズな近似を行うには、周波数変換法またはウィンドウ処理法を使用することを検討してください。

ウィンドウ処理法

ウィンドウ処理法は、理想的なインパルス応答にウィンドウ関数を乗算して対応するフィルターを生成し、理想的なインパルス応答にテーパーを適用します。周波数サンプリング法と同様、ウィンドウ処理法は必要な周波数応答を近似するフィルターの周波数応答を作成します。しかしウィンドウ処理法は、周波数サンプリング法よりも良い結果を出す傾向があります。

ツールボックスには、ウィンドウをベースにしたフィルター設計に 2 つの関数 fwind1 と fwind2 が用意されています。fwind1 は、ユーザーが指定した 1 つまたは 2 つの 1 次元ウィンドウから作成された 2 次元のウィンドウを使うことにより、2 次元フィルターを設計します。fwind2 は、指定した 2 次元ウィンドウを直接使用して 2 次元フィルターを設計します。

fwind1 は、使用する 2 次元ウィンドウを作成する方法として、2 つの異なる方法をサポートしています。

• 1 つの 1 次元ウィンドウを変換し、回転と同様の処理を使用して、ほぼ環状に対称な 2 次元ウィンドウを作成します。

• 外積を計算することにより、 2 つの 1 次元ウィンドウから四角形の分離可能なウィンドウを作成します。

次の例は、fwind1 を使用して、必要な周波数応答 Hd から 11 行 11 列のフィルターを作成します。この例では Signal Processing Toolbox の関数 hamming を使用して、1 次元ウィンドウを作成します。fwind1 はこれを 2 次元のウィンドウに拡張します。

Hd = zeros(11,11); Hd(4:8,4:8) = 1;
[f1,f2] = freqspace(11,'meshgrid');
mesh(f1,f2,Hd), axis([-1 1 -1 1 0 1.2]), colormap(jet(64))
h = fwind1(Hd,hamming(11));
figure, freqz2(h,[32 32]), axis([-1 1 -1 1 0 1.2])

必要な 2 次元周波数応答 (左) と実際の 2 次元周波数応答 (右)

必要な周波数応答行列の作成

フィルター設計関数 fsamp2、fwind1、および fwind2 は、必要な周波数応答の大きさの行列をベースにフィルターを作成します。周波数応答は、さまざまな入力周波数に対するフィルターのゲインを記述する数学関数です。

必要な周波数応答行列を作成するには、関数 freqspace を使用します。freqspace は、任意の大きさの応答に対して等間隔の正しい周波数値を返します。freqspace によって返されるもの以外の周波数点を使用して必要な周波数応答行列を作成すると、非線形位相のような予期せぬ結果が生じます。

たとえば、0.5 のカットオフを持つ円の理想的なローパス周波数応答を作成するには、以下を使用します。

[f1,f2] = freqspace(25,"meshgrid");
Hd = zeros(25,25); d = sqrt(f1.^2 + f2.^2) < 0.5;
Hd(d) = 1;
mesh(f1,f2,Hd)

理想的な円のローパス周波数応答

周波数応答に対して、fsamp2、fwind1、fwind2 によって作成されるフィルターは実数になることに注意してください。この結果は、ほとんどのイメージ処理アプリケーションに適しています。一般的にこれを達成するために必要な周波数応答は、周波数の原点 (f1 = 0, f2 = 0) に関して対称になります。

フィルターの周波数応答の計算

関数 freqz2 は、2 次元フィルターに対する周波数応答を計算します。出力引数を設定しないと、freqz2 は、周波数応答のメッシュ プロットを作成します。たとえば、次の FIR フィルターについて考えてみましょう。

h =[0.1667    0.6667    0.1667
    0.6667   -3.3333    0.6667
    0.1667    0.6667    0.1667];

次のコマンドは、h の 64 x 64 点の周波数応答を計算し、プロットします。

freqz2(h)

2 次元フィルターの周波数応答

周波数応答行列 H と周波数点ベクトル f1 および f2 を得るには、出力引数を使用します。

[H,f1,f2] = freqz2(h);

freqz2 は、1.0 の値がサンプリング周波数の 1/2、あるいは π ラジアンに対応するように周波数 f1 と f2 を正規化します。

簡単な m 行 n 列の応答に対しては、上に示すように、freqz2 は 2 次元の高速フーリエ変換関数 fft2 を使用します。また任意の周波数点ベクトルも指定できますが、この場合、freqz2 は時間がかかるアルゴリズムを使用します。

高速フーリエ変換、および線形フィルター処理、フィルター設計への応用の詳細については、フーリエ変換を参照してください。

参考

freqspace | freqz2 | fsamp2 | ftrans2 | fwind1 | fwind2

×

例を開く

この例の変更されたバージョンがあります。編集された方の例を開きますか?

いいえ、変更されたバージョンを上書きします。 はい