%張潤寧論文中多普勒補償問題, 巴克碼（加噪聲）
%改變fd的值可以改變多普勒頻率的大小
%改變SNR的值可以改變信噪比的大小
clear all;
T=0.2*10.^-6;                         %子脈沖寬度T
B=1/T;                                %子脈沖帶寬B=1/T
fd=0.05*B                              %多普勒頻率fd    

SNR=20                                 %信噪比SNR 
Pnoise=0.01;                          %噪聲功率Pnoise
Psignal=Pnoise*10.^(SNR/10);          %信號功率Psignal
%Psignal=1;
%Pnoise=0

sequence=[1 1 1 1 1 -1 -1 1 1 -1 1 -1 1];

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

%%%%%%%%%%%%%%%%%
%加多普勒信息和噪聲
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);          	%全景匹配濾波器

MFOut3=conv(s3,matchfilter);
[PSL_A,ISL_A]=PSLandISL(MFOut3)
mm_A=max(abs(MFOut3))

MFOut1=conv(SignalR,matchfilter);
[PSL_1,ISL_1]=PSLandISL(MFOut1)
mm_1=max(abs(MFOut1))

figure(1)
 subplot(2,2,3)
 plot(20*log10(0.01+abs(MFOut1)/max(abs(MFOut1))));
 title('C(補償前)S/N=6dB,fd/B=0.0')
 y=length(MFOut1);
 axis([0,y,-40,5]);
 grid on
 subplot(2,2,4)
 plot(20*log10(0.01+abs(MFOut3/max(abs(MFOut3)))));
 title('D(補償后)S/N=6dB,fd/B=0.0')
 y=length(MFOut3);
 axis([0,y,-40,5]);
 grid on
 %subplot(2,2,3)
 %plot(s)
 %subplot(2,2,4)
 %plot(s3)

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