% Program to show the signal building components in terms of sinusoidals 
%
% By using a fixed sampling frequency, we need to make sure that all the 
% frequency components should be less than the folding frequency fs/2.
echo on;
fs=1000; % fs>fN=2fB 
T=1/fs; % Sampling interval 
%
% frequency components
f1=100; 
f2=320;
f3=180;
T0=1/20; % f0=maximum common factor(100,320,180)=20 Hz is the fundamental frequency.
         % using factor() function to find all the factors first. 
t=0:T:T0-T; % colect the samples covers one period
              % of the continuous time signal.
x1=cos(2*pi*f1*t);
x2=cos(2*pi*f2*t);
x3=cos(2*pi*f3*t);
% plot 3 signal components
figure; % creat a figure and make 3 plots (3 rows and 1 column)in one figure 
subplot(3,1,1); plot(t,x1,'r'); title('signal component 1 (100Hz)');
subplot(3,1,2); plot(t,x2,'b'); title('signal component 2 (320Hz)');
subplot(3,1,3); plot(t,x3,'g'); title('signal component 3 (180Hz)');
% 
pause % press any key to resume
% Let us build 3 signals
%
build1=3*x1+0.5*x2+x3; % recipe: three cups of x1, a half cup of x2 puls one cup of x3
%
figure;
subplot(4,1,1); plot(t,3*x1,'r'); title('amplified signal component 1 (100Hz)'); 
subplot(4,1,2); plot(t,0.5*x2,'b'); title('attenuated signal component 2 (320Hz)');
subplot(4,1,3); plot(t,x3,'g'); title('signal component 3 (180Hz)');
subplot(4,1,4); plot(t,build1,'k'); title('new signal build1');
% 
pause % press any key to resume
% change recipe
x1shift5=[x1(6:end),x1(1:5)]; % time shift 5*T, phase shift 2*pi*f1*5*T
x2shift8=[x2(9:end),x2(1:8)];
build2=3*x1shift5+0.5*x2shift8+x3;
figure;
subplot(4,1,1); plot(t,3*x1shift5,'r'); title('amplified and shifted signal component 1 (100Hz)');
subplot(4,1,2); plot(t,0.5*x2shift8,'b'); title('attenuated and shifted signal component 2 (320Hz)');
subplot(4,1,3); plot(t,x3,'g'); title('signal component 3 (180Hz)');
subplot(4,1,4); plot(t,build2,'k'); title('new signal build2');
%
pause % press any key to resume
% change the recipe again
build3=x1shift5+5*x2shift8+x3/3;
figure;
subplot(4,1,1); plot(t,x1shift5,'r'); title('shifted signal component 1 (100Hz)');
subplot(4,1,2); plot(t,5*x2shift8,'b'); title('amplified and shifted signal component 2 (320Hz)');
subplot(4,1,3); plot(t,x3/3,'g'); title('attenuated signal component 3 (180Hz)');
subplot(4,1,4); plot(t,build3,'k'); title('new signal build3');
% 
pause % press any key to resume
%
% Fourier series analysis
N=length(t); % length of data
fr=fs/N; % frequency resolution of 
f=0:fr:fs-fr;
fsx1=fft(x1)/N;
ax1=abs(fsx1);
zindx=find(ax1>10^-10);
px1=zeros(1,N);% only calculate phase angle for the samples with nonzero amplitude
px1(zindx)=angle(fsx1(zindx)); 
fsx2=fft(x2)/N;
ax2=abs(fsx2);
zindx=find(ax2>10^-10);
px2=zeros(1,N);% only calculate phase angle for the samples with nonzero amplitude
px2(zindx)=angle(fsx2(zindx)); 
fsx3=fft(x3)/N;
ax3=abs(fsx3);
zindx=find(ax3>10^-10);
px3=zeros(1,N);% only calculate phase angle for the samples with nonzero amplitude
px3(zindx)=angle(fsx3(zindx)); 
figure;
subplot(5,2,1);stem(f,ax1); title('Magnitude of FS of component 1 (100Hz)');axis tight;
subplot(5,2,2);stem(f,px1); title('Phase of FS of component 1 (100Hz)');axis tight;
subplot(5,2,3);stem(f,ax2); title('Magnitude of FS of component 2 (320Hz)');axis tight;
subplot(5,2,4);stem(f,px2); title('Phase of FS of component 2 (320Hz)');axis tight;
subplot(5,2,5);stem(f,ax3); title('Magnitude of FS of component 3 (180Hz)');axis tight;
subplot(5,2,6);stem(f,px3); title('Phase of FS of component 3 (180Hz)');axis tight;
%
pause % press any key to resume
% 
% for the amplified and shifted signal
fsx1shift5=fft(3*x1shift5)/N;
ax1shift5=abs(fsx1shift5);
zindx=find(ax1shift5>10^-10);
px1shift5=zeros(1,N);% only calculate phase angle for the samples with nonzero amplitude
px1shift5(zindx)=angle(fsx1shift5(zindx)); 
subplot(5,2,7);stem(f,ax1shift5); title('Magnitude of FS of amplified and shifted component 1 (100Hz)');axis tight;
subplot(5,2,8);stem(f,px1shift5); title('Phase of FS of amplified and shifted component 1 (100Hz)');axis tight;
%
fsx2shift8=fft(5*x2shift8)/N;
ax2shift8=abs(fsx2shift8);
zindx=find(ax2shift8>10^-10);
px2shift8=zeros(1,N);% only calculate phase angle for the samples with nonzero amplitude
px2shift8(zindx)=angle(fsx2shift8(zindx)); 
subplot(5,2,9);stem(f,ax2shift8); title('Magnitude of FS of amplified and shifted component 2 (320Hz)');axis tight;
subplot(5,2,10);stem(f,px2shift8); title('Phase of FS of amplified and shifted component 2 (320Hz)');axis tight;
%
pause % press any key to resume
%
% Fourier Series analysis of the combined signals
fsbuild1=fft(build1)/N;
abuild1=abs(fsbuild1);
zindx=find(abuild1>10^-10);
pbuild1=zeros(1,N);% only calculate phase angle for the samples with nonzero amplitude
pbuild1(zindx)=angle(fsbuild1(zindx)); 
figure;
subplot(3,2,1);stem(f,abuild1); title('Magnitude of FS of signal build1');axis tight;
subplot(3,2,2);stem(f,pbuild1); title('Phase of FS of signal build1');axis tight;
%
fsbuild2=fft(build2)/N;
abuild2=abs(fsbuild2);
zindx=find(abuild2>10^-10);
pbuild2=zeros(1,N);% only calculate phase angle for the samples with nonzero amplitude
pbuild2(zindx)=angle(fsbuild2(zindx)); 
subplot(3,2,3);stem(f,abuild2); title('Magnitude of FS of signal build2');axis tight;
subplot(3,2,4);stem(f,pbuild2); title('Phase of FS of signal build2');axis tight;
%
fsbuild3=fft(build3)/N;
abuild3=abs(fsbuild3);
zindx=find(abuild3>10^-10);
pbuild3=zeros(1,N);% only calculate phase angle for the samples with nonzero amplitude
pbuild3(zindx)=angle(fsbuild3(zindx)); 
subplot(3,2,5);stem(f,abuild3); title('Magnitude of FS of signal build3');axis tight;
subplot(3,2,6);stem(f,pbuild3); title('Phase of FS of signal build3');axis tight;

%end of file