# EmbVision(CQ版) チュートリアル（２）

フィルタリングと周波数特性

新潟大学 村松　正吾，高橋　勇希

Copyright (c), All rights reserved, 2014-2025, Shogo MURAMATSU and Yuki TAKAHASHI

## Contents

- 一次元信号の周波数特性
- 一次元信号のフィルタリング
- 一次元フィルタの周波数応答
- 二次元信号の周波数特性
- 二次元信号のフィルタリング
- 二次元フィルタの周波数応答
- 演習課題

Part1 | | メニュー | | Part3

概要

本演習では、MATLABにて一次元信号のフィルタリングと周波数解析法、 画像情報のフィルタリングと周波数解析法について簡単に学ぶ。

準備として、開いている全ての Figure を close 関数で 閉じておく。

    close all

## 一次元信号の周波数特性

まず、予め用意されているオーディオデータ chirp を load 関数を利用して読み出して準備する。

    load chirp

オーディオデータは変数 y に倍精度実数型配列として保持される。 なお、標本化周波数は変数 Fs に保持される。

    whos y Fs

      Name          Size             Bytes  Class     Attributes

      Fs            1x1                  8  double
      y         13129x1             105032  double

次に、予め用意されているオーディオデータ gong を読み出し、 length 関数を利用して長さを揃え、 chirp のデータと混合する。

    c = y; % 変数 c に代入
    load gong
    g = y; % 変数 g に代入

    x = g(1:length(c))+c; % 長さを調整して混合

変数 x には混合したオーディオデータが保持される。 plot 関数で確認しよう。

    plot(x)

オーディオ再生環境があれば audioplayer 関数を利用して、オーディオ再生も可能である。

まず、オーディオ再生オブジェクトを生成する。

    player = audioplayer(x,Fs);
    whos player

オーディオ再生は、 play メソッドに オブジェクト player を指定して実行される。

    play(player)

spectrogram 関数を利用することで、 短時間フーリエ変換を利用したオーディオデータの周波数解析を実行できる。

- 窓長: 256

として実行しよう。

    figure(1)
    spectrogram(x,256)
    caxis([ -70 10 ])
    colorbar

上記の操作により周波数特性（スペクトログラム）が表示される。 横軸は正規化周波数( を 1 と正規化)、 縦軸は標本インデックス( 秒単位)である。

なお、値の大きさ[dB] が分かり易いように caxis 関数と colorbar 関数を併用した。

view 関数で視点を変えてみよう。

    view([-15 30])

zlim 関数を利用して、 Z軸の座標を　-70 から 10 に調整する。

    zlim([ -70 10 ])

[ トップ ]

## 一次元信号のフィルタリング

次に、オーディオデータ x に、線形フィルタ処理

を施してみよう。ここで、

- : フィルタ入力
- : フィルタ出力
- : フィルタ係数（インパルス応答）
- : フィルタ次数

とする。

MATLAB では線形フィルタ処理に conv 関数を利用できる。

フィルタ係数

として線形フィルタ処理を実行してみよう。

    h = [ 1 1 1 ]/3;
    figure(2)
    impz(h)

    y = conv(h,x);

変数 y にはフィルタ処理結果が保持されている。 出力 y の周波数特性（スペクトログラム）を確認しよう。

    figure(3)
    spectrogram(y,256);
    caxis([ -70 10 ])
    colorbar

    view([ -15 30 ])
    zlim([ -70 10 ])

入力 x と出力 y のスペクトログラムを比較してみて欲しい。 どのようなことに気が付くだろうか？

- 大よそ、正規化周波数0.4以上の高周波成分が減衰している。
- 特に、0.667付近の減衰が大きい。

ということに注意して観察して欲しい。

なお、処理結果をオーディオ再生により確認する場合は、

    player = audioplayer(y,Fs);
    play(player)

と実行すればよい。

[ トップ ]

## 一次元フィルタの周波数応答

線形フィルタによる周波数特性の変化は、 フィルタの周波数応答により確認できる。

何故ならば、時間領域での畳込み演算は

のように周波数（DTFT: 離散時間フーリエ変換）領域では 積演算に対応するためである。ここで、

- : 入力 の周波数特性
- : 出力 の周波数特性
- : フィルタ係数（インパルス応答） の周波数応答

である。

フィルタ係数 の周波数応答は freqz 関数 により確認できる。

    figure(4)
    freqz(h)

振幅応答を確認すると、正規化周波数 0.3 から　0.4 の間から高域に渡り -6 [dB] 以上の減衰特性をもち、特に 0.6 から 0.7 の間で大きく減衰する特性が 確認できる。

なお、

より、

- で、
- で、

となることが確認できる。

[ トップ ]

## 二次元信号の周波数特性

では、画像データに対するフィルタリングに話を進めよう。 先に示したオーディオデータと同様に線形フィルタ処理を施すことができる。

まず、 開いている全ての Figure を close 関数で で閉じた後、予め用意されている画像データ cameraman を 読み込んで表示しよう。

    close all
    figure(1)
    X = imread('cameraman.tif');
    imshow(X)

画像データは変数 X に符号なし整数８ビット型配列として保持されるので、 これを倍精度実数型に変換しよう。

    X = im2double(X);
    whos X

      Name        Size              Bytes  Class     Attributes

      X         256x256            524288  double

fft2 関数を利用することで、 画像データの二次元の周波数解析を実行できる。 画像サイズに合わせて

- 二次元FFT点数:

と設定し二次元周波数解析を実行しよう。

    F = fft2(X,256,256);
    whos F

      Name        Size               Bytes  Class     Attributes

      F         256x256            1048576  double    complex

変数 F に二次元離散フーリエ変換(DFT)係数が保持される。 なお、複素数として結果が得られるため、絶対値をとって 振幅特性を求めてみよう。

    Fmag = abs(F);
    whos Fmag

      Name        Size              Bytes  Class     Attributes

      Fmag      256x256            524288  double

変数 Fmag には、振幅特性として実数配列が保持される。 surface 関数を利用して、 特性を可視化しよう。

surface プロットを調整するためのハンドルオブジェクトとして 変数 hsrfc を用意しておこう。

    figure(2)
    hsrfc = surface(20*log10(Fmag));
    hsrfc.EdgeColor = 'none';

ここでは、デシベルに換算している点に注意する。

グラフが見やすいようカラーバーを設置する。

    colorbar
    caxis([ -20 80 ])

中心が直流となるように fftshift 関数を用いて 配列をシフトする。

    hsrfc.ZData = fftshift(hsrfc.ZData); % Z軸のシフト
    hsrfc.CData = fftshift(hsrfc.CData); % カラー軸のシフト

正規化周波数となるよう座標の調整を行う。

    fstep = 1/256; % 周波数標本点の間隔
    hsrfc.XData = -0.5:fstep:0.5-fstep;
    hsrfc.YData = -0.5:fstep:0.5-fstep;
    xlabel('F_x (\times\pi rad/sample)')
    ylabel('F_y (\times\pi rad/sample)')

視点を変える。

    view([ -15 30 ])
    zlim([ -20 80 ])

[ トップ ]

## 二次元信号のフィルタリング

次に、画像データ X に、二次元線形フィルタ処理

を施してみよう。ここで、

- : フィルタ入力
- : フィルタ出力
- : フィルタ係数（インパルス応答）
- : 垂直フィルタ次数
- : 水平フィルタ次数

とする。

MATLAB では二次元線形フィルタ処理に conv2 関数もしくは imfilter 関数を利用できる。

フィルタ係数 として配列

を利用して線形フィルタ処理を実行してみよう。

    H = [ 1 1 1 ; 0 0 0 ; -1 -1 -1 ];
    Y = conv2(H,X);

変数 Y にはフィルタ処理結果が保持されている。 ただし、

    size(Y)

    ans =

       258   258

のように、もとの配列 X よりもサイズが縦横 2 画素づつ増加している。 これは、サイズ の画像にサイズ の線形フィルタをかけるとその出力のサイズが

に増加する性質による。

上下左右、1画素ずつ削って、入力画像 X のサイズに出力画像 Y のサイズを 調整しよう。 end 関数を利用した配列インデックスの 最大値指定を利用すると便利である。

    Y = Y(2:end-1,2:end-1);

また、出力画像 Y は、演算の結果、

    min(Y(:))

    ans =

       -2.7882

のように、負の値を含む。このため画像として表示する際には工夫が必要である。

実数型画像の場合、imshow 関数は、画素値が 0 から 1 にスケールされていると 仮定して表示を行うので、負の値が0.5 以下、正の値が0.5以上となるよう Y の値を調整する。

    figure(3)
    imshow(Y+0.5)

先のフィルタは、垂直方向の微分 の離散近似を出力する。

出力 Y の周波数特性を確認しよう。

    figure(4)
    F = fft2(Y,256,256);
    Fmag = abs(F);
    hsrfc = surface(20*log10(Fmag));
    hsrfc.EdgeColor = 'none';
    colorbar
    caxis([ -20 80 ])
    hsrfc.ZData = fftshift(hsrfc.ZData);
    hsrfc.CData = fftshift(hsrfc.CData);
    hsrfc.XData = -0.5:fstep:0.5-fstep;
    hsrfc.YData = -0.5:fstep:0.5-fstep;
    xlabel('F_x (\times\pi rad/sample)')
    ylabel('F_y (\times\pi rad/sample)')
    view([ -15 30 ])
    zlim([ -20 80 ])

入力 X と出力 Y の周波数特性を比較してみて欲しい。 どのようなことに気が付くだろうか？

- 直流におけるピークがなくなる。
- 水平方向の高域に減衰が見られる。

ということに注意して観察して欲しい。

[ トップ ]

## 二次元フィルタの周波数応答

二次元線形フィルタによる周波数特性の変化は、 一次元の場合と同様に、フィルタの周波数応答により確認できる。

何故ならば、空間領域での畳込み演算は

のように周波数（DSFT: 離散空間フーリエ変換）領域では 積演算に対応するためである。ここで、

- : 入力 の周波数特性
- : 出力 の周波数特性
- : フィルタ係数（インパルス応答） の周波数応答

である。

フィルタ係数 の周波数応答は freqz2 関数により確認できる。

    figure(5)
    freqz2(H)

振幅応答を確認すると、直流と水平方向の高周波数成分に対する減衰特性をもち、 垂直方向については帯域通過特性をもつことが確認できる。

なお、

より、

- で、
- で、
- で、
- で、

となることが確認できる。

[ トップ ]

## 演習課題

課題2-1. 水平微分フィルタ

フィルタ係数 として配列

（水平微分の離散近似フィルタ）を用意し、 画像ファイル cameraman.tif のグレースケール画像に対して 線形フィルタ処理を施し、処理結果を画像ファイル cameramangradx.tif に保存せよ。 （負の値を考慮して、値0.5 によりかさ上げすること）。

また、フィルタの周波数特性をグラフで確認せよ。

(処理例）

課題2-2. 勾配の大きさと偏角

垂直微分フィルタ出力 と 水平微分フィルタ出力 から、

- 勾配の大きさ :
- 勾配の方向　 :

を計算し、処理結果をそれぞれ画像ファイル cameramangradmag.tif と cameramangradang.tif に保存せよ。ただし、勾配の方向については、 値の範囲を 0 から 1 に換算すること。

(処理例）

Part1 | | メニュー | | トップ | | Part3