global dt t df N f
close all
clear Eb_N0 Pe
N=2^14; 	 %采樣點數
L=64;   		     
M=N/L   	 %碼元數
Rb=2;		 %碼速率是2Mb/s
Ts=1/Rb;       %碼元間隔
dt=Ts/L;        %時域采樣間隔
df=1/(N*dt);   %頻域采樣間隔
T=N*dt ;        %截短時間
Bs=N*df/2   ;  %系統帶寬
Na=4;           %示波器掃描寬度為4個碼元
f0=5;

alpha=0.5
t=linspace(-T/2,T/2,N);
f=linspace(-Bs,Bs,N)+eps;
g1=sin(pi*t/Ts)./(pi*t/Ts);
g2=cos(alpha*pi*t/Ts)./(1-(2*alpha*t/Ts).^2);
g=g1.*g2;		 %升余弦脈沖波形
G=t2f(g);
GG=abs(G);
GR=sqrt(GG);
GT=GR;

EPIN=zeros(1,N)+eps;
EPOUT=zeros(1,N)+eps;
EP=zeros(1,N)+eps;

for loop1=1:10
 Eb_N0(loop1)=(loop1-1)  ;
 eb_n0(loop1)=10^(Eb_N0(loop1)/10);%分貝值變為真值
 Eb=1;	                     %每碼元能量為1
 n0=Eb/eb_n0(loop1);         %信道噪聲譜密度 
 noiseP=n0*Bs; 			     %噪聲功率 
 n_err=0;  

for j1=1:3
 a=sign(randn(1,M));   %發送M個碼元,序列為a  －1，＋1
 
  for ii=1:L
        anrz(ii+[0:M-1]*L)=a;
    end
 
 imp=zeros(1,N);      %產生沖激序列
 imp(L/2:L:N)=a/dt;
 IMP=t2f(imp);
 S=IMP.*GT;         %經過發送成型濾波器后的輸入信號的傅氏變換
 PIN=S.*conj(S)/T;
 EPIN=(EPIN*(j1-1)+PIN+eps)/j1;   %經過發送成型濾波器后的輸入信號的功率譜
 logpin=10*log10(EPIN+eps);
 s=f2t(S);
 ss=s.*cos(2*pi*f0*t);
 SS=t2f(s);
 P=SS.*conj(S)/T;
 EP=(EP*(j1-1)+P+eps)/j1; 
 logp=10*log10(EP+eps);
 
 n_ch=sqrt(noiseP)*randn(size(t)); %信道噪聲
 nr=real(f2t(t2f(n_ch.*cos(2*pi*f0*t)).*GR)); %輸出的噪聲實部!!!!!!!!!!!!!
 
 R=t2f(ss.*cos(2*pi*f0*t)).*GR  ; %接收到的信號  !!!!!!!!!!!!!!!!!!!!!!!!
 POUT=R.*conj(R)/T;
 EPOUT=(EPOUT*(j1-1)+POUT+eps)/j1; 
 logpout=10*log10(EPOUT+eps);
 r=f2t(R);
 
 rey=real(r)+nr;
 REY=t2f(rey); 
 yy=rey(L/2:L:N);
 aa=sign(yy);
 n_err=n_err+length(find(aa~=a)) ; 
end
Pe(loop1)=n_err/(M*j1);	%誤碼率 %M碼元數   
for ii=1:L
        aanrz(ii+[0:M-1]*L)=aa;
    end
if loop1==10
figure(1)
hold on
 tt=[0:dt:Na*64*dt];
 for jj=1:Na*32:N-Na*32
  plot(tt,rey(jj:jj+Na*32));
 end
grid 
axis([0,2,-1.5,1.5]);
 xlabel('t (us)')
 ylabel('s(t) (V)');
 title('接收端抽樣前的眼圖  ');
 set(1,'position',[10,10,300,200])
end
end

figure(2)
plot(f,logpin,'r')
xlabel('f(MHZ)');
ylabel('Ps(f)(dB)');
title('發送信號的功率譜');
set(2,'position',[250,10,300,200]);
axis([-5,5,-100,50]);

figure(3)
hold on
 tt=[0:dt:Na*32*dt];
 for jj=1:Na*32:N-Na*32
  plot(tt,s(jj:jj+Na*32));
 end
grid 
axis([0,2,-3,3]);
 xlabel('t (us)')
 ylabel('s(t) (V)');
 title('基帶信號的眼圖  ');
set(3,'position',[550,10,300,200]);

figure(4)
plot(f,logp,'r')
xlabel('f(MHZ)');
ylabel('Ps(f)(dB)');
title('頻帶信號的功率譜');
set(4,'position',[800,10,300,200]);
axis([-10,10,-100,50]);

 figure(5) 
  semilogy(Eb_N0,Pe,'g');       %y軸為對數坐標
eb_n0=10.^(Eb_N0/10);	%還原為真值 
  hold on
  semilogy(Eb_N0,0.5*erfc(sqrt(eb_n0)),'red');
  axis([0,9,1e-4,1])
 title('pe曲線');
  xlabel('Eb/N0  dB')
  ylabel('Pe')   
set(5,'position',[10,250,300,200])  ;

  figure(6)
  plot(t,f2t(S),'b');
  axis([-5,5,-2.5,2.5]);
  grid;
  title('發送基帶信號波形');
  xlabel('t  (uV)')
  ylabel('s(t)') 
  set(6,'position',[250,250,300,200])  ;

  
  figure(7)
  plot(f,abs(S),'b');
  axis([-3,3,0,40]);
  grid;
  title('發送基帶信號頻譜');
  xlabel('f (MHz)')
  ylabel('S(f)') 
  set(7,'position',[550,250,300,200])  ;

  
  figure(8)
  plot(t,ss,'b');
  axis([-5,5,-2,2]);
  grid;
  title('發送信號波形（調制后）');
  xlabel('t  (us)')
  ylabel('ss(t)') 
  set(8,'position',[800,250,300,200])  ;

  
  figure(9)
  plot(f,abs(t2f(ss)),'b');
  axis([-8,8,0,40]);
  grid;
  title('發送信號頻譜（調制后）');
  xlabel('f (MHz)')
  ylabel('SS(f)') ;
  set(9,'position',[10,500,300,200])  ;

  
  figure(10)
  plot(t,rey,'b');
  axis([-5,5,-2.5,2.5]);
  grid;
  title('接收信號波形');
  xlabel('t  (us)')
  ylabel('y(t)') 
  set(10,'position',[250,500,300,200])  ;

  figure(11)
    plot(t,anrz,'g');
    axis([-10,10,-1.5,1.5]);
   xlabel('t(us)');
ylabel('a(V)') ;
title('輸入碼形');
  set(11,'position',[550,500,300,200])  ;

figure(12)
    plot(t,aanrz,'g');
    axis([-10,10,-1.5,1.5]);
   xlabel('t(us)');
ylabel('aa(V)') ;
title('輸出端判決得到的碼形');
  set(12,'position',[800,500,300,200])  ;