vtb3_3.m

在MATLAB中开发震动相关仿真时需要用到的

function [ap,bp]=vtb3_3(dat,t,n)
%VTB3_3  Fourier series approximation to a function.
%[a,b]=VTB3_3(dat,t,n) returns Fourier coefficients of a function
%  The coefficients are numerical approximations of the true
%  coefficients.
%    dat is a vector of data representing the function
%    t   is the corresponding time vector
%    n   is the desired number of terms to use in the Fourier series
% Intermediary plots show the impact of each successive term 
% on the total seres if no output arguments are specified.
% EXAMPLE: Manually performs the steps of the command vtb3_3(5)
% f=[ -1:.04:.96 1:-.04:-.96]'+1;
% t=(0:length(f)-1)/length(f)';
% plot(t,f)
% [a,b]=vtb3_3(f,t,5);
% vtb3_3(f,t,5)
%
% VTB3_3(N) displays the N term Fourier approximation to a 
% triangular input.  The approximation is plotted versus time 
% normalized by the period of the wave. 
%
% VTB3_3 displays the 5 term Fourier approximation to a 
% triangular input.  The approximation is plotted versus time 
% normalized by the period of the wave. 
%
% Note that these results are only an approximation, and the quality
% depends on the number of points used, and the proper selection of begining
% and end points. 


% Copyright Joseph C. Slater, Dec 1996
% Revised 02/29/00 - Now can run with no arguments.
% Revised 11/11/98 - Example changed to match default function 
%                    (Example 3.3.1)
% Revised 12/10/97 - Improved location of legend to avoid covering up data.
%                    Disclaimer on quality

if nargin==0
	vtb3_3(5)
  else
if nargin==1
  n=dat;
  tau1=0:.01:.5;
  Ftr1=(4*tau1-1);
  tau2=.51:.01:.99;
  Ftr2=3-4*tau2;
  t=[tau1 tau2]';
  dat=[Ftr1 Ftr2]';
end

if size(dat,1)==1
  dat=dat';
end

if size(t,1)==1
  t=t';
end

len=length(dat)/2;
plot(dat),grid on
fs=(fft(dat))/len;
fs(1:10);
a0=fs(1);

a=[a0; real(fs(2:length(fs/2)))];
b=-imag(fs(2:length(fs/2)));
len=len*2;
dt=2*pi/len;
tp=(0:dt:2*pi-dt)';
% fudge 12/10/97
%a=a+(t(2)-t(1))/2;


datapprox=a(1)/2+zeros(size(dat));
plot(t,dat,t,datapprox)%
grid on
aa=version;ll=length(aa);
grid on
context=['Contribution of terms n=' num2str(i-1)];
legend('Function','New Approximation')%
if nargout==0
  pause
end
%delete(h1)
%delete(h2)

for i=2:n+1
%  a(i)
%  b(i-1)
  newdat=a(i)*cos(tp*(i-1))+b(i-1)*sin(tp*(i-1));
  datapprox=datapprox+newdat;
  if nargout==0
    legend off
    plot(t,dat,t,datapprox,'o',t,datapprox-newdat,'x',t,newdat)
    grid on
%    pause
    aa=version;ll=length(aa);
    if ~strcmp(aa(ll-2:ll),'3.5')
%     figure
%     h1=uicontrol('style','pushbutton','units','normal','position',[.91 .95... 
%     .075 .05],'string','Print','callback','print');
%     h2=uicontrol('style','pushbutton','units','normal','position',[.91 .89...
%     .075 .05],'string','Close','callback','delete(gcf)');
    end
    grid on
    context=['Contribution of terms n=' num2str(i-1)];
    legend('Function','New Approximation','Old Approximation',context,0),pause
%    delete(h1)
%	delete(h2)
  end
end

legend off
plot(t,dat,t,datapprox),grid on

legend('Function','Approximation')
aa=version;ll=length(aa);
%if ~strcmp(aa(ll-2:ll),'3.5')
 %figure
% uicontrol('style','pushbutton','units','normal','position',[.91 .95 .075...
% .05],'string','Print','callback','print')
% uicontrol('style','pushbutton','units','normal','position',[.91 .89 .075...
% .05],'string','Close','callback','delete(gcf)')
%end

if nargout~=0
  ap=a(1:n+1);bp=b(1:n);
end
end