www.pudn.com > GaoJie.rar > gaojie.m

%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)