clc;
clear;
fc=10.7e6;%10.7e6;%載頻
fs=0.408e6*2;
df1=1000;%多普勒
df2=2000;%多普勒
df3=3000;%多普勒
m=14.5;td=m/fs; %時延6
N=8192*2;%fs;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%白噪聲過濾低通波器 作為調制信號
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
y1 = randn(1,N); %wgn(1,fs,0);

ffty1=fft(y1);
cofft=ffty1.*conj(ffty1);% plot(abs(cofft));

f2=3.5e5;
% fs=0.408e6*2;
b = fir1(127,f2/fs);
 scope=(1:length(b))*fs/length(b);
% figure(gcf+1);
figure(1)
plot(scope,20*log10(abs(fft(b))));
modsig=conv(y1,b);
% plot(abs(fft(modsig)));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NN=N/8;%4*1024;%512*2*2*2%0.25s  data length
t=(0:1/fs:(N-1)/fs);%采樣時間序列
deviation=4;fm=5000;
% ref=cos(2*pi*fc*t+deviation*sin(2*pi*fm*t));% 參考信號
ref=cos(2*pi*fc*t+deviation*cumsum(modsig(1:length(t))));% 參考信號

ref=conv(ref,b);  %%%%%%%%%%%%%%%%%%%%%
% scope=(-2:0.1:2)  %觀察幅值的分布狀況
% figure(gcf+1);    %
% hist(ref,scope);  % 察看其柱狀圖
% grid on;

% scope=(1:NN)*fs/NN;
% figure(1);plot(scope,20*log10(abs(fft(ref(1:NN)))));

% taget=cos(2*pi*(fc+df)*(t-td)+deviation*sin(2*pi*fm*(t-td)));% 目標信號
taget1=0.0001*cos(2*pi*(fc+df1)*(t)+deviation*cumsum(modsig(1:length(t))));% 目標信號
taget2=0.0001*cos(2*pi*(fc+df2)*(t)+deviation*cumsum(modsig(1:length(t))));% 目標信號
taget3=0.0001*cos(2*pi*(fc+df3)*(t)+deviation*cumsum(modsig(1:length(t))));% 目標信號
taget1=conv(taget1,b);
n1=10;n2=20;n3=30;
delay1=[zeros(1,n1) 1];taget1=conv(taget1,delay1);
delay2=[zeros(1,n2) 1];taget2=conv(taget2,delay2);
delay3=[zeros(1,n3) 1];taget3=conv(taget3,delay3);
taget=taget1(1:N);%   +taget2(1:N)+taget3(1:N);

% snr=200;
% taget=awgn(taget,snr,'measured') ;

ref=ref(1:N);
sig=taget(1:N);
% load ref.dat;
% load sig.dat;
len=N;
L=100;
w=zeros(1,L);
% sig=cos(2*pi*20*t);
% ref=sig+0.2*randn(size(sig));
% 
% tem=sig;
% sig=ref;
% ref=tem;

%figure(gcf+1)
sig=sig+circshift(ref,[0 10])+circshift(ref,[0 200])+0.2*circshift(ref,[0 300])+0.05*circshift(ref,[0 400])+0.005*circshift(ref,[0 100])+0.001*circshift(ref,[0 200])++0.0001*circshift(ref,[0 500])+0.00001*circshift(ref,[0 1000]);
figure(2)
subplot(311);plot(sig(100:end));title('回波信號');
subplot(312);plot(ref(100:end));title('參考信號');


% u=1/abs(ref*ref');
u=0.0001
for k=1:(len-L+1)
    outp(k+L-1)=w*ref(k:k+L-1).';
    err(k+L-1)=sig(k+L-1)-outp(k+L-1);
    w=w+2*u*err(k+L-1)*conj(ref(k:k+L-1));     
end;
% % 變步長法 
% beita=0.999
% gama=0.01
%  MUmax=(1/(ref*ref'/N*L))*0.4;
%  mu=MUmax;
%     for k=1:(len-L+1)
%     outp(k+L-1)=w*ref(k:k+L-1).';
%     err(k+L-1)=sig(k+L-1)-outp(k+L-1);
%     w=w+2*mu*err(k+L-1)*conj(ref(k:k+L-1));
%     mu=beita*mu+(1-beita)*gama*MUmax;
%    end;


tem=sig';
tem1=ref'; 
subplot(313);plot(err(100:end));title('自適應對消后的信號');%(fftshift(abs(fft(err(100:end)))));
errsig=err(100:end)+ref(100:end);
%subplot(414);plot(errsig);
figure(5);
scope=(1:length(sig(1000:end)))*fs/length(sig(1000:end))-fs/2;
subplot(121);plot(scope,20*log10(fftshift(abs(fft(sig(1000:end))))));grid on;title('回波信號頻譜');
scope=(1:length(err(8000:end)))*fs/length(err(8000:end))-fs/2;
subplot(122);plot(scope,20*log10(fftshift(abs(fft(err(8000:end))))));grid on;title('自適應對消后信號頻譜');

scope=(1:length(taget(1000:end)))*fs/length(taget(1000:end))-fs/2;
figure(6);plot(scope,taget(1000:end));grid on;title('理想回波信號');

scope=(1:length(taget(1000:end)))*fs/length(taget(1000:end))-fs/2;
figure(7);plot(scope,20*log10(fftshift(abs(fft(taget(1000:end))))));grid on;title('理想回波信號頻譜');