%BPSK的識別率分析，一個碼元對應一個周期正弦波（Fs＝2000）
clear
figure(64)
hold on;
for iiii=1:20 %信噪比循環
k=0;
for ii=1:10  %仿真100次
    iiii
    ii
df=0.25;
fs=20;%20MHz
fc=2;
rb=0.2;
ts=1/fs;
m=256;%m是碼元數，自己隨便定義
n=40;%n是每個碼元取點數，即對應的正弦波4個周期取多少個點
a=sj(m,n);
%sj(m,n)的解釋參見sj.m,sj(m,20)的作用相當于上面的for循環
t=0:1:m*n-1;
Fc1=sin(8*pi*t/n);%碼元0發射的載波
Fc2=sin(8*pi*t/n+pi);%碼元1發射的載波
b=1-a;
u0=a.*Fc2+b.*Fc1;

%產生白噪聲
snr=iiii;%S/N
snr_lin=10^(snr/10);%線性信噪比
signal_power=(norm(u0)^2)/length(u0);
noise_power=signal_power/snr_lin;
noise_std=sqrt(noise_power);
noise=noise_std*randn(1,length(t));
u1=u0+noise;%加白噪聲
uls=length(u1);
maxu1=max(u1);
u=u1/maxu1;

%求瞬時平方包絡
v=blpf(u);
%均值
miu=mean(v);
%方差
delt=var(v);
%求R參數
r=delt/(miu^2);


%信號平方
sqrtu=u1.*u1;

%求功率譜
nfft=1024;
u2=Rx_set(u,nfft);  %這個100是取了100個自相關系數a
u3=fftshift(abs(fft(u2)));%求功率譜
ff=-1*fs/2:fs/nfft:fs/2;

%求AWGN因子
e1=mean(u3);
j=0;
  for i=1:nfft+1
      if u3(i)>e1
          j=j+1;
          ee(j)=u3(i);
      end
  end
 e2=mean(ee);  
 j=0;
  for i=1:nfft+1
      if u3(i)>e2
          j=j+1;
          eee(j)=u3(i);
      end
  end
e3=mean(eee); 
a1=(e2-e1)/e1;
a2=(e3-e2)/e1;

%求平方功率譜及共振峰個數
sqrtu2=Rx_set(sqrtu,nfft);  %這個nfft是取了nfft個自相關系數a
sqrtu3=fftshift(abs(fft(sqrtu2)));%求功率譜
ls=length(sqrtu3);
lsl=round((fs/2+2*fc-1)/fs*ls);
lsr=round((fs/2+2*fc+1)/fs*ls);
c=20;
gzf=0;
for iii=lsl:lsr
    if sqrtu3(iii)>20
        gzf=1;
    end
end



%頻譜分析
[M,m,df1]=fftseq(a,ts,df);
M=M*ts;
f=[0:df1:df1*(length(m)-1)]-fs/2;
[U,u,df1]=fftseq(u,ts,df);
U=U*ts;

[SQRTU,sqrtu1,df1]=fftseq(sqrtu,ts,df);
SQRTU=SQRTU*ts;



%識別
if r<0.55
    if a2<10.5 & a2>0.7
        if a1<=2.2
            if gzf==1
                k=k+1;
            end

        end
    end
end

end %一次信噪比循環里的仿真循環結束
rtio(iiii)=(k/ii)*100;
end %信噪比循環結束
iiii=1:20;
plot(iiii,rtio,'*')
grid;
hold off;
title('BPSK 識別率(仿真100次)')