四种通信序列

一、Costas序列（常用在跳频中）

在数学中，Costas 数组可以在几何上视为 n 个点的集合，每个点位于 n×n 方块平铺中一个正方形的中心，使得每一行或每一列仅包含一个点，并且每对点之间的 n(n − 1)/2 个位移向量都是不同的。这导致了一个理想的“图钉”自模糊函数，使得这些数组在声纳和雷达等应用中非常有用。可以由对数Welch方法构建。[4]

下图是N=4的costas矩阵

这个代码老师没给，仿真时直接上网查输进去的QAQ，但是给了一本书《Radar Signals》(Nadav Leavnopn, Eli Mozaeeson)（这本书年龄比我爷爷都大……），本人用了，一堆报错，真服了……以下是修改改后的代码[10]

function positions = costas(N)
% 修复后的函数，修正本原元判断逻辑

positions = [];
if isprime(N+1) && rem(N+1,2)==1  % Type 1
    p = N+1;
    primitive_roots = [];
    for e = 2:p-1
        powers = mod(e.^(1:p-1), p);
        if all(powers(1:p-2) ~= 1) && powers(p-1) == 1
            primitive_roots = [primitive_roots; e];
        end
    end
    if ~isempty(primitive_roots)
        positions = zeros(length(primitive_roots), p-1);
        for i = 1:length(primitive_roots)
            e = primitive_roots(i);
            positions(i, :) = mod(e.^(1:p-1), p);
        end
    end
    
elseif isprime(N+2) && rem(N+2,2)==1  % Type 2
    p = N+2;
    primitive_roots = [];
    for e = 2:p-1
        powers = mod(e.^(1:p-1), p);
        if all(powers(1:p-2) ~= 1) && powers(p-1) == 1
            primitive_roots = [primitive_roots; e];
        end
    end
    if ~isempty(primitive_roots)
        positions = zeros(length(primitive_roots), N);
        for i = 1:length(primitive_roots)
            e = primitive_roots(i);
            freq_sequence = mod(e.^(2:p-1), p) - 1;
            positions(i, :) = freq_sequence;
        end
    end
    
elseif isprime(N+3) && rem(N+3,2)==1  % Type 3 (检查2是否为GF(p)的本原元)
    p = N+3;
    e = 2;
    powers = mod(e.^(1:p-1), p);
    if all(powers(1:p-2) ~= 1) && powers(p-1) == 1
        positions = mod(e.^(3:p-1), p) - 2;
    end
end
end

函数中使用了 Welch 构造方法，这是一种常用的构造 Costas 数组的方法。
函数中还涉及到有限域（GF(p)）的概念，这是数学和密码学中的一个概念，用于确定哪些元素是原始元素。

1	positions = costas(4);

典型的应用案例：

N	构造类型	应用领域	示例序列（部分）	来源
4	Welch1	通信跳频	`[2,4,3,1]` 或 `[3,4,2,1]`	[7] [8]
10	Welch1	雷达抗干扰	`[2,4,8,5,10,9,7,3,6,1]`	[8]
29	未明确	高频声纳	`[3,21,23,...,25]`（29元素序列）	[8]
15	Welch2	密码学中的密钥序列生成	基于GF(17)生成	[9]

接下来介绍的都是相位编码

二、m序列

m序列是最大长度序列的简称（Maximum length sequence，MLS），就是线性以为寄存器通过线性反馈产生的最长序列，规定级数为r时，产生的码序列为
MLS的实际应用包括测量脉冲响应（例如，房间混响或来自海洋拖曳源的到达时间 [1]）,此外还可作为伪随机序列用于直接序列扩频和跳频传输的数字通信系统中。[2]

1.m序列的生成

假设生成多项式是,则移位器如图所示，若移位器系数为1，则之前的线路与异或器相连，如果为0则不相连。

MLS 是通过最大线性反馈移位寄存器生成的。下图显示了一个具有长度为 4 的移位寄存器的 MLS 生成系统。它可以用以下递归关系表示：

其中 n 是时间索引， + 表示模2 加法。对于比特值 0 = FALSE 或 1 = TRUE，这相当于异或操作。

此处有介绍的视频：通信原理4.6.1m序列_哔哩哔哩_bilibili

代码：

function seq = m_seq(oct,gen)

% function seq = m_seq(oct)
% Generates an m-sequence using generator given by oct.
% Reference:  Dilip V. Sarwate and Michael B. Pursley, 1980. pp. 599. Fig.1.


s = min(find(gen));
gen = gen(s+1:end);
n = size(gen, 2);
N = 2^n-1;
gen = fliplr(gen);

seq = zeros(1, n); seq(n) = 1;
for i=1:N-n
        next_bit = mod(sum(seq(i:i+n-1)&gen), 2);
        seq = [seq, next_bit];
end

m = 5 ; l = 2^m-1 ; %阶数：5  
oct1 = 45 ;
gen1  = [1 0 0 1 0 1] ;%生成多项式
mseq1 = m_seq(oct1,gen1);%生成m序列

2.m序列属性

0和1的数量大致相同
MLS 的线性自相关近似于δ函数[3]

3. 脉冲响应的提取

如果系统的冲激响应是 h[n]，而 MLS 是 s[n]，那么

对两边与s[n]进行相关处理

$hi_{sy}=h[n]*hi{ss}$

由于性质2，当长度足够长时， $hi_{ss}$ 近似为一个冲击函数，所以s[n]与y[n]的互相关就是脉冲响应

$h[n]=hi_{sy}$

https://readit.site/a/XYlZG/Maximum_length_sequence

三、Gold序列

Gold序列常用于CDMA[5]和GPS卫星导航[6]

Gold codes have bounded small cross-correlations within a set, which is useful when multiple devices are broadcasting in the same frequency range.

引文所说，在一个集合内，互相关性能很小，因此在同意同一频带多用户同通信时很有用

Gold序列的互相关极大值

1.Gold序列的生成

gold序列是由两个M序列异或得到

https://readit.site/a/popNM/Gold_code

function seqs = gold(oct1,gen)
% 此函数用于生成gold序列，输入参数为oct1：和gen：生成多项式系数
% function seqs = gold(oct1)
% This function generates N+2=2^n+1 Gold sequences when (n mod 4 != 0);
% and N+1=2^n Gold-like sequences when (n mod 4 == 0).
% oct1 is the generator polynomial in oct form.
% reference: Dilip V. Sarwate and Michael B. Pursley, 1980.

% some equavalent generators.
% m=3, 13
% m=4, 23
% m=5, 45
% m=6, 103
% m=7, 211

u = m_seq(oct1,gen);
N = size(u, 2);
n = log2(N+1);
t = 1 + 2^(floor(n/2+1));
v = u(mod(0:t:N*t-1, N)+1);

if mod(n, 4) ~=0
        seqs = [u; v];
        for i = 1:N
                v = rshift(v);
                seqs = [seqs; xor(u, v)];
        end
        fprintf(1, 'Gold sequences\n');
else
        seqs = [u];
        v1 = u(mod(0:t:N*t-1, N)+2);
        v2 = u(mod(0:t:N*t-1, N)+3);
        for i = 1:N/3
                seqs = [seqs; xor(u, v); xor(u, v1); xor(u, v2)];      
                v = rshift(v); v1 = rshift(v1); v2 = rshift(v2);
        end
        fprintf(1, 'Gold-like sequences\n');
end


function seq = m_seq(oct,gen)

% function seq = m_seq(oct)
% Generates an m-sequence using generator given by oct.
% Reference:  Dilip V. Sarwate and Michael B. Pursley, 1980. pp. 599. Fig.1.


s = min(find(gen));
gen = gen(s+1:end);
n = size(gen, 2);
N = 2^n-1;
gen = fliplr(gen);

seq = zeros(1, n); seq(n) = 1;
for i=1:N-n
        next_bit = mod(sum(seq(i:i+n-1)&gen), 2);
        seq = [seq, next_bit];
end

function y = rshift(x)
y = [x(:, end), x(:, 1:end-1)];

m = 5 ; l = 2^m-1 ; %阶数：5  
oct1 = 45 ;
gen  = [1 0 0 1 0 1] ;%生成多项式
seqs = gold(oct1,gen);%生成gold序列


gold1 = seqs(1,:);
gold2 = seqs(15,:);
gold3 = seqs(30,:);

四、Kasami序列

Kasami 序列是长度为 2 N −1 的二进制序列，其中 N 是偶整数。Kasami 序列具有良好的交叉相关值，接近 Welch 下界。Kasami 序列分为两类——小集和大集。

1kasami序列的生成：

选择两个不同的M序列： Kasami序列通常由两个长度相同的M序列（最大长度序列）生成。第一个M序列的生成通常使用一个基准的反馈多项式，而第二个M序列则通过将第一个序列进行适当的位移来获得。
通过异或操作生成Kasami序列：通过将这两个M序列进行异或操作得到Kasami序列。对于一个给定的序列，可以通过改变异或操作的位置（即通过移位）来生成多个Kasami序列。

eg. 假设我们要生成一个长度为15的Kasami序列，使用两个4阶M序列。

步骤1：生成两个M序列

假设两个M序列的生成多项式分别为：

P1 = x^4 + x + 1，对应的反馈结构生成M序列 M1。
P2 = x^4 + x^3 + 1，对应的反馈结构生成M序列 M2

步骤2：生成Kasami序列

选择一个移位值 d（例如，d = 3），然后将第二个M序列 M2 进行移位并与 M1 异或，生成Kasami序列。就获得了长度为15的kasami序列

function seqs = lkasami(oct,gen)

% function seqs = lkasami(oct)
% This function generates Large Kasami sequences.
% oct generate an m-sequence of period 2^n-1
% oct should be in oct form. n should be even.
% The first (N+1)^(0.5) sequences are the small set Kasami sequences. 
% reference: Dilip V. Sarwate and Michael B. Pursley, 1980.

% Qinghua Zhao
% Aug. 2001 at UCSD

%1.先生成m序列
seqs = [];
u = m_seq(oct,gen);
N = size(u, 2); 
n = log2(N+1);
if mod(n, 2) ~= 0
	fprintf(1, 'No Kasami sequence exist for generator of this order.\n');
	beep;
	return;
end 

t = 1 + 2^(n/2+1); % no need to take floor since it will be an interger.
v = u(mod(0:t:N*t-1, N)+1);

% generate w
N_2 = 2^(n/2)+1;
w = u(1:N_2:end);
if sum(w)==0 % if all zero sequence, then pick another one. It won't be all one sequence since it can only be an m-sequence. 
	w = u(2:N_2:end);
end
w = ones(N_2, 1)*w;
w = reshape(w', 1, N);

fprintf(1, 'First %d sequences are small set Kasami\n', N_2-1); 

if mod(n, 4) ~= 0
	% generate gold sequences 
	G_uv = [u; v];
	for i = 1:N
	        G_uv = [G_uv; xor(u, v)];
        	v = rshift(v);
	end

	% generate kasami sequences
	% the first 2^(n/2) sequences are the small set kasami sequences.
	for k = 1:size(G_uv, 1)
		seqs = [seqs; G_uv(k,:)];
		for i = 1:N_2-2
        		seqs = [seqs; xor(G_uv(k,:), w)];
        		w = rshift(w);
		end
	end
else
	% generate gold like sequences
	v1 = u(mod(0:t:N*t-1, N)+2);
	v2 = u(mod(0:t:N*t-1, N)+3);
	Ht_u = u;
	for i = 1:N/3
 		Ht_u = [Ht_u; xor(u, v); xor(u, v1); xor(u, v2)];	
        	v = rshift(v); v1 = rshift(v1); v2 = rshift(v2);
	end

	% generate kasami sequences
	% the first 2^(n/2) sequences are the small set kasami sequences.
	for k = 1:size(Ht_u, 1)
		seqs = [seqs; Ht_u(k,:)];
		for i = 1:N_2-2
        		seqs = [seqs; xor(Ht_u(k,:), w)];
        		w = rshift(w);
		end
	end
	for i = 1:(N_2-2)/3
		seqs = [seqs; xor(v, w); xor(v1, w); xor(v2, w)];
        	w = rshift(w);
	end


end

% 该函数用于生成m序列，基于给定的生成器oct。
% 参考文献：Dilip V. Sarwate和Michael B. Pursley，1980年。第599页。图1。

function seq = m_seq(oct,gen)

% gen = oct2gen(oct);     % oct2gen may be removed from future version
%1.查找非零索引，并将序列从该位置之后的位置提取出来
s = min(find(gen));
gen = gen(s+1:end);

%2.计算生成器的位数和序列的长度
n = size(gen, 2);
N = 2^n-1;

%3.将生成器反转
gen = fliplr(gen);

seq = zeros(1, n); seq(n) = 1;
for i=1:N-n
        next_bit = mod(sum(seq(i:i+n-1)&gen), 2);
        seq = [seq, next_bit];
end

return

function y = rshift(x)
y = [x(:, end), x(:, 1:end-1)];
return


m=4; 
oct = 23;
gen = [1 0 0 1 1];

seqs = lkasami(oct,gen);%生成gold序列
kasami1 = seqs(6,:);
kasami2 = seqs(15,:);

五、其它编码

Phase-Coded Pulse: Barker码
Chirplike Phase Codes:
- Frank Code
- P1,P2和Px Codes
- Zadoff-Chu Code
- Golomb Polyphase Codes
Huffman Code

附上常用的本征多项式：

参考文章：

[1] Gemba, Kay L.; Vazquez, Heriberto J.; Fialkowski, Joseph; Edelmann, Geoffrey F.; Dzieciuch, Matthew A.; Hodgkiss, William S. (October 2021). “A performance comparison between m-sequences and linear frequency-modulated sweeps for the estimation of travel-time with a moving source”. The Journal of the Acoustical Society of America. 150 (4): 2613–2623. Bibcode:2021ASAJ…150.2613G. doi:10.1121/10.0006656. PMID 34717519. S2CID 240355915.

[2] Buracas GT, Boynton GM (July 2002). “Efficient design of event-related fMRI experiments using M-sequences”. NeuroImage. 16 (3 Pt 1): 801–13. doi:10.1006/nimg.2002.1116. PMID 12169264. S2CID 7433120.

[3] Jacobsen, Finn; Juhl, Peter Moller (2013-06-04). Fundamentals of General Linear Acoustics. John Wiley & Sons. ISBN 978-1118636176. A maximum-length sequence is a binary sequence whose circular autocorrelation (except for a small DC-error) is a delta function.

[4] Costas (1965); Gilbert (1965); An independent discovery of Costas arrays, Aaron Sterling, October 9, 2011.

[5] George, Maria; Hamid, Mujtaba; Miller, Andy (2001-01-10). “Gold Code Generators in Virtex Devices” (PDF). Virtex Series, Virtex-II Series, and Spartan-II family (Application note). 1.1. Xilinx. XAPP217. Archived from the original (PDF) on 2008-07-05. (9 pages)

[6] “Transmitted GPS Signals”. The GPS System. kowoma GPS. 2009-04-19. Archived from the original on 2012-08-04.

[7] GOLOMB S W, TAYLOR H. Constructions and properties of Costas arrays[J/OL]. Proceedings of the IEEE, 1984: 1143-1163. https://doi.org/10.1109/proc.1984.12994. DOI:10.1109/proc.1984.12994.

[8] https://blogs.mathworks.com/steve/2007/02/15/costas-arrays/

[9]KEN TAYLOR,SCOTT RICKARD,KONSTANTINOS DRAKAKIS．Costas Arrays: Survey, Standardization, and MATLAB Toolbox[J]．ACM transactions on mathematical software,2011,37(4):p.38:59-38:89

[10] Levanon, Nadav, E. Mozeson, R. Signals and Chaim Levanon. “RADAR SIGNALS.” (2013).