M序列与伪随机二值图

本人没研究过算法，只是偶尔写程序接触到了一点，做个记录。
M序列是伪随机二进制序列，英文maximum length sequence (MLS)就是最长序列的意思。
它的产生是通过线性反馈移位寄存器(LFSR，linear feedback shift register)，也就是多项式+一个移位反馈。逻辑上算是比较简单的。
下面是一个4级的LFSR。

因为是伪随机，所以实际上是有周期的。m级LFSR序列最长周期为2^m-1，因此要获得较好的随机性，多项式的次数要越高越好。
如何衡量生成序列的随机性有多好呢？有个简单的可视化方法，就是伪随机二值图。
用Python写了一个简单程序，作图如下。可以看到m=8的时候伪随机序列还可以明显看到周期现象，到m=16时已经和np.random.randint(2, size=(100, 100))生成的伪随机序列看不出差异了。

附上Python小程序，多项式公式参考liquidsdr.org

import numpy as np
import matplotlib.pyplot as plt

def mseq_m4(seed=1, length=100):
    """ M-sequence with m=4
    """
    m_state = seed
    sequence = []

    for _ in range(length):
        bit = m_state & 0x1
        sequence.append(bit)
        feedback = ((m_state >> 0) ^ (m_state >> 3) ^ (m_state >> 4)) & 0x1
        m_state = (m_state >> 1) | (feedback << 3)

    return sequence

def mseq_m8(seed=1, length=100):
    """ M-sequence with m=8
    """
    m_state = seed
    sequence = []

    for _ in range(length):
        bit = m_state & 0x1
        sequence.append(bit)
        feedback = ((m_state >> 0) ^ (m_state >> 4) ^ (m_state >> 5) ^ (m_state >> 6) ^ (m_state >> 8)) & 0x1
        m_state = (m_state >> 1) | (feedback << 7)

    return sequence

def mseq_m16(seed=1, length=100):
    """ M-sequence with m=16
    """
    m_state = seed
    sequence = []

    for _ in range(length):
        bit = m_state & 0x1
        sequence.append(bit)
        feedback = ((m_state >> 0) ^ (m_state >> 4) ^ (m_state >> 13) ^ (m_state >> 15) ^ (m_state >> 16)) & 0x1
        m_state = (m_state >> 1) | (feedback << 15)

    return sequence


if __name__ == '__main__':
    col, row = 100, 100
    mseq_m4 = mseq_m4(1, col * row)
    binary_m4 = np.array(mseq_m4).reshape(col, row)

    mseq_m8 = mseq_m8(1, col * row)
    binary_m8 = np.array(mseq_m8).reshape(col, row)

    mseq_m16 = mseq_m16(1, col * row)
    binary_m16 = np.array(mseq_m16).reshape(col, row)

    np.random.seed(0)
    binary_np = np.random.randint(2, size=(col, row))

    # 创建一个包含多个子图的图像窗口
    fig, axs = plt.subplots(1, 4, figsize=(20, 6))

    # 在每个子图中绘制一个二进制数组
    axs[0].imshow(binary_m4, cmap='gray', interpolation='nearest')
    axs[0].set_title('M-sequence, m=4')

    axs[1].imshow(binary_m8, cmap='gray', interpolation='nearest')
    axs[1].set_title('M-sequence, m=8')

    axs[2].imshow(binary_m16, cmap='gray', interpolation='nearest')
    axs[2].set_title('M-sequence, m=16')

    axs[3].imshow(binary_np, cmap='gray', interpolation='nearest')
    axs[3].set_title('np.random')

    plt.show()