fseries.m

很多matlab的源代码

function [v,c,e,h1,mm]=fseries(x,t0,mm,ty,pw1,pw2)
% FSERIES Fourier series of periodic signal.
%
%	[V,C,E,H,M] = FSERIES(x,T,k,w,p1,p2) Fourier Series of signal x
%	T=[P td] where P = period(0,P), td = shift from origin [Default:td=0]
%	x is a row of n (n=odd) values over (0,P) INCLUDING BOTH ENDPOINTS
%	OR x is a string function of t over (0,P) e.g. '2*t' or 'tri(t-1)'
%	k = ARRAY of harmonics for intermediate reconstruct [Def: [1 3 5 32]]
%	w = smoothing window e.g. 'hamming' for reconstruction  [Def. None]
%	p1,p2 = parameters for some windows. Type HELP WINDOW for info.
%	V is a max(k)x7 array whose rows contain:
%	V = [k a(k) b(k) c(k) ang pwr cum.pwr] up to highest harmonic in k.
%	C = [c1 c2] returns convergence   E = [e1 e2] returns envelope
%	Use c1/(k^c2) e1*([sin(k*e2)/(k*e2)]^c2) to overplot.
%	H = one period reconstructions of x(t) at each harmonic in k.
%	H(1,:) is time. Smoothed reconstruction is last row if window is used.
%	M = harmonics at which x reconstructed. If window was used, M = [0 M].
%	FSERIES PLOTS RESULTS. If ANY ELEMENT of k is 0, plots are suppressed.
%
%	FSERIES (with no input arguments) invokes the following example:
%
%	% Find FS for a rect pulse with width=1, P=4 and a hamming window
%         >>fx='(t>=0 & t<1)+.5*(t==1)';fseries(fx,4,[1 3 25],'hamming')


% ADSP Toolbox: Version 2.0 
% For use with "Analog and Digital Signal Processing", 2nd Ed.
% Published by PWS Publishing Co.
%
% Ashok Ambardar, EE Dept. MTU, Houghton, MI 49931, USA
% http://www.ee.mtu/faculty/akambard.html
% e-mail: [e-mail]
% Copyright (c) 1998


%	USAGE: Use defaults & a full rectified sine: >>fseries('sin(pi*t)',1);
%	USAGE: For array x, use midpoint values  at discontinuities as follows: 
%	At ALL discontinuities over 0<=t<=P OR WITHIN the open interval 0<t<P. 
%	USAGE: For rect pulse with P=4, width=1:  x='(t>=0 & t<1)+.5*(t==1)'  

if nargin==0,help fseries,disp('Strike a key to see results of last example')
pause,fx='(t>=0 & t<1)+.5*(t==1)';fseries(fx,4,[1 3 25],'hamming');
return,end

vx=matverch;
%if exist('version')==5,mlv=4;else,mlv=3;end
tl=0;td=0;
narg=nargin;
if narg<3,mm=[1 3 5 32];end,
if narg<2,t0=1;end
vr=version;st=0;n=1024;
%if mlv==3,if vr(1)=='S',n=512;end,end
th=abs(t0(1));T=th-tl;
if prod(size(t0))==2,td=rem(t0(2),T);end
if isstr(x),st=1;t=tl:T/n:th;x=eval(x);else,n=length(x)-1;t=tl:T/n:th;end
x(1)=0.5*(x(1)+x(n+1));x(n+1)=x(1);
pl=max(mm);i=find(mm==0);
if ~isempty(i),pl=0;mm(i)=[];end
ln=length(mm);
if ln==0,error('Array of harmonics not specified'),return,end
mm=sort(mm);k=mm(ln);m=k+1;km=m;
if n<2*(km-1),error('# of harmonics must be < n/2'),return,end

j=sqrt(-1);F=fft(x(1:n))/n;d=F(1:km);d=d.*exp(-j*(0:km-1)*2*pi*td/T);
d1=real(d);d2=imag(d);d=d1.*(abs(d1)>eps)+j*d2.*(abs(d2)>eps);
c=abs([d(1) 2*d(2:km)]);p=c.*c;p=[p(1) p(2:km)/2];pc=cumsum(p);
a=real([d(1) 2*d(2:km)]);b=-2*imag(d);%%tk=180*atan2(-b,a)/pi;
tk=180*angle(d)/pi;
hk=0:km-1;i=find(tk==180);[sm,sn]=size(i);t1=-180*ones(sm,sn);tk(i)=t1;
v0=[hk',a',b',c',tk',p',pc'];v=v0(1:m,:);
if st==1,x=x(1:4:n+1);t=t(1:4:n+1);end

o=200;
%if mlv==3,if vr(1)=='S',o=100;end,end
u=tl:T/o:th;d=2*pi*u/T;
h=zeros(ln,o+1);
y=c(1)*ones(1,o+1);
m1=mm;j0=0;
if narg==4,sc=windo(ty,2*k+1);end,
if narg==5,sc=windo(ty,2*k+1,pw1);end
if narg==6,sc=windo(ty,2*k+1,pw1,pw2);end
a(1)=[];b(1)=[];
if narg>3,sc=sc(k+2:2*k+1);y1=y;end,
yh=h;
yh(1,:)=y;
for jj=1:k,
sm=(a(jj)*cos(jj*d)+b(jj)*sin(jj*d));y=y+sm;
if jj<=7,yh(jj+1,:)=sm;end
if jj==m1(1),j0=j0+1;h(j0,:)=y;m1(1)=[];end,
if narg>3,y1=y1+sm*sc(jj);end
end,

if narg>3,h=[h;y1];mm=[0 mm];end,
h1=[u;h];sz=min(k,7);sz=num2str(sz);

if km>16

p=0;tol=5e-3;p1=0;y=c(2:km);ya=y/y(1);z=find(ya<tol);s=length(z);
if s>1,w=z(1);if any((z/w)-(1:s))~=0,z=find(ya<tol/20);end,end
s=length(z);i=1:km-1;i(z)=[];i(1:4)=[];y1=y(i);l=0;j=length(find(i<16));
c1=fliplr(sort(y1));c=[];e=[];
if y1==c1,
p=i(j-3);r=i(j);a=y(p);b=y(r);
s1=(r/p).^(1:3);d=abs(a/b-s1);
[z1,g1]=min(d);q=b*(r^g1);c=[q g1];l=1;
end

p=0;
if s>1,w=z(1);if rem(z,w)==0*(1:s),pa=w/2;p=fix(pa);r=p+w;p1=1;end,end
if p>0,
a=y(p);b=y(r);g=[1:3];wa=pi/w;qa=a./(sinc(p*wa/pi).^g);
qb=abs(b./(sinc(r*wa/pi).^g));df=abs(qa-qb)./qa;[q0,i]=min(df);
if q0<=tol,q=qa(i)*(sinc(pa*wa/pi)^i)*(pa^i);c=[q i];l=1;
y=y(1:w);y=-diff(y);if all(y>.0005),e=[qa(i) wa];l=2;end,end
end

end

if pl>0,disp('STRIKE A KEY BETWEEN PLOTS'),pause(2)
ds=min(17,m);vs=v0(1:ds,:);kk=pl;dp=pl+1;
vp=v0(1:dp,:);kp=vp(:,1);
if l==2,z=0:.05:kk;d=e(1)*(abs(sinc(z*e(2)/pi)).^c(2));end
if l>0,z1=.5:.05:kk;d1=c(1)./(z1.^c(2));end
subplot,hold off,
if l==2,plot(z,d,'-.'),hold on,end,
dtplot(kp,vp(:,4),'o'),grid,title('magnitude vs k'),hold on
ax=axis;
if l>0,axis(ax);plot(z1,d1,':'),axis(ax);end,pause,hold off

if vx<4, eval('clg');else,eval('clf');end

ax=[ax(1) ax(2) -180 180];axis(ax);
dtplot(kp,vp(:,5),'o'),axis(ax);axesn;
grid,title('phase vs k'),pause,hold off
plot(kp,vp(:,7),'*',kp,vp(:,6),'o'),hold on,dtplot(kp,vp(:,6),'.'),grid
title('power(o) & cumulative power(*) vs k'),pause,hold off

plot(u,yh),grid,title(['dc & first ' sz ' harmonics']),pause
T0=t+td;
plot(u,h,':'),hold on,plot(T0,x),grid


title('One period reconstructions(:) of x(t)(-)'),hold off,pause
zl=num2str(k);ha=h(ln,:);
plot(u,ha),grid,title(['Final reconstruction to ' zl ' harmonics']),pause
if mm(1)==0,hs=h(ln+1,:);
plot(u,ha,'-',u,hs,'-.'),grid
title(['Actual(-) and smoothed(-.) reconstruction to ' zl ' harmonics']),pause
end
v1=['    Harmonic   a(k)      b(k)      Mag.      Ang.     Power    CUM.PWR'];
disp(' '),disp(v1),disp(vs)
if l>0,disp(['convergence is 1/(k^' num2str(c(2)) ')']),end
end