%張潤寧論文中多普勒補償問題,對于多個散射中心的目標，M序列
%2002/4/2，目標為多個散射中心構成（面目標？）
clear all;
T=0.2*10.^-6;                                                           %子脈沖寬度T
B=1/T;                                                                  %子脈沖帶寬B=1/T
fd=0.05*B;                                                              %多普勒頻率fd    
lim_amp=0;                                                              %lim_amp==1為限幅脈壓
SNR=10;                                                                 %信噪比SNR 
Pnoise=0.01;                                                            %噪聲功率Pnoise
Psignal=Pnoise*10.^(SNR/10);                                            %信號功率Psignal
%Psignal=1;
%Pnoise=0;

x=[0,0,1,0,1,1];  
y=[6,5];						                                        %反饋公式中系數(shù)
k=length(x);					
ms=x;
for i=k+1:2^k-1
   ms(i)=xor(ms(i-y(1)),ms(i-y(2))) ;
end;
ms=ms*-2+1;
LengthSeq=length(ms);                                                    %序列長度P=SeqLength 
matchfilter=ms(LengthSeq:-1:1);          	                             %全景匹配濾波器
m=1:LengthSeq;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%多個散射中心的情況信號回波的疊加
DotNum=4;
Amplitude=[6,8,10,8];
TimeDelay=[0,10,15,25];

MaxLeng=(max(TimeDelay)+LengthSeq+1);
U_matrix=zeros(DotNum,MaxLeng);
 for m=1:DotNum
   U_matrix(m,TimeDelay(m)+1:TimeDelay(m)+LengthSeq)=ms;
 end
 
MFOut3=[];
for kk=1:20                                                              %改變kk的取值范圍可以做非相干積累          
    rand('state',sum(100*clock));
    s=sqrt(Psignal)*Amplitude*U_matrix.*exp(j*2*pi*fd*((1:length(U_matrix))*T-TimeDelay(1)*T)+sqrt(Pnoise)*rand(1,MaxLeng));
%s=Amplitude*U_matrix;
if lim_amp==1
	s(find(s>0))=1;
	s(find(s<0))=-1;
end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%
s1=s.*(s);                                                               %張潤寧文章中的
%s1=s.*conj(s);                                                          %改進的。
%%%%%%%%                                                              

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%用兩次內插的方法來實現(xiàn)分頻處理
SeqLength=length(U_matrix);
 u=s1(1:(SeqLength+1)/2);
 for k=0:(SeqLength-1)/2
     %k
     EvenN=2*k;
     OddN=2*k+1;
     v(EvenN+1)=u(EvenN/2+1);
     if (OddN<SeqLength)
     v(OddN+1)=(u((OddN-1)/2+1)+u((OddN+1)/2+1))/2;
     end
 end
 s2=v;
 %s2為分頻后的結果
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 s3=s.*conj(s2);
 MFOut33(kk,:)=conv(s3,matchfilter);
 MFOut11(kk,:)=conv(s,matchfilter);

 kk=kk+1;
end
MFOut3=sum(abs(MFOut33),1);
MFOut1=sum(abs(MFOut11),1);

%MFOut1=conv(s,matchfilter);

figure(2)
subplot(2,2,1)
plot(abs(MFOut11(1,:))/max(abs(MFOut11(1,:))))
title('（A）fd=0.05B匹配脈壓輸出')
axis([0,length(MFOut11),0,1.0])
grid;zoom
subplot(2,2,2)
plot(abs(MFOut33(1,:))/max(abs(MFOut33(1,:))));
title('（B）fd=0.05B多普勒補償后脈壓輸出 ')
axis([0,length(MFOut33),0,1.0])
grid;zoom


subplot(2,2,3)
plot(abs(MFOut1)/max(abs(MFOut1)))
title('（C）fd=0.05B匹配脈壓30次非相干積累')
axis([0,length(MFOut1),0,1.0])
grid;zoom
subplot(2,2,4)
plot(abs(MFOut3)/max(abs(MFOut3)));
title('（D）fd=0.05B多普勒補償后30次非相干積累 ')
axis([0,length(MFOut3),0,1.0])
grid;zoom
%subplot(2,2,3)
%plot(abs)