%QPSK simulation with Gray coding
%Run from editor debug(F5)
%JC 2/16/07
%The purpose of this m-file is to show a baseband simulated version of QPSK with
%Gray coding which may give valid results (still trying to figure out if it is correct) 
%when compared to theoritical analysis.
%The simulation assumes a perfect system. I have kluged this together from
%various programs that I have seen on the internet and hope it may be
%somewhat useful to others to play with. I have provided comments and notes for review. 
clear all;
close all;

Fd=44;          %符号率
Fs=88;          %采样率
Td=1/Fd;        %符号间隔
Ts=1/Fs;        %采样间隔Ts=T/2=Ti;
%Ti=Td;          %内插滤波周期

coeff_5 = firrcos(28,22,0.25,88,'rolloff','sqrt');
dc_rig = zeros(1,29);       %寄存器用来进行接收端根升余弦滤波
c1 = 1/128;         %Ki
c2 = 1/(2048+128);  %Kp
c3 = 1/2048;        % 针对-0.02情况。鉴相器的鉴相增益系数和数控振荡器的压控灵敏度的乘积
c_1 = inline('1/6*u^3-1/6*u');
c_2 = inline('-1/2*u^3+1/2*u^2+u');
c_3 = inline('1/2*u^3-u^2-1/2*u+1');
c_4 = inline('-1/6*u^3+1/2*u^2-1/3*u');

% num = [0.01 1];
% den = [1 1.01 1];
% [Af,Bf,Cf,Df] = tf2ss(num,den);

%randn('state',0);%keeps bits the same on reruns
nr_data_bits=320000;% 0's and 1's, keep even number-Takes ~1 minute for a run of 1 million
%64000 allows bits and complex values to be shown in array editor
nr_symbols=nr_data_bits/2;
b_data = (randn(1, nr_data_bits) > .5);%random 0's and 1's
b = [b_data];
% Map the bits to be transmitted into QPSK symbols using Gray coding. The
% resulting QPSK symbol is complex-valued, where one of the two bits in each
% QPSK symbol affects the real part (I channel) of the symbol and the other
% bit the imaginary part (Q channel). Each part is subsequently
% modulated to form the complex-valued QPSK symbol.
%
% The Gray mapping resulting from the two branches are shown where
% one symbol error corresponds to one bit error going counterclockwise.
% imaginary part (Q channel)
%         ^
%         |
%  10 x   |   x 00   (odd bit, even bit)
%         |
%  -------+------->  real part (I channel)
%         |
%  11 x   |   x 01
%         |
% Input:
%   b = bits {0, 1} to be mapped into QPSK symbols
%
% Output:
%   d = complex-valued QPSK symbols  0.70711 + 0.70711i, etc


d=zeros(1,length(b)/2);
%definition of the QPSK symbols using Gray coding.
for n=1:length(b)/2
    p=b(2*n);
    imp=b(2*n-1);
    if (imp==0)&(p==0)
        d(n)=exp(j*pi/4);%45 degrees
    end
    if (imp==1)&(p==0)
        d(n)=exp(j*3*pi/4);%135 degrees
    end
    if (imp==1)&(p==1)
        d(n)=exp(j*5*pi/4);%225 degrees
    end
    if (imp==0)&(p==1)
        d(n)=exp(j*7*pi/4);%315 degrees
    end
end


pn_15 = [1 1 1 1 0 1 0 1 1 0 0 1 0 0 0];
pn_15_c=pn_15+1i*pn_15;
qpsk=cat(2,pn_15_c,d);

pn_regist = zeros(1,15);    %帧前导相关寄存器

%插零上采样,成型滤波
Ep_insert_sig = zeros(1,length(qpsk)*2);  % insert 8 zero
Ep_insert_sig(:,1:2:end) = qpsk;
Filter_sig(1,:) = conv(Ep_insert_sig(1,:),coeff_5);         % low_filter (fc = 11,fs =88);低通滤波


%加噪声
SNR=5;
sigma = sqrt(10.0^(-SNR/10.0));
sigma=sigma/2;  
sig_tx_channel=(real(Filter_sig)+sigma.*randn(size(Filter_sig))) +1i.*(imag(Filter_sig)+sigma*randn(size(Filter_sig)));

%重采样
y_rx_sig = resample(sig_tx_channel,8798,8800);


y_rx_rrc = conv(y_rx_sig,coeff_5);
%进行定时同步

%这里Ti=T
%NCO操作的时钟为Ts
%内插操作的时钟为Ti
%NCO控制字更新的时钟为Ts?Ti?T
%定时差错检测的时钟为T,每个符号周期计算一次
    
i_idx=1;           %当前内插位置                 采样时钟Ts
j_idx=1;           %当前输出符号索引              选通时钟T=2*Ti=4*Ts
k_idx=1;           %内插滤波器输出索引            内插时钟Ti=2*Ts  
%uk=0;                 %初始化小数延迟
%strobe=0;             %用于判决当前采样点为最佳采样或是中间采样    
tnco_out(i_idx)=0.9;         %NCO寄存器的值
Wn(j_idx)=0.5;             %NCO控制字的值Wn通过重复T/Ts即可获得Wm

    
ut=[];                 %小数采样延迟输出
symbol_gad=[];         %两倍符号内插序列输出
time_error=[];         %定时偏差输出
symbols_recovery=[];
    
%omega=0.0005;  %环路滤波归一化带宽，fs归一化为1
                          %调制方式                    QPSK           16QAM
c1=6.9151e-4;         %比例环路滤波增益          6.9151e-4      1.14534e-4     
c2=2.3474e-7;         %积分环路滤波增益          2.3474e-7      4.04e-8   
    
while i_idx<length(y_rx_rrc)-2
        
      %NCO更新
      n_temp=tnco_out(i_idx)-Wn(j_idx);
      tnco_out(i_idx+1)=mod(n_temp,1);
        
      if n_temp<-10^-10
           
          strobe=mod(k_idx,2);            
          ut(k_idx)=2* tnco_out(i_idx);
            
          % 进行内插滤波
          % 当tnco_out(mk+1)<0时，说明mk为内插点，需要将x(mk-1),x(mk),x(mk+1),x(mk+2)进行滤波            
          symbol_gad(k_idx)=farrow(y_rx_rrc(i_idx-1:i_idx+2),ut(k_idx));
            
          if strobe==0              
              symbols_recovery(j_idx)=symbol_gad(k_idx);
              if k_idx>2
                time_error(j_idx)=real(symbol_gad(k_idx-1)*(symbol_gad(k_idx)-symbol_gad(k_idx-2)));
              else
                time_error(j_idx)=real(symbol_gad(k_idx-1)*(symbol_gad(k_idx)));   
              end
              %更新NCO控制字
              if j_idx>1
                Wn(j_idx+1)=Wn(j_idx)+c1*(time_error(j_idx)-time_error(j_idx-1))+c2*time_error(j_idx);
              else
                Wn(j_idx+1)=Wn(j_idx)+c1*(time_error(j_idx))+c2*time_error(j_idx);
              end
              %Wn(j_idx)=Wn(j_idx-1)+(Kp+Ki)*time_error(j_idx)+Ki*time_error(j_idx-1);
              
              %符号选通，应该输出符号                
              j_idx=j_idx+1;
          end
          
          %内插滤波输出符号前移
          k_idx=k_idx+1;  
      end
      
      %采样位置前移
      i_idx=i_idx+1;
end    
disp(['Time recovery is finished.']);

% frame_sync_flag=0;
% bits_recv=[];
% for rdx=1:length(symbols_recovery)
%     if frame_sync_flag==1
%         %帧已经同步上
%         r=symbols_recovery(rdx);
%         bhat=[real(r)<0;imag(r)<0];%detector 
%         bits_recv=cat(2,bits_recv,bhat);
%     else
%         pn_regist(:,1:end-1) = pn_regist(:,2:end);
%         pn_regist(:,end) = symbols_recovery(rdx);
%         pn_mul_result = real(pn_regist(1,:)*pn_15_c');
%         if pn_mul_result>8
%             frame_sync_flag=1;
%         end
%     end
% end

%比较误比特率
