clc;
clear all
t=1:300;
s1=sin(4*pi*t/300);%有用信號
C=1;  %噪音強度
N=3000; %某點的觀測次數
for t=1:300
   z=C*randn(1,N)+sin(4*pi*t/300);%在某個點觀察了N次的測量結果
xx(1)=0;
Q=0.001;
R=1;
p(1)=0.2;
for k=2:1:N
xs(k)=xx(k-1);
ps(k)=p(k-1)+Q;
K(k)=ps(k)/(ps(k)+R);
xx(k)=xs(k)+K(k)*(z(k)-xs(k));
p(k)=(1-K(k))*ps(k);
end
x(t)=xx(k);
end
 
plot(s1);
pause;
hold on 
plot(x,'k')