| 查看: 813 | 回复: 0 | |||
[求助]
算法错误，求原因和解决办法
|
|
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  1for j=1 1+1ss=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+1Sp = 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 |
回复此楼» 猜你喜欢
2026国自然函评费到账
已经有17人回复
-
让我中一个面上吧！
已经有15人回复
-
2026年的国家社科基金项目通讯评审的新规则与新动向、新挑战
已经有5人回复
-
跳槽后在研项目怎么办？
已经有10人回复
-
什么时候开奖？
已经有10人回复
-
今天放榜吗？
已经有16人回复
-
93BebMhtakh前后11位开头都是大写
已经有8人回复
-
范进中举一文的中心思想
已经有4人回复
-
只有每年这种时候来逛逛小木虫
已经有26人回复
-
filecode，4个jtjc了
已经有16人回复
找到一些相关的精华帖子，希望有用哦~
关于晶体结构解析中的A、B类错误的解决办法
已经有4人回复
- LARS算法问题求助
已经有6人回复
- 关于迭代算法的收敛问题，求高人指点！
已经有3人回复
- 急！！！新手求助phylip 软件具体操作以及错误解决方法
已经有4人回复
- vasp.5.2编译出错，跪求解决办法！
已经有10人回复
- bruker topspin 难题解决方法
已经有23人回复
- 发现生活、工作中的问题，并提出使用人工智能方法去解决问题的方案
已经有13人回复
- 2个A类错误，求解决办法
已经有19人回复
- 苯基硼酸酯过柱子消失的问题，求解决方法
已经有9人回复
- 【求助】gaussian计算错误L301，请高手赐教
已经有14人回复
- 【交流】常见A、B类错误及解决办法大收集
已经有61人回复
- 【求助】DLS中数据算法问题
已经有16人回复
- 【求助】遗传算法求解eil51.tsp问题，离最优解还很远。怎么办？
已经有26人回复
- 【交流】icp测试过程中碰到的问题及解决方法
已经有5人回复
- 【求助】求助B 类错误的解决方法
已经有5人回复
- 【资源】液相色谱常见问题及处理方法
已经有43人回复
- 【原创】GAUSSION计算常见错误及解决方案
已经有32人回复
点击这里搜索更多相关资源科研从小木虫开始，人人为我，我为人人
40