clear all;
close all;
set(0,'defaultAxesFontSize',20)
set(0,'defaultTextFontSize',20)
set(0,'defaultAxesFontName','Helvetica')
set(0,'defaultTextFontName','Helvetica')


# 問2 振幅1V, 時間幅 1/4 sec, 周期 1 secの方形パルス列の、三角関数型フーリエ級数
# の片側スペクトル成分を図示せよ。
A=1 #volt
tau = 1/4 #sec
T = 1 #sec
t= -2*T: T/50: 2*T;

for i=1:17
  N=ceil(1.5^i); #N次高調波まで計算
  hold on;

  C=zeros(1,N+1);
  f=zeros(1,N+1);
  for n = 0:N
    f(n+1)=n/T;#周波数
    if (n == 0)
      C(n+1) = A*tau/T; #=C0
      v=C(n+1);
    else
      x= n*pi*tau/T;
      C(n+1) = 2*C(0+1)*sin(x)/x; #=Cn
      v += C(n+1)*cos(2*pi*n*t/T);
    endif
  endfor

  subplot(2,1,1);
    plot (t, v,'linewidth', 2);
    xlabel('t [sec]')#,'fontsize',20)
    ylabel('Amplitude [V]')#,'fontsize',20)
    s = "N=";
    title(strcat( s, num2str(N)));

  subplot(2,1,2);
    if N<21
      stem(f(1:N+1),C(1:N+1));
      #close all;
    else  
      stem(f(1:21),C(1:21),'linewidth', 2);
    endif
    xlabel('f [Hz]')#,'fontsize',20)
    ylabel('Amplitude [V]')#,'fontsize',20)
    
  refresh;
endfor

print eq_1_22.pdf