% Fourier series and Fourier transforms
%
% @Kefu Xue, Ph.D., May 2001
%
% consider a signal
% x(t)=3cos(20*t)+sin(8t-1)-2cos(32t+2)
% find the period of the signal 
%
wf=4; %rad/sec
Tp=2*pi/wf;
Np=32; % number of samples per period.
deltaT=Tp/Np;
t=0:deltaT:(Np-1)*deltaT;
x=3*cos(20*t)+sin(8*t-1)-2*cos(32*t+2);
% find the Fourier series coefficients
Dn=fft(x)/Np; % divide Np for the normalization. 
zindx=find(abs(Dn)<10^-7); 
% find index of the Dn's with very small amplitude.
Dn(zindx)=0; 
% set them to zero so they won't creat large false phase value.
n=0:Np-1;
subplot(2,1,1);stem(n,abs(Dn));title('magnitude');
subplot(2,1,2);stem(n,angle(Dn));title('phase');
%
pause
%
Dnshift=fftshift(Dn); % fftshift() creats correct index value.
n=-Np/2:Np/2-1;
% the coefficients Dn are now correctely labeled 
figure;subplot(2,1,1);stem(n,abs(Dnshift));title('magnitude after shift');
subplot(2,1,2);stem(n,angle(Dnshift));title('phase after shift');
%
pause
%
% consider a sawtooth signal
w0=2*pi/10; % fundamental frequency
Np=16;
deltaT=10/Np;
t=0:deltaT:(3*Np-1)*deltaT;
D0=5; % D0
%
% using only D1 and D10
n=1:10;
Dn=j*5/pi./n; %Dn
% signal synthesis
xsyn1=D0+2*abs(Dn(1))*cos(w0*t+angle(Dn(1)));
xsyn3=xsyn1+2*abs(Dn(2))*cos(2*w0*t+angle(Dn(2)))+2*abs(Dn(3))*cos(3*w0*t+angle(Dn(3)));
xsyn5=xsyn3+2*abs(Dn(4))*cos(4*w0*t+angle(Dn(4)))+2*abs(Dn(5))*cos(5*w0*t+angle(Dn(5)));
xsyn7=xsyn5+2*abs(Dn(6))*cos(6*w0*t+angle(Dn(6)))+2*abs(Dn(7))*cos(7*w0*t+angle(Dn(7)));
xsyn10=xsyn7+2*abs(Dn(8))*cos(8*w0*t+angle(Dn(8))) ...
        +2*abs(Dn(9))*cos(9*w0*t+angle(Dn(9)))+2*abs(Dn(10))*cos(10*w0*t+angle(Dn(10)));
figure;plot(t,xsyn1,'g',t,xsyn3,'r',t,xsyn5,'k',t,xsyn10,'b');

%
pause

close all