%	fig_filter1.m
%	illustrate filtering of a square wave by an RC circuit

f0 = 0.5;				% 0.5kHz square wave
k = -9:9;				% plot these Fourier coefficients
ck = 0.5 * (k == 0);			% DC term
ko = min(k):2:max(k);			% odd k's
ck(1:2:length(k)) = 1 ./ (j*ko*pi);	% ck's for odd k

t = linspace(-0.5, 4.3, 1001);

x = mod(t, 2) < 1;			% square wave trick
clf, subplot(321), plot(t, x)
xlabel('t [msec]'), ylabel('x(t)'), title('Square wave')
axis([range(t)' -0.1 1.1])

subplot(322), stem(k*f0, abs(ck), 'filled')
xlabel('Frequency [kHz]'), ylabel('|c_k|')
title('Square wave magnitude spectrum')
axis([range(k*f0)' 0 0.52])

rc = 1;
rc = 2;
rc = 0.1;
subplot(323), plot(t, exp(-t/rc)/rc .* (t >= 0))
xlabel('t [msec]'), ylabel('h(t)'), title('RC impulse response')
axis([range(t)' -0.1 1.1/rc])
text(2, 0.7/rc, sprintf('RC=%g msec', rc))

f = linspace(min(k)*f0, max(k)*f0, 101);
H = 1 ./ (1+j*2*pi*f*rc);		% frequency response
subplot(324), plot(f, abs(H))
xlabel('Frequency [kHz]'), ylabel('|H(j2\pi f)|')
title('RC filter magnitude response')
axis([range(k*f0)' 0 1])

f = k * f0;				% frequencies of harmonics
H = 1 ./ (1+j*2*pi*f*rc);		% H(j*k*omega_0)
subplot(326), stem(k*f0, abs(ck .* H), 'filled')
xlabel('Frequency [kHz]'), ylabel('|c_k| H(j \omega_0 k)')
title('Output signal magnitude spectrum')
axis([range(k*f0)' 0 0.52])

ck = ck .* H;		% coefficients after filtering
y = zeros(size(t));
for ik=1:length(k)	% synthesize output signal using Fourier series
	y = y + ck(ik) * exp(j*2*pi*f0*k(ik)*t);
end
subplot(325), plot(t, Real(y))
xlabel('t [msec]'), ylabel('y(t)'), title('Filtered square wave')
axis([range(t)' -0.1 1.1])

%print('fig_filter1', '-deps')