%張潤寧論文中多普勒補償問題,【2.^6-1=63】位m－序列（加噪聲）
%改變fd的值可以改變多普勒頻率的大小
%改變SNR的值可以改變信噪比的大小
clear all;
T=0.2*10.^-6;                         %子脈沖寬度T
B=1/T;                                %子脈沖帶寬B=1/T
fd=0*B                              %多普勒頻率fd    

SNR=0                                 %信噪比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];						                                       %反饋公式中系數
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;
sequence=ms;

SeqLength=length(sequence);          %序列長度P=SeqLength 

MFOut3=zeros(1,SeqLength+SeqLength-1);
MFOut1=zeros(1,SeqLength+SeqLength-1);;
%for kkk=1:1
%%%%%%%%%%%%%%%%%
%加多普勒信息和噪聲
Vector=1:SeqLength;
%SignalR=BPCSeq.*exp(1i*2*pi*fd*T*Vector);                                                    %回波信號加多普勒信息
rand('state',sum(100*clock));
SignalR=sqrt(Psignal)*sequence.*exp(1i*2*pi*fd*T*Vector)+sqrt(Pnoise)*randn(1,SeqLength);            %SignalR_N為加噪聲和多普勒信息的回波信號
%%%%%%%%%%%%%%%%%

s1=SignalR.*SignalR;

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%用兩次內插的方法來實現分頻處理
 u=s1(1:(SeqLength+1)/2);
 for k=0:(SeqLength-1)/2
     %k
     Even=2*k;
     Odd=2*k+1;
     v(Even+1)=u(Even/2+1);
     if (Odd<SeqLength)
     v(Odd+1)=(u((Odd-1)/2+1)+u((Odd+1)/2+1))/2;
     end
 end
 s2=v;
 %s2為分頻后的結果
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 s3=SignalR.*conj(s2);

matchfilter=sequence(SeqLength:-1:1);          	%全景匹配濾波器

MFOut31=conv(s3,matchfilter);

%mm_A=max(abs(MFOut3))
MFOut3=MFOut3+MFOut31;

MFOut11=conv(SignalR,matchfilter);

%mm_1=max(abs(MFOut1))
MFOut1=MFOut1+MFOut11;
%end
[PSL_A,ISL_A]=PSLandISL(MFOut3)
[PSL_1,ISL_1]=PSLandISL(MFOut1)
figure(1)
   subplot(2,2,1)
 plot(20*log10(0.01+abs(MFOut1)/max(abs(MFOut1))));
 axis([0,length(MFOut1),-40,5]);
 title('A(補償前)Noise=0,fd/B=0')
 xlabel('距離單元')
 ylabel('幅度（dB）');
 grid on

 

%匹配濾波器的輸出基本是和論文中的圖一致的（P26）
%但是補償后的脈壓輸出有很大的不同
%問題估計存在于插值實現分頻的計算中
%經過檢查，不是插值的問題
%完成，開始錯在公式s3=s.*conj(s2),寫成s3=s2.*conj(s2)了
%如何在信號中加噪聲？？x=s+n? 解決了見公式[s=s+Pnoise*rand(1,SeqLength)]
%2002/3/25增加了計算SNR損失的函數PSLandISLandMSSLandLSNRandLPG