fzzyfft.m

这本书是matlab经典书籍

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;							%使总采样点是2的整数倍 
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);						%给出n=0,1,...,N-1上的DFT数据值 
cn=W/N; 						%据式(4.13.2.2-2)计算n=0,1,...,N-1上的FS系数 
%
z_cn=find(abs(cn)<1.0e-10);%寻找有限字长运算而产生（原应为0）的"小"复数
cn(z_cn)=zeros(length(z_cn),1); %强制那些"小"复数为0						<23>
cn_SH=fftshift(cn);			%据式(4.13.2.2-3)计算							<24>
									%n=-N/2,...,-1,0,1,...,(N/2)-1上的FS系数
C=[cn_SH cn_SH(1)];			%形成关于0对称的（N+1）个FS系数
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;				%产生总点数为（N+1）关于0对称的序列 
fn=n1*f; 						%关于0对称的频率分度序列