%16QAM系统的调制、解调的仿真，信源为单极性不归零码(NRZ)

clear all  %清空工作空间的变量

close all  %关闭所有图形窗口

global dt df t f N T %定义全局变量

k=13;

N=2^k;    %总的取样点数

L=32;                %L为每个码元的取样点数

M=N/L;                %M码元总数

Rb=3;                %码元速率Rb=2Mb/s

Ts=1/Rb;                %码元宽度Ts

dt=Ts/L;        %时域取样间隔

df=1/(N*dt);%频域取样间隔

T=N*dt;     %时间截断长度

Bs=N*df/2;  %带宽

Na=4;      %示波器扫描宽度

alpha=0.5;   %升余弦滚降系数

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);

EP=zeros(size(f));     %输入信号功率谱

EPa=zeros(size(f));    %调制后信号功率谱

EPo=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=0.5;

n0=Eb/(eb_n0(loop1)); %信道的噪声谱密度

sita=n0*Bs; %信道中噪声功率

n_err=0;     %误码计数

for loop2=1:5

a=round(rand(1,M));   %产生随机码

for i0=1:L;

s(i0+[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,

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);           %调制后信号的傅氏变换

Pp=abs(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));

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);

Po=Se.*conj(Se)/T;

EPo=(EPo*(loop1-1)+Po)/loop1; %输出信号功率谱

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(EPo+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('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)

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('调制后信号功率谱密度图');

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

title('I路解调并滤波后信号波形眼图');

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')

title('误码率曲线')

figure(4)

plot(I,Q,'.');

axis([-4,4,-4,4])

xlabel('I ');

ylabel('Q ');

title('16QAM星座图');

figure(5)

plot(f,Pp);

xlabel('w');

ylabel('V/HZ');

title('已调制信号频谱图');