goldcode.m

张贤达老师 通信信号处理4－9章的程序 这本书是张老师的经典之作。

function gold_code = goldcode(r)
%  GOLDCODE	Gold Code generation
%    GOLD_CODE = GOLDCODE(R) produces 2^R + 1 Gold code, including 2 m-sequences.
%    The code length is 2^R - 1.

gold_code	= [];
%load table				% get the table of characteristic polynomial for m-sequences
table =['#2#7,*#3#13,*#4#23,*#5#45,*75,67,#6#103,*147,155,#7#211,*217,235,367,277,325,203,*313,345,#8#435,551,747,453,545,537,703,543,#9#1021,*1131,1461,1423,1055,1167,1541,1333,1605,1751,1743,1617,1553,1157,#10#2011,*2415,3771,2157,3515,2773,2033,2443,2461,3023,3543,2745,2431,3177,#11#4005,*4445,4215,4055,6015,7413,4143,4563,4053,5023,5623,4577,6233,6673,#13#20033,23261,24623,23517,30741,21643,30171,21277,27777,35051,34723,34047,32535,31425,#14#42103,43333,51761,40503,77141,62677,44103,45145,76303,64457,57231,64167,60153,55753,#15#100003,*102043,110013,102067,104307,100317,177775,103451,110075,102061,114725,103251,100021,*100201,*#16#210013,234313,233303,307107,307527,306357,201735,272201,242413,270155,302157,210205,305667,236107,#17#400011,*400017,400431,525251,410117,400731,411335,444257,600013,403555,525327,411077,400041,*400101,*#18#1000201,*1000247,1002241,1002441,1100045,1000407,003011,1,020121,1101005,1000077,1001361,1001567,1001727,1002777,#19#2000047,2000641,2001441,2000107,2000077,2000157,2000175,2000257,2000677,2000737,2001557,2001637,2005775,2006677,#20#4000011,*4001051,4004515,6000031,4442235,#21#10000005,10040205,10020045,10040315,10000635,10103075,10050335,10002135,17000075,#22#20000003,*20001043,22222223,25200127,20401207,20430607,20070217,#23#40000041,*40404041,40000063,40010061,50000241,402201,40006341,40405463,40103271,41224445,40435651,#24#100000207,125245661,113763063,#25#200000011,*200000017,204000051,200010031,200402017,252001251,201014171,204204057,200005535,209014731,#26#400000107,430216473,402365755,426225667,510664323,473167545,411335571,#27#1000000047,1001007071,1020024171,1102210617,1250025757,1020560103,1112225171,1035530241,#28#2000000011,*2104210431,2000025051,2020006031,2002502115,2001601071,#29#4000000005,*4004004005,4000010205,4010000045,4400000045,4002200115,4001040115,4004204435,4100060435,4040003075,4004064275,#30#10040000007,10104264207,10115131333,11362212703,10343244533,#31#20000000011,*20000000017,20000020411,21042104211,20010010017,20005000251,20004100071,20202040217,20000200435,20060140231,21042107357,#32#40020000007,40460216667,40035532523,42003247143,41760427607,#33#100000020001,*100020024001,104000420001,100020224401,111100021111,100000031463,104020466001,100601431001,#34#201000000007,201472024107,377000007527,225213433257,227712240037,251132516577,211636220473,200000140003,#35#400000000005,*#36#1000000004001,*#37#2000000012005,#38#4000000000143,#39#10000000000021,*#40#20000012000005,#61#200000000000000000047,#89#400000000000000000000000000151,'];

 

% to find the characteristic polynomial
num	= ['#', num2str(r), '#'];
index	= findstr(table, num);
if isempty(index)
   error('Not found in table');
else
   index	= index + length(num);
   result	= table(index : length(table));
   e_index	= findstr(result, '#');
   if ~isempty(e_index)
      result	= result(1 : e_index - 1);
   end
end

Code_len	= pow2(r) - 1;

num	= [];
count	= 0;
record	= [];
for index = 1 : length(result)
   if ~isempty(findstr('0123456789', result(index)))		% is a digital character
      num	= [num, result(index)];
   elseif ~isempty(num)
      count	= count + 1;
      dec_n	= base2dec(num, 8);
      record(count)	= dec_n;				% record decimal form of characteristic polynomial
      num	= [];								% clear number recording buffer
   end
end
if ~isempty(num)								% process for the last number
   count	= count + 1;
   dec_n	= base2dec(num, 8);
   record(count)	= dec_n;
end

if count < 2								% no enough table item
   disp('No enough polynomial in table');
   return;
end

% to choose optimal pair randomly
char_poly	= zeros(r, count);
m_seq	= zeros(Code_len, count);
for k = 1 : count
   num	= dec2bin(record(k));
   dec_n	= (num == '1');
   dec_n	= reshape(dec_n(1 : r), [r, 1]);
   char_poly(:, k)	= dec_n;
   m_seq(:, k)	= maxseq(dec_n);
end

r_choose	= floor(rand(1) * count + 1);
r_choose	= r_choose - (r_choose > count);			% to ensure to be less than count
indices	= floor((r + 2) / 2);
coe_op_pair	= pow2(indices) + 1;
found	= 0;				% condition of looking for optimal pair
for k = r_choose : r_choose + count - 2
   if k > count
      k0	= k - count;
   else
      k0	= k;
   end
   
   for l = k + 1 : r_choose + count - 1
      if l > count
         l0	= l - count;
      else
         l0	= l;
      end
      
      % computing cross-covariance
      result	= covar_bin(m_seq(:, k0), m_seq(:, l0));
      if result <= coe_op_pair
         found	= 1;
         break;
      end
   end
   
   if found == 1
      break;
   end
end

if found < 1
   disp('No optimal pair available');
   return;
end

count	= Code_len + 2;
gold_code	= zeros(Code_len, count);
gold_code(:, 1)	= m_seq(:, k0);
gold_code(:, 2)	= m_seq(:, l0);
parent1	= m_seq(:, k0);
parent2	= m_seq(:, l0);

for index	= 3 : count
   result	= bitxor(parent1, parent2);
%   result	= bitand(result, 1);
   gold_code(:, index)	= result;
   parent1	= [parent1(2 : end); parent1(1)];
end
gold_code = sign(gold_code-0.5);

function sequence = maxseq(polynomial)
%	MAXSEQ	Maximum-Length Sequence generator
%		SEQUENCE = MAXSEQ(POLYNOMIAL) generates the m-sequence corresponding to
%		the characeristic polynomial POLYNOMIAL. All arrays are assumed to be in
%		the vector form.

polynomial	= polynomial(:);			% to turn into vector form
n	= length(polynomial);
vector	= zeros(n, 1);
vector(1)	= 1;			% initial state [an,...,a1,a0] = vector = e1
sequence	= vector;

% vector(1) corresponds to a(n), polynomial(1) corresponds to feedback of a(n)
for k = 1 : pow2(n) - n - 1
   feed_back	= polynomial' * vector;
   feed_back	= bitand(feed_back, 1);			% be same as mod(2) operation
   sequence	= [sequence; feed_back];
   vector	= [vector(2 : n); feed_back];
end

function covar_m_seq = covar_bin(m_seq1, m_seq2)
%  COVAR_BIN	Covariance compution of binary data
%		COVAR_M_SEQ = COVAR_BIN(M_SEQ1, M_SEQ2) gives the maximum covarian of the 
%		two m_sequences M_SEQ1 and M_SEQ2.

N	= length(m_seq1);
covar_m_seq	= 0;
for k = 1 : N
   covar1	= bitxor(m_seq1, m_seq2);
   num	= sum(covar1);
   num	= abs(N - num * 2);
   if num > covar_m_seq
      covar_m_seq	= num;
   end
   m_seq1	= [m_seq1(2 : end); m_seq1(1)];
end