Scilabでのプログラミング

質問

Scilabを使って、FIRフィルタ（移動平均）のグラフを書きたいと思っています。
http://www.heg.co.jp/dspnyuumon/dsp1-3.htm
上記のサイトにある、移動平均のブロック図（図3）のグラフ（つまり図4のグラフ）を作りたいのです。
サイトには
「xに現在サンプルした値が入っており、yにその結果が入る。
x1は現在から1サンプル前の値、x2は現在から2サンプル前の値、同様にx3
y = x;
y = y + x1;
y = y + x2;
y = y + x3; /* 4つの値を足し合わせる */
y = y / 4; /* 4つの平均をとる */
x3 = x2; /* x3に1つ前の値であるx2を代入する */
x2 = x1;
x1 = x;
この処理を1サンプル取り込むごとに実行する。」

という流れが書いてありますが、実際にこれをScilabで作るならどのように作るればよいのでしょうか？
c言語ならfor文で回して作れそうなのですが、Scilabでfor文を入れるとうまくいきませんでした。

上記のサイトの図4のようなグラフをScilabで作りたいという事です。
プログラミングの得意な方、何とかお願いします。

回答No.1

sinisorsa

図４のグラフだけ描きたいのであれば、フィルタの計算は必要ありません。
このフィルタの周波数特性の計算行えばよいのです。

このフィルタの伝達関数は
　　H(z)=a0+a1・z^(-1)+a2・z^(-2)+a3・z^(-3)
で与えられます。
規格化角周波数をωとしたとき、ｚ＝e^(jω)とおけば、
周波数特性は
　　H(e(jω))=a0+a1・e^(-jω)+a2・e^(-2jω)+a3・e^(-3jω)
a0=a1=a2=a3=0.25として、複素数のままScilabで計算して
プロットすればよいのですが、このFIRフィルタは線形位相条件を
満たしますので、
　　H(e(jω))=0.5・e^(-jω3/2){cos(ω/2)+cos(3ω/2)}
と整理されます。

図４と同様に、標本化周波数を10000Hzとし、周波数を0から5000Hzまで
（基本帯域）の振幅特性を求めましょう。

Freq=0:10:5000;
SamplingFreq=10000;
W = 2*%pi*Freq/SamplingFreq;
H=0.5*(cos(W/2)+cos(W*3/2));
plot2d(Freq,20*log10(abs(H)),rect=[0,-80,5000,0]);

複素数バージョンは、
Freq=0:10:5000;
SamplingFreq=10000;
W = 2*%pi*Freq/SamplingFreq;
Zinv = exp(-%i*W);
H=0.25*(1+Zinv+Zinv^2+Zinv^3);
plot2d(Freq,20*log10(abs(H)),rect=[0,-80,5000,0]);

もうひとつの方法は、インパルス応答を求め、これの離散フーリエ変換
を行っても求められる。
インパルス応答は、
h(n)=0 （n<0のとき）
h(0)=h(1)=h(2)=h(3)=0.25
h(n)=0 （n>4のとき）
です。
離散フーリエ変換はScilabではfftという関数です。
以上　参考まで