function [t,y,S,A,B]=fzzyquad(a,T,Nf,K)
% 用數值積分求Fourier級數展開系數A,B和近似波形S 。如調用本指令時，沒有任何
% 輸出宗量，那么將自動隨積分進程動態地畫出原波形和各階近似波形。
%	a		是被展開函數的時間區間的左端。
%	T		是被展開函數的周期。
%	Nf		是所需展開的最高諧波階次
%	K		繪波形圖所采用的數據長度
%	t,y		是被展開的時間函數數據。
%	S		是一個矩陣。它每行的長度與t,y相同。
%				它的第（n+1)行是所有n次及更低次諧波迭加生成的近似波形。
%	A,B		是展開系數向量。A(n+1),B(n+1)分別存放cos，sin項n次諧波的展開系數。
%				A(1)是直流項
if nargin<4;K=100;end
if (nargin<3|isempty(Nf));Nf=15;end
k=1:K;
t=a+k*T/length(k);
y=time_fun(t,T);
S=zeros(Nf+1,length(k));
a0=mean(y); 
S(1,:)=a0;
A=zeros(1,Nf+1);B=zeros(1,Nf+1);
A(1)=a0;
tolzzy=[0.001,0.018];
for n = 1:Nf
   A(n+1) = quad8('cos_y',a,a+T,tolzzy,[],n,T)/T*2;
   B(n+1) = quad8('sin_y',a,a+T,tolzzy,[],n,T)/T*2;
   S(n+1,:)=S(n,:)+A(n+1)*cos(2*n*pi*t/T)+B(n+1)*sin(2*n*pi*t/T);
   if nargout==0
      plot(t,y,'b');hold on
      plot(t,S(n+1,:),':r');
      aa=axis;dx=aa(2)-aa(1);dy=aa(4)-aa(3);
      text(aa(1)+0.1*dx,aa(3)+.75*dy,'n = ')
      text(aa(1)+0.2*dx,aa(3)+.75*dy,num2str(n));
      pause(1);hold off
   end
end 
%------------------------------
function y=time_fun(t,T)
% 這是需要用戶自己編寫的待展開的時間函數 y(t+k*T)=y(t)
% 編寫程序時，一定要注意使表達式適于數組運算
%	t		是時間數組
%	T		是周期
y=zeros(size(t));ii=find(t>=0.5 & t<=1.5);
y(ii)=ones(size(ii)).*(t(ii)-0.5);y(y==1.0)=0.5;