Sawtooth Example Harman Ex8.2

% EX8_2.M Plot the Fourier series of the function f(t) in Example 8.1
%   f(t)=0  -pi < t < 0         Sawtooth pulse train
%   f(t)=t   0  < t < pi
%
% Plot f(t) for 5 and 20 terms in the series
clear
t =[-pi:.031:pi];               % Time points for plotting
sizet=size(t);
fn = pi/4*(ones(sizet));        % Fourier approximation at each t
yplt=zeros(sizet);              %  for plot of f(t)
% 5 terms
for n=1:5
 fn=fn+ (1/pi)*(-2*cos((2*n-1)*t)/(2*n-1)^2)-((-1)^n*sin(n*t)/n); 
end					
%
for k=1:length(t)               % Create f(t) 
  if t(k) < 0 
   yplt(k)=0;
  else
   yplt(k)=t(k);
  end 
end
clf                             % Clear any figures 
subplot(2,1,1),plot(t,fn,t,yplt,'--');
xlabel('t')
ylabel('f(t)'),grid
title('Fourier series approximation to f(t)  Example 8.1 and 8.2')
legend(['N=',num2str(n)],'f(t)') % Annotate the graph
% Add 15 more terms
for n=6:20
 fn=fn+ (1/pi)*(-2*cos((2*n-1)*t)/(2*n-1)^2)-((-1)^n*sin(n*t)/n); 
end
subplot(2,1,2),plot(t,fn,t,yplt,'--');
xlabel('t')
ylabel('f(t)')
legend(['N=',num2str(n)],'f(t)')
%
% Modify the program to compute an arbitrary number of terms 
%  in the series (i.e. input n). Plot the graph for many terms 
%  and notice the overshoot at the ends of the interval no matter
%  how many terms are taken.
%    (This is explained in the text as the Gibbs phenomenon.)

------------------------------





Fourier Series of pulses

%CODE from M-file fourierSeriestest1a.m   A square wave of period T=2*pi
% Show fundamental, 3 terms, 5 terms, 1,3,5 and then 55 odd terms to N=101
clear, close, clf, clc  % Clear Workspace, figures, command window
nmax = 5; %  Iterations for the Fourier series for first four plots
t = linspace(-2*pi,2*pi,1000); % Domain of plot using 1000 points
% Plot of function f(x)- the square wave from -2pi to 2pi (Two periods)
for i = 1:length(t)
if t(i) <= -pi | (t(i) > 0 & t(i) <= pi)
f(i) = 1;
end
if (t(i) > -pi & t(i) <= 0) | t(i) > pi
f(i) = -1;
end
end
figure(1),plot(t,f),title('Press a key for next Odd Harmonic, T=2\pi')
xlabel('Time in Seconds')
hold on    % Plot over current graph
% Calculate and plot first 3 partial sums as a key is pressed
A0 = 0; y = A0; % No dc term for the square wave 
for n=1:2:nmax
An = 0;             % Odd Function (The A terms could be dropped!)
Bn = 4/(n*pi);      % Sine coefficients
y = y + An*cos(n*t) + Bn*sin(n*t);  % The time series
plot(t,y),hold on
fprintf('Press a Key'),pause
end       
xlabel('Time in Seconds'),ylabel('y'),title('Fourier Series of Square Wave')
legend('f','n = 1','n = 3','n = 5')  
hold off
% Do (101-1)/2 +1  55 odd terms Make separate figure
A0=0, y101 = A0;
for n1=1:2:101
An = 0;
Bn = 4/(n1*pi);
y101 = y101 + An*cos(n1*t) + Bn*sin(n1*t);
end 
fprintf('Press a Key'),pause
figure(2)
plot(t,y101),grid,xlabel('Time in Seconds'),ylabel('y101')
title('Fourier Series of Square Wave Train for n=1,3,5,...101')

-------------------------------------

% Complex Fourier series and Spectrum
% Description: Plots the truncated Fourier Series and specturm of an even
% rectangular pulse train,amplitude A,period is 2*pi seconds,w=1 r/s,tau=pi. 
% Here f(t)=A*tau/T +(2*A/pi)(Tau/T)*SUM(sinc(n*w0*tau/2)*cos(n*wo*t)
% A*tau/T = 1*2pi/pi=2.
clear, clf, clc;      % clear all variables & Figures & Command Window                     
A=2,tau=pi, T=2*pi    % Define Amplitude, pulse width, period
N = 21;               % Select arbitrary summation limit (use N odd)
w0 = 2*pi/T;          % w0=1, f=1/2pi, T=2pi fundamental frequency (rad/s)
t = -pi:0.01:3*pi;   % declare time values- Plot two cycles
% Compute yce, the Fourier Series in complex exponential form
% The complex values for n=1,2,... are 1/2 the trig values a_n
c0 = A*tau/T;        % dc bias Here c0=2*(1/2)=1 
yce = c0*ones(size(t));% initialize yce to c0-This adds the dc value to all
%  NOTE:   sinc= sin(pi*x)/(pi*x) function using MATLAB. So use w0/pi
for n = -N:2:N,
cn =(A)*(tau/T)*sinc(n*(w0/pi)*tau/2) ;    % 1/(j*n*pi)
yce = yce + (cn*exp(j*n*w0*t));  % Fourier Series computation
end
figure(1),plot(t,yce);			  % plot truncated exponential FS
xlabel('t (seconds) Period is 2*pi'); ylabel('y(t)');
ttle = ['Truncated Exponential Fourier Series with N = ',...
         num2str(N)];  % Convert N from number to ASCII
title(ttle);           % Put on a title
% Draw the amplitude spectrum from exponential Fourier Series
figure(2)                          % put next plots on figure 2
subplot(2,1,1)
stem(0,c0);                        % plot c0 at nwo = 0
hold;
for n = -N:2:N,                    % loop over series index n
    cn =(A)*(tau/T)*sinc(n*(w0/pi)*tau/2) ;    % 1/(j*n*pi)
  %cn = (1/2)*((4*A)/(j*n*pi));     % Fourier Series Coefficient
  stem(n*w0,abs(cn))               % plot |cn| vs nwo
end
for n = -N+1:2:N-1,                % loop over even series index n 
  cn = 0;                          % Fourier Series Coefficient
  stem(n*w0,abs(cn));              % plot |cn| vs nwo
end
xlabel('w (rad/s)'),ylabel('|cn|')
ttle = ['Amplitude Spectrum with N = ',num2str(N)];
title(ttle),grid, hold;

% Draw the phase spectrum from exponential Fourier Series
subplot(2,1,2)
stem(0,angle(c0)*180/pi);          % plot angle of c0 at nwo = 0
hold;
for n = -N:2:N,                    % loop over odd series index n 
  cn = 2/(j*n*w0);                 % Fourier Series Coefficient
  stem(n*w0,angle(cn)*180/pi);     % plot |cn| vs nwo
end
for n = -N+1:2:N-1,                % loop over even series index n 
  cn = 0;                          % Fourier Series Coefficient
  stem(n*w0,angle(cn)*180/pi);     % plot |cn| vs nwo
end
xlabel('w (rad/s)'),ylabel('angle(cn) (degrees)')
ttle = [' Phase Spectrum with N = ',num2str(N)];
title(ttle);
grid;
[X,Y] = ginput(5)  % Pick off the values
%      n*w0     Amplitude 
%       0       1.0
%X =  1.0365  0.6377    (2/pi)   % w in rad/sec and amplitude
%    3.1094   0.2126    (2/3pi)  % Trig coefficients are 2*Cn
%    5.0672   0.1108
%    7.1401   0.0808
%    8.9827   0.0689