function [A,B,C,fn,t,w]=fzzyfft(T,M,Nf)
% 利用FFT，計算[0,T}區間上定義的時間波形的Fourier級數展開系數A,B和頻譜C,fn 。
%	T		時間波形周期
%	M		用作2的冪次
%	Nf		輸出諧波的階次，決定A,B的長度為(Nf+1)。Nf不要超2^(M-1)。
%	A,B		分別是Fourier級數中cos,sin展開項的系數。A(1)是直流量。
%	C		是定義在[-fs/2,fs/2]上的頻譜
%	t,w		是原時間波形數據對
if (nargin<2 | isempty(M));M=8;end										%		<9>
if nargin<3;Nf=6;end														%		<10>
N=2^M;
f=1/T;
w0=2*pi*f; 
dt=T/N;
n=0:1:(N-1);
t=n*dt;
w=time_fun(t,T);															%		<17>
%
W=fft(w);
cn=W/N;
%
z_cn=find(abs(cn)<1.0e-10);
cn(z_cn)=zeros(length(z_cn),1); 										%		<23>
cn_SH=fftshift(cn);														%		<24>
C=[cn_SH cn_SH(1)];
A(1)=C(N/2+1);
A(2:N/2+1)=2*real(C((N/2+2):end));
B(2:N/2+1)=-2*imag(C((N/2+2):end));
if Nf>N/2;error(['第三輸入宗量 Nf 應小于 ' int2str(N/2-1)]);end
A(Nf+2:end)=[];
B(Nf+2:end)=[];
n1=-N/2:1:N/2; 
fn=n1*f;