%信源為單極性不歸零碼（NRZ），取樣點數(shù)自己定義.......16QAM


global dt df t f N T %定義全局變量
close all  %關閉以前的應用窗口
clear Eb_N0 Pe

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

N=2^k;    %總的取樣點數(shù)
L=32;		%L為每個碼元的取樣點數(shù)
M=N/L;		%M碼元總數(shù)	
Rb=2;		%碼元速率Rb=2Mb/s
Ts=1/Rb;		%碼元寬度Ts
dt=Ts/L;	%時域取樣間隔
df=1/(N*dt);%頻域取樣間隔
T=N*dt;     %時間截斷長度
Bs=N*df/2;  %帶寬
Na=2;      %示波器掃描寬度為2個碼元
alpha=0;   %升余弦滾降系數(shù)

t=linspace(-T/2,T/2,N); %頻域橫坐標
f=linspace(-Bs,Bs,N)+eps;%時域橫坐標

%升余弦
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])
                %設定窗口位置及大小 
figure(2)
 set(2,'Position',[515,250,500,400])
                %設定窗口位置及大小 
                
figure(3)
 set(3,'Position',[50,10,500,200])
                %設定窗口位置及大小 

figure(4)
 set(4,'Position',[600,10,400,300])
                %設定窗口位置及大小                 
A=1;
fc=1/Ts;
m=A*cos(2*pi*fc*t);
n=-A*sin(2*pi*fc*t);      % lianggezaibo
EP=zeros(size(f));    
EPa=zeros(size(f)); 
 
for loop1=1:15    %誤碼率曲線
 Eb_N0(loop1)=2*(loop1-1);  %Eb/N0 in dB
 %loop1=1;
 %Eb_N0(loop1)=30;  %Eb/N0 in dB
 eb_n0(loop1)=10^(Eb_N0(loop1)/10);
 Eb=1;	
 n0=Eb/(eb_n0(loop1)); %信道的噪聲譜密度   
 sita=n0*Bs; %信道中噪聲功率
 n_err=0;     %誤碼計數(shù)

 for loop2=1:5  
 a=round(rand(1,M));   %產生隨機碼
 
for l3=1:L;  
s(l3+[0:M-1]*L)=a;   %將輸出信號變成雙極性不歸零碼
end
S=t2f(s);
P=S.*conj(S)/Ts;
EP=(EP*(loop1-1)+P)/loop1; %發(fā)送功率譜


 i=zeros(1,M/2);         %串并轉換
 q=zeros(1,M/2);  
 I=zeros(1,N);         %串并轉換
 Q=zeros(1,N); 
 
 i=a(1:2:M-1);
 q=a(2:2:M);
 ii=zeros(1,M/4);         %串并轉換
 qq=zeros(1,M/4); 
 yii=zeros(1,M/2);         %串并轉換
 yqq=zeros(1,M/2); 
 y=zeros(1,M); 
 for l1=1:2:(M/2-1);
        if i(l1)-i(l1+1)==1,
             ii((l1+1)/2)=1;
        elseif  i(l1)-i(l1+1)==-1
                ii((l1+1)/2)=-1;
        elseif  i(l1)==0
               ii((l1+1)/2)=-3;
           else  ii((l1+1)/2)=3;
        end
     
       if q(l1)-q(l1+1)==1,
             qq((l1+1)/2)=1;
        elseif  q(l1)-q(l1+1)==-1
                qq((l1+1)/2)=-1;
        elseif  q(l1)==0
               qq((l1+1)/2)=-3;
           else  qq((l1+1)/2)=3;
        end
 end

for l2=1:4*Ts/dt;   
   I(l2+[0:M/4-1]*4*L)=ii;
   Q(l2+[0:M/4-1]*4*L)=qq;
end;

Ia=I.*m;    %兩路分別調制
Qa=Q.*n;

st=Ia+Qa;  %并串轉換--矩形星座
ST=t2f(st);           %信源信號的傅式變換
Pa=ST.*conj(ST)/T;     
EPa=(EP*(loop1-1)+Pa)/loop1; %發(fā)送功率譜

n_ch=sqrt(sita)*randn(size(t));  %信道噪聲

st=st+n_ch; %加噪聲

Iaa=st.*m;  %對兩路信號分別解調
Qaa=st.*n;

Ir=real(f2t(t2f(Iaa).*GR)); %分別通過濾波器
Qr=real(f2t(t2f(Qaa).*GR)); 

yi=Ir(2*L:4*L:N);  %抽樣判決
yq=Qr(2*L:4*L:N);

for l2=1:2:M/2-1
    k1=abs(yi((l2+1)/2)+1);
    k2=abs(yi((l2+1)/2)+0.5);
    k3=abs(yi((l2+1)/2)-0.5);
    k4=abs(yi((l2+1)/2)-1);
    k=[k1,k2,k3,k4];
    
    p1=abs(yq((l2+1)/2)+1);
    p2=abs(yq((l2+1)/2)+0.5);
    p3=abs(yq((l2+1)/2)-0.5);
    p4=abs(yq((l2+1)/2)-1);
    p=[p1,p2,p3,p4];
    
    if k1==min(k)
        yii(l2)=0;yii(l2+1)=0;
    elseif k2==min(k)
        yii(l2)=0;yii(l2+1)=1;
    elseif k3==min(k)
        yii(l2)=1;yii(l2+1)=0;
    elseif k4==min(k)
        yii(l2)=1;yii(l2+1)=1;
    end
    
    
     if p1==min(p)
        yqq(l2)=0;yqq(l2+1)=0;
     elseif p2==min(p)
        yqq(l2)=0;yqq(l2+1)=1;
     elseif p3==min(p)
        yqq(l2)=1;yqq(l2+1)=0;
     elseif p4==min(p)
         yqq(l2)=1;yqq(l2+1)=1;
     end
    
 end
 y(1:2:M-1)=yii;  %并串轉換
 y(2:2:M)=yqq;
 
for l4=1:L;  
e(l4+[0:M-1]*L)=y;   %將輸出信號變成雙極性不歸零碼
end  
Se=t2f(e);
Pa=Se.*conj(Se)/T;
EPa=(EPa*(loop1-1)+Pa)/loop1; %經調制后的PSK信號功率譜
n_err=n_err+length(find(y~=a));  
end
Pe(loop1)=n_err/(M*loop2);

end

figure(1)
 subplot(2,2,1)
 plot(t,s,'r');
 axis([-T/64,T/64,min(s)-0.5,max(s)+0.5])
 xlabel('t (ms)');
 ylabel('s (V/KHz)');
 title('信源信號波形');

 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('P(f) (dBm/MHz)');
 title('信源信號功率譜密度圖形');
 
 subplot(2,2,3)
 plot(t,e,'b');
 axis([-T/64,T/64,min(e)-0.5,max(e)+0.5])
 xlabel('t (ms)');
 ylabel('e  (V/KHz)');
 title('輸出信號濾波波形');
 
 subplot(2,2,4)
 bb=30+10*log10(EPa+eps);
 plot(f,bb,'g')
 axis([-20,20,-20,max(bb)+0.5])
 grid
 xlabel('f(Mhz)');
 ylabel('Po(f) (dBm/MHz)');
 title('輸出信號功率譜密度圖形');
 
 
 figure(2)
 subplot(2,3,1)
 plot(t,I,'r');
 axis([-T/16,T/16,min(I)-0.5,max(I)+0.5])
 xlabel('t (ms)');
 ylabel('I  (V/KHz)');
 title('I路信號波形');
 
 subplot(2,3,2)
 plot(t,st,'b');
 axis([-T/16,T/16,min(st)-0.5,max(st)+0.5])
 xlabel('t (ms)');
 ylabel('st  (V/KHz)');
 title('調制信號波形');
 
 subplot(2,3,3)
 plot(t,Iaa,'b');
 axis([-T/16,T/16,min(Iaa)-0.5,max(Iaa)+0.5])
 xlabel('t (ms)');
 ylabel('Iaa  (V/KHz)');
 title('解調信號波形');
 
 subplot(2,3,4)
 plot(t,Q,'b');
 axis([-T/16,T/16,min(Q)-0.5,max(Q)+0.5])
 xlabel('t (ms)');
 ylabel('Q  (V/KHz)');
 title('Q路信號波形');
 
 subplot(2,3,5)
 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('QAM信號功率譜密度圖形');
 
 
 subplot(2,3,6)
 plot(t,Ir,'b');
 axis([-T/16,T/16,min(Ir)-0.5,max(Ir)+0.5])
 xlabel('t (ms)');
 ylabel('Ir  (V/KHz)');
 title('I路解調信號濾波波形');
  
 
  
 figure(3)
 subplot(1,2,1)
  tt=[0:dt:Na*L*dt];
  hold on
  for jj=1:Na*L:N-Na*L
    plot(tt,Ir(jj:jj+Na*L)); %接收信號眼圖
  end
 
 subplot(1,2,2) 
  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,14,1e-4,1])
  xlabel('Eb/N0')
  ylabel('Pe')
  
  
 figure(4)
 plot(I,Q,'.');
 axis([-4,4,-4,4])
 xlabel('I ');
 ylabel('Q ');
 title('16QAM星座圖');