5次のバタワースフィルタを作って実験しました。
T = 0.01; //サンプリング時間Tastinにつかう wa = 50; //実カットオフ周波数:バタワースの設計時 function [F]=Hz(z); //F= (1+z^-1)/(1+2*z^-1); //F= (1+2*(z^-1)+(z^-2))/(429+802*(z^-1)+373*(z^-2)); //F=1/(z^2+sqrt(2)*z+1); //F=1/((z^3)+2*(z^2)+2*z+1); //F=1/((z^4)+2.61*(z^3)+3.41*(z^2)+2.61*z+1); F=1/((z^5)+3.246*(z^4)+5.236*(z^3)+5.236*(z^2)+3.246*z+1); //5次バタワースフィルタ //F=(1-z^-1)/0.01; endfunction function [F]=Hz1(z); //カットオフ周波数変更 F=Hz(z/wa); endfunction function [s]=Tas(z); //s=(2/0.1)(z-1)/(z+1); k=1/2; //z=(1+k*s)/(1-k*s); s=2/T*(z-1)/(z+1); endfunction //w1 = 0:2*%pi/1000:2*%pi; w1 = -1*%pi:2*%pi/1000:%pi; //周波数軸 //u(1,101); for i = 1:1001 //k = Hz(exp(%i*w1(1,i))); //z //k = Hz1(%i*w1(1,i)*5); //s k = Hz1(Tas(exp( %i*w1(1,i)))); //双一次変換後 wd = 2*atan(wa*T/2); y(i)=abs( k ); //振幅 y2(i)=atan(imag(k),real(k)); //位相 end disp(wd); disp(%pi/wd); disp('プリワープ後 カットオフ周波数'); //y(1001)=0; //グラフの軸をきれいにするため figure(1); plot(w1,y'); figure(2); plot(w1,y2'); z = poly(0,"z"); k= Hz1(Tas(z)); r =denom(k); m=roots(r); disp(m); disp(k);
こんなScilabのプログラムで差分方程式を導出し、エクセルで実験しました(IIRはconvolでできないので他に方法がなかった)。
まず問題としては、実カットオフ周波数wcが低かったり、サンプリング時間Tsが小さすぎて Ts*wc が小さすぎると差分方程式の係数がかなり大きく、極がかなり単位円のそばになって不安定になる危険があるということです。
とりあえず、wa = 50rad/s , Ts =0.01s でフィルタを作ると
2 3 4 5
1 + 5z + 10z + 10z + 5z + z
---------------------------------------------------------------
2 3 4 5
-456.336 + 3012z - 8091.36z + 11100.16z - 7820.304z + 2287.84z
となり、フィルタをかけると(元データは前回の記事)
となります。wa*Ts = 0.5=2PI/12.6 であり、カットオフもは700/80 = 8.75 とずれはありますがまあ、理論通りカットされています。
フィルタ後のデータのピークの大きさはあまり参考にならないですね。
前回のFIRフィルタは40次も使ったのにもかかわらず、今回は10次のフィルタでここまで除去できるうえに、サイドローブの影響が出ないのは便利なフィルタです。
微分器にローパスフィルタを入れた不完全微分器の実験をしようと思います。
コメント
コメントを投稿