%　移動平均の周波数特性
%  2005.08.10　小川作成　Octave用

clear;

dT = 100e-6; 	%[sec]　 サンプリング周期
n_ave = 100;	%[回]    移動平均の個数 (n_ave >= 2)

cal_method = 3; %計算方法 0:Z変換を用いて移動平均を記述した後，双一次変換で連続時間系に
                %           変換して計算←サンプリング短く平均回数多いときはハングしてしまう．
                %         3:Z変換を用いた方法を指定した次数で近似した方法低周波領域ではほぼ一致．
                %         4:0の方法と3の方法の比較．

n_aft = 10;  %[次]    近似後の次数

% 1 + sdT
one_plus_sdT = [dT 1];

% 1 - sdT
one_minus_sdT = [-dT 1];

% (1 + sdT)^(n-1)
one_plus_sdT_n_1 = [1];
for ii = 1:(n_ave-1)
one_plus_sdT_n_1 = conv(one_plus_sdT_n_1,one_plus_sdT);
end

% (1 + sdT)^n
one_plus_sdT_n = conv(one_plus_sdT_n_1,one_plus_sdT);

% (1 - sdT)^n
one_minus_sdT_n = [1];
for ii = 1:n_ave
one_minus_sdT_n = conv(one_minus_sdT_n,one_minus_sdT);
end

%--------------------------------------------------------
if (cal_method == 0)||(cal_method == 4)
num_mave0 = (one_plus_sdT_n - one_minus_sdT_n)/n_ave;
den_mave0 = [one_plus_sdT_n - conv(one_minus_sdT,one_plus_sdT_n_1)];
sys_mave0 = tf2sys(num_mave0,den_mave0,0);
end;

if (cal_method == 0)
bode(sys_mave0);
end;

%----------------------------------------------------------
if (cal_method == 3)||(cal_method == 4)
% (1 + sdT)^nの近似
one_plus_sdT_n2 = 1;
for ii = 1:n_aft
one_plus_sdT_n2 = [ prod((n_ave-ii+1):n_ave)*(dT^ii)/prod(1:ii) one_plus_sdT_n2];
end;

% (1 - sdT)^nの近似
one_minus_sdT_n2 = 1;
for ii = 1:n_aft
one_minus_sdT_n2 = [ prod((n_ave-ii+1):n_ave)*((-dT)^ii)/prod(1:ii) one_minus_sdT_n2];
end;

% (1 + sdT)^(n-1)の近似
one_plus_sdT_n_1_2 = 1;
for ii = 1:n_aft
one_plus_sdT_n_1_2 = [ prod((n_ave-1-ii+1):(n_ave-1))*(dT^ii)/prod(1:ii) one_plus_sdT_n_1_2];
end;

num_mave3 = (one_plus_sdT_n2 - one_minus_sdT_n2)/n_ave;
num_mave3 = [0 num_mave3];
den_mave3 = [0 one_plus_sdT_n2] - conv(one_minus_sdT,one_plus_sdT_n_1_2);

sys_mave3 = tf2sys(num_mave3,den_mave3,0);

end;

if (cal_method == 3)
%    bode(sys_mave3);
    [Gain3,Phase3,w3]=bode(sys_mave3);
       
    subplot(2,1,1);
    grid('on')
    ylabel('Gain[dB]')

    %範囲指定
%    axis([ 0.01 10 -20 0],'auto');
    axis([ 1 1000 -40 0]);

    semilogx(w3/(2*pi()),20*log10(Gain3),'-;;');

    subplot(2,1,2);

    xlabel('[Hz]')
    ylabel('Phase[deg]')
    grid('on')

    %範囲指定
%    axis([ 0.01 10 -180 180],'auto');
    axis([ 1 1000 -180 180]);

    semilogx(w3/(2*pi()),Phase3,'-;;');
    
end;

%----------------------------------------------------------
if (cal_method == 4)
[Gain0,Phase0,w0]=bode(sys_mave0);
[Gain3,Phase3,w3]=bode(sys_mave3);

%hold on
semilogx(w0*2*pi(),20*log10(Gain0));
hold on

semilogx(w3*2*pi(),20*log10(Gain3));
hold off
end;



