www.pudn.com > CPM4TR.rar > cpm_sim.m
% clc
clear all
close all
T=1;
L=3;
h=2/7;
m=8; % 八进制
p = 2/h; % 累计相位的状态个数
OVS=24;
OVR=24;
EbN0 = 10;
SNR = EbN0 - 10 * log10(OVR) + 10 * log10(log2(m));
SNR = 10;
% SNR = 1;
block_num = 10; %发送信息的桢数
block_size = 26; %发送信息的桢大小 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_size+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*dmod_sig1(index1:index1+OVR-1)'; %'
offset_sig = offset_sig * dfe_sig^-0.5;
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 correlator32 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');