はじめに

先日、授業でMATLABを使って簡単な信号処理を行う演習があったのですが、公式以外で参考になる記事があまり見当たらなかったので備忘録がてらまとめておこうと思います。この記事では、ローパスフィルタのMATLAB実装のコードを載せています。

開発環境

Windows 10（64bit）
MATLAB R2015b（64bit）

ローパスフィルタ

ローパスフィルタとは、特定の周波数以上の信号の成分を低減させるためのフィルタです。詳細はWikipediaを見てください。
このフィルタは一定以上の周波数成分をカットすることができるのです。逆に言うと、どの周波数以上をカットするのかこちら側で決める必要があります。
どうやって決めるかですが、私の知る方法では（他にも良い方法があるかもしれません）信号をフーリエ変換してパワースペクトルのプロットから決定することが多いように思います。他には、信号の特徴（この信号はある周波数帯域に存在するとか）が分かっているときはその情報を使ったりします。以下、パワースペクトルをプロットする自作関数のコードです。

power_spectrum.m

%% フーリエ変換　→　パワースペクトル表示
% signal: 入力信号
% Fs    : サンプリング周波数
% TITLE : プロットに表示するタイトル
function powerSpectrum(signal, Fs, TITLE)
% 信号の長さを取得
L = length(signal);
% L点の離散フーリエ変換（目的次第で第2引数の入力を変更）
Y = fft(signal,L);
% パワーの計算
P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);

f = Fs*(0:(L/2))/L;

% パワースペクトルをプロット
figure('Name',TITLE,'NumberTitle','off')
plot(f,P1)
T = ['Signal-Sided Amplitude Spectrum of audio ' TITLE];
title(T)
xlabel('f (Hz)')
ylabel('|P1(f)|')
end

この関数を用いて、適当な音声信号をプロットすると次の左の図のようになります。

処理前	処理後

これを見ると明らかに5000Hzにノイズが入力されているのが分かります。これは作為的に作ったデータなので綺麗に5000Hzにだけノイズが入ってますが、普通はこんなプロットになることはないです。
さて、これでカットする必要のある周波数が分かったので、さっそくフィルタをかけていきます。以下、カットオフ周波数を指定してフィルタをかけ、処理後のデータを返してくれる自作関数のコードです。フィルタの性能も表示されます。ここでは、MATLABの標準関数で実装したので、butter関数を使用していますが、DSP System Toolboxが入っている場合はfdesign.lowpassという関数が既にあるので、こちらを使った方が次数でフィルタ設計するよりも狙ったフィルタ設計ができると思います。

lowpass.m

%% ローパスフィルタ
% 出力
% ret   : フィルタ処理後の信号データ
% 入力
% dataIn: 入力信号
% Fs    : サンプリング周期
% Fc    : カットオフ周波数
% order : フィルタの次数（大きいほど理想のフィルタに近づく）
function ret = lowpass(dataIn,Fs,Fc,order)
% フィルタ設計
% 
[b,a] = butter(order,Fc/Fs*2,'low');
freqz(b,a)
ret = filtfilt(b,a,dataIn);
end

フィルタをかけると処理後の結果が得られ、5000Hzの成分がちゃんとカットされていることが確認できます。音声信号の場合は、sound関数で音を聴いて確認もできますね。ちゃんとキーンとする高周波の音は除去されてました。
今回作成した自作の関数を用いたサンプルコードも載せておきます。適当に4000Hzをカットオフ周波数としています。

audio_filter.m

%% 音声信号処理
% filename : 音声ファイルの名前
function audio_filter(filename)
close all;

[audio_data,Fs] = audioread(filename);

% パワースペクトル表示
powerSpectrum(audio_data,Fs,filename);

% ローパスフィルタ
filtered_audio_data = lowpass(audio_data,Fs,4000,100);
sound(filtered_audio_data,Fs)

powerSpectrum(filtered_audio_data,Fs_sine,['filtered' filename]);

end

butter関数を用いたハイパスフィルタとバターワースフィルタの設計もlowpass.mの引数パラメータを少し変えるだけでできるので、参考にしてください。