%calculation of Fourier Coef. for Example #3

clear
syms L x n
%Compute c0

Int1=(1/L)*3*cos(15*x);

c0=int(Int1,x,0,L)

%Compute cn

Int2=(2/L)*3*cos(15*x)*cos(2*pi*n/L*x);

cn=int(Int2,x,0,L);

cn=simplify(cn);

%Compute dn

Int3=(2/L)*3*cos(15*x)*sin(2*pi*n/L*x);

dn=int(Int3,x,0,L);

dn=simplify(dn);

%Further Simplify using sin(pi*n)=0 and cos(pi*n)=((-1)^n)

C=char(cn);

C=strrep(C,'sin(pi*n)','0');

C=strrep(C,'cos(pi*n)','((-1)^n)');

cn=sym(C);

cn=simplify(cn);

C=char(cn);

C=strrep(C,'(-1)^(2*n)','1');

C=strrep(C,'(-1)^(1+2*n)','(-1)');

cn=sym(C);

cn=simplify(cn)

D=char(dn);

D=strrep(D,'sin(pi*n)','0');

D=strrep(D,'cos(pi*n)','((-1)^n)');

dn=sym(D);

dn=simplify(dn);

D=char(dn);

D=strrep(D,'(-1)^(2*n)','1');

D=strrep(D,'(-1)^(1+2*n)','(-1)');

dn=sym(D);

dn=simplify(dn)

%Create inline functions for the coefficients

cnf=inline(char(vectorize(cn)),'n','L');

dnf=inline(char(vectorize(dn)),'n','L');

%Compute the first 20 coefficients and set L=1

c0=subs(c0,L,1)

cn=subs(cn,L,1)

dn=subs(dn,L,1)

M=20,L=1

n=1:M;

an(n)=cnf(n,L);

bn(n)=dnf(n,L);

%Plot the series and the function on [0, 1]

figure(1),clf reset

FSPlot(L,c0,an,bn,0,1)

hold on

x=linspace(0,1);

plot(x,3*cos(15*x),'r')

%Plot the series and the function on [0, 3]

figure(2),clf reset

FSPlot(L,c0,an,bn,0,3)

hold on

x=linspace(0,3,300);

plot(x,3*cos(15*x),'r')

%Plot the coefficients

figure(3), clf reset

axes('pos',[.1 .35 .8 .6])

plot(n,an,'r*',n,bn,'b+')

axis([0 M -2.1 2.1])

legend('Cosine Coefs.','Sine Coefs.')

xlabel('frequency index'), ylabel('amplitude')

p=get(gca,'position')

axes('pos',[p(1) .15 p(3) .01],'xlim',[0 M]*2*pi/L)

xlabel('radians/unit time')

%Plot the square root of the sum of the squares of cn and dn

figure(4), clf reset

axes('pos',[.1 .35 .8 .6])

plot(n,sqrt(an.^2+bn.^2),'b*')

xlabel('frequency index') , ylabel('amplitude')

p=get(gca,'position')

axes('pos',[p(1) .15 p(3) .01],'xlim',[0 M]*2*pi/L)

xlabel('radians/unit time')

%calculation of Fourier Coef. for Example #4

clear

syms L x n a

%Compute c0

Int1=(1/L);

c0=int(Int1,x,a,L)

c0=subs(c0,L,1);

c0=subs(c0,a,2/3)

%Compute cn

Int2=(2/L)*cos(2*pi*n/L*x);

cn=int(Int2,x,a,L)

cn=simplify(cn);

%Compute dn

Int3=(2/L)*sin(2*pi*n/L*x);

dn=int(Int3,x,a,L)

dn=simplify(dn);

%Further Simplify using sin(pi*n)=0 and cos(pi*n)=((-1)^n)

C=char(cn);

C=strrep(C,'sin(pi*n)','0');

C=strrep(C,'cos(pi*n)','((-1)^n)');

cn=sym(C);

cn=simplify(cn);

C=char(cn);

C=strrep(C,'(-1)^(2*n)','1');

C=strrep(C,'(-1)^(1+2*n)','(-1)');

cn=sym(C);

cn=simplify(cn)

D=char(dn);

D=strrep(D,'sin(pi*n)','0');

D=strrep(D,'cos(pi*n)','((-1)^n)');

dn=sym(D);

dn=simplify(dn);

D=char(dn);

D=strrep(D,'(-1)^(2*n)','1');

D=strrep(D,'(-1)^(1+2*n)','(-1)');

dn=sym(D);

dn=simplify(dn)

%Substitute in the values for L and a

cn=subs(cn,L,1);

cn=subs(cn,a,2/3)

dn=subs(dn,L,1);

dn=subs(dn,a,2/3)

%Create inline functions for the coefficients

cnf=inline(char(vectorize(cn)))

dnf=inline(char(vectorize(dn)))

%Compute the first 20 coefficients

M=20, a=2/3

n=1:M;

an(n)=cnf(n);

bn(n)=dnf(n);

%Plot the series and the function on [0, 1]

figure(1),clf reset

FSPlot(1,c0,an,bn,0,1)

hold on

x=linspace(0,1);

y=zeros(1,length(x)); index=find(x>=a);y(index)=1;

plot(x,y,'r')

%Plot the series and the function on [0, 3]

figure(2),clf reset

FSPlot(1,c0,an,bn,0,3)

%Plot the coefficients

figure(3), clf reset

plot(n,an,'r*',n,bn,'b+')

legend('Cosine Coefs.','Sine Coefs.')

%Plot the square root of the sum of the squares of cn and dn

figure(4), clf reset

plot(n,sqrt(an.^2+bn.^2),'b*')

%calculation of Fourier Coef. for Example #5

clear

syms L x n a b

%Compute c0

Int1=(1/L*(a+b*x));

c0=int(Int1,x,0,L)

c0=subs(c0,L,1);

c0=subs(c0,a,2);

c0=subs(c0,b,1/4)

%Compute cn

Int2=(2/L)*(a+b*x)*cos(2*pi*n/L*x);

cn=int(Int2,x,0,L)

cn=simplify(cn);

%Compute dn

Int3=(2/L)*(a+b*x)*sin(2*pi*n/L*x);

dn=int(Int3,x,0,L)

dn=simplify(dn);

%Further Simplify using sin(pi*n)=0 and cos(pi*n)=((-1)^n)

C=char(cn);

C=strrep(C,'sin(pi*n)','0');

C=strrep(C,'cos(pi*n)','((-1)^n)');

cn=sym(C);

cn=simplify(cn);

C=char(cn);

C=strrep(C,'(-1)^(2*n)','1');

C=strrep(C,'(-1)^(1+2*n)','(-1)');

cn=sym(C);

cn=simplify(cn)

D=char(dn);

D=strrep(D,'sin(pi*n)','0');

D=strrep(D,'cos(pi*n)','((-1)^n)');

dn=sym(D);

dn=simplify(dn);

D=char(dn);

D=strrep(D,'(-1)^(2*n)','1');

D=strrep(D,'(-1)^(1+2*n)','(-1)');

dn=sym(D);

dn=simplify(dn)

%Substitute in the values for L=1, a=2, b=0.25

cn=subs(cn,L,1);

cn=subs(cn,a,2);

cn=subs(cn,b,1/4)

dn=subs(dn,L,1);

dn=subs(dn,a,2);

dn=subs(dn,b,1/4)

%Create inline functions for the coefficients

dnf=inline(char(vectorize(dn)))

%Compute the first 20 coefficients

M=20, A=2, B=1/4, L=1

n=1:M;

an(n)=zeros(1,M);

bn(n)=dnf(n);

%Plot the series and the function on [0, 1]

figure(1),clf reset

FSPlot(1,c0,an,bn,0,1)

hold on

x=linspace(0,1);

plot(x,A+B*x,'r')

%Plot the series and the function on [0, 3]

figure(2),clf reset

FSPlot(1,c0,an,bn,0,3)

%Plot the coefficients

figure(3), clf reset

axes('pos',[.1 .35 .8 .6])

plot(n,an,'r*',n,bn,'b+')

axis([0 M -0.1 0.1])

legend('Cosine Coefs.','Sine Coefs.')

xlabel('frequency index') , ylabel('amplitude')

p=get(gca,'position')

axes('pos',[p(1) .15 p(3) .01],'xlim',[0 M]*2*pi/L)

xlabel('radians/unit time')

%Plot the square root of the sum of the squares of cn and dn

figure(4), clf reset

axes('pos',[.1 .35 .8 .6])

plot(n,sqrt(an.^2+bn.^2),'b*')

xlabel('frequency index') , ylabel('amplitude')

p=get(gca,'position')

axes('pos',[p(1) .15 p(3) .01],'xlim',[0 M]*2*pi/L)

xlabel('radians/unit time')

%calculation of Fourier Coef. for Example #6

clear

syms L x n

%Compute c0

Int1=(1/L)*(3*cos(15*x)*16*x^2*(x-L)^2);

c0=int(Int1,x,0,L);

c0=subs(c0,L,1)

%Compute cn

Int2=(2/L)*(3*cos(15*x)*16*x^2*(x-L)^2)*cos(2*pi*n/L*x);

cn=int(Int2,x,0,L);

cn=subs(cn,L,1);

cn=simplify(cn);

%Compute dn

Int3=(2/L)*(3*cos(15*x)*16*x^2*(x-L)^2)*sin(2*pi*n/L*x);

dn=int(Int3,x,0,L);

dn=subs(dn,L,1);

dn=simplify(dn);

%Further Simplify using sin(pi*n)=0 and cos(pi*n)=((-1)^n)

C=char(cn);

C=strrep(C,'sin(pi*n)','0');

C=strrep(C,'cos(pi*n)','((-1)^n)');

cn=sym(C);

cn=simplify(cn);

C=char(cn);

C=strrep(C,'(-1)^(2*n)','1');

C=strrep(C,'(-1)^(1+2*n)','(-1)');

cn=sym(C);

cn=simplify(cn)

D=char(dn);

D=strrep(D,'sin(pi*n)','0');

D=strrep(D,'cos(pi*n)','((-1)^n)');

dn=sym(D);

dn=simplify(dn);

D=char(dn);

D=strrep(D,'(-1)^(2*n)','1');

D=strrep(D,'(-1)^(1+2*n)','(-1)');

dn=sym(D);

dn=simplify(dn)

%Create inline functions for the coefficients

cnf=inline(char(vectorize(cn)))

dnf=inline(char(vectorize(dn)))

%Compute the first 20 coefficients

M=20,L=1

n=1:M;

an(n)=cnf(n);

bn(n)=dnf(n);

%Plot the series and the function on [0, 1]

figure(1),clf reset

FSPlot(1,c0,an,bn,0,1)

hold on

x=linspace(0,1);

plot(x,(3*cos(15*x)*16.*x.^2.*(x-L).^2),'r')

plot(x,(3*cos(15*x)),'k')

%Plot the coefficients

figure(3), clf reset

axes('pos',[.1 .35 .8 .6])

plot(n,an,'r*',n,bn,'b+')

axis([0 M -2.1 2.1])

legend('Cosine Coefs.','Sine Coefs.')

xlabel('frequency index') , ylabel('amplitude')

p=get(gca,'position')

axes('pos',[p(1) .15 p(3) .01],'xlim',[0 M]*2*pi/L)

xlabel('radians/unit time')

%Plot the square root of the sum of the squares of cn and dn

figure(4), clf reset

axes('pos',[.1 .35 .8 .6])

plot(n,sqrt(an.^2+bn.^2),'b*')

xlabel('frequency index') , ylabel('amplitude')

p=get(gca,'position')

axes('pos',[p(1) .15 p(3) .01],'xlim',[0 M]*2*pi/L)

xlabel('radians/unit time')

%Computations for Example #7

%Numerical Computation of the first M Fourier Coeffs

%of the function f on the interval [0, L]

clear,

figure(1), clf reset,

ind=0; col=['k','b','r','g'];sy='*';

for L=[1,2,4,8]

M=ceil(30*L/(2*pi)); ind=ind+1;

f='3*cos(15*t)'

IntC0=inline(f,'t')

IC=[f,'.*cos(2*pi*n*t/L)']

IntC=inline(IC,'t','n','L')

ID=[f,'.*sin(2*pi*n*t/L)']

IntD=inline(ID,'t','n','L')

for i=1:M

c(i)=2/L*quadl(IntC,0,L,1.e-6,[],i,L);

d(i)=2/L*quadl(IntD,0,L,1.e-6,[],i,L);

end

c0=1/L*quad(IntC0,0,L);

%Plot the square root of the sum of the squares of cn and dn

m=0:M; rr=2*pi/L*m;

pp=[abs(c0),sqrt(c.^2+d.^2)];

str=[sy,col(ind)];

plot(rr,pp,str,rr,pp,col(ind))

hold on

axis([0,30,0,3])

end

xlabel('radians/unit time')

ylabel('amplitude')

text(20,2.6,'Figure 14.')

figure(2), clf reset

ph=[0,atan2(d,c)];

plot(rr,ph,'*b',rr,ph,'b')

xlabel('radians/unit time')

ylabel('phase angle in radians')

axis([0,30,-pi,pi])

grid on

text(6,2.5,'Figure 15.')

%Computation for Example #8

%Computes an approximate Fourier Series for the

%PL approximation to the data on the interval [0,L]

clear

N=100;L=8; %N=number of terms in the series

f=inline('3*cos(15*x)'), %Test function

n=71; x=linspace(0,L,n);%n=number of samples

y=(f(x)+6*(rand(1,n)-.5));

[c0,cn,dn]=FSPLBasisF(N,L,x(2));

c0=y(1)*c0;cn=y(1)*cn;dn=y(1)*dn;

for i=2:n-1

[cc0,ccn,dcn]=FSPLBasis(N,L,x(i-1),x(i),x(i+1));

c0=y(i)*cc0+c0;cn=cn+y(i)*ccn;dn=dn+y(i)*dcn;

end

[cc0,ccn,dcn]=FSPLBasisL(N,L,x(n-1));

c0=c0+y(n)*cc0;cn=cn+y(n)*ccn;dn=dn+y(n)*dcn;

figure(1), clf reset

FSPlot(L,c0,cn,dn,0,L)

grid on

hold on

plot(x,y,'ro','markersize',7,'linewidth',2)

xx=linspace(0,L,500);

plot(xx,f(xx),'k','linewidth',2)

axis([0,L,-6,6])

text(5,3.5,'Figure 20.')

c0

M=ceil(30*L/(2*pi));

nn=0:M;

figure(2), clf reset

mag=abs(cn+dn*sqrt(-1));

plot(2*pi/L*nn,[abs(c0),mag(1:M)],'k')

grid on

axis([0,35,0,3])

xlabel('radians/unit time')

ylabel('amplitude')

text(20,1.5,'Figure 21.')

figure(3), clf reset

ph=[0,atan2(dn,cn)];

rr=2*pi/L*nn;

plot(rr,ph(1:M+1),'*b',rr,ph(1:M+1),'b')

xlabel('radians/unit time')

ylabel('phase angle in radians')

axis([0,30,-pi,pi])

grid on

text(6,2.5,'Figure 22.')

z=fft(y);

figure(4), clf reset

plot(2*pi/L*(0:n-1),[abs(z(1))/n,2*abs(z(2:n))/n],'k')

grid on

axis([0,35,0,3])

%Example #9

clear
x=[1 -2 3 6 0 -1];
X=fft(x);
xx=ifft(X);
R=[(1:length(x))',x',X',xx']

%Example #10

clear, clf reset
x=[1 -2 3 6 0, -1]; L=1;
N=length(x);
Data=[(0:N)*L/N;x, x(1)]
plot(Data(1,:),Data(2,:),'k*','markersize',10,'linewidth',2)
hold on
X=fft(x), M=floor((N+1)/2)
c0=X(1)/N, cn=2*real(X(2:M))/N,dn=-2*imag(X(2:M))/N
if M==N/2
cn(M)=X(M+1)/N, dn(M)=0
end
FSplot(L,c0,cn,dn,0,L)
grid on

%calculation of Fourier Coef. for PL Basis Functions

syms a b c L x n w
%Calculation of cn and dn for a basis function Theta defined on a,b,c
Int1=(2/L)*(x-a)/(b-a)*cos(w*x);
Int2=(2/L)*(c-x)/(c-b)*cos(w*x);
cn=int(Int1,x,a,b)+int(Int2,x,b,c);
cn=factor(cn);
cn=subs(cn,'w','2*pi*n/L')
Int3=(2/L)*(x-a)/(b-a)*sin(w*x);
Int4=(2/L)*(c-x)/(c-b)*sin(w*x);
dn=int(Int3,x,a,b)+int(Int4,x,b,c);
dn=factor(dn);
dn=subs(dn,'w','2*pi*n/L')
%Calculation of the special case of the first basis function
cnfirst=int(Int2,x,0,c);
%cnfirst=simplify(cnfirst);
cnfirst=subs(cnfirst,'b',0);
cnfirst=subs(cnfirst,'w','2*pi*n/L')
%Calculation of the special case of the last basis function
dnfirst=int(Int4,x,0,c);
%dnfirst=simplify(dnfirst);
dnfirst=subs(dnfirst,'b',0);
dnfirst=subs(dnfirst,'w','2*pi*n/L')
cnlast=int(Int1,x,a,L);
%cnlast=simplify(cnlast);
cnlast=subs(cnlast,'w','2*pi*n/L');
cnlast=subs(cnlast,'b',sym('L'));
cnlast=simplify(cnlast);
dnlast=int(Int3,x,a,L);
%dnlast=simplify(dnlast);
dnlast=subs(dnlast,'w','2*pi*n/L');
dnlast=subs(dnlast,'b',sym('L'));
dnlast=simplify(dnlast);
%Further Simplify using sin(pi*n)=0 and cos(pi*n)=((-1)^n)
C=char(cnlast);
C=strrep(C,'sin(2*pi*n)','0');
C=strrep(C,'cos(2*pi*n)','1');
cnlast=sym(C);
cnlast=simplify(cnlast)
D=char(dnlast);
D=strrep(D,'sin(2*pi*n)','0');
D=strrep(D,'cos(2*pi*n)','1');
dnlast=sym(D);
dnlast=simplify(dnlast)

function FSPlot(L,c0,cn,dn,a,b)
%Plots the truncated Fourier Series having coefficients
%c0, cn and dn defined on the interval [0, L] and
%ploted on the interval [a, b]
%c0, cn, dn must be row vectors containing the coefficients
M=length(cn);
x=linspace(a,b,15*M);
n=1:M;
y=c0+cn*cos(2*pi*n'*x/L)+dn*sin(2*pi*n'*x/L);
plot(x,y)