function output=lr_lmmse_estimation(input,pilot_inter,pilot_sequence,pilot_num,trms,t_max,snr,cp);
%trms為多經信道的平均延時，t_max為最大延時,此處所有的時間都是已經對采樣間隔做了歸一化后的結果
beta=17/9;
[N,NL]=size(input);
Rhh=zeros(N,N);
for k=1:N
    for l=1:N
        Rhh(k,l)=(1-exp((-1)*t_max*((1/trms)+j*2*pi*(k-l)/N)))./(trms*(1-exp((-1)*t_max/trms))*((1/trms)+j*2*pi*(k-l)/N));
    end
end
%以下過程是對Rhh矩陣進行特征值分解，并按特征值大小排隊，選最大的16個特征值
[U,D]=eig(Rhh);%U為滿秩酉陣，D是以特征值為主對角線的對角陣
dlamda=diag(D);%取D對角元素構成列向量
[dlamda_sort,IN]=sort(dlamda);%按升序排列各特征值，dlamda_sort為排序結果，IN為各元素在原向量中的位置索引
for k=1:N
    lamda_sort(k)=dlamda_sort(N-k+1);%按照降序排列各特征值
    IN_new(k)=IN(N-k+1);
end
%以下按照IN_new順序排列U的各列
U_new=zeros(N,N);
for k=1:N
    U_new(:,k)=U(:,IN_new(k));
end
%以下只取前cp個特征值構成新的對角陣
delta=zeros(N,1);
for k=1:N
    if k<=cp
       delta(k)=lamda_sort(k)/(lamda_sort(k)+beta/snr);
    else
        delta(k)=0;
    end
end
D_new=diag(delta);
output=zeros(N,NL-pilot_num);
i=1;
count=0;
while i<=NL
    Hi=input(:,i)./pilot_sequence;
    Hlr=U_new*D_new*(U_new')*Hi;
    count=count+1;
    if  count*pilot_inter<=(NL-pilot_num)
        for p=((count-1)*pilot_inter+1):count*pilot_inter
            output(:,p)=input(:,(i+p-(count-1)*pilot_inter))./Hlr;
        end
    else
        for p=((count-1)*pilot_inter+1):(NL-pilot_num)
            output(:,p)=input(:,(i+p-(count-1)*pilot_inter))./Hlr;
        end
    end
    i=i+pilot_inter+1;
  end