% 基于LMS的自適應時延估計
%  Lmstde_a
close all; clear all; clc
L=5;           %取5個信號周期;L=6,7,8,9,10
fs=1000;       %抽樣頻率
t=(0:fs)/fs;   %抽樣時間間隔
s1=sin(2*pi*L*t); %正弦信號 
m=length(t);  %m=1001

N=50;
delta=0.001;
J=10;
if J<0
    d=[zeros(1,-J) s1]; 
    x=[s1 zeros(1,-J)];
else
    x=[zeros(1,J) s1];
    d=[s1 zeros(1,J)];
end
M=length(x); 
t1=(0:M-1)/fs;
y=zeros(1,M); %輸出初始化
e=zeros(1,M);
h=zeros(1,2*N+1); 
for n=N:M-N-1
    x1=x(n-N+1:n+N+1);
    y(n)=h*x1';                   %計算n時刻的輸出
    e(n)=d(n)-y(n);               %計算n時刻的誤差
    h=h+delta*e(n)*x1;            %調整濾波器的系數
    hh(n,:)=h;
end
figure(1)
subplot(311);
plot(t1,d,'r',t1,x,'b');
legend('參考信號','輸入信號'); grid;
title('s信號經過一段延時后的信號');
subplot(312);
plot(t1,d,'r',t1,y,'b'); 
legend('參考信號','輸出信號'); grid;
title('經自適應濾波器處理后的信號');
subplot(313);
plot(t1,e); grid;
title('誤差信號');

for k=N : M-N-1
    h=hh(k,:);
    [hmax,hloc]=max(h);
    Hmax(k)=hmax; Hloc(k)=hloc;
end
HL=length(Hmax);

[HHM,HHL]=max(Hmax);
%figure(2);
%plot(hh(HHL,:));
delay=(Hloc(HHL)-N-1);  %延遲的樣點數
delaytime=delay/fs;     %計算出延遲的時間
fprintf('延遲樣點數=%4d   延遲時間=%5.6f(秒)\n',delay,delaytime);