923307293 发表于 2020-10-27 17:19:10

MATLAB中prony分析法的问题

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

x=0:0.0001:1;   %原信号
y=V22(3:end);   %提取
% h=fopen('ws.txt','w');

p=6   %p为模型阶数
sf=1/0.0001    %sf为采样频率
=exprony_ma(y,p,sf)
%调用prony算法函数
function =exprony_ma(y,p,sf)
nm=2*p;                        %%% p为模型的阶数
n=fix(length(y)/2);
dt=1/sf;                        %%% 由采样频率计算时间间隔
h=y(1:2*n);
L=length(h);
M=L/2;
for k=1:nm
    x1(:,k)=h(k:M-1+k);
end
for k=1:M
    x2(k,:)=-h(nm+k);
end
B=x1\x2;         %% 最小二乘法求解prony 多项式系数
B(nm+1)=1;
B1=B(nm+1:-1:1);
V=roots(B1);      %% 求特征方程的根
F1=abs(log(V))/(2*pi*dt);%% 初步计算频率
D1=log(abs(V))/(2*pi*dt);%% 初步计算衰减因子
for k=0:(2*n-1)
    Va(k+1,:)=conj(V').^k;%% 范得蒙矩阵 Va
end
S1=2*(inv(conj(Va')*Va)*conj(Va')*h);
=sort(F1);      %%%频率F1从小到达排序
m=0;
for k=1:nm-1
    if F2(k)~=F2(k+1)   %%去掉其中的共轭项
      continue;
    end
    m=m+1;
    ii=I(k);         %% 提取有效的元素序数
    F(m)=F1(ii);
    D(m)=D1(ii);
    A(m)=abs(S1(ii));
    theta(m)=angle(S1(ii))*180/pi;
end
end

LJT1994 发表于 2020-12-15 20:10:14

请问楼主,是不是Plot(x,y),横轴是表示频率的变量,纵轴是幅值的变量,就可以了哇。。plot(A(m),F1)。我觉得是这样修改的,你也是这样绘制的嘛
页: [1]
查看完整版本: MATLAB中prony分析法的问题