prony谱线估计算法已修正

clear all;
close all;
clc;
x1=0.05:0.05:10;
dt=0.05;%抽样时间
fs=1/dt;%抽样频率

wn=randn（1200）;
% x=（10*x1）.*cos（2*pi*5*x1+pi/2）+（25*x1）.*cos（2*pi*2.5*x1+pi/2）+wn;
x=（5*x1）.*cos（2*pi*4*x1+pi/2）+（10*x1）.*cos（2*pi*4.5*x1+pi/2）+wn;

N=length（x）-20;%N为了满足后面r（ij）的计算，所以减少N的取值
figure（1）;
plot（x）;
grid on;
pe=6;%生设置阶数为6，记住pe和p不同，p是要算的，pe是先给出的
%书上公式有问题，加法x维度超出，减法维度刚好不超出，所以只能将加法的维度拉回，减法的维度不管
for i=2:pe+1
    for j=1:pe+1
        for k=pe:N-1
            r1（ij）= [x（k+j-2）+x（k-j+2）]*[conj（x（k+i-2））+conj（x（k-i+2））];
        end
    end
end

Re=r1（2:pe+1:）;%构建Re矩阵

%求解有效秩p
[U S V]=svd（Re）;%奇异值分解
for k=1:pe
    b=S（kk）/S（11）;
    if（b<0.05）
        break
    end
end 
p=k;

sp=zeros（p+1p+1）;%初始化sp矩阵p+1维
for j=1:p
    for i=1:（pe+1-p）
        sp=（S（jj）.^2）*V（i:i+pj）*V（i:i+pj）‘+sp;%V（1:71）是第1行到第7行的第一列的数据，窗是6行1列   
    end
end

inv_sp=inv（sp）;     %求矩阵sp的逆  
%根据（4.4.58），利用整体最小二乘法求出待求参数
for i=1:p+1       
    a_compute（1i）=inv_sp（i1）/inv_sp（11）;    %