gardner 算法matlab實現

% 仿真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));