clear;%an example to use Fourier series to approximate a periodic square pulse sequence in one time periodA = 1;%magnitudeTs = 5e-4;%sample interval of the continuous time signalTp = 1;     F0 = 1/Tp;%period and frequency of the periodic signaltau = 0.025;%width of the square pulset = -Tp/2:Ts:Tp/2;%generate the time vectort_len = length(t);%obtain the length of the time vectormid_point = (t_len+1)/2;%obtain the mid point indexfor i = 1:t_len    t_now = (i-mid_point)*Ts;%obtain the corresponding time of the element t(i);    %generate the square pulse    if (abs(t_now)>tau/2)        xt(i) = 0;    else        xt(i) = A;    endend%x_app is the approximation of xt using only the DC component.%since C_0 = A*tau/Tp, we have the following.x_app(1:t_len) = A*tau/Tp;%x_app1 is the approximation of xt using the DC component and 1 sinusoid.x_app1 = x_app;K_max = 1;for k = 1:K_max    coef(k) = A*tau/Tp*sin(pi*tau*k*F0)/(pi*tau*(k*F0));    x_app1 = x_app1+2*coef(k)*cos(2*pi*k*F0*t);end%x_app2 is the approximation of xt using the DC component and 10 sinusoid.x_app2 = x_app;K_max = 10;for k = 1:K_max    coef(k) = A*tau/Tp*sin(pi*tau*k*F0)/(pi*tau*(k*F0));    x_app2 = x_app2+2*coef(k)*cos(2*pi*k*F0*t);end%x_app3 is the approximation of xt using the DC component and 100 sinusoid.x_app3 = x_app;K_max = 100;for k = 1:K_max    coef(k) = A*tau/Tp*sin(pi*tau*k*F0)/(pi*tau*(k*F0));    x_app3 = x_app3+2*coef(k)*cos(2*pi*k*F0*t);end%x_app4 is the approximation of xt using the DC component and 200 sinusoid.x_app4 = x_app;K_max = 200;for k = 1:K_max    coef(k) = A*tau/Tp*sin(pi*tau*k*F0)/(pi*tau*(k*F0));    x_app4 = x_app4+2*coef(k)*cos(2*pi*k*F0*t);end%the above code can be further simplified by using "function", however, a%simple copy and paste of code can get the same result.%plot the signalfigure(1);plot(t, xt);hold on;    grid on;plot(t, x_app,'r');plot(t, x_app1,'g');plot(t, x_app2,'m');plot(t, x_app3,'c');plot(t, x_app4,'k');ylabel('magnitude');xlabel('t (sec)');legend('square sequence', 'DC appx', '1-sinusoid appx', '10-sinusoid appx', '100-sinusoid appx', '200-sinusoid appx');%compute the error signalerr = xt-x_app;err1 = xt-x_app1;err2 = xt-x_app2;err3 = xt-x_app3;err4 = xt-x_app4;%plot the error signalfigure(2);plot(t, err,'r');hold on;    grid on;plot(t, err1,'g');plot(t, err2,'m');plot(t, err3,'c');plot(t, err4,'k');ylabel('magnitude');xlabel('t (sec)');legend('err in DC appx', 'err in 1-sinusoid appx', 'err in 10-sinusoid appx', 'err in 100-sinusoid appx', '200-sinusoid appx');%compute the error powerP_err = err*err'*Ts/Tp;P_err1 = err1*err1'*Ts/Tp;P_err2 = err2*err2'*Ts/Tp;P_err3 = err3*err3'*Ts/Tp;P_err4 = err4*err4'*Ts/Tp;%show the error powerP_errP_err1P_err2P_err3P_err4%plot the Fourier series of the square pulse, you can change the parameters%and see how the plot changes.figure(3);stem(-F0*K_max:F0:F0*K_max, [wrev(coef), A*tau/Tp, coef],'k');ylabel('C_k');xlabel('f (Hz)');hold on;    grid on;legend('Fourier coefficients of the square pulse',4);xlim([-50, 50]);