% MATLAB demonstration of a square wave filtered by a 1st-order lowpass
% filter, H(w) = 1/(1 + j*w).
% The input signal to the filter is f(t) shown on page 196 of the Lathi
% textbook.
%
% ELEC 320, Rich Kozick,

% First, get Fourier series for input signal:
N = input('Number of harmonics: ');

clear C H
t = -3*pi:6/1000:3*pi;  % Time values to evaluate function: -3 to 3, steps of 6/1000
T = 2*pi;		  % Period
w0 = 2*pi/T;	  % Fundamental frequency
wc = 1;           % Cutoff frequency of filter

C0 = 0.5;
fN = C0*ones(1,length(t));   % dc component

H0 = 1/(1 + j*0);   % Gain of system at 0 freq.
yN = C0*H0*ones(1,length(t));

% Compute truncated trig. Fourier series

for n=1:N
  if (rem(n,2) == 0)	% Even harmonic
    C(n) = 0;
    theta(n) = 0;
  else			% Odd harmonic
    C(n) = 2 / (pi*n);
    if (sin(n*pi/2) < 0)
        theta(n) = pi;
    else
        theta(n) = 0;
    end
  end
  fN = fN + C(n)*cos(n*w0*t + theta(n));
  
  H(n) = 1/(1 + j*(n*w0/wc));  % Gain of system at nth harmonic
  yN = yN +C(n)*abs(H(n))*cos(n*w0*t + +theta(n) + angle(H(n)));
end

figure(1)
subplot(211)
% stem((0:N)*w0, [C0, C], 'b', (0:N)*w0, [C0, C] .* abs([H0, H]), 'r')
stem((0:N)*w0, [C0, C], 'b')
hold on
stem((0:N)*w0, [C0, C] .* abs([H0, H]), 'r*')
hold off
xlabel('FREQUENCY (rad/sec)')
% legend('INPUT', 'OUTPUT')
title('MAGNITUDE SPECTRA')
axis([0 N 0 0.8])
subplot(212)
plot(t,fN, 'b', t, yN, 'r')
legend('INPUT', 'OUTPUT')
title(['FILTERED SQUARE WAVE W/ HARMONICS <= ',num2str(N)])
xlabel('TIME (sec)')

