关于prony算法
最近一直在搞prony算法，资料是张贤达的《现代信号处理》和束红春的《电力工程信号处理应用》，这两本书都有关于prony扩展算法的内容，甚至有相关步骤，我也根据这些内容自己编了程序。为了测试 在信号x=160*sin(2*pi*f1*t+pi/5)+150*exp(-3*t).*sin(2*pi*f2*t+pi/4);取时窗10ms 100个数据点，但是分离不出这两信号，后面尽管修改程序自己设定介数p ，增加数据窗长度还是失败，而且有数据稳定性也有问题。
但是用mathwork网站上下的pronytool工具箱还是能比较容易的分离出来的，于是我查看pronytool的原程序，发现工具包里边的算法和上边两本书上的算法是不一样的，好像是根据零极点和滤波脉冲响应什么的，我实在不是太懂。
所以在这里诚心向各位请教
%% 数据准备
clear;
clc;
format long
% load('1000kV示范工程线路','t','vX0043a')
% x = vX0043a(201:500);
% t = t(201:500);
f1 = 49;
f2 = 51;
t=0.0001:0.0001:0.01;
x = 160*sin(2*pi*f1*t+pi/5)+150*exp(-3*t).*sin(2*pi*f2*t+pi/4);
dt = 0.0001;
N = length(x);
pe = floor(N/2);
%% 构造样本矩阵
Re=zeros(pe+1,pe+1);
for i = 2:pe+1
for j = 1:pe+1
for n = pe:N-1
Re(i,j) = Re(i,j)+x(n-j+2)*x(n-i+2);
end
end
end
Re(1,:) = [];
%% SVD_TLS确定介数p及a
= svd(Re); %%%%%%奇异值分解
% 求p值
% 计算全部奇异值平方和
sum_all = 0;
for i = 1:pe
sum_all = sum_all+S(i,i)^2;
end
% 归一化比值Ak/A 求p值
sum_k = 0;
k = 1;
while 1 - sqrt(sum_k/sum_all) > 0.0000000000000000001 & i<=pe
sum_k = sum_k+S(k,k)^2; %%%%%%%计算k个奇异值平方和
k = k+1;
end
p = k-1;
% 求Sp部分
Sp=zeros(p+1,p+1); %%%%%%s生成（p+1）X（p+1）维矩阵Sp
for j = 1:p
for i = 1:(pe+1-p)
Sp = Sp+S(j)^2*V(i:i+p,j)*V(i:i+p,j)';
end
end
% SS = zeros(pe,pe+1);
% for i = 1:p
% SS(i,i) = S(i,i);
% end
% B = U*SS*V;
% 求Sp逆矩阵
inv_Sp=inv(Sp);
if isinf(inv_Sp(1,1)) == 1
inv_Sp = pinv(Sp);
end
% 求a
a=inv_Sp(2:p+1,1)/inv_Sp(1,1);
%% 求z
y=;
z=roots(y);
%% 求x的近似值x_j
%求前p近似值等于测量值 x_j(1:p)
x_j=zeros(N,1);
for i = 1:p
x_j(i)=x(i);
end
%求x的N-p+1个近似值 x_j(p+1:N)
for n = p+1:N
for i = 1:p
x_j(n)=x_j(n)-a(i)*x_j(n-i);
end
end
%% 画图 x、x_j
hold on;
plot(t,x,'k');
plot(t,x_j,'r');
hold off;
%% 求取 b=inv(H)*Z'*x_j
% 求取N X p维vandermode矩阵Z
Z=zeros(N,p);
for i=1:N
Z(i,:)=z'.^(i-1);
end
%求取H
H=zeros(p,p);
for i=1:p
for j=1:p
m=(conj(z(i))*z(j));
H(i,j)=(m^N-1)/(m-1);
end
end
% 求取b
b=inv(H)*Z'*x';
%% 计算振幅Amp 频率Fre 衰减因子Damp 相位 Pha
for i = 1:p
Amp(i) = abs(b(i));
Fre(i) = atan(imag(z(i)/real(z(i))))/(2*pi*dt);
Damp(i) = log(abs(z(i)))*dt;
Pha(i) = atan(imag(b(i)/real(b(i))));
end 楼主加油，我们都看好你哦。 这个。。。太复杂。额额。不懂哟 有人做过这个的话帮个忙啊 好专业，看不懂帮顶 我的学长似乎在做这个啊！他可能会，帮你请教一下 我的学长似乎在做这个啊！他可能会，帮你请教一下 回复 6# sjx0323
真的吗 那真是太谢谢了 可以联系我吗？qq345749437 回复 7# xian2006
他今天不在，明天来了帮你请教一下！ 回复 7# xian2006
看来这个问题解决不了，我似乎发现你貌似就是我的学长啊！呵呵 我帮学长顶顶
页:
[1]
2
