DTFTのScilabプログラム

始めに

こんにちは、ニートです。
ScilabにはDFT（FFT）をする機能は標準で付いていますが、DTFT（z変換）する機能は付いていないと思われます。そこで、Scilabを用いたDTFTを行うプログラムを紹介します。
ディジタルフィルタの設計には必要となります。

DTFTとは

DTFTとは、z変換した離散伝達関数の周波数特性を見るものです。
例えば、直近5サンプルを平均化するシステムの離散伝達関数は

(1+z^-1+z^-2+z^-3+z^-4)/5

となります。このシステムに離散信号入力を通すと、離散信号の周波数特性に離散システムをDTFTした結果を掛けたものが、出力離散信号の周波数特性となります。
このシステムをDTFTすると、

これは周波数特性的にはLPFに近い特性です。正規価格周波数0.5rad程度でカットオフするシステムと考えられます。

差分方程式の周波数特性確認プログラム

T=0.01;                                 //サンプリング時間
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=z/((z^5)+3.246*(z^4)+5.236*(z^3)+5.236*(z^2)+3.246*z+1);
    //F=(1-(z^-1));
    //F=z
    //F=(1+4*z^-1+5*z^-2-5*z^-4-4*z^-5-z^-6)/(22.878-78.20*z^-1+111.00*z^-2-80.91*z^-3+30.12*z^-4-4.563*z^-5)          //近似微分器
    //F=(1-2*z^-1+z^-2)/T^2;
    F=(1+z^-1+z^-2+z^-3+z^-4)/5
endfunction
    

w1 = -1*%pi:2*%pi/1000:%pi;              //周波数軸

for i = 1:1001
   k = Hz(exp(%i*w1(1,i)));               //z = exp(j*w)を代入し応答を見る
   //k = Hz(exp(%i*w1(1,i)*T));           //z = exp(j*w*T)の場合
   //k = Hz(%i*w1(1,i));                  //s = jw の場合
   y(i)=abs( k );                         //振幅
   y2(i)=atan(imag(k),real(k));           //位相
end

figure(1);
plot(w1,y');

figure(2);
plot(w1,y2');

差分方程式F(z)に z=exp(j*w)を代入し、各wに対し値を取ればDTFTできたことになります。

離散信号のDTFT

インパルス応答をフィルタとして入力信号とたたみこんで使用する場合、インパルス応答のDTFTの結果がわかれば、フィルタの性能がわかります。インパルス応答をFFTしてもインパルス応答が繰り返おこると仮定したものとなるので正しいフィルタの性能評価とはならないので注意が必要です。
Scilabのプログラムは次の方法でうまくいきます。

function[F] =DTFT(w);
    F=0;
    for m=1:n
         F = F + h(m)*exp(-1*%i*(m-1)*w);   //DTFT、重ね合わせ
    end
endfunction
    
//w1 = 0:2*%pi/1000:2*%pi;
w1 = -1*%pi:2*%pi/1000:%pi;              //周波数軸

//u(1,101);
for i = 1:1001
   k = DTFT(w1(1,i));
   y(i)=abs( k );                         //振幅
   y2(i)=atan(imag(k),real(k));           //位相
end

//y(1001)=0;                              //グラフの軸をきれいにするため
figure(1);
plot(w1,y');

figure(2);
plot(w1,y2');
xgrid(1);

これを実行する前に、配列h[n]に対象のデータを保存します。nはデータの個数です。
下のh[]はFIR型LPF となるようにsinwt/t　で作りました。

clear;
n=41              //データ個数  : -1 がフィルタ次数
h=[-1.94988E-17
-0.016753152
1.94988E-17
0.018724111
-1.94988E-17
-0.021220659
1.94988E-17
0.024485376
-1.94988E-17
-0.028937262
1.94988E-17
0.035367765
-1.94988E-17
-0.045472841
1.94988E-17
0.063661977
-1.94988E-17
-0.106103295
1.94988E-17
0.318309886
0.5
0.318309886
1.94988E-17
-0.106103295
-1.94988E-17
0.063661977
1.94988E-17
-0.045472841
-1.94988E-17
0.035367765
1.94988E-17
-0.028937262
-1.94988E-17
0.024485376
1.94988E-17
-0.021220659
-1.94988E-17
0.018724111
1.94988E-17
-0.016753152
-1.94988E-17];

結果は

位相特性は悪い気がしますが、LPFとなっています。
ちなみに今回使ったh[]をDFTすると:plot(fft(h));とすると

DTFTとはまた違った結果となります。どういった関係性があるかはよくわからないですが・・・違った結果になるのはわかっています