MATLAB中prony分析法的问题

923307293

现在我根据论文搞了一个prony的程序，但是只有输出的四个量，我想把导入的数据画成横轴为频率，纵轴为幅值的频谱图，应该怎么改

1. x=0:0.0001:1;   %原信号
   
2. y=V22(3:end);   %提取
   
3. % h=fopen('ws.txt','w');
   
4. 
   
5. p=6     %p为模型阶数
   
6. sf=1/0.0001    %sf为采样频率
   
7. [F,D,A,theta]=exprony_ma(y,p,sf)
   
8. %调用prony算法函数
   
9. function [F,D,A,theta]=exprony_ma(y,p,sf)
   
10. nm=2*p;                          %%% p为模型的阶数
    
11. n=fix(length(y)/2);
    
12. dt=1/sf;                        %%% 由采样频率计算时间间隔
    
13. h=y(1:2*n);
    
14. L=length(h);
    
15. M=L/2;
    
16. for k=1:nm
    
17.     x1(:,k)=h(k:M-1+k);
    
18. end
    
19. for k=1:M
    
20.     x2(k,:)=-h(nm+k);
    
21. end
    
22. B=x1\x2;           %% 最小二乘法求解prony 多项式系数
    
23. B(nm+1)=1;
    
24. B1=B(nm+1:-1:1);
    
25. V=roots(B1);        %% 求特征方程的根
    
26. F1=abs(log(V))/(2*pi*dt);  %% 初步计算频率
    
27. D1=log(abs(V))/(2*pi*dt);  %% 初步计算衰减因子
    
28. for k=0:(2*n-1)
    
29.     Va(k+1,:)=conj(V').^k;  %% 范得蒙矩阵 Va
    
30. end
    
31. S1=2*(inv(conj(Va')*Va)*conj(Va')*h);
    
32. [F2,I]=sort(F1);        %%%频率F1从小到达排序
    
33. m=0;
    
34. for k=1:nm-1
    
35.     if F2(k)~=F2(k+1)   %%去掉其中的共轭项
    
36.         continue;
    
37.     end
    
38.     m=m+1;
    
39.     ii=I(k);           %% 提取有效的元素序数
    
40.     F(m)=F1(ii);
    
41.     D(m)=D1(ii);
    
42.     A(m)=abs(S1(ii));
    
43.     theta(m)=angle(S1(ii))*180/pi;
    
44. end
    
45. end

复制代码

LJT1994

请问楼主，是不是Plot(x,y)，横轴是表示频率的变量，纵轴是幅值的变量，就可以了哇。。plot(A(m),F1)。我觉得是这样修改的，你也是这样绘制的嘛