Scilabを用いたDTFT

FFT(DFT)の機能を持つソフトはいくつか有りますが、フィルタ学ぶためにはインパルス、差分方程式をDTFT(z変換)しないと正しい結果には結びつかないでしょう。矩形窓をDTFTした場合とDFTした場合では結果が異なるのと同じです(DFTは繰り返しを仮定しているためそうなる)。

Scilabを用いてコードを作りました。
まず、DTFTは
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');

こうするとうまくいきました。警告が出ますが、結果は合っていました。
clear;
n=5              //データ個数
h=[0.25
0.25
0
-0.25
-0.25
]
この形で別ファイルにインパルス応答を入れて使います。

差分方程式H(z)からの応答は
function [F]=Hz(z);
    F= (1+z^-1)/(1+2*z^-1);
   
endfunction
    
//w1 = 0:2*%pi/1000:2*%pi;
w1 = -1*%pi:2*%pi/1000:%pi;              //周波数軸
T = 0.01;                                //サンプリングタイム

//u(1,101);
for i = 1:1001
   k = Hz(exp(%i*w1(1,i)*T));
   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');
これでうまくいきました。



コメント