clear
syms t  n   k  x 
T=4;
Nf=input('請輸入信號合成的諧波次數：');
%選擇信號類型，確定出信號一個周期內的符號表達式
q=input('請選擇信號代碼（1-周期矩形波信號，2-周期三角波信號，3-周期鋸齒波信號）：');
p=(q~=1)&(q~=2)&(q~=3);
while p==1;
      q=input('請選擇信號代碼（1-周期矩形波信號，2-周期三角波信號，3-周期鋸齒波信號）：');
      p=(q~=1)&(q~=2)&(q~=3);
end
if q==1
    x=sym('Heaviside(t+1)-Heaviside(t-1)');
end
if q==2
    x=sym('(t+2)/4*(Heaviside(t+2)-Heaviside(t-2))');
end
if q==3
    x=sym('(-abs(t)+1)*(Heaviside(t+1)-Heaviside(t-1))');
end
% 下面的程序段計算傅里葉級數的系數
A0=int(x/T,t,-T/2,T/2);                   %求出三角函數展開系數A0
As=int(2*x*cos(2*pi*n*t/T)/T,t,-T/2,T/2); %求出三角函數展開系數As
Bs=int(2*x*sin(2*pi*n*t/T)/T,t,-T/2,T/2); %求出三角函數展開系數Bs
A_sym(1)=double(vpa(A0));                %獲取串數組A0所對應的ASC2碼數值數組
for k=1:Nf
A_sym(k+1)=double(vpa(subs(As,n,k))); %獲取串數組A所對應的ASC2碼數值數組
    B_sym(k+1)=double(vpa(subs(Bs,n,k))); %獲取串數組B所對應的ASC2碼數值數組
end
disp('顯示an列表：第1元素是直流項，其后元素依次是1,2,3...次諧波cos項展開系數');
a_n=A_sym %輸出cn:第1元素是直流項，其后元素依次是1,2,3...次諧波cos項展開系數
disp('顯示bn列表：第2,3,4,...元素依次是1,2,3...次諧波sin項展開系數');
b_n=B_sym  %輸出dn:第2,3,4,...元素依次是1,2,3...次諧波sin項展開系數
% 下段程序計算原周期信號的時域波形
t=-8:0.01:8;
f=a_n(1);
for l=1:Nf;
    f=f+a_n(l+1).*cos(2*l*pi*t/T)+b_n(l+1).*sin(2*l*pi*t/T); 
end;
if q==1
    y=0;
   for r=-2:2;
      y=y+u(t+r*4+1)-u(t+r*4-1);
   end;
end
if q==2
    y=0;
   for r=-2:2;
      y=y+(t+r*4+2)/4.*(u(t+r*4+2)-u(t+r*4-2));
   end;
end
if q==3
    y=0;
   for r=-2:2;
      y=y+(-abs(t+r*4)+1).*(u(t+r*4+1)-u(t+r*4-1));
   end;
end
% 下段程序繪制原周期信號和逼近信號的波形圖
figure(1);
clf
plot(t,y,'r:');hold on;
plot(t,f);
title('非正弦周期信號形波的合成')
axis([-8,8,-0.4,1.3]);
% 下段程序繪制信號的頻譜圖
figure(2);
n=0:Nf;
c_n=sqrt(a_n.^2+b_n.^2);
subplot(2,1,1);
stem(n,c_n,'.r');
title('非正弦周期信號的幅度譜|cn|');
fi=-atan(b_n./a_n);
subplot(2,1,2);
stem(n,fi,'.k');
title('非正弦周期信號的相位譜');