% The script file pwopro03 gives examples of Fourier series
% expansions of periodic functions. The first example is a discontinuous
% skew symmetric function, the second a continuous function with a discontinuous
% derivative.

l=1.00; N=100; dl=l/N; c=1.00; wcl=c/l;
A=-2*l+dl:dl:2*l; A=A';
BA=wcl*dl:wcl*dl:c-wcl*dl;BA=BA';BB=0; BC=-c+wcl*dl:wcl*dl:c-wcl*dl;
BC=BC'; BD=0; BE=-c+wcl*dl:wcl*dl:0; BE=BE';
B=[BA;BB;BC;BD;BE];
plot(A,B,'k.',A,zeros(size(A)),'k')
title('The periodic function with mean values at discontinuities')
pause
close

BAC=wcl*dl:wcl*dl:c; BAC=BAC'; BCC=-c+wcl*dl:wcl*dl:c; BCC=BCC';
BEC=-c+wcl*dl:wcl*dl:0; BEC=BEC';
BBC=[BAC;BCC;BEC];
plot(A,BBC,'k.',A,zeros(size(A)),'k')
title('The periodic discontinous function')
pause
close

B1=2*c/pi*sin(pi/l*A); % Sums of FS terms calculated from formulae
B2=B1-c/pi*sin(2*pi/l*A);
B3=B2+2*c/pi/3*sin(3*pi/l*A);
B4=B3-2*c/pi/4*sin(4*pi/l*A);
B5=B4+2*c/pi/5*sin(5*pi/l*A);
plot(A,B,'k.',A,zeros(size(A)),'k',A,B1,A,B2,A,B3,A,B4,A,B5)
title('The function and comparison by consecutive sums of FS terms') 
pause
close

AP=A(2*N:4*N); BP=B(2*N:4*N);
plot(AP,BP,'k.')
title('The period represented by Fourier series')  
pause
close
FBP=fft(BP);
AK1=real(FBP(2:30)/2/N);
BK1=imag(FBP(2:30)/2/N);
AD=1:30-1;
plot(AD,AK1,'r*',AD,BK1,'k.')
title('The first Fourier coefficients from MATlAB procedure')
pause 
close
BPI=ifft(FBP);
plot(AP,BPI,'k.')
title('The inverse transformation')
pause 
close


BF=wcl*dl:wcl*dl:c;BF=BF'; BG=c-wcl*dl:-wcl*dl:0;BG=BG';
BH=[BF;BG;BF;BG];
C0=c/2*ones(size(A));
C1=C0-4*c/pi/pi*cos(pi/l*A);
C2=C1-4*c/pi/pi/9*cos(pi/l*A);
C3=C2-4*c/pi/pi/25*cos(pi/l*A);
plot(A,BH,'k.',A,zeros(size(A)),'k',A,C0,A,C1,A,C2,A,C3)
title('The continous symmetric periodic function and FS approximations')
pause
close
COFA=[c,-4/pi/pi,-4/pi/pi/9,-4/pi/pi/25,-4/pi/pi/49]';

BHS=BH(2*N:4*N); AHS=A(2*N:4*N);
plot(AHS,BHS,'k.')
title('The period represented by Fourier series')
pause 
close
FBHS=fft(BHS);
AK2=real(FBHS(1:20)/N);
BK2=imag(FBHS(1:20)/N);
AD2=AP(1:20);
AD2=0:20-1;
plot(AD2,AK2,'r*',AD2,BK2,'k.')
title('The first coefficients in FS by Matlab')
pause 
close
COFN=[AK2(1),AK2(2),AK2(4),AK2(6),AK2(8)]'
[COFA,COFN]
BI=0:0.2:1; BJ=0.8:-0.2:0; BK=[BI,BJ]; BK=BK'; 
AK=0:0.2:2; AK=AK';
plot(AK,BK,'k*')
title('The 10 discrete points,last omitted,approximation by 10 terms') 
pause
close
AKR=AK(1:10);
WA=zeros(10,10);
for r=1:6;
   WA(:,r)=cos(2*pi/2*(r-1)*AKR);
end
for r=7:10;
   WA(:,r)=sin(2*pi/2*(r-6)*AKR);
end
BKR=BK(1:10);
FCD=inv(WA)*BKR; 
plot(AKR,FCD,'r o',AKR, zeros(10,1),'k')
title('The coefficients found from the approximation') 
pause
close

COF=[AK2(1)/2;AK2(2);0;AK2(4);0;AK2(6)];
COP=[COF,FCD(1:6)]   






