fourier_series.m

matlab 源程序

clear;
%an example to use Fourier series to approximate a periodic square pulse sequence in one time period

A = 1;%magnitude
Ts = 5e-4;%sample interval of the continuous time signal
Tp = 1;     F0 = 1/Tp;%period and frequency of the periodic signal
tau = 0.025;%width of the square pulse

t = -Tp/2:Ts:Tp/2;%generate the time vector
t_len = length(t);%obtain the length of the time vector
mid_point = (t_len+1)/2;%obtain the mid point index
for 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;
    end
end

%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 signal
figure(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 signal
err = xt-x_app;
err1 = xt-x_app1;
err2 = xt-x_app2;
err3 = xt-x_app3;
err4 = xt-x_app4;

%plot the error signal
figure(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 power
P_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 power
P_err
P_err1
P_err2
P_err3
P_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]);