%16QAM：信源為單極性不歸零碼（NRZ），取樣點數自己定義
%信源為單極性不歸零碼（NRZ），取樣點數自己定
義.......16QAM
global dt df t f N T %定義全局變量
close all %關閉以前的應用窗口
clear Eb_N0 Pe
k=input('取樣點數＝2^k, k=[13]'); %采樣點數自定
義，默認為2^13
if k==[], k=13; end
N=2^k; %總的取樣點數
L=32; %L 為每個碼元的取樣點數
M=N/L; %M 碼元總數
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; %升余弦滾降系數
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)
通信原理（軟件）實驗報告 2005.4 page 26 of 37
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); %載波
EP=zeros(size(f));
EPa=zeros(size(f));
for loop1=1:15 %誤碼率曲線
Eb_N0(loop1)=2*(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; %誤碼計數
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; %發送功率譜
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,%10=>1
ii((l1+1)/2)=1;
elseif i(l1)-i(l1+1)==-1 %01=-1
ii((l1+1)/2)=-1;
elseif i(l1)==0 %00=-3
ii((l1+1)/2)=-3;
else ii((l1+1)/2)=3; %11=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; %s1
Q(l2+[0:M/4-1]*4*L)=qq; %s2
end;
Ia=I.*m; %兩路分別調制s1*f1
Qa=Q.*n; %s2*f2
st=Ia+Qa; %并串轉換--矩形星座
ST=t2f(st); %信源信號的傅式變換
Pa=ST.*conj(ST)/T;
EPa=(EP*(loop1-1)+Pa)/loop1; %發送功率譜
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)); %幅度變成+-1，+-0。5
yi=Ir(2*L:4*L:N); %抽樣判決
yq=Qr(2*L:4*L:N);%以下部分參考別人的了~~因
為一開始想取樣然后round 一下，可是后來發現不
通信原理（軟件）實驗報告 2005.4 page 27 of 37
行
for l2=1:2:M/2-1 %1=10,0.5=11,-0.5=01,-1=00
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; %經調制后的信號
功率譜
n_err=n_err+length(find(y~=a));
end
Pe(loop1)=n_err/(M*loop2);
end
figure(1)
subplot(2,3,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,3,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,3,3)
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,4)
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,3,5)
bb=30+10*log10(EPa+eps);
plot(f,bb,'g')
axis([-20,20,-20,max(bb)+0.5])
grid
通信原理（軟件）實驗報告 2005.4 page 28 of 37
xlabel('f(Mhz)');
ylabel('Po(f) (dBm/MHz)');
title('輸出信號功率譜密度圖形');
subplot(2,3,6)
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 信號功率譜密度圖形');
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,Iaa,'b');
axis([-T/16,T/16,min(Iaa)-0.5,max(Iaa)+0.5])
xlabel('t (ms)');
ylabel('Iaa (V/KHz)');
title('I 路信號解調信號波形');
subplot(2,3,3)
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 路濾波后解調信號濾波波形');
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)
plot(t,Qaa,'b');
axis([-T/16,T/16,min(Qaa)-0.5,max(Qaa)+0.5])
xlabel('t (ms)');
ylabel('Qaa (V/KHz)');
title('Q 路信號解調信號波形');
subplot(2,3,6)
plot(t,Qr,'b');
axis([-T/16,T/16,min(Qr)-0.5,max(Qr)+0.5])
xlabel('t (ms)');
ylabel('Qr (V/KHz)');
title('Q 路濾波后解調信號濾波波形');
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 星座圖');