% Sampling Sinusoids, Periods, Fourier Series:
%
% @ Kefu Xue, Ph.D., June 2006
%
f0=20; % fundamental frequency from Maximum Common Factor: 40pi
T0=1/f0; % shortest period of the signal
Omegas1=450*pi; % selected sampling frequency 1
T1=2*pi/Omegas1; % sampling 40 samples over a period of the signal
N1=45; % plot 1 period
L5=20;L3=12;L2=8; % number of periods covered for each frequency components
t1=[0:N1-1]*T1; % time sampling vector
s11=2*cos(200*pi*t1-1);
s12=sin(120*pi*t1+pi/3);
s13=-0.5*cos(80*pi*t1+2);
x1=s11+s12+s13;
numprd1=N1/(T0/T1);
% plot the signals 
figure;subplot(4,1,1);stem(t1,s11,'r');axis tight; 
title(['The first sinusoid (Omega5) shown with ' num2str(L5) ' periods.']);
subplot(4,1,2);stem(t1,s12,'b');axis tight; 
title(['The second sinusoid (Omega3) shown with ' num2str(L3) ' periods.']);
subplot(4,1,3);stem(t1,s12,'g');axis tight; 
title(['The second sinusoid (Omega2) shown with ' num2str(L2) ' periods.']);
subplot(4,1,4);plot(t1,x1,'b');axis tight; 
title(['One period covers ' num2str(numprd1) ' periods of original signal.']);
xlabel('time in second');
%
% Different sampling frequency
%
Omegas2=500*pi; % selected sampling frequency 2
T2=2*pi/Omegas2; % sampling 40 samples over a period of the signal
N2=25; % plot 1 period
L25=10;L23=6;L22=4; % number of periods covered for each frequency components
t2=[0:N2-1]*T2; % time sampling vector
s21=2*cos(200*pi*t2-1);
s22=sin(120*pi*t2+pi/3);
s23=-0.5*cos(80*pi*t2+2);
x2=s21+s22+s23;
numprd2=N2/(T0/T2);
% plot the signals 
figure;subplot(4,1,1);stem(t2,s21,'r');axis tight; 
title(['The first sinusoid (Omega5) shown with ' num2str(L25) ' periods.']);
subplot(4,1,2);stem(t2,s22,'b');axis tight; 
title(['The second sinusoid (Omega3) shown with ' num2str(L23) ' periods.']);
subplot(4,1,3);stem(t2,s22,'g');axis tight; 
title(['The second sinusoid (Omega2) shown with ' num2str(L22) ' periods.']);
subplot(4,1,4);plot(t2,x2,'b');axis tight; 
title(['One period covers ' num2str(numprd2) ' periods of original signal.']);
xlabel('time in second');
%
% DFS
%
% Sampled signal 1
%
X1=fft(x1); % perform Discrete-Fourier Transform
C1=X1/N1; % Fourier Series Coefficients
[lin,lindex]=find(abs(C1)>10^-3);
f1=[0:N1-1]*Omegas1/N1/pi; % form the frequency axis
linephase=zeros(size(C1));
linephase(lindex)=angle(C1(lindex));
figure;subplot(2,1,1);stem(f1,abs(C1));axis tight;
title('Magnitude Line Spectrum');ylabel('Magnitude');
subplot(2,1,2);stem(f1,linephase);axis tight;
title('Phase Line Spectrum');ylabel('Phase (radian)');
xlabel('frequency (pi radian/seconds)');
%
% Sampled singal 2
%
X2=fft(x2);
C2=X2/N2; % Fourier Series Coefficients
[lin,lindex2]=find(abs(C2)>10^-3);
f2=[0:N2-1]*Omegas2/N2/pi; % form the frequency axis
linephase2=zeros(size(C2));
linephase2(lindex2)=angle(C2(lindex2));
figure;subplot(2,1,1);stem(f2,abs(C2));axis tight;
title('Magnitude Line Spectrum');ylabel('Magnitude');
subplot(2,1,2);stem(f2,linephase2);axis tight;
title('Phase Line Spectrum');ylabel('Phase (radian)');
xlabel('frequency (pi radian/seconds)');
%
% display the Coefficients
%
disp('The Fourier Series Coefficients (in the order of): C2, C3, C5, C5, C3, C2');
disp(['Magnitude: ' num2str(abs(C1(lindex)))]);
disp(['Phase (radian): ' num2str(linephase(lindex))]);
disp(['Line Spectrum in Matlab index : ' num2str(lindex)]);
disp(['Line Spectrum in Math index : ' num2str(lindex-1)]);
disp(['Line Spectrum in frequency (pi radian): ' num2str((lindex(1:3)-1)*Omegas1/N1/pi) ... 
    '   ' num2str((lindex(4:end)-1-N1)*Omegas1/N1/pi)]);