prony算法问题

123

以下为prony算法matlab程序。
拟合实际曲线时用的算法如下图。
但是这不是相当于拿原先的数据复制粘贴吗？

%打开指定文件，并对信号进行Prony分析计算

1. 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)
   
2.     format long;
   
3.     load input3.txt;
   
4.     x_in=input3;%xlsread('shiyishi.xls');
   
5.     x = x_in;
   
6.     cpu=cputime;
   
7.          N=size(x,1);
   
8.         %取N/2的整数部分为初始的P
   
9.         P=floor(N/2);
   
10. 
    
11.         %去直流环节
    
12.         x_Sum = 0;
    
13.         for i=1:N
    
14.             x_Sum = x_Sum + x(i);
    
15.         end
    
16.         x_av = x_Sum / N;
    
17.     if x_av > 10E-10
    
18.             for i=1:N
    
19.                 x(i) = x(i) - x_av;
    
20.             end
    
21.     end
    
22.         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
23. 
    
24.         %滤波环节
    
25.         %fL=1;
    
26.     fL=2;
    
27.         if fL>1
    
28.             for i=fL+1:N
    
29.                 x(i-fL)=0;
    
30.                 for j=1:fL
    
31.                     x(i-fL) = x(i-fL)+(1/fL)*x(i-j+1);
    
32.                 end
    
33.             end
    
34.         end
    
35.         N=N-fL;
    
36.         tt=0:1:N-1;
    
37.         
    
38.         %P要求为偶数
    
39.         if mod(P, 2) ~= 0
    
40.            P = P - 1;
    
41.         end
    
42.         
    
43.         P1=P;
    
44.         if mod(P, 2) == 0
    
45.            % Generate R，生成样本矩阵
    
46.            for i=1:P1
    
47.                   for j=1:P1+1
    
48.                       ss=0;
    
49.                       for k=P1:N-1
    
50.                         %前向预测误差
    
51. %                         ss=ss+x(k+2-j,1)*x(k+1-i,1);
    
52.                         %后向预测误差
    
53.                         ss=ss+x(k+1-P1+i)*x(k+1-P1-1+j);
    
54.                         %同时考虑双向误差
    
55. %                         ss=ss+x(k+2-j,1)*x(k+1-i,1)+x(k+1-P1+i)*x(k+1-P1-1+j);              
    
56.                       end   
    
57.                       R(i,j)=ss;
    
58.                   end
    
59.            end
    
60. 
    
61. %               divide R into R1 and R2, and get A; R1*A=R2;
    
62.            for i=1:P
    
63.               for j=1:P
    
64.                  R1(i,j)=R(i,j+1);
    
65.               end
    
66.            end
    
67.            for i=1:P
    
68.                    R2(i)=R(i,1);
    
69.            end
    
70.         
    
71.            %添加的代码，使用SVD－TLS算法%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
72.            %首先进行SVD分解
    
73.            [u,s,v]=svd(R1);
    
74.            %根据奇异值确定实际的阶数M
    
75.            M=0;
    
76.            %计算全部奇异值的均方和
    
77.            Sum_SVD=0;
    
78.            for i=1:P
    
79.               Sum_SVD = Sum_SVD+s(i,i)^2;
    
80.            end
    
81.            cur_sum = 0;
    
82.            i=1;
    
83.            while 1 - sqrt(cur_sum/Sum_SVD) > 0.000001 & i<=P
    
84.               cur_sum = cur_sum + s(i,i)^2;
    
85.               M = M + 1;
    
86.               i = i + 1;
    
87.            end
    
88.            %输出预测的阶数M
    
89.            M;
    
90. 
    
91.            %生成Sp矩阵
    
92.            Sp=zeros(M+1, M+1);
    
93.            for j=1:M
    
94.               for i=1:(P-1)-M+1
    
95.                  Sp = Sp+s(j,j)^2*v(i:i+M,j)*v(i:i+M,j)';
    
96.               end
    
97.            end
    
98.         
    
99.            %根据公式生成最佳最小二乘解
    
100.        Sp_inv=inv(Sp);
     
101.        if isinf(Sp_inv(1,1)) == 1
     
102.            Sp_inv = pinv(Sp);
     
103.        end
     
104. 
     
105.            a_SVD = Sp_inv(1+1:M+1, 1)/Sp_inv(1,1);
     
106.            P = M;
     
107.            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     
108.            % Get Zi
     
109.            cof_SVD(1)=1;
     
110.            for i=1:M
     
111.               cof_SVD(i+1)=a_SVD(i);
     
112.            end
     
113.            Z=roots(cof_SVD);
     
114.            
     
115.        % Get y,y should be very close to x
     
116.            for i=1:P
     
117.               y(i)=x(i);
     
118.            end
     
119.         
     
120.            for i=P+1:N
     
121.               y(i)=0;
     
122.            end
     
123.            for i=P+1:N
     
124.               for j=1:P
     
125.                  y(i)=y(i)-a_SVD(j)*x(i-j);
     
126.               end
     
127.            end
     
128.             
     
129.            % Get B
     
130.            for i=1:N
     
131.               for j=1:P
     
132.                  Fy(i,j)=Z(j)^(i-1);
     
133.                end
     
134.            end
     
135.            Fz=Fy';
     
136.            FH=Fy'*Fy;
     
137.            Fn=inv(FH);
     
138.            B = inv(Fy'*Fy)*Fy'*y';
     
139.         
     
140.            % Calculate Amplitude, Phasor, Frequency and damp
     
141.        dt=0.01;
     
142.            for i=1:P
     
143.               Amp(i)=abs(B(i));
     
144.               Pha(i)=atan(imag(B(i))/real(B(i)));
     
145.               Fre(i)=atan(imag(Z(i))/real(Z(i)))/(2*pi*dt);
     
146.               damp(i)=log(abs(Z(i)))/dt;
     
147.            end
     
148. 
     
149.        %调试用的代码
     
150.        if isnan(Amp(1)) == 1
     
151.            error('幅值的求解出现错误！');
     
152.        end
     
153. 
     
154.            % Get xc,verify if xc is nearly equal to x
     
155.            for i=1:P
     
156.                  if(real(B(i))>=0.0)
     
157.                      bth(i)=Pha(i);
     
158.                  else
     
159.                      bth(i)=pi+Pha(i);
     
160.                  end
     
161.            end
     
162.        %对Pha重新幅值
     
163.        for i=1:P
     
164.            if Pha(i) < 0
     
165.                Pha(i) = Pha(i) + 2*pi;
     
166.            end
     
167.        end
     
168. 
     
169.            for i=1:N
     
170.                  xc(i)=0.0;
     
171.                  for j=1:P
     
172.                      xc(i)=xc(i)+Amp(j)*exp(damp(j)*(i-1)*dt)*cos(2*pi*Fre(j)*(i-1)*dt+bth(j));
     
173.                  end
     
174.            end
     
175. 
     
176. %       if showfigure == 1
     
177. %              %显示特征根的位置
     
178. %              figure;
     
179. %              plot(damp, 2*pi*Fre, 'r*');
     
180. %       end
     
181. 
     
182.            xj=xc';
     
183.            er=0;
     
184.            all = 0;
     
185.            for i=1:N         
     
186.                   er=er+(x(i)-xc(i))^2;
     
187.                   all = all+x(i)^2;
     
188.            end
     
189.            SNR=-20*log(sqrt(er/all));
     
190.            
     
191.            % Calculate Prony spectrum
     
192.            %频谱的取值范围为0－5Hz
     
193.            f=0:0.01:4.99;
     
194.            for j=1:size(f,2)
     
195.                sptr(j)=0;
     
196.                sptr1(j)=0;
     
197.                sptr2(j)=0;
     
198.                angl(j)=0;
     
199.                tmpangle = 0;
     
200.                for k=1:P
     
201.                   sptr1(j)=sptr1(j)+Amp(k)*cos(bth(k))*(2*damp(k)/(damp(k)^2+(2*pi*(f(j)-Fre(k)))^2));
     
202.                   sptr2(j)=sptr2(j)+Amp(k)*sin(bth(k))*(2*damp(k)/(damp(k)^2+(2*pi*(f(j)-Fre(k)))^2));
     
203.                end
     
204.                sptr(j)=sqrt(sptr1(j)^2+sptr2(j)^2);
     
205.                tmpangle = atan2(sptr1(j),sptr2(j));
     
206.                angl(j) = tmpangle*360/pi;
     
207.            end
     
208.        %对幅频响应进行归一化，并且寻找主导频率模式
     
209.        %寻找sptr的最大值
     
210.        max_sptr=0;
     
211.        main_f=0;
     
212.        for j=1:size(f,2)
     
213.            if sptr(j) > max_sptr
     
214.                max_sptr = sptr(j);
     
215.                %主导频率出现在最大幅频响应位置
     
216.                if f(j) ~= 0;
     
217.                    main_f = f(j);
     
218.                else
     
219.                    max_sptr = 0;
     
220.                end
     
221.            end
     
222.        end
     
223.        main_f;
     
224.        %找出真正计算的频率值
     
225.        f_err = 10;
     
226.        f_main = 0;
     
227.        for i=1:size(Amp, 2)
     
228.            if abs(main_f - Fre(i)) < f_err
     
229.                f_main = Fre(i);
     
230.                damp_main = damp(i);
     
231.                f_err = abs(main_f - Fre(i));
     
232.            end
     
233.        end
     
234.        main_f = f_main;
     
235.        main_damp = damp_main;
     
236.        %进行归一化
     
237.        for j=1:size(f,2)
     
238.            sptr(j)=sptr(j)/max_sptr;
     
239.            Amp_Response(j) = sptr(j);
     
240.        end
     
241. 
     
242.            for i=1:size(f,2)
     
243.                if angl(i) > 0
     
244.                    angl(i)=angl(i)-360;
     
245.                end
     
246.            end
     
247.            showfigure = 1;
     
248.        if showfigure == 1
     
249.               %在一个图上显示时域拟和曲线和Prony频谱分布
     
250.               figure;
     
251.               subplot(2,1,1);
     
252.               %plot(tt,x(1:N),'b', tt,y, 'r', tt,xc, '*b');
     
253. 
     
254.           plot(tt,x(1:N),'r', tt,xc, '*b');
     
255.           subplot(2,1,2);
     
256.               plot(f,sptr);
     
257.               %subplot(2,2,3);
     
258.               %plot(f,sptr);
     
259.               %subplot(2,2,4);
     
260.               %plot(f,angl);
     
261.       end
     
262. 
     
263.           cpu=cputime-cpu;
     
264.       format short;
     
265.         end
     
266. 
     

复制代码

特斯拉之魂

老哥，可以交流下吗，我qq257589181。感激不尽

123

有什么问题可以在这里说，我懂得就会讲啦