PRONY算法最后日志

**

PRONY算法最后日志

**

• 1.优点

（1）采样时间较短（相对于FFT），即可求得出信号中的低频含量；
（2）计算出的频率之间的间隔（分辨率）与采样时间无关，分辨率很高；
（3）计算结果中含有衰减系数，可视为对动态信号有一定的对应性。

• 2.缺点

（1）模型矩阵阶数决定可检测的频率个数，即矩阵阶数-1≥频率个数。信号含有频率必须有限次，微小噪声信号会对计算结果产生影响（对频率低于噪声信号频率的计算影响较大），矩阵阶数必须足够大。
（2）模型矩阵阶数决定计算矩阵，即计算矩阵＝4倍的模型矩阵阶数。计算中间占用数据空间大，接近指数增加。
（3）算法中包含广义逆矩阵运算、复数矩阵幂运算（4倍模型矩阵阶数幂）、复数矩阵逆矩阵运算、高阶方程求根运算。计算量大，计算时间长。
（4）除广义逆矩阵运算外可采用单精度模式而不影响最后计算结果（可节省90%的时间），广义逆矩阵运算中含有大量乘法，同时必须采用双精度模式，无法提高运算效率。

• 3.结论：

（1）prony算法不适合底层硬件运算，仅仅适合上位机（64位）中运用；
（2）Prony算法仅仅适合频率个数有限次的频谱分析。

PRONY算法（matlab可运行.m文件）

M=16; % 谐波个数
fs=256; % 采样频率
t=0.25;
N=fs*t; % 数据长度
n=0:N-1; % 数据索引
f=[50,6,11,17,30,41,48,1,61,71,75,80,88,105];
A=[70,10,10,10,10,10,10,10,10,10,10,10,10,10];
Q=[pi/3,pi/3,pi/3,pi/3,pi/3,pi/3,pi/3,pi/3,pi/3,pi/3,pi/3,pi/3,pi/3,pi/3];
k=[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0];
s=zeros(1,N); % 初始化
for i=1:M % 产生谐波信号
s=s+exp(k(i))*A(i)*cos(2*pi*f(i)*n./fs+Q(i));
end
p=16; % p 为模型阶数
nm=2*p;
n=fix(length(s)/2);
dt=1/fs;
h=s(1,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=pinv(x1)*x2;
B(nm+1)=1;
B1=B(nm+1:-1:1);
V=roots(B1);
F1=abs(log(V))/(2*pi*dt);
for k=0:(2*n-1)
    Va(k+1,:)=conj(V').^k;
end
S1=2*(inv(conj(Va')*Va)*conj(Va')*h);
[F2,I]=sort(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);
    A1(m)=abs(S1(ii));
    theta(m)=angle(S1(ii))*180/pi;
end
subplot(3,1,1); plot(s);
xlabel('点数');ylabel('幅值');grid on;
subplot(3,1,2); stem(F,A1);
xlabel('频率(Hz)');ylabel('幅值');grid on;   
subplot(3,1,3); stem(F,theta);
xlabel('频率(Hz)');ylabel('相位');grid on; 

如有问题的话，可关注订阅号，会一一解答问题。