function [x,y,zx,zy,XXE,YYE,XERB,YERB,XSTD,YSTD]=IMM(T,mont,d)
%輸入參數：觀測噪聲方差d,采樣間隔T,蒙特卡洛次數mont
%返回值：
totalTime=800;  
%真實航跡產生
[x,y]=realTrack(T,totalTime);
%隨機初始化
randn('state',sum(100*clock));
D=d*d;
N=ceil(totalTime/T);
for n=1:mont;
    %觀測數據產生
    for i=1:N
        zx(i)=x(i)+d*randn(1);
        zy(i)=y(i)+d*randn(1);
    end
    %濾波(IMM)
    [XXE,YYE]=immkf(zx,zy,T,D);   
    XXE(1)=zx(1);XXE(2)=zx(1);
    YYE(1)=zy(1);YYE(2)=zy(1);
    %誤差矩陣
    XER(1:N,n)=x(1:N)-(XXE(1:N))';
    YER(1:N,n)=y(1:N)-(YYE(1:N))';
end
%end for mont

%濾波誤差的均值
XERB=mean(XER,2);
YERB=mean(YER,2);

%濾波誤差的標準差
XSTD=std(XER,1,2);
YSTD=std(YER,1,2);

% %作圖
% figure(1)
% plot(x,y,'r');hold on;
% plot(zx,zy,'g');hold on;
% plot(XXE,YYE,'b');
% axis([1500 5000 -5000 11000]);
% title('IMM算法仿真');
% legend('真實軌跡','觀測數據','濾波估計');
% 
% figure(2)
% subplot(2,2,1)
% plot(XERB)
% axis([0 500 -50 50])
% xlabel('觀測次數')
% ylabel('X方向濾波誤差均值')
% subplot(2,2,2)
% plot(YERB)
% axis([0 500 -50 50])
% xlabel('觀測次數')
% ylabel('Y方向濾波誤差均值')
% subplot(2,2,3)
% plot(XSTD)
% axis([0 500 0 150])
% xlabel('觀測次數')
% ylabel('X方向濾波誤差標準值')
% subplot(2,2,4)
% plot(YSTD)
% axis([0 500 0 150])
% xlabel('觀測次數')
% ylabel('Y方向濾波誤差標準值')