prony 算法

%打开指定文件，并对信号进行Pon分析计算
function [M Amp Fre damp Pha main_f main_damp x xc y Amp_Response er all N dt]=x_svdprony（x_in dt fL showfigure）
    format long;
    x_in=[0.1999 0.1953 0.1898 0.1839 0.1768 0.1727 0.1684 0.1665 0.1664 0.1664 0.1663];
    x = x_in‘
    cpu=cputime;
	N=size（x1）;
	%取N/2的整数部分为初始的P
	P=floor（N/2）;
    
	%去直流环节
	x_Sum = 0;
	for i=1:N
	    x_Sum = x_Sum + x（i）;
	end 
	x_av = x_Sum / N;
    if x_av > 10E-10
    	for i=1:N
	        x（i） = x（i） - x_av;
    	end
    end
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

	%滤波环节
	%fL=1;
    fL=2;
	if fL>1
	    for i=fL+1:N
	        x（i-fL）=0;
	        for j=1:fL
	            x（i-fL） = x（i-fL）+（1/fL）*x（i-j+1）;
	        end
	    end
	end
	N=N-fL;
	tt=0:1:N-1;
	
	%P要求为偶数
	if mod（P 2） ~= 0
	   P = P - 1;
	end
	
	P1=P;
	if mod（P 2） == 0
	   % Generate R，生成样本矩阵
	   for i=1:P1
	          for j=1:P1+1
	              ss=0;
	              for k=P1:N-1
	                %前向预测误差
	                ss=ss+x（k+2-j1）*x（k+1-i1）;
	                %后向预测误差
	                %ss=ss+x（k+1-P1+i）*x（k+1-P1-1+j）;
	                %同时考虑双向误差
	                %ss=ss+x（k+2-j1）*x（k+1-i1）+x（k+1-P1+i）*x（k+1-P1-1+j）;              
	              end   
	              R（ij）=ss;
	          end
	   end

   	   % divide R into R1 and R2 and get A; R1*A=R2;
	   for i=1:P
	      for j=1:P
	         R1（ij）=R（ij+1）;
	      end
	   end
	   for i=1:P
     	      R2（i）=R（i1）;
	   end
	
	   %添加的代码，使用SVD－TLS算法%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	   %首先进行SVD分解
	   [usv]=svd（R1）;
	   %根据奇异值确定实际的阶数M
	   M=0;
	   %计算全部奇异值的均方和
	   Sum_SVD=0;
	   for i=1:P
	      Sum_SVD = Sum_SVD+s（ii）^2;
	   end
	   cur_sum = 0;
	   i=1;
	   while 1 - sqrt（cur_sum/Sum_SVD） > 0.000001 & i<=P
	      cur_sum = cur_sum + s（ii）^2;
	      M = M + 1;
	      i = i + 1;
	   end
	   %输出预测的阶数M
	   M;

	   %生成Sp矩阵
	   Sp=zeros（M+1 M+1）;
	   for j=1:M 
	      for i=1:（P-1）-M+1
	         Sp = Sp+s（jj）^2*v（i:i+Mj）*v（i:i+Mj）‘;
	      end
	   end
	
	   %根据公式生成最佳最小二乘解
       Sp_inv=inv（Sp）;
       if isinf（Sp_inv（11）） == 1
           Sp_inv = pinv（Sp）;
       end
       
	   a_SVD = Sp_inv（1+1:M+1 1）/Sp_inv（11）;
	   P = M;
	   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	   % Get Zi
	   cof_SVD（1）=1;
	   for i=1:M
	      cof_SVD（i+1）=a_SVD（i）;
	   end
	   Z=roots（cof_SVD）;
	   
       % Get yy should be very close to x
	   for i=1:P
	      y（i）=x（i）;
	   end
	
	   for i=P+1:N
	      y（i）=0;
	   end
	   for i=P+1:N
	      for j=1:P
	         y（i）=y（i）-a_SVD（j）*x（i-j）;
	      end
	   end
	    
	   % Get B
	   for i=1:N
	      for j=1:P
	         Fy（ij）=Z（j）^（i-1）;
	       end
	   end
	   Fz=Fy‘;
	   FH=Fy‘*Fy;
	   Fn=inv（FH）;
	   B = inv（Fy‘*Fy）*Fy‘*y‘;
	
	   % Calculate Amplitude Phasor Frequency and damp
       dt=0.01;
	   for i=1:P
	      Amp（i）=abs（B