FIRフィルタに関して

参考文献及びWEB

フィルタとは？

フィルタは必要な信号を取り出すものです。

ローパスフィルタ(low pass filter,LPF)@低い周波数を通過させる
ハイパスフィルタ(high pass filter,HPF)@高い周波数を通過させる
バンドパスフィルタ(band pass filter,BPF)@ある帯域を通過させる
バンドエリミネーションフィルタ(band elimination filter,BEF)@ある帯域以外を通過させる

振幅特性とは？

フィルタに入力した信号の振幅が、何倍になって出力されるかを表現したもの。

通過域では1倍 -> 信号はそのまま通過する
阻止域では0倍 -> 信号は除去される

通過域と阻止域との境界にあたる周波数を、遮断周波数 f_c と呼びます。

FIRフィルタとは？

ディジタルフィルタには大きく分けて2種類のフィルタがあります。

・FIRフィルタ: FIR(Finite Impulse Response 有限インパルス応答)システムを応用したフィルタ
・IIRフィルタ: IIR(Infinite Impulse Response 無限インパルス応答)システムを応用したフィルタ

FIRフィルタの形状

FIRフィルタの形状は、畳込みについてで示された形状そのものです。

$y[n]= \sum_{k=0}^{N-1} (h[k] x[n-k])$

= h[0]x[n-0] + h[1]x[n-1] + h[2]x[n-3] + ... + h[N-1]x[n-(N-1)]

ここで、
フィルタの係数(タップ数)は個ある場合、 h[N] で表されます。

また、フィルタを計算するためにはフィルタ係数分だけの過去の入力値を保存しなければなりません。
(現在の入力値と過去の入力値を含めてフィルタ係数分です。)

そこで、
x[N] 個の入力バッファが必要となります。

FIRの意味は有限時間内にインパルス応答が終了するということです。

畳込みについてで説明したとおり、フィルタ係数はある入力に対する出力の予想を示しています。
ですからフィルタの場合はある入力をいれたとき、いらない信号は除去するような係数を求めること
になります。

係数の個数分のデータが入力として与えられたとき、出力は正常に出力されます。
係数の個数分より少ないデータが入力として与えられたとき、それは過度状態と等価と考えられます。

例えば、係数3で、データが1つしかない場合、
現在は

(1) 現在の予測と
(2) 1つ前の過去からの現在の予測と
(3) 2つ前の過去からの現在の予測

からなりたっているため、

現在の予測しかできず、
1つ前及び2つ前の予測が定まらないためこの状態を過度状態と呼びます。

またフィルタ係数が h[3] の3個だとして、
データが x[10] の10個だとすると

y[11]= h[0]*x[11] + h[1]*x[10] + h[2]*x[9];
     = h[0]*0     + h[1]*0     + h[2]*x[9];
     = h[2]*x[9];

またフィルタ係数が h[5] の5個だとして、
データが x[10] の10個だとすると

y[13]= h[0]*x[13] + h[1]*x[12] + h[2]*x[11] + h[3]*x[10] + h[4]*x[9];
     = h[0]*0     + h[1]*0     + h[2]*0     + h[3]*0     + h[4]*x[9];
     = h[4]*x[9];

つまり最後のデータから係数N-1個余計に出力することができます。この場合は、

係数が3つの場合は、x[9] + 2つ分 = x[11]
係数が5つの場合は、x[9] + 4つ分 = x[13]

となります。

まとめると
(1)係数個分入力データを入れたときから正しい出力がでる。
(2)最後のデータから「係数-1」分余計に出力して出力が終わる。

をフィルタのタップ数
N-1 をフィルタの次数
と呼びます。

直流除去フィルタを設計してみる

MatlabでFilter Design & Analysis Toolを起動してみる

Matlabで

>fdatool

とすると画面のようなFilter Design & Analysis Toolが現れる

周波数-群遅延特性を見る

出力が入力より125サンプル遅れることが分かる

125 x 1/250 = 0.5[s] 遅れる。

係数をヘッダとして出力する

ターゲット→Cヘッダを作成
なぜか251個の係数ができる

◆ヘッダ
Head Code

FIRフィルタを作成し、フィルタ特性どおりの出力かを確認する

filterTest.zip

コードを書く

■ヘッダ

class FIRFilter{
private:
 double *hn;
 double *x;
 int    coff_size;
public:
 FIRFilter(double *coff,int coff_size);
 double doFilter(double in);
};

■コード

FIRFilter::FIRFilter(double* coff,int coff_size){
 x = new double[coff_size];
 hn = new double[coff_size];
 this->coff_size = coff_size;
  
 for(int i=0;i<coff_size;i++){
   hn[i]=coff[i];
   x[i]=0;
 }
}

double FIRFilter::doFilter(double in){
 double y=0;
 x[0]=in;
	
 for(int i=0;i<coff_size;i++){
  y=y+hn[i]*x[i];
 }

 for(int i=coff_size-1;i>0;i--){
  x[i]=x[i-1];
 }
 return y;
}

■使用方法

FIRFilter fir(COFF,COFF_SIZE);
for(int i=0;i<10000;i++){
 double time=i*0.004;
 double data=sin.nextWave();
 
 sprintf_s(str,"%f %f\n",time,data);
 m_OFStrm1 << str;
		
 double output=fir.doFilter(data);
 sprintf_s(str,"%f %f\n",time,output);
 m_OFStrm2 << str;
}

結果の検討

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  check program
%                 programming by embedded.samurai
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
echo off
clear all
close all

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 生データを表示する
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1);
subplot(2,1,1);
%//250Hz
dt=0.004;
load 'result1.txt'
[row,column]=size(result1)

time=result1(:,1);
data=result1(:,2);

% 生データを表示する
plot(time,data,'r-','linewidth',2);
hold on;

xhani=500*dt; 

axis([0 xhani -1 1]);

h=gca;
set(h,'LineWidth',2,...,
'FontSize',15,'FontName','Century');

grid on; 

% キャプション
h=legend('正弦波',...,
'3',...,
2);
set(h,'FontSize',15,'FontName','MSゴシック'); 

xlabel('Time[s]','Fontsize',15,'FontName','Century');
ylabel('Amplitude','Fontsize',15,'FontName','Century')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(2,1,2);

fftdata=fft(data);
power=abs(fftdata).^2;

freqdata=[0:row-1]'/dt/row

plot(freqdata,power,'b-','linewidth',2);
hold on;

axis([0 50 -Inf Inf]);

h=gca;
set(h,'LineWidth',2,...,
'FontSize',15,'FontName','Century');

grid on;

% キャプション
h=legend('正弦波の周波数',...,
'2',...,
1);
set(h,'FontSize',15,'FontName','MSゴシック');

xlabel('Frequency[Hz]','Fontsize',15,'FontName','Century');
ylabel('Amplitude','Fontsize',15,'FontName','Century')

print -djpeg sin_wave.jpg 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(2);
subplot(2,1,1);
%//250Hz
dt=0.004;
load 'result2.txt'
[row,column]=size(result1) 

time2=result2(:,1);
data2=result2(:,2);

% 生データを表示する
plot(time2,data2,'r-','linewidth',2);
hold on;

xhani=500*dt;

axis([0 xhani -1 1]);

h=gca;
set(h,'LineWidth',2,...,
'FontSize',15,'FontName','Century');

grid on;

% キャプション
h=legend('正弦波(フィルタ通過後)',...,
'3',...,
4);
set(h,'FontSize',15,'FontName','MSゴシック'); 

xlabel('Time[s]','Fontsize',15,'FontName','Century');
ylabel('Amplitude','Fontsize',15,'FontName','Century')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(2,1,2);

fftdata2=fft(data2);
power2=abs(fftdata2).^2;

freqdata=[0:row-1]'/dt/row

plot(freqdata,power2,'b-','linewidth',2);
hold on;

axis([0 50 -Inf Inf]);

h=gca;
set(h,'LineWidth',2,...,
'FontSize',15,'FontName','Century');

grid on;

% キャプション
h=legend('正弦波の周波数(フィルタ通過後)',...,
'2',...,
1);
set(h,'FontSize',15,'FontName','MSゴシック'); 

xlabel('Frequency[Hz]','Fontsize',15,'FontName','Century');
ylabel('Amplitude','Fontsize',15,'FontName','Century') 

print -djpeg sin_wave_filter.jpg

■20Hzの正弦波をフィルタに通してみる

出力が群遅延分 250 x (1/250) = 0.5[s]後から始まっている。
20Hzの正弦波はそのまま出力されている。

■1Hzの正弦波をフィルタに通してみる

出力が群遅延分 250 x (1/250) = 0.5[s]後から始まっている。
1Hzの正弦波は減衰している。

周波数と位相、群遅延の関係を求めてみる

計測では、入力と出力との間に群遅延分の遅れがあった。
この点をもう少し考えてみる。

ラジアンと度の関係をおさらいする

$2 \pi$ ラジアンが360度なので、

1ラジアンは、
$2\pi : 360 = 1 : x$
$x = \frac{360}{2\pi}$ 度

1度は、
$x = \frac{2\pi}{360}$ ラジアン

■ 経過時間からラジアンを求める

正弦波の1周期は $2\pi$ ラジアン

$y= A \sin(\frac{2\pi}{T} \cdot t +\theta)$

ので、時間とラジアンの関係は、

$x\[\rm{rad}\]= \frac{2\pi}{T} \cdot t[\rm{s}]$

さて、図から
周波数20.01953[Hz]のときは、位相 -2883.516度なので、

正弦波形は、

$y= \sin(2\pi\cdot20.01953 \cdot t + \frac{2\pi}{360} \cdot -2883.516)$
$y= \sin(40.0391\pi t + \frac{2\pi}{360} \cdot -2883.516)$
$y= \sin(40.0391\pi( t + \frac{2}{360} \cdot \frac{-2883.516}{40.0391}))$
$y= \sin(40.0391\pi( t - 0.4001))$

つまり0.4001秒遅れている。

周波数79.84924[Hz]のときは、位相 -13652.86度なので、

$y= \sin(2\pi\cdot79.84924 \cdot t + \frac{2\pi}{360} \cdot -13652.86)$
$y= \sin(159.6985\pi t + \frac{2\pi}{360} \cdot -13652.86)$
$y= \sin(159.6985\pi( t + \frac{2}{360} \cdot \frac{-13652.86}{159.6985}))$
$y= \sin(159.6985\pi( t - 0.4750))$

線形なので同じ0.4001になるはずだが、ずれる。→計算精度の問題?。

■ 演算精度を変えてやってみる。

Matlabでformat long eとして演算

$y= \sin(2\pi\cdot20.01953 \cdot t + \frac{2\pi}{360} \cdot -2883.516)$ のとき

>> (2/360)*(2883.516/(2*20.01953))
ans = 4.000976379898363e-001

$y= \sin(2\pi\cdot79.84924 \cdot t + \frac{2\pi}{360} \cdot -13652.86)$ のとき

>> (2/360)*(13652.86/(2*79.84924))
ans = 4.749526872279701e-001

なので計算精度の問題ではない。

■ 直線はどんな式になるのか調べる。

y = ax + b とすると

$-2883.516 = a \cdot 20.01953 + b$
$-13652.86 = a \cdot 79.84924 + b$

59.8297a = -1.076934400000000e+004
a=-1.799999665717863e+002=-179.99

そして
$b = -2883.516 - a \cdot 20.01953$
$b = -2883.516 - (-179.99) \cdot 20.01953 = 719.7992$

よって y = -179.99x + 719.7992 となる。

あまり(719.7992)があるために、両者の位相差が一致しなくなる。

あまりは考えない

y = -180x のようにあまりがないと
20Hzのとき-3600度の位相遅れ
60Hzのとき-10800度の位相遅れ
80Hzのとき-14400度の位相遅れ

(a)20Hzのとき
$y= \sin(2\pi\cdot20 \cdot t + \frac{2\pi}{360} \cdot -3600)$
$y= \sin(40\pi t + \frac{2\pi}{360} \cdot -3600))$
$y= \sin(40\pi( t + \frac{2}{360} \cdot \frac{-3600}{40}))$
$y= \sin(40\pi( t - 0.5))$

つまり0.5秒遅れている。

(b)60Hzのとき
$y= \sin(2\pi\cdot60 \cdot t + \frac{2\pi}{360} \cdot -10800)$
$y= \sin(120\pi t + \frac{2\pi}{360} \cdot -10800))$
$y= \sin(120\pi( t + \frac{2}{360} \cdot \frac{-10800}{120}))$
$y= \sin(120\pi( t - 0.5))$

つまり0.5秒遅れている。

(c)80Hzのとき
$y= \sin(2\pi\cdot80 \cdot t + \frac{2\pi}{360} \cdot -14400)$
$y= \sin(160\pi t + \frac{2\pi}{360} \cdot -14400))$
$y= \sin(160\pi( t + \frac{2}{360} \cdot \frac{-14400}{160}))$
$y= \sin(160\pi( t - 0.5))$

つまり0.5秒遅れている。

位相を考えるとときは、このようにまず傾きをもとめて求める周波数の位相差を求め
その位相差を計算に使わないと計算が合わなくなる。

群遅延を考えてみると

群遅延の定義は、 $\frac{d\phi}{d\omega}$ なので、
周波数が変化したときに位相がどのくらい変化するかの傾きを表している。
つまり周波数が分かればそのときの位相差が決定する。

図では位相差は時間で表されている。
つまり、20,40,60...[Hz]のとき500[ms]の遅れが生じると出力されている。

群遅延は傾きなので、本来は、
-180[度/Hz]一定の直線になる。

しかしユーザは、周波数が変化したとき出力波形が入力波形よりどのくらい遅れるか知りたいので、

以下のような計算をして位相遅れを計算してグラフに表示していることに注意！

============= 以下の計算をしたのちグラフに表示している ===================
(a) 20Hzのとき -180[度/Hz] * 20[Hz] = -3600度の位相遅れ
(b) 60Hzのとき -180[度/Hz] * 60[Hz] = -10800度の位相遅れ
(c) 80Hzのとき -180[度/Hz] * 80[Hz] = -14400度の位相遅れ