周波数フィルタの作成

脳波・脳磁図の波形データが読み込めたとしても、生波形には通常生体信号以外の雑音も含まれていて、周波数フィルタによる信号処理が必要となります。 MATLABを使えば簡単にフィルタをかけることができます。 sptool関数を使えばGUIでフィルタが作れますが、まずは自分でフィルタをかけれるようになりましょう。平成24年9月20日の広大霞キャンパス基礎棟2Fで行ったMATLABセミナーの資料の一部です。参加人数は7人くらいでした。

1～50Hzのチャープ信号の作成
チャープ信号の表示
無限インパルス応答（IIR）・楕円フィルタの作成
IIRフィルタの周波数応答の表示
IIRフィルタ前後の波形表示
IIRフィルタの時間周波数変化表示
有限インパルス応答（FIR）・最小二乗線形フィルタ
FIRフィルタの周波数応答の表示
FIRフィルタ前後の波形表示
FIRフィルタの時間周波数変化表示

1～50Hzのチャープ信号の作成

chirp関数を使って1～50Hzに時間とともに変調するチャープ信号を作成します。

sfreq=600;% サンプリング周波数
t=linspace(0,5,5*sfreq);% 5秒
fstart=1;fend=50;% １～50Hz
y=chirp(t,fstart,t(end),fend);% チャープ信号

チャープ信号の表示

figure;
subplot(2,1,1);% 波形表示
plot(t,y);xlabel('sec');ylabel('Amp');
win=2^8;% window
subplot(2,1,2);% 時間周波数変化
[S,F,T,P]=spectrogram(y,win,win/2,win,sfreq);
imagesc(T,F,P);
set(gca,'YDir','normal','ylim',[fstart,fend]);
xlabel('sec');ylabel('Hz');

無限インパルス応答（IIR）・楕円フィルタの作成

ellip関数を使って楕円関数を作成します。 ellipord関数を使えば最低次数を出力するらしいのですが、不慣れなので次数は決め打ちで３にしてます。

bandpass=[8,14]; % 通過周波数帯域
n=3;Rp=0.1;Rs=60;% 次数は3
[b,a]=ellip(n,Rp,Rs,bandpass*2/sfreq);

IIRフィルタの周波数応答の表示

freqz関数を使います。

figure;freqz(b,a,1000,sfreq);
set(findobj(gcf,'type','axes'),'xlim',[0,30]);

IIRフィルタ前後の波形表示

フィルタ処理後は次数+α分位相がズレるのでfiltfilt関数でズレを補正します。

figure;
for n=1:3
    subplot(3,1,n);
    switch n
        case 1;yy=y;txt='元波形';
        case 2;yy=filter(b,a,y);txt='フィルタ後波形';
        case 3;yy=filtfilt(b,a,y);txt='位相補正後フィルタ波形';
    end;
    plot(t,yy);title(txt);
    xlabel('time');ylabel('amp');
end;
set(findobj(gcf,'type','axes'),'ylim',[-2,2]);

IIRフィルタの時間周波数変化表示

figure;
win=2^8;
for n=1:3
    subplot(3,1,n);
    switch n
        case 1;yy=y;txt='元データ';
        case 2;yy=filter(b,a,y);txt='フィルタ処理後';
        case 3;yy=filtfilt(b,a,y);txt='位相補正＋フィルタ処理後';
    end;
    [S,F,T,P]=spectrogram(yy,win,win/2,win,sfreq);
    if n==1
        clim=[min(P(:)),max(P(:))];
    end;
    imagesc(T,F,P);title(txt);
    xlabel('sec');ylabel('Hz');
end;
h=findobj(gcf,'type','axes');
set(h,'YDir','normal','ylim',[fstart,fend],'clim',clim);

有限インパルス応答（FIR）・最小二乗線形フィルタ

FIRもいろいろあるんですが、firls関数を使ってみました。 FIRは分母が1なのでIIRよりも次数を増やさないといいフィルタになりません。

bandpass=[8,14]; % 通過周波数帯域
w=0.05*[-1,1];
f=[0,bandpass(1)+w,bandpass(2)+w,sfreq/2]/(sfreq/2);
a=1;% FIRなので分母は1
n=128;% 次数は128
bp=[0,0,1,1,0,0];
b=firls(n,f,bp);

FIRフィルタの周波数応答の表示