% Demonstration of frequency response of system for pulse train input

% Ideal input signal as a function of time
tv = -5:0.05:5;
xv_ideal = zeros(size(tv));
xv_ideal(find(abs(tv)<=0.5)) = 1;
xv_ideal(find((abs(tv)>=1.5) & (abs(tv)<=2.5))) = 1;
xv_ideal(find((abs(tv)>=3.5) & (abs(tv)<=4.5))) = 1;

% Fourier series approximation to input as a function of time
N = 40;
w0 = pi;
xv = zeros(size(tv));
Xv = zeros(2*N+1,1);
for k=-N:N
  ckx = 0;
  if mod(k,4)==1, ckx = 1/(k*pi); end
  if mod(k,4)==3, ckx = -1/(k*pi); end
  if k==0, ckx = 0.5; end
  xv = xv + ckx*exp(j*w0*k*tv);
  Xv(k+N+1) = ckx;
end
xv = real(xv);  % should be real anyway

% Output Fourier series as a function of time
% Fourier series approximation as a function of time
w0 = pi;  R = 1;  C = 1;
xva = zeros(size(tv));
yva = zeros(size(tv));
Yv = zeros(2*N+1,1);
for k=0:N
  
  % Positive frequency component
  ckx = 0;
  if mod(k,4)==1, ckx = 1/(k*pi); end
  if mod(k,4)==3, ckx = -1/(k*pi); end
  if k==0; ckx = 0.5; end
  xva = xva + ckx*exp(j*w0*k*tv);
  Hw = 1/(1 + j*w0*k*R*C);
  ckv = ckx*Hw;
  yva = yva + ckv*exp(j*w0*k*tv);
  Yv(k+N+1) = ckv;
  
  % Negative frequency component
  if k~=0  
    cmkx = 0;
    if mod(-k,4)==1, cmkx = 1/(-k*pi); end
    if mod(-k,4)==3, cmkx = -1/(-k*pi); end
    xva = xva + cmkx*exp(j*w0*-k*tv);
    Hw = 1/(1 + j*w0*-k*R*C);
    cmkv = cmkx*Hw;
    yva = yva + cmkv*exp(-j*w0*k*tv);
    Yv(-k+N+1) = cmkv;
  end
  xva = real(xva);  yva = real(yva);  % should be real anyway
  
end

% Plot signals
fh = figure;
subplot(2,1,1);
plot(tv,xv_ideal,'r',tv,xva,'g',tv,yva,'b');  axis([-5 5 -0.2 1.2]);
legend('Actual input', 'Approximate input', 'Approximate output');

subplot(2,1,2);
shx = stem(-N:N,abs(Xv),'g.');
hold on;
shx = stem(-N+0.1:N+0.1,abs(Yv),'b.');
hold off;

return;
