24小时热门版块排行榜

返回列表

172547204

铁虫 (初入文坛)

应助: 0 (幼儿园)
金币: 10
散金: 10
帖子: 7
在线: 4.5小时
虫号: 2016540
注册: 2012-09-20
性别: MM

[求助] 算法错误，求原因和解决办法

Sample Text：??? Attempted to access Amp(2); index out of bounds because numel(Amp)=1.

Error in ==> x_svdprony at 263
      xn(i)=xn(i)+Amp(2)*exp(damp(2)*(i-1)*dt)*cos(2*pi*Fre(2)*(i-1)*dt+bth(2));

Sample Text
%打开指定文件，并对信号进行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=[32.32
29.29
27.56
25.95
24.44
23.02
21.68
20.41
19.23
18.11
17.05
16.06
15.12
14.24
13.41
12.63
11.89
11.2
10.55
9.93
9.36
8.81
8.3
7.81
7.36
6.93
6.53
6.15
5.79
5.45
5.13
4.83
4.55
4.29
4.04
3.8
3.58
3.37
3.17
2.99
2.81
2.65
2.5
2.35
2.21
2.08
1.96
1.85
1.74
1.64
1.54
1.45
1.37
1.29
1.21
1.14
1.07
1.01
0.95
0.9
0.84
0.79
0.75

];
x = x_in;
cpu=cputime;
N=size(x,1);%返回的是数组的行数也就是有多少个输入数据
%取N/2的整数部分为初始的P
P=floor(N/2);%四舍五入取整数

%去直流环节
x_Sum = 0;
for i=N-5: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

1
for j=1

1+1
               ss=0;
               for k=P1:N-1
                  %前向预测误差
                  ss=ss+x(k+2-j,1)*x(k+1-i,1);
                  %后向预测误差
                  %ss=ss+x(k+1-P1+i)*x(k+1-P1-1+j);
                  %同时考虑双向误差
                  %ss=ss+x(k+2-j,1)*x(k+1-i,1)+x(k+1-P1+i)*x(k+1-P1-1+j);
               end
               R(i,j)=ss;
            end
   end

      % divide R into R1 and R2, and get A; R1*A=R2;
   for i=1

for j=1

         R1(i,j)=R(i,j+1);
      end
   end
   for i=1

           R2(i)=R(i,1);
   end

   %添加的代码，使用SVD－TLS算法%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %首先进行SVD分解
   [u,s,v]=svd(R1);
   %根据奇异值确定实际的阶数M
   M=0;
   %计算全部奇异值的均方和
   Sum_SVD=0;
   for i=1

      Sum_SVD = Sum_SVD+s(i,i)^2;
   end
   cur_sum = 0;
   i=1;
   while 1 - sqrt(cur_sum/Sum_SVD) > 0.000001 & i<=P
      cur_sum = cur_sum + s(i,i)^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(j,j)^2*v(i:i+M,j)*v(i:i+M,j)';
      end
   end

   %根据公式生成最佳最小二乘解
   Sp_inv=inv(Sp);
   if isinf(Sp_inv(1,1)) == 1
         Sp_inv = pinv(Sp);
   end

   a_SVD = Sp_inv(1+1:M+1, 1)/Sp_inv(1,1);
   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 y,y should be very close to x
   for i=1

      y(i)=x(i);
   end

   for i=P+1:N
      y(i)=0;
   end
   for i=P+1:N
      for j=1

         y(i)=y(i)-a_SVD(j)*x(i-j);
      end
   end

   % Get B
   for i=1:N
      for j=1

         Fy(i,j)=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

      Amp(i)=abs(B(i));
      Pha(i)=atan(imag(B(i))/real(B(i)));
      Fre(i)=atan(imag(Z(i))/real(Z(i)))/(2*pi*dt);
      damp(i)=log(abs(Z(i)))/dt;
   end

   %调试用的代码
   if isnan(Amp(1)) == 1
         error('幅值的求解出现错误！');
   end

   % Get xc,verify if xc is nearly equal to x
   for i=1

         if(real(B(i))>=0.0)
               bth(i)=Pha(i);
         else
               bth(i)=pi+Pha(i);
         end
   end
   %对Pha重新幅值
   for i=1

         if Pha(i) < 0
            Pha(i) = Pha(i) + 2*pi;
         end
   end

   for i=1:N
         xc(i)=0.0;
         for j=1

               xc(i)=xc(i)+Amp(j)*exp(damp(j)*(i-1)*dt)*cos(2*pi*Fre(j)*(i-1)*dt+bth(j));
         end
   end
%.........。。。自己编的代码。。。。。。。。
  xm(i)=0.0;
  for i=1:N
      xm(i)=xm(i)+Amp(1)*exp(damp(1)*(i-1)*dt)*cos(2*pi*Fre(1)*(i-1)*dt+bth(1));
  end
xn(i)=0.0;
  for i=1:N
      xn(i)=xn(i)+Amp(2)*exp(damp(2)*(i-1)*dt)*cos(2*pi*Fre(2)*(i-1)*dt+bth(2));
end
% xk(i)=0.0;
%for i=1:N
%    xk(i)=xk(i)+Amp(3)*exp(damp(3)*(i-1)*dt)*cos(2*pi*Fre(3)*(i-1)*dt+bth(3));
%end
   figure;
  % subplot(2,1,1);
      % plot(tt,x(1:N),'b', tt,y, 'r', tt,xc, '*b');
      %plot(tt,xm,'r', tt,xn, 'b',tt,xc,'g',tt,x(1:N), '*b',tt,xk,'m');
      plot(tt,xm,'r', tt,xn, 'b',tt,xc,'g',tt,x(1:N), '*b');%二阶画图命令
      %plot(tt,xc,'r',tt,x(1:N), '*b');
      ylabel('电压（V）');
      legend('拟合的第一条指数分量','拟合的第二条指数分量','拟合曲线','实际测量值');
      % legend('拟合曲线','实际测量值');
%。。。。。。。。。。自己编的代码。。。。。。。。
%    if showfigure == 1
%    %显示特征根的位置
%       figure;
%       plot(damp, 2*pi*Fre, 'r*');
%    end

   xj=xc';
   er=0;
   all = 0;
   for i=1:N
            er=er+(x(i)-xc(i))^2;
            all = all+x(i)^2;
   end
   SNR=-20*log(sqrt(er/all));

   % Calculate Prony spectrum
   %频谱的取值范围为0－5Hz
   f=0:0.01:4.99;
   for j=1:size(f,2)
         sptr(j)=0;
         sptr1(j)=0;
         sptr2(j)=0;
         angl(j)=0;
         tmpangle = 0;
         for k=1

            sptr1(j)=sptr1(j)+Amp(k)*cos(bth(k))*(2*damp(k)/(damp(k)^2+(2*pi*(f(j)-Fre(k)))^2));
            sptr2(j)=sptr2(j)+Amp(k)*sin(bth(k))*(2*damp(k)/(damp(k)^2+(2*pi*(f(j)-Fre(k)))^2));
         end
         sptr(j)=sqrt(sptr1(j)^2+sptr2(j)^2);
         tmpangle = atan2(sptr1(j),sptr2(j));
         angl(j) = tmpangle*360/pi;
   end
   %对幅频响应进行归一化，并且寻找主导频率模式
   %寻找sptr的最大值
   max_sptr=0;
   main_f=0;
   for j=1:size(f,2)
         if sptr(j) > max_sptr
            max_sptr = sptr(j);
            %主导频率出现在最大幅频响应位置
            if f(j) ~= 0;
               main_f = f(j);
            else
               max_sptr = 0;
            end
         end
   end
   main_f;
   %找出真正计算的频率值
   f_err = 10;
   f_main = 0;
   for i=1:size(Amp, 2)
         if abs(main_f - Fre(i)) < f_err
            f_main = Fre(i);
            damp_main = damp(i);
            f_err = abs(main_f - Fre(i));
         end
   end
   main_f = f_main;
   main_damp = damp_main;
   %进行归一化
   for j=1:size(f,2)
         sptr(j)=sptr(j)/max_sptr;
         Amp_Response(j) = sptr(j);
   end

   for i=1:size(f,2)
         if angl(i) > 0
            angl(i)=angl(i)-360;
         end
   end
   showfigure = 1;
      if showfigure == 1
%    %显示特征根的位置
%       figure;
%       plot(damp, 2*pi*Fre, 'r*');
%    end
   % if showfigure == 1
      %在一个图上显示时域拟和曲线和Prony频谱分布
      %  figure;
      % subplot(2,1,1);
      %plot(tt,x(1:N),'b', tt,y, 'r', tt,xc, '*b');
      % plot(tt,x(1:N),'r', tt,xc, '*b');
%修改的       subplot(2,1,2);
%修改的       plot(f,sptr);
      %subplot(2,2,3);
      %plot(f,sptr);
      %subplot(2,2,4);
      % plot(f,angl);
   end

   cpu=cputime-cpu;
   format short;
end

end

回复此楼