cpm_sim.asv

cpm调制解调实现,关联长度为3

% clc
clear all
close all

T=1;
L=3;
h=2/7;
m=8;            % 八进制
p = 2/h;        % 累计相位的状态个数
OVS=40;
OVR=40;
EbN0 = 10;
SNR = EbN0 - 10 * log10(OVR) + 10 * log10(log2(m));
SNR = 10;
% SNR = 1;

block_num = 3;      %发送信息的桢数
block_size = 96;     %发送信息的桢大小   block_num*block_size为发送的信息bit数
frm_num = block_num;    
frm_size = block_size + 4;  %每桢加一个前导码和三个归零码

t=[0:T/OVS:L*T];
gt=( 1-cos( 2*pi*t/(L*T) ) )/(2*L*T);

for ii=1:L*OVS
    QT(ii)=sum(gt(1:ii).*(T/OVS));
end
rand('state',1);
source_bit = 2*(randint(1,block_num*block_size,m)) - (m-1);     %生成发送信息

% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% 对信息序列成桢，相当于每一条发的序列，在每一桢加三个symbol使状态归零
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
framed_bit = zeros(1,block_num*(block_size+4));
tail = [0,-(m-1),-(m-1)];
prem = -(m-1);
v = [];
for ii = 1:block_num
    v = source_bit( (ii-1)*block_size+1:ii*block_size );
    tail(1) = mod(p - mod(sum((v+m-1)/2),p),m);
    tail(1) = tail(1)*2 - (m-1);
    framed_bit((ii-1)*frm_size+1:ii*frm_size) = [prem,source_bit((ii-1)*block_size+1:ii*block_size),tail];
end;

% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
%                                       MODULATE
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
init_state = [0,0,0];
phas = zeros(1,length(framed_bit));
acc_ph = zeros(1,frm_num*frm_size*OVS);
mod_bit = [-(m-1),-(m-1),framed_bit];	% init state = (0,-7,-7) 使第一桢从初始状态（phase=0,ak2=-7,ak1=-7）开始调制

for jj = 3:length(mod_bit)
		kk = jj - 2;
	if(kk>1)
		phas(kk) = mod(phas(kk-1) + pi*h*mod_bit(jj-3),2*pi);
	end;
    acc_ph((kk-1)*OVS+1:kk*OVS) = mod(  2*pi*h*( mod_bit(jj-2)*QT(2*OVS+1:3*OVS) + mod_bit(jj-1)*QT(OVS+1:2*OVS) + mod_bit(jj)*QT(1:OVS) ) + phas(kk),2*pi );
end;
mod_sig = exp(j*acc_ph);
sp_sig = mod_sig(1:frm_size*OVS);                       % 已知信息,用来估计频偏
base_sig = mod_sig(1:OVS);    % 已知信息,用来估计没跳得相偏
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
%     PLOT ACC_PH
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
% x = 1./OVS:1./OVS:length(mod_bit)-2;
% plot(x,acc_ph);
% ylim([0,2*pi]);

% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
%											HOPPING   
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 

%为每一跳加上一个随机相位
trans_sig1 = zeros(1,length(mod_sig));
% rand('state',3);
for ii = 1:frm_num
	index1 = (ii-1)*frm_size*OVS + 1;
	index2 = ii*frm_size*OVS;
	trans_sig1(index1:index2) = mod_sig(index1:index2)*exp(j*2*pi*rand(1));
    trans_sig1(index1:index2) = mod_sig(index1:index2)*exp(j*2);
end;
% trans_sig1 = mod_sig;


% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
%											CHANNEL
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
noisy_sig1 = awgn(trans_sig1, SNR, 0);
noisy_sig1 = trans_sig1;    %没有噪声的情况  

% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
% %                                        DEMODULATE
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
    % 晶体频偏带来的误差 dw*t
hop_period = 0.005;
resu = hop_period/(frm_size*OVS);
dw = 100;
t = [resu:resu:hop_period];
fe = zeros(1,length(t));
for i=1:length(t)
    fe(i) = mod(2*pi*dw*t(i), 2*pi);
end;
% plot(t,fe);

rec_sig1 = zeros(1,length(noisy_sig1));
for ii = 1:frm_num
    for jj = 1:frm_size*OVS
        index1 = (ii-1)*frm_size*OVS + jj;
        rec_sig1(index1) = noisy_sig1(index1)*exp(j*fe(jj));
    end;
end;

	% 把同步引导码作为已知信息，估计本地的晶体频偏
fe_sig = zeros(1,frm_size);
for ii = 1: frm_size
	index = (ii-1)*OVR;
	fe_sig(ii) = rec_sig1(index+1:index+OVR) * sp_sig(index+1:index+OVR)';
end;

dfe_sig = fe_sig(2:frm_size) * fe_sig(1:frm_size-1)';
dfe_sig = dfe_sig/abs(dfe_sig);
angle(dfe_sig)
abs( angle(dfe_sig) - 2*pi*dw*hop_period/frm_size)/(2*pi*dw*hop_period/frm_size)
 
dmod_sig1 = zeros(1,length(rec_sig1));

    % 补偿频偏估计 
for ii = 1:frm_num
    acc_dfe_sig = 1;
    for jj = 1:frm_size
        index = ((ii-1)*frm_num+jj-1)*OVS;
        acc_dfe_sig = acc_dfe_sig * dfe_sig;
        dmod_sig1(index+1:index+OVS) = rec_sig1(index+1:index+OVS) * conj(acc_dfe_sig);
    end;
end;

% dmod_sig1 = noisy_sig1;
    %每个hopping第一个符号做相位估计
for ii = 1:frm_num
	index1 = (ii-1)*frm_size*OVR + 1;
	index2 = ii*frm_size*OVR;
	offset_sig = base_sig*rec_sig1(index1:index1+OVR-1)';		%'
    offset_sig = offset_sig * conj(
    angle(offset_sig)
    dmod_sig1(index1:index2) = dmod_sig1(index1:index2) * offset_sig;
end;

% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
% state_from是记录状态转移的矩阵 correlator是记录各种状态下的相关器，跟m,h,L,OVR有关，都是多维矩阵，前三个index分别是(phase,ak2,ak1)
% 两个load分别load m=8,h=2/7,OVR=8的correlator和state_from文件，如果改变了OVR，需要使用gen_mtchd_fltr重新计算correlator，一般如果
% 没有改变m和h，不要重新计算state_from，否则比较慢。
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
load state_from state_from;             % state_from save how the states transfer
load correlator40 correlator;             % correlator h = 2/7 m = 8 Ns = 8
% correlator = gen_mtchd_fltr(OVR, m, h, QT);      % generate mtchd_fltr(m,m,m,Ns)
% save correlator32 correlator
% [state_from_input, state_from, to_state_output, to_state, phase_state] = gen_trellis(m, h, L);     % generate state trans table()

rec_sym_num = length(dmod_sig1)/OVR;
if( frm_num*frm_size ~= rec_sym_num)
    disp('transmit or receive error!');
end;

branch_metric = zeros(p,m,m,m);								%每个状态（p,m,m）出去8条路径(m)
acc_metric_left = zeros(p,m,m);								%累计的metric
acc_metric_right = zeros(p,m,m);
surv_path = ones(rec_sym_num/frm_num,p,m,m);					%记录winning route
sub_acc_every_state = zeros(1,m);							%一个中间保存metric的存储空间
dmod_bit = zeros(1,rec_sym_num);							%调制出来的信息，其中每桢还包含了三个归零码

for frm = 1:frm_num
        acc_metric_left = zeros(p,m,m);
        acc_metric_left(1,1,1) = OVR*2;                             % 把初始状态的metric设置为最大，表示开始符号肯定是从这个状态出发
		for sym = 1:rec_sym_num/frm_num							% for each symbol
    			% left trevills states	%对每个状态的m条出发路径分别计算metric	
		    for phk = 1:p							% for each state 
		    for ak2 = 1:m							% for each state 
		    for ak1 = 1:m							% for each state 
		        for ak = 1:m			% for each incoming branch
                    temp = dmod_sig1((frm-1)*frm_size*OVR+(sym-1)*OVR+1 : (frm-1)*frm_size*OVR+sym*OVR ) * conj( reshape(correlator(ak2,ak1,ak,:),OVR,1) ); 		% 复数相关
		            branch_metric(phk,ak2,ak1,ak) = real(temp * exp(-j*pi*h*(phk-1)));								% 得到Zk
%                     sprintf('phk=%d ak2=%d ak1=%d ak=%d,val=%f',phk,ak2,ak1,ak,real(temp * exp(-j*pi*h*(phk-1))))
                end;
		    end;end;end;
			 	
			 	% right trevills states 			%对每个状态的m条到达路径选最大的metric，并累计acc_metric_right上
		    for phk = 1:p									% for each state
		    for ak2 = 1:m									% for each state
		    for ak1 = 1:m									% for each state
		
		        for ii = 1:m									% for each incoming branch
		            v = reshape( state_from(phk,ak2,ak1,ii,:),1,3);						% which state from
		            sub_acc_every_state(ii) = acc_metric_left(v(1),v(2),v(3)) + branch_metric(v(1),v(2),v(3),ak1);
		        end;
		
		        [acc_metric_right(phk,ak2,ak1),surv_path(sym,phk,ak2,ak1)] = max(sub_acc_every_state);	% update acc_metric
		
		    end;end;end;			
				acc_metric_left = acc_metric_right;
		end;
		
		maximum = 0;					%由于是归零码，每桢的最后一个symbol从状态（p,1,1)里挑选最大的一个metric
		for ii = 1:p					
		    if(acc_metric_right(ii,1,1) > maximum)
				v = [ii,1,1];
				maximum = acc_metric_right(ii,1,1);
			end;
		end;
		
			
		for ii = frm_size:-1:1		%retrieve the route
		    dmod_bit((frm-1)*frm_size+ii) = 2*(v(3)-1) - (m-1);
		    v = reshape( state_from(v(1),v(2),v(3),surv_path(ii,v(1),v(2),v(3)),:),1,3);
		end;
    sprintf('%d frames demodualted!',frm)
 end;	%	end of (frm = 1:frm_num)

%reshape(acc_metric_right,p*m*m,1);
[err_num,err_rate] = symerr(framed_bit,dmod_bit);	% 计算误码率
sprintf('SER:%d',err_rate)

%%debug  Locate error-symbol position
jj = 0;
v = [];
for ii = 1:length(framed_bit)
    if framed_bit(ii) ~= dmod_bit(ii)
        jj = jj+1;
        v(jj) = ii;
    end;
end;
disp('done');