基于prony算法的参数辨识算法的仿真的详细版欢迎订阅本博：

https://blog.csdn.net/ccsss22/article/details/115358232

1.问题描述：

 建立如下被测信号：



被测信号中包含四个振荡模态，在数据窗宽度同样为10s的前提下，利用不同的采样频率做普罗尼计算。 

2.部分程序：

 

function X = func_Prony(Signal,dt);

s     = Signal; 

L     = length(s(1:length(Signal))); 

Order = ceil(L/2); 

R     = []; 

K1    = 0;   

K2    = 0; 

%扩展矩阵 

while K1 <= Order 

    K2  = 1; 

    Re = []; 

    while K2 <= Order 

        u  = Order - K2 +1; 

        v  = L - K2; 

        m  = Order - K1 +1; 

        l  = L - K1; 

        r  = sum(s(u:v).*conj(s(m:l))); 

        Re = [Re,r]; 

        K2  = K2 + 1; 

    end 

    R = [R,Re']; 

    K1 = K1 + 1; 

end 

%计算阶数 

Order = func_Order(R); 

%计算相关参数

K2 = Order-1; 

Re = R(2:end,2:K2+1); 

b  = R(:,1); 

b  = b(2:end); 

a  = pinv(Re)*(-b); 

R1 = R(1,:); 

R1 = R1(1:K2+1); 

a1 = [1 a']; 

Ep = sum(R1.*a1); 

P  = [1 a']; 

z  = roots(P); 

%估计序列X 

ks    = 1:K2; 

X(ks) = s(ks); 

for Order = K2+1 : L 

    Lij      = 1:K2; 

    X(Order) = sum(-a'.*s(Order-Lij)); 

    Order    = Order+1; 

end 

Zh = []; 

for mm = 0:L-1 

    Zh = [Zh,z.^mm]; 

end 

Z    = Zh'; 

Z    = conj(Z); 

Zhh  = Z'; 

b    = (Zhh*Z)^(-1)*Zhh*X'; 

%最后得到的四个参数值

A     = abs(b)

f     = angle(z)/2/pi/0.001

a     = log(abs(z))/dt

theta = angle(b)/2/pi/dt



 

3.仿真结论：

       注意，这里论文中你所给的那个公式，貌似有点小错误，这里我们使用了两组公式进行计算，一组是你所提供的公式，一组是我们给的测试数据。

       仿真结果如下所示：



A-27-6

 

        这里的思想主要是广义旁瓣相消器是LCMV的一种等效的结构。其主要涉及到的公式有：



这个降秩矩阵S的计算公式，和原论文相同。



首先介绍一下，为什么性能会提高，

这个，有的时候也说不准，这么做一定会提示性能，我们多半都是在理论上说得通的前提下，通过仿真去分析性能，然后来完善理论，这里我弄完这个部分之后，提升性能的原因有二：

一：一个是采用GSC算法进行降序之后，获得的降秩性能更优，具体可看GSC的相关理论；

二：第二是，这里采用Wq补偿的方式，在实际中，由于直接使用W降秩矩阵的时候，由于迭代得到的值，肯定会存在迭代步进的设置而导致最后迭代误差的影响，那么这里，通过Wq进行补偿，会在一定程度上弥补这个缺陷。

        我们通过仿真可知，最后的误码率性能提升结果如下所示：



 

对于同样的误码率值，大概提升了0.3db，

对于同样的误码率值，大概提升了0.3db，

另外当SNR<0的时候，两个性能比较接近。但是放大之后仔细看，改进后的效果还是有略微的提升。

1.问题描述：

 建立如下被测信号：



被测信号中包含四个振荡模态，在数据窗宽度同样为10s的前提下，利用不同的采样频率做普罗尼计算。根据公式的基本表达式：



所以最后，进行估计得到的参数为：



    我们的程序设计中，按照论文中的说明进行设计，对于非连续时间的函数方程式：



普罗尼算法步骤：





3）特征解求取：



由式递推求解









 

 

2.部分程序：

 

function X = func_Prony(Signal,dt);

s     = Signal; 

L     = length(s(1:length(Signal))); 

Order = ceil(L/2); 

R     = []; 

K1    = 0;   

K2    = 0; 

%扩展矩阵 

while K1 <= Order 

    K2  = 1; 

    Re = []; 

    while K2 <= Order 

        u  = Order - K2 +1; 

        v  = L - K2; 

        m  = Order - K1 +1; 

        l  = L - K1; 

        r  = sum(s(u:v).*conj(s(m:l))); 

        Re = [Re,r]; 

        K2  = K2 + 1; 

    end 

    R = [R,Re']; 

    K1 = K1 + 1; 

end 

%计算阶数 

Order = func_Order(R); 

%计算相关参数

K2 = Order-1; 

Re = R(2:end,2:K2+1); 

b  = R(:,1); 

b  = b(2:end); 

a  = pinv(Re)*(-b); 

R1 = R(1,:); 

R1 = R1(1:K2+1); 

a1 = [1 a']; 

Ep = sum(R1.*a1); 

P  = [1 a']; 

z  = roots(P); 

%估计序列X 

ks    = 1:K2; 

X(ks) = s(ks); 

for Order = K2+1 : L 

    Lij      = 1:K2; 

    X(Order) = sum(-a'.*s(Order-Lij)); 

    Order    = Order+1; 

end 

Zh = []; 

for mm = 0:L-1 

    Zh = [Zh,z.^mm]; 

end 

Z    = Zh'; 

Z    = conj(Z); 

Zhh  = Z'; 

b    = (Zhh*Z)^(-1)*Zhh*X'; 

%最后得到的四个参数值

A     = abs(b)

f     = angle(z)/2/pi/0.001

a     = log(abs(z))/dt

theta = angle(b)/2/pi/dt



 

3.仿真结论：

       注意，这里论文中你所给的那个公式，貌似有点小错误，这里我们使用了两组公式进行计算，一组是你所提供的公式，一组是我们给的测试数据。

       仿真结果如下所示：



 

注意，这里我们还用了自己的测试数据进行测试：



 

参数估计结果如下所示：


			

			
			

			
4.9379
			
			

			
4.9379
			
			

			
2.8537
			
			

			
2.8537
			
		

			

			
			

			
-15.0000
			
			

			
15.0000
			
			

			
-120.0000
			
			

			
120.0000
			
		

			

			
			

			
-0.0125
			
			

			
-0.0125
			
			

			
-0.0500
			
			

			
-0.0500
			
		

			

			
			

			
2.1250
			
			

			
-2.1250
			
			

			
0.7500
			
			

			
-0.7500
			
		
 

A-27-6