自适应滑窗prony非平稳间谐波谐波参数测量程序

%自适应滑窗PRONY分析算法 tr=[0:0.001:1];fs=1000;l=length(tr); xtest=[10*sin(2*pi*90*tr(1:250)+pi/3)+2*sin(2*pi*60*tr(1:250)+pi/4),4*sin(2*pi*10*tr(251:l))+6*sin(2*pi*30*tr(251:l))]; p=4; np=fs/50;%十个工频周期的采样点数 eth=0.01;%eth---误差阈值 timerecord=l;%总记录点数 tmin=0.2;%最小分析窗长为10个基波周期，即0.2s n0=0;nf=10*np-1; t=0;myerr=10; F=[];%存储每个滑窗的测量频率 Am=[];T=[];tt=[];ff=[];aa=[];coutwin=0;f=0;am=0;ph=0; while nf<timerecord flag=1; coutwin=coutwin+1; while myerr>eth & flag<=np & nf<timerecord [am,f,ph]=myprony(xtest(n0+1:nf+1),p,fs); xcal=zeros(1,nf-n0+1); in=find(f>0); for i=1:length(in); xcal=xcal+2*am(i)*sin(2*pi*f(i)/fs*[n0:nf]+pi/2-ph(i)); end myerr=err(xtest(n0+1:nf+1),xcal); flag=flag+1; nf=nf+1;%增加一个工频周期 end t=[n0:nf]'/fs; F=[F,f(in)];Am=[Am,am(in)];T=[T,{t}]; n0=nf+1;nf=n0+np-1;%滑动到下个窗 myerr=10; end for i=1:length(T) tt=[tt,T{i}']; ff=[ff,repmat(F(:,i),1,length(T{i}))]; aa=[aa,repmat(Am(:,i),1,length(T{i}))]; end plot(tt,ff);xlabel('t/s');ylabel('f/Hz'); fdim=[0:500];tdim=tt; [X,Y]=meshgrid(tdim,fdim); Z=zeros(length(tr),length(fdim)); for i=1:length(tdim) for j=1:length(fdim) if all(fix([fdim(j),fdim(j)]'-ff(:,i)))==0 for k=1:2 if abs(fdim(j)-ff(k,i))<0.1 Z(i,j)=aa(k,i); end end end end end % figure;surf(X,Y,Z'); pcolor(X,Y,Z');colormap jet, shading interp; contour(x,Y,Z',8,'k');