t=[0:.002:2];
N=200;
n = [0:N];
j = sqrt(-1);
%
% Square Wave
cn = 2*ones(1,length(n))./(j*pi*n);
cn(1:2:N+1) = 0; % Set even terms to zero
%
% Periodic pulses
%
delta = .1;
cn = exp(-j*pi*n*delta).*sin(pi*n*delta)./(pi*n);
cn(1) = delta;
%
% Half-wave rectified sinusoid
%
cn(1) = 1/pi;
cn(2) = -j*1/4;
cn(3:N+1) = -(1/pi)*ones(1,N-1)./([2:N].^2-1);
cn(4:2:N+1) = 0;
%
% Full-wave rectified sinusoid
%
cn(1) = 2/pi;
cn(2) = 0;
cn(3:N+1) = -(2/pi)*ones(1,N-1)./([2:N].^2-1);
cn(4:2:N+1) = 0;
%
% Form and plot Fourier series representations
%
x = cumsum([cn(1)*ones(1,length(t));zeros(N-1,length(t))]...
   +(cn(2:N+1).'*ones(1,length(t)).*exp(j*2*pi*[1:N]'*t))... % Positive integers
   +(cn(2:N+1)'*ones(1,length(t)).*exp(-j*2*pi*[1:N]'*t))); % Negative integers
figure(1);plot(t,x(N,:));set(gca,'fontsize',18);
% figure(1);plot(t,x);set(gca,'fontsize',18);
%
% Filter by modifying the Fourier series
%
f = n; % Assume period = 1;
%
% First order lowpasss filter
%
RC = 2;
H=ones(1,length(f))./(j*2*pi*f*RC+1);
figure(2)
subplot(211);
loglog(f,abs(H));set(gca,'fontsize',18);grid
subplot(212);
semilogx(f,angle(H)*180/pi);set(gca,'fontsize',18);grid
axis([0 N -90 0])
%
% First order highpass 
%
RC = .05;
H=j*2*pi*f*RC./(j*2*pi*f*RC+1);
figure(2)
subplot(211);
loglog(f,abs(H));set(gca,'fontsize',18);grid
subplot(212);
semilogx(f,angle(H)*180/pi)grid
axis([0 N 0 90])
%
% Filter the signal x(t)
%
dn = cn.*H; % Multiply Fourier coefficients by transfer function
y = cumsum([dn(1)*ones(1,length(t));zeros(N-1,length(t))]...
   +(dn(2:N+1).'*ones(1,length(t)).*exp(j*2*pi*[1:N]'*t))... % Positive integers
   +(dn(2:N+1)'*ones(1,length(t)).*exp(-j*2*pi*[1:N]'*t))); % Negative integers
figure(3);
subplot(211);plot(t,x(N,:));set(gca,'fontsize',18);
subplot(212);plot(t,y(N,:));set(gca,'fontsize',18);a=axis; a(3)=0;axis(a)
% figure(3);
% subplot(211);plot(t,x);subplot(212);plot(t,y);axis([0 2 -2 2])
% set(gca,'fontsize',18);
