%區分調頻信號和相位編碼信號

%設計者：宋虎
%設計時間：2005-9-24

%參考文獻：

%**********線性調頻信號******************************
Fs=100e6; %采樣頻率100MHz
deltt=1/Fs; %采樣間隔
tao=7e-6; %線性調頻信號時寬 ??！也可以設置成可調參數！
SNR=3; %信噪比 ??！也可以設置成可選參數！
N=fix(tao/deltt); %采樣點數N
n=0:N-1;
t=n*deltt; %時間序列
f0=74.5e6; %起始頻率，可以設置成可調！
%      f0=75e6+(rand(1)-0.5)*10e6;
 Bw=20e6;  %帶寬10MHz
 k=Bw/tao;  %調頻系數k=帶寬(Hz)/T時寬（s);
 %phase=(2*rand(1)-1)*pi; %相位-pi~pi隨機
phase=1;
signal_LFM=exp(i*(2*pi*f0.*t+pi*k*(t.^2)+phase));
%*****************************************************
%***********二相編碼信號**********************
%首先產生仿真用的二相編碼信號 s(t)=a*exp(j*[2*pi*fc*t+pi*dt+phase0]);
%有關參數設置同參考資料p28

N=1024;      %采樣點數
phase0=1;   %初相弧度
%fc=40e6;     %載頻40MHZ
fc=90e6;
Nc=16;       %碼的長度16
%SNR=10;       %信噪比10dB
T=N*deltt;
Tc=T/Nc;   %碼元寬度
%numMC=20;  %仿真次數20
%*************************************
mayuan=[0 0 0 0 1 0 0 1 1 0 1 0 1 1 0 0]; %調制碼
n=0:N-1;
%生成二相編碼信號，關鍵在于相位在0,pi之間變化
BPSKphase=zeros(1,N);                        
for ll=0:N-1
    mayuan_number=fix(ll*deltt/Tc);
  if (mayuan_number<0)
      mayuan_number=0;
  elseif (mayuan_number>Nc -1) 
      mayuan_number=Nc -1;
  end
BPSKphase(ll+1)=mayuan(1,mayuan_number+1)*pi;
end
%*********生成二相編碼信號+噪聲
%for ii=1:numMC
%noise=randn(1,N)+i*randn(1,N);  %I,Q支路的噪聲，高斯分布，合成(復噪聲）功率2
%signal_A=sqrt(2)*10^(SNR/20);  %信號幅度
signal_BPSK=exp(i*(2*pi*fc*n*deltt+BPSKphase(n+1)+phase0)); %二相編碼信號
%signal_BPSK_noise=signal_BPSK+noise; %信號+噪聲




%計算相位突變概率
phase_LFM=angle(signal_LFM)*180/pi;
phase_BPSK=angle(signal_BPSK)*180/pi;


phase_LFM_length=length(phase_LFM);
phase_BPSK_length=length(phase_BPSK);

sum_absphase_LFM=0;
sum_absphase_BPSK=0;

for ii=1:phase_LFM_length-1
    phaseLFM_cha=phase_LFM(ii+1)-phase_LFM(ii);
   if (phaseLFM_cha>180) 
       phaseLFM_cha=phaseLFM_cha-360;
   end 
   if (phaseLFM_cha<-180)
       phaseLFM_cha=phaseLFM_cha+360;
   end
    if abs(phaseLFM_cha)>90
    sum_absphase_LFM=sum_absphase_LFM+1;
    end
end

for kk=1:phase_BPSK_length-1
    phaseBPSK_cha=phase_BPSK(ii+1)-phase_BPSK(ii);
   
    if abs(phaseBPSK_cha)>90
    sum_absphase_BPSK=sum_absphase_BPSK+1;
    end
end

%計算譜的平坦度,在帶寬內求
%找出帶寬

fft_LFM=abs(fft(signal_LFM));

figure(1)
plot(fft_LFM)


fft_BPSK=abs(fft(signal_BPSK));

figure(2)
plot(fft_BPSK)

%找到譜的最大值，并前后搜索，如果有連續三個值小于最大值的1/2，則認為搜索結束，便得到帶寬

[max_fft_LFM max_location_LFM]=max(fft_LFM);%最大值

[max_fft_BPSK max_location_BPSK]=max(fft_BPSK);%最大值

gate_LFMfft=max_fft_LFM/2;
gate_BPSKfft=max_fft_BPSK/2;

for ii=1:length(fft_LFM)-3 % 前后搜索
    if fft_LFM(ii)<gate_LFMfft && fft_LFM(ii+1)<gate_LFMfft && fft_LFM(ii+2)<gate_LFMfft && fft_LFM(ii+3)>gate_LFMfft
        
       f_start_LFM=ii+3;
    end
    if  fft_LFM(ii)>gate_LFMfft && fft_LFM(ii+1)<gate_LFMfft && fft_LFM(ii+2)<gate_LFMfft && fft_LFM(ii+3)<gate_LFMfft
    
       f_end_LFM=ii;
    end
end

for ii=1:length(fft_BPSK)-3  % 前后搜索
    if fft_BPSK(ii)<gate_BPSKfft && fft_BPSK(ii+1)>gate_BPSKfft % && fft_BPSK(ii+2)>gate_BPSKfft %&& fft_BPSK(ii+3)>gate_BPSKfft
        
       f_start_BPSK=ii+1;
    end
    if  fft_BPSK(ii)>gate_BPSKfft && fft_BPSK(ii+1)<gate_BPSKfft %&& fft_BPSK(ii+2)<gate_BPSKfft %&& fft_BPSK(ii+3)<gate_BPSKfft
    
       f_end_BPSK=ii;
    end
end
%%%%?????????BPSK信號的帶寬計算有問題？  f_start_BPSK f_end_BPSK  都是同一個值

%于是就可以得到帶寬內的譜值
fft_LFM_inBw=fft_LFM(f_start_LFM:f_end_LFM);
fft_LFM_inBw_dB=20*log10(fft_LFM_inBw);

fft_BPSK_inBw=fft_BPSK(f_start_BPSK:f_end_BPSK);
fft_BPSK_inBw_dB=20*log10(fft_BPSK_inBw);

fft_LFM_inBw_dB_var=var(fft_LFM_inBw_dB);
fft_BPSK_inBw_dB_var=var(fft_BPSK_inBw_dB);