%========================================================================
%               周期信號的傅里葉級數   程序名：CTFS
%   本程序根據三角形式傅里葉級數的系數的計算公式，采用符號計算法計算出信號
%   的傅里葉級數的系數，然后再根據你輸入的諧波次數合成（逼近）原周期信號。
% 	Nf  	合成信號的諧波的次數,由實驗者從鍵盤輸入，諧波次數可以不限，但
%   建議取50以內，因為次數太大，計算機的運算量越大，所需計算時間越長。
%	A_sym	第1元素是直流項，其后元素依次是1,2,3...次諧波cos項展開系數
%	B_sym	第2,3,4,...元素依次是1,2,3...次諧波sin項展開系數
%   tao     脈沖寬度。周期矩形信號的寬度為0.5s，鋸齒波為4s，三角波為2s。
%   T=4     信號周期為4s
%========================================================================
clear, close all
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('2*(t)/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
Fn=int(x*exp(-j*2*pi*n*t/T)/T,t,-T/2,T/2);
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碼數值數組
    F_sym(k+1)=double(vpa(subs(Fn,n,k)));
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項展開系數
Fn=F_sym
Fn(1)=a_n(1);%A0%(1);
%===========================================================================
%                 下段程序計算原周期信號的時域波形
%===========================================================================
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.*(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('Approximation of aperiodic signal')
%axis([-8,8,-1.2,1.2]);
%====================================
%      下段程序繪制信號的頻譜圖
%====================================
figure(2);
n=0:Nf;
c_n=sqrt(a_n.^2+b_n.^2);
subplot(2,1,1);
stem(n,abs(Fn),'.r');
%stem(n,c_n,'.r');
title('Amplitude spectra of aperiodic signal |cn|');
%title('Magnitude spectra of aperiodic signal  Fn');
fi=-atan(abs(imag(Fn))./(abs(real(Fn))+eps));
for L=1:length(n);
    
if b_n(L)==0
    if a_n(L)<=-0.0001
        fi(L)=pi;
    else
        fi(L)=0;
    end
end
end
subplot(2,1,2);
stem(n,fi,'.k');
title('Phase spectra of aperiodic signal');
xlabel('index n');
%====================================