% ===============================================================
% Fourier Analysis Demonstration
% DRB 31-Mar-2010
% ===============================================================

% ---------------------------------------------------------------
% SYMBOLIC ANALYSIS
% ---------------------------------------------------------------
syms t n real
N = 20;                 % number of coefficients to compute
t0 = 0;                 % start of fundamental period
t1 = pi;                % end of fundamental period
x = exp(-t/2);          % signal over fundamental period

T0 = t1-t0;             % fundamental period
w0 = (2*pi)/T0;         % fundamental frequency (rad/sec)

% compute a0 coefficient
a0 = int((1/T0)*x,t,t0,t1);

% compute symbolic expression for a and b coefficients
an = int((2/T0)*x*cos(n*w0*t),t,t0,t1);
bn = int((2/T0)*x*sin(n*w0*t),t,t0,t1);

% ---------------------------------------------------------------
% NUMERIC ANALYSIS
% ---------------------------------------------------------------
% compute numeric values for a and b coefficients
a = zeros(1,N);
b = zeros(1,N);
for n = 1:N,
    a(n) = int((2/T0)*x*cos(n*w0*t),t,t0,t1);
    b(n) = int((2/T0)*x*sin(n*w0*t),t,t0,t1);
end

% plot reconstructed signal
t = linspace(t0-T0,t1+T0,1000);     % time vector
z = zeros(1,length(t));             % signal
z = z+a0;                           % add in a0 coefficient
for n = 1:N,
    z = z+a(n)*cos(w0*n*t)+b(n)*sin(w0*n*t);
end
subplot(3,1,1);
plot(t,z);
grid on;
xlabel('time');
ylabel('signal');

% plot Fourier magnitude spectrum
C = sqrt(a.^2+b.^2);
subplot(3,1,2);
stem([0:N],[abs(a0) C]);
xlabel('n');
ylabel('C_n');

% plot Fourier phase spectrum
theta = atan2(-b,a);
subplot(3,1,3);
stem([0:N],[0 theta]);
xlabel('n');
ylabel('\theta_n');






