global dt t f df N T
close all
clear Eb_N0 Pe

N=2^10; %采樣點數
L=32;    %每碼元的采樣點數
M=N/L   %碼元數
Rb=2;		%碼速率是2Mb/s
Ts=1/Rb;        %碼元間隔
dt=Ts/L;        %時域采樣間隔
df=1/(N*dt)     %頻域采樣間隔
T=N*dt          %截短時間
Bs=N*df/2       %系統帶寬
alpha=0.5 	%滾降系數
fc=15
Na=4;           %示波器掃描寬度為4個碼元
t=linspace(-T/2,T/2,N); %時域橫坐標
f=linspace(-Bs,Bs,N)+eps;   %頻域橫坐標

carry1=sqrt(2)*cos(2*pi*fc*t);  %同相載波
carry2=-sqrt(2)*sin(2*pi*fc*t); %正交載波

figure(1)
set(1,'Position',[50,100,400,600])

figure(2)
set(2,'Position',[500,100,400,600])
figure(4)
set(4,'Position',[500,100,400,600])

%升余弦
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;

for loop1=1:10
    Eb_N0(loop1)=(loop1-1)  %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;     %誤碼計數
    EP=zeros(size(f))+eps;
    EP1=zeros(size(f))+eps;
    for loop2=1:3            %循環3次求Pe取平均值
        a=round(rand(1,M));
        b1=a(1:2:M);         %b1為同相支路序列，b2為正交支路序列
        c1=zeros(1,0.25*M);
        for i=1:0.25*M
            if b1(2*i-1)==1&b1(2*i)==0
                c1(i)=1;
            elseif b1(2*i-1)==1&b1(2*i)==1
                c1(i)=3;
            elseif b1(2*i-1)==0&b1(2*i)==0
                c1(i)=-1;
            elseif b1(2*i-1)==0&b1(2*i)==1
                c1(i)=-3;
            end
        end
        b2=a(2:2:M);
        c2=zeros(1,0.25*M);
        for j=1:0.25*M
            if b2(2*j-1)==1&b2(2*j)==0
                c2(j)=1;
            elseif b2(2*j-1)==1&b2(2*j)==1
                c2(j)=3;
            elseif b2(2*j-1)==0&b2(2*j)==0
                c2(j)=-1;
            elseif b2(2*j-1)==0&b2(2*j)==1
                c2(j)=-3;
            end
        end
        
        imp1=zeros(1,N);    %產生沖激序列(同相支路)imp1
        imp1(2*L:4*L:N)=c1/dt;
        imp2=zeros(1,N);    %產生沖激序列(正交支路)imp2
        imp2(2*L:4*L:N)=c2/dt;
        IMP1=t2f(imp1);
        IMP2=t2f(imp2);
        si=real(f2t(IMP1.*GT));
        S=t2f(si);
        sq=real(f2t(IMP2.*GT));
        ssi=si.*carry1;   %發送信號同相支路
        ssq=sq.*carry2;   %發送信號正交支路
        ss=ssi+ssq;
        n0=sqrt(sita)*randn(size(t)); %信道噪聲
        d=ss+n0;       %接收信號
        r1=d.*carry1;
        r2=d.*carry2;
        Y1=t2f(r1).*GR;
        Y2=t2f(r2).*GR;
        y1=real(f2t(Y1)); 
        y2=real(f2t(Y2)); 
        aa1=y1(2*L:4*L:N);       %同相支路取樣
        aaa1=zeros(1,0.5*M);
        for ii=1:0.25*M
            if aa1(ii)>=2
                aaa1(2*ii-1)=1;aaa1(2*ii)=1;
            elseif aa1(ii)>=0
                aaa1(2*ii-1)=1;aaa1(2*ii)=0;
            elseif aa1(ii)>=-2
                aaa1(2*ii-1)=0;aaa1(2*ii)=0;
            elseif aa1(ii)<-2
                aaa1(2*ii-1)=0;aaa1(2*ii)=1;
            end 
        end
        aa2=y2(2*L:4*L:N);       %正交支路取樣
        aaa2=zeros(1,0.5*M);
        for jj=1:0.25*M
            if aa2(jj)>=2
                aaa2(2*jj-1)=1;aaa2(2*jj)=1;
            elseif aa2(jj)>=0
                aaa2(2*jj-1)=1;aaa2(2*jj)=0;
            elseif aa2(jj)>=-2
                aaa2(2*jj-1)=0;aaa2(2*jj)=0;
            elseif aa2(jj)<-2
                aaa2(2*jj-1)=0;aaa2(2*jj)=1;
            end 
        end
        aa=zeros(1,M);
        aa(1:2:M)=aaa1;
        aa(2:2:M)=aaa2;
        P1=S.*conj(S)/T;      %發送根號升余弦信號的功率譜
        EP1=(EP1*(loop2-1)+P1+eps)/loop2;
           
        P=Y1.*conj(Y1)/T;      %升余弦信號的功率譜
        EP=(EP*(loop2-1)+P+eps)/loop2;
        n_err=n_err+length(find(aa~=a))  
    end
    Pe(loop1)=n_err/(M*loop2)      %平均誤碼率＝誤碼總數/（碼元數×循環次數）
    figure(3) 
    semilogy(Eb_N0,Pe,'g');
    eb_n0=10.^(Eb_N0/10);
    hold on
    semilogy(Eb_N0,0.5*erfc(sqrt(eb_n0)));
    axis([0,9,1e-4,1])
    xlabel('Eb/N0')
    ylabel('Pe')
    legend('16QAM誤比特率曲線','2PSK誤比特率曲線')
end

figure(1)							%調制解調信號的波形
subplot(2,1,1)					%發送的調制信號
plot(t/Ts,ssi);
hold on
xlabel('t/Ts');
ylabel('信號幅度');
axis([-5,5,-5,5]);				%只取20個碼元進行觀察
title('16QAM同相支路信號');
plot(t/Ts,si*sqrt(2),'g');
plot(t/Ts,-si*sqrt(2),'g');


subplot(2,1,2)					%接收的解調信號
plot(t/Ts,y1/sqrt(2));
xlabel('t/Ts');
ylabel('信號幅度');
axis([-5,5,-3,3]);				%只取20個碼元進行觀察
title('同相支路解調信號波形）');

figure(2)
subplot(2,1,1);
plot(f,30+10*log10(EP1),'g');
axis([-3,+3,-50,50])
xlabel('f (MHz)')
ylabel('Ps(f)  (dBm/MHz)')
title('調制信號功率譜密度')


subplot(2,1,2);    %畫眼圖
hold on
tt=[0:dt:Na*L*dt];
for jj=1:Na*L:N-Na*L
    plot(tt,si(jj:jj+Na*L)); 
    xlabel('t in us')
    ylabel('y(t) in V')
    title('調制信號眼圖')
end
figure(4)
subplot(2,1,1);
plot(f,30+10*log10(EP),'g');
axis([-3,+3,-50,50])
xlabel('f (MHz)')
ylabel('Ps(f)  (dBm/MHz)')
title('解調信號功率譜密度')


subplot(2,1,2);    %畫眼圖
hold on
tt=[0:dt:Na*L*dt];
for jj=1:Na*L:N-Na*L
    plot(tt,y1(jj:jj+Na*L)); 
    xlabel('t in us')
    ylabel('y(t) in V')
    title('解信號眼圖')
end

figure(5)
subplot(2,1,1)
stem(c1,c2,'filled');
grid on
axis([-5,5,-5,5]);

subplot(2,1,2)
stem(aa1,aa2,'filled');
grid on
axis([-5,5,-5,5]);