%信源為單極性歸零碼，取樣點(diǎn)數(shù)自己定義 

global dt df t f N T %定義全局變量
close all  %關(guān)閉以前的應(yīng)用窗口
clear Eb_N0 Pe

k=input('取樣點(diǎn)數(shù)＝2^k, k=[13]'); %采樣點(diǎn)數(shù)自定義，默認(rèn)為2^13
if k==[], k=13; end

%Rt=input('占空比[0：1](默認(rèn)0.5,為1則是NRZ碼)');
%if Rt==[],Rt=0.5;
%end;
Rt=0.5;   %占空比是0.5
N=2^k;    %總的取樣點(diǎn)數(shù)
L=32;		%L為每個(gè)碼元的取樣點(diǎn)數(shù)
M=N/L;		%M碼元總數(shù)	
Rb=2;		%碼元速率Rb=2Mb/s
Ts=1/Rb;		%碼元寬度Ts
dt=Ts/L;	%時(shí)域取樣間隔
df=1/(N*dt);%頻域取樣間隔
T=N*dt;     %時(shí)間截?cái)嚅L度
Bs=N*df/2;  %帶寬
Na=2;      %示波器掃描寬度為2個(gè)碼元
alpha=0.5;   %升余弦滾降系數(shù)
%f=[-Bs+df/2:df:Bs];   %頻域橫坐標(biāo)
%t=[-T/2+dt/2:dt:T/2]; %時(shí)域橫坐標(biāo)
t=linspace(-T/2,T/2,N); %頻域橫坐標(biāo)
f=linspace(-Bs,Bs,N)+eps;%時(shí)域橫坐標(biāo)


%。。。。。。。。。。。。。。。。。升余弦
hr1=sin(pi*t/Ts)./(pi*t/Ts);
hr2=cos(alpha*pi*t/Ts)./(1-(2*alpha*t/Ts).^2);
hr=hr1.*hr2;
HR=abs(t2f(hr));
GT=sqrt(HR);
GR=GT;

figure(1)
 set(1,'Position',[10,250,500,400])
                %設(shè)定窗口位置及大小 
figure(2)
 set(2,'Position',[515,250,500,400])
                %設(shè)定窗口位置及大小 
                
figure(3)
 set(3,'Position',[360,10,300,200])
                %設(shè)定窗口位置及大小 
%。。。。。。。。。。。。。。。。。功率譜矩陣，初值為零
EP=zeros(size(f)); 
EPa=zeros(size(f)); 
EPr=zeros(size(f)); 
EPo=zeros(size(f)); 
 
for loop1=1:20
    Eb_N0(loop1)=3*(loop1-1) ;   %Eb/N0 in dB，橫坐標(biāo)
    eb_n0(loop1)=10^(Eb_N0(loop1)/10);%求其功率
    Eb=0.1;
    n0=Eb/eb_n0(loop1);  %信道的噪聲譜密度
    sita=n0*Bs;    %信道中噪聲功率
    n_err=0;     %誤碼計(jì)數(shù)
for loop2=1:10  
%a=2*round(rand(1,M));   %產(chǎn)生隨機(jī)碼
a=round(rand(1,M));   %產(chǎn)生隨機(jī)碼
s=zeros(1,N);
ss=zeros(1,N);
for ii=1:Rt*Ts/dt;   %Rt為占空比，dt為采樣間隔，
                     %Rt*Ts/dt為每個(gè)碼元中的非零值點(diǎn)數(shù)
s(ii+[0:M-1]*L)=a;   %信源變成單極性歸零碼
end
ss(L/2:L:N)=a/dt;
SS=t2f(ss); 
S=t2f(s);           %信源信號(hào)的傅式變換
P=S.*conj(S)/T;     
EP=(EP*(loop1-1)+P)/loop1; %信源功率譜


s1=real(f2t(SS.*GT)); %經(jīng)發(fā)送濾波后的信源信號(hào)
S1=t2f(s1);
Pa=S1.*conj(S1)/T;
EPa=(EPa*(loop1-1)+Pa)/loop1;%信道前信號(hào)功率譜


n_ch=sqrt(sita)*randn(size(t));  %信道噪聲
nr=real(f2t(t2f(n_ch).*GR));    %輸出噪聲

dtt=Ts/2;
C=abs(1-0.5*exp(-j*2*pi*f*dtt));%雙徑信道
sr=real(f2t(S.*HR.*C))+nr; %接收信號(hào)
%sr=real(f2t(SS.*HR))+nr; %接收信號(hào)
SR=t2f(sr);
Pr=SR.*conj(SR)/T;
EPr=(EPr*(loop1-1)+Pr)/ loop1;%接收信號(hào)功率譜
y=sr(L/2:L:N);       %采樣
aa=sign(y-0.05);    %判決
%aa=sign(y-0.125);    %判決
aa=(aa+1)/2;         %判決輸出信號(hào)
%aa=aa+1;
n_err=n_err+length(find(aa~=a)); %誤碼率計(jì)數(shù)
end
Pe(loop1)=n_err/(M*loop2);   %誤碼率
e=zeros(1,N);
for jj=1:Rt*Ts/dt;   %Rt為占空比，dt為采樣間隔，
                     %Rt*Ts/dt為每個(gè)碼元中的非零值點(diǎn)數(shù)
e(jj+[0:M-1]*L)=aa; %將判決輸出信號(hào)轉(zhuǎn)化為半占空比
end
E=t2f(e);
Po=E.*conj(E)/T;
EPo=(EPo*(loop1-1)+Po)/loop1;%輸出信號(hào)功率譜
end
figure(1)
 subplot(2,2,1)
 plot(t,s,'r');
 axis([-T/16,T/16,min(s)-0.5,max(s)+0.5])
 xlabel('t (ms)');
 ylabel('s  (V/KHz)');
 title('單極性歸零信號(hào)波形');
 
 subplot(2,2,2)
 bb=30+10*log10(EP+eps);
 plot(f,bb,'g')
 axis([-20,20,-20,max(bb)+0.5])
 grid
 xlabel('f(Mhz)');
 ylabel('Ps(f) (dBm/MHz)');
 title('信源功率譜密度圖形');
 
 subplot(2,2,3)
 plot(t,e,'r');
 axis([-T/16,T/16,min(e)-0.5,max(e)+0.5])
 xlabel('t (ms)');
 ylabel('e  (V/KHz)');
 title('輸出信號(hào)波形');
 
 subplot(2,2,4)
 bb=30+10*log10(EPo+eps);
 plot(f,bb,'g')
 axis([-20,20,-20,max(bb)+0.5])
 grid
 xlabel('f(Mhz)');
 ylabel('Po(f) (dBm/MHz)');
 title('輸出信號(hào)功率譜密度');
 
 figure(2)
 subplot(2,3,1)
 plot(t,s1,'b');
 axis([-T/16,T/16,min(s1)-0.5,max(s1)+0.5])
 xlabel('t (ms)');
 ylabel('s1  (V/KHz)');
 title('信道前信號(hào)波形');
 
 subplot(2,3,2)
  tt=[0:dt:Na*L*dt];
  hold on
  for jj=1:Na*L:N-Na*L
    plot(tt,s1(jj:jj+Na*L)); %信道前信號(hào)眼圖
  end
  title('信道前升余弦滾降眼圖')
  
 subplot(2,3,3)
 bb=30+10*log10(EPa+eps);
 plot(f,bb,'g')
 axis([-20,20,-20,max(bb)+0.5])
 grid
 xlabel('f(Mhz)');
 ylabel('Pa(f) (dBm/MHz)');
 title('信道前信號(hào)功率譜密度圖形');
 
subplot(2,3,4)
 plot(t,sr,'b');
 axis([-T/16,T/16,min(sr)-0.5,max(sr)+0.5])
 xlabel('t (ms)');
 ylabel('sr  (V/KHz)');
 title('接收信號(hào)波形');
 
  subplot(2,3,5)
  tt=[0:dt:Na*L*dt];
  hold on
  for jj=1:Na*L:N-Na*L
    plot(tt,sr(jj:jj+Na*L)); %接收信號(hào)眼圖
  end
  title('信道后升余弦滾降眼圖')
  
 subplot(2,3,6)
 bb=30+10*log10(EPr+eps);
 plot(f,bb,'g')
 axis([-20,20,-20,max(bb)+0.5])
 grid
 xlabel('f(Mhz)');
 ylabel('Pr(f) (dBm/MHz)');
 title('接收信號(hào)功率譜密度圖形');
 
  figure(3) 
  semilogy(Eb_N0,Pe,'g');
  eb_n0=10.^(Eb_N0/10);
  hold on
  plot(Eb_N0,0.5*erfc(sqrt(eb_n0/2)));
  axis([0,19,1e-4,1])
  xlabel('Eb/N0')
  ylabel('Pe')
  title('誤碼率曲線')