% 仿真4比特原始数据与星座图的编码映射过程；

% 完成16QAM信号的调制解调；

% 基带信号符号速率 ps ＝1Mbps；

% 成形滤波器的滚降因子 a=0.8；

% 载波信号频率fc=2MHz ；

% 采样频率 Fs=8MHz ；

% 绘制16QAM信号的频谱及时域波形；

% 采用相干解调法仿真其解调过程；

% 绘制解调前后的基带信号时域波形；

% 将原始基带数据、QAM已调数据、滤波器系数写入相应的文本文件中。

clc;

close all;

ps=1*10^6;   %码速率为1MHz

a=0.8;       %成形滤波器系数

Fs=8*10^6;   %采样速率

fc=2*10^6;   %载波频率

N=4000;      %仿真数据的长度

t=0:1/Fs:(N*Fs/ps-1)/Fs;%产生长度为N,频率为fs的时间序列

s=randint(N,1,16);       %产生随机16进制数据作为原始数据

Bs=dec2bin(s,4);         %将十进制数据转换成4比特的二进制数据

%对Bs的高2比特进行差分编码

%取高2比特分别存放在A，B变量中

A=s>7;

B=(s-A*8)>3;

%差分编码后的数据存放在C，D中

C=zeros(N,1);D=zeros(N,1);

for i=2:N

    C(i)=mod(((~mod(A(i)+B(i),2))&mod(A(i)+C(i-1),2)) + (mod(A(i)+B(i),2)&mod(A(i)+D(i-1),2)),2);

    D(i)=mod(((~mod(A(i)+B(i),2))&mod(B(i)+D(i-1),2)) + (mod(A(i)+B(i),2)&mod(B(i)+C(i-1),2)),2);

end

%差分编码后的高2比特数据与原数据低2比特合成映射前的数据DBs

DBs=C*8+D*4+s-A*8-B*4;

%完成调制前的正交同相支路数据映射

I=zeros(1,N);Q=zeros(1,N);

for i=1:N

    switch DBs(i)

        case 0, I(i)=3; Q(i)=3;

        case 1, I(i)=1; Q(i)=3;

        case 2, I(i)=3; Q(i)=1;

        case 3, I(i)=1; Q(i)=1;

        case 4, I(i)=-3;Q(i)=3;

        case 5, I(i)=-3;Q(i)=1;

        case 6, I(i)=-1;Q(i)=3;

        case 7, I(i)=-1;Q(i)=1;

        case 8, I(i)=3; Q(i)=-3;

        case 9, I(i)=3; Q(i)=-1;

        case 10,I(i)=1; Q(i)=-3;

        case 11,I(i)=1; Q(i)=-1;

        case 12,I(i)=-3;Q(i)=-3;

        case 13,I(i)=-1;Q(i)=-3;

        case 14,I(i)=-3;Q(i)=-1;

        otherwise,I(i)=-1;Q(i)=-1;

    end

end

%对编码数据以Fs频率采样

Ads_i=upsample(I,Fs/ps);

Ads_q=upsample(Q,Fs/ps);

%加噪声

% SNR=30;

% Ads_i=awgn(Ads_i,SNR);

% Ads_q=awgn(Ads_q,SNR);

%设计平方根升余弦滤波器

n_T=[-2 2];

rate=Fs/ps;

T=1;

Shape_b = rcosfir(a,n_T,rate,T,'sqrt');

%对采样后的数据进行升余弦滤波;

rcos_Ads_i=filter(Shape_b,1,Ads_i);

rcos_Ads_q=filter(Shape_b,1,Ads_q);

%产生同相正交两路载频信号

f0_i=cos(2*pi*fc*t);

f0_q=sin(2*pi*fc*t);       

%产生16QAM已调信号

qam16=rcos_Ads_i.*f0_i+rcos_Ads_q.*f0_q;

figure;

plot(qam16);

srdata = rcos_Ads_i + rcos_Ads_q * 1i;

scatterplot(srdata(length(srdata)*0.9:2:length(srdata)));

%%%%%　仿真输入测试的PSK基带数据　%%%

bitstream=randint(1,N,2);

psk2=pskmod(bitstream,2);

xI=zeros(1,Ns);

xQ=zeros(1,Ns);

xI(1:8:8*N)=real(psk2);%8倍插值

xQ(1:8:8*N)=imag(psk2);

%截短后的根升余弦匹配滤波器

h1=rcosfir(0.8,[-8,8],4,1,'sqrt');

hw=kaiser(65,3.97);

hh=h1.*hw.';

aI1=conv(xI,hh);

bQ1=conv(xQ,hh);

L=length(aI1);

%仿真输入数据

% aI=[aI1(22:2:L),0,0];%2倍抽取

% bQ=[bQ1(22:2:L),0,0];

%

% scatterplot(psk2(length(psk2)*0.5:length(psk2)));

aI=real(srdata(202:2:length(srdata)));%2倍抽取？为何先8倍插值，再2倍抽取？

bQ=imag(srdata(202:2:length(srdata)));

ma=max(abs(aI));mb=max(abs(bQ));

m=max(ma,mb);

aI=aI/m;bQ=bQ/m;

%

N=floor(length(aI)/4);

Ns=4*N;  %总的采样点数

bt=0.001;

c2=2^(-14);

c1=2^(-6);

% c1=8/3*bt;

% c2=32/9*bt*bt;

i=3;    %用来表示Ts的时间序号,指示n,n_temp,nco,

w=[0.5,zeros(1,N-1)];  %环路滤波器输出寄存器，初值设为0.5

n=[0.7 zeros(1,Ns-1)]; %NCO寄存器，初值设为0.9

n_temp=[n(1),zeros(1,Ns-1)];

u=[0.6,zeros(1,2*N-1)];%NCO输出的定时分数间隔寄存器，初值设为0.6

yI=zeros(1,2*N);       %I路内插后的输出数据

yQ=zeros(1,2*N);       %Q路内插后的输出数据

time_error=zeros(1,N); %Gardner提取的时钟误差寄存器

ik=time_error;

qk=time_error;

k=1;    %用来表示Ti时间序号,指示u,yI,yQ

ms=1;   %用来指示T的时间序号,用来指示a,b以及w

strobe=zeros(1,Ns);

ns=length(aI)-2;

while(i<ns)

    n_temp(i+1)=n(i)-w(ms);

%      n_temp(i+1)=n(i)-0.5;

    if(n_temp(i+1)>0)

        n(i+1)=n_temp(i+1);

    else

        n(i+1)=mod(n_temp(i+1),1);

        %内插滤波器模块

        FI1=0.5*aI(i+2)-0.5*aI(i+1)-0.5*aI(i)+0.5*aI(i-1);

        FI2=1.5*aI(i+1)-0.5*aI(i+2)-0.5*aI(i)-0.5*aI(i-1);

        FI3=aI(i);

        yI(k)=(FI1*u(k)+FI2)*u(k)+FI3;

        FQ1=0.5*bQ(i+2)-0.5*bQ(i+1)-0.5*bQ(i)+0.5*bQ(i-1);

        FQ2=1.5*bQ(i+1)-0.5*bQ(i+2)-0.5*bQ(i)-0.5*bQ(i-1);

        FQ3=bQ(i);

        yQ(k)=(FQ1*u(k)+FQ2)*u(k)+FQ3;

        strobe(k)=mod(k,2);

        %时钟误差提取模块，采用的是Gardner算法

        if(strobe(k)==0)

            %取出插值数据

            ik(ms)=yI(k);

            qk(ms)=yQ(k);

            %每个数据符号计算一次时钟误差

            if(k>2)

               Ia=(yI(k)+yI(k-2))/2;

               Qa=(yQ(k)+yQ(k-2))/2;

               time_error(ms)=[yI(k-1)-Ia ]  *(yI(k)-yI(k-2))+[yQ(k-1)-Qa]  *(yQ(k)-yQ(k-2));

            else

                time_error(ms)=(yI(k-1)*yI(k)+yQ(k-1)*yQ(k));

            end

            %环路滤波器,每个数据符号计算一次环路滤波器输出

            if(ms>1)

                w(ms+1)=w(ms)+c1*(time_error(ms)-time_error(ms-1))+c2*time_error(ms-1);

                %w(ms+1)=w(ms)+c1*(time_error(ms)-time_error(ms-1));

            else

                w(ms+1)=w(ms)+c1*time_error(ms)+c2*time_error(ms);

            end

            ms=ms+1;

        end

          k=k+1;

%         u(k)=n(i)/w(ms);

          u(k) = 0.5* n(i);%%近视运算

    end

    i=i+1;

end

%

figure;

subplot(311);plot(u);xlabel('运算点数');ylabel('分数间隔');

subplot(312);plot(time_error);xlabel('运算点数');ylabel('定时误差');

subplot(313);plot(w);xlabel('运算点数');ylabel('环路滤波器输出');

%

%

%

iq=ik+qk*sqrt(-1);

L=length(iq)

% off=6;

scatterplot(iq(1:end));

scatterplot(iq(L*0.5:end));

scatterplot(iq(L*0.9:end));