clear all
close all
clc
 
 
a=.5;
db=1;
L=801;
N=800;
 %for ff=0:10
U=0.0005;%正則化因子
aa=1/500;
 
for f=0:39
d=50;
v=1400;
t0=0.3;
 
fm=25;
A=1;
FS0=500;
tx=sqrt(t0*t0+((f*d)*(f*d))/(v*v));
 
[t,s]=ricker(FS0,fm,tx,A);
sigpower=10*log10(sum(s.^2)/length(s));
disp(sigpower);
x=awgn(s,-5,'measured','db');
w=x-s;
disp('初始信噪比');
snr1=10*log10((sum(s.^2)/length(s))/(sum(w.^2)/length(w)));disp(snr1)
 
 
bs=10^(db/10);%?????????ù??±???
%w=randn(1,L);%?ú?ú?ù????0?????ú????
%s=zeros(1,L);
%x=zeros(1,L);
r_ss=zeros(1,L);% r_ss??s???à????????r_ss?¨1????????r_ss(0),??????
r_vv=zeros(1,L);
r_xx=zeros(1,L);
r_xs=zeros(1,L);
%s(1)=w(1);%s(0)=0
%for i=2:L
    %s(i)=a*s(i-1)+w(i);
    %end
for i=1:L
    r_ss(1)=r_ss(1)+s(i)*s(i);
end
r_ss(1)=r_ss(1)/L;
r_vv(1)=r_ss(1)/bs;% r_vv(1)??°×???ù??????????
%v=sqrt(r_vv(1))*randn(1,L);%???ú?ú?ù????0??????????r_vv(0)??°×???ù
%x(1:L)=s(1:L)+v(1:L); % x(n)????????????????????×é
r_xx(1)=r_ss(1)+r_vv(1);% r_xx??x(n)??×??à??????
r_xs(1)=r_ss(1);% r_xs??x(n)??s(n)?????à??????
for i=2:N
    %r_ss(i)=a^(1-i)*r_ss(1);
     r_ss(i)=a^(i-1)*r_ss(1);
    r_vv(i)=0;
    r_xx(i)=r_ss(i)+r_vv(i);
    r_xs(i)=r_ss(i);
end
R_xs=zeros(N,1);
R_xx=zeros(N,N);
for i=1:N
    k=1;
    for j=i:N
        R_xx(i,j)=r_xx(k);
        k=k+1;
    end
   if i>1
       k=2;
      for j=i-1:1
          R_xx(i,j)=r_xx(k);
          k=k+1;
      end
  end
  R_xs(i)=r_xs(i);
end

R_xx_aa= circshift(R_xx,[0,-2*aa*500]);
R_xxaa= circshift(R_xx,[0,2*aa*500]);
R_xs_aa= circshift(R_xs,[0,-2*aa*500]);
R_xsaa= circshift(R_xs,[0,2*aa*500]);

RR=(50+2*U/aa)*R_xx-(U/aa)*R_xx_aa-(U/aa)*R_xxaa;
RX=(50+2*U/aa)*R_xs-(U/aa)*R_xs_aa-(U/aa)*R_xsaa;

hopt=(inv(RR)*RX);
%E=r_ss(1)-(conj(R_xs))'*hopt
%Emin=abs(r_ss(1)-(conj(R_xs))'*hopt)
 
y=conv(hopt',x);
snr1=10*log10(((sum(y(1:801)).^2)/length(s))/(sum((y(1:801)-s).^2)/length(s)));disp(snr1)
er=0;
for i=1:L
    er=er+((s(i)-1*y(i)).^2)/L;% s z 
end
disp(er)
hold on;
figure(1);
plot(t,s+f+1,'b');
axis([0.2,1.6,0,41]);
hold on;
%figure(2);
%plot(t,w+f+1,'k');
%axis([0.2,1.6,0,41]);
%hold on;
figure(3);
plot(t,x+f+1,'-m');
axis([0.2,1.6,0,41]);
hold on;
 

figure(4);
plot(0:1/500:1299*0.002,y+f+1,'-r');
axis([0.2,1.6,0,41]);
hold on;
 
xlabel('n'),ylabel('觀測數據x(n)');
%grid on;
title('觀測數據x的波形');
end
%if ff>10
%    clear figure(ff-10);
%end
 %end
%figure(5);
%plot(R_xx)