Multitone test signal generation (and Rudin-Shapiro sequence) in Python

readme.md

Reproduction of multitone test signals from Boyd - Multitone signals with low crest factor

• The zero-phase method just produces impulses, which are not very useful because the energy is so concentrated in time, it's easy to distort, get bad SNR, etc.
• The Rudin-Shapiro method produces noise-like tones, probably best for most things because it's white-like at all points in the signal.
• The Newman method produces sweep-like tones, which are good for anything else that would be tested using a sweep.
  • Playing this on my computer causes my cat to freak out, try to burrow into my computer and bite my arm hard, refusing to let go. None of the other signals has this effect.

Figure 1 zero-phase.png

Figure 4 Rudin-Shapiro.png

Figure 6 Newman.png

multitone.py

"""
Multitone test signal generation.

Created on Fri Sep 28 12:11:26 2012

Reproduction of signals from
http://www.stanford.edu/~boyd/papers/pdf/multitone_low_crest.pdf

These have exactly flat frequency spectra at the FFT frequency bin centers.
"""
from numpy import (pi, zeros, empty, arange, sqrt, exp, log10, amax,
                   linspace, mean, absolute)
from numpy.random import random
from numpy.fft import rfft, rfftfreq, irfft


# Originally from matplotlib.mlab:
def rms_flat(a):
    """
    Return the root mean square of all the elements of *a*, flattened out.
    """
    return sqrt(mean(absolute(a)**2))


def rudinshapiro(N):
    """
    Return first N terms of Rudin-Shapiro sequence.

    https://en.wikipedia.org/wiki/Rudin-Shapiro_sequence

    Confirmed correct output to N = 10000:
    https://oeis.org/A020985/b020985.txt
    """
    def hamming(x):
        """
        Hamming weight of a binary sequence.

        http://stackoverflow.com/a/407758/125507
        """
        return bin(x).count('1')

    out = empty(N, dtype=int)
    for n in range(N):
        b = hamming(n << 1 & n)
        a = (-1)**b
        out[n] = a

    return out


def multitone(type, fs, length, N=None, N0=1):
    """
    Return a multitone signal of given type.

    fs = sampling rate
    length = length in samples
    N = number of tones
    N0 = starting tone
    """
    if N is None:
        N = length//2

    if type == 'impulse':
        """
        Zero-phase method
        Figure 1
        Worst crest factor possible (18 dB for N=32)
        """
        phase = zeros(N)

    elif type == 'random':
        """
        Random phase method
        Noise, with crest factor on order of sqrt(log(N)) (5 dB for N=32)
        """
        phase = random(N)*2*pi

    elif 'rudin' in type:
        """
        Rudin-Shapiro sequence method
        Figure 4
        Noise-like, with crest factor under 6 dB when N is a power of 2.
        """
        phase = -pi*rudinshapiro(N)
        phase[phase == -pi] = 0

    elif 'newman' in type:
        """
        Newman method
        Figure 6
        Sweep-like, with crest factor about 4.6 dB at higher N
        """
        k = arange(N) + N0
        phase = (pi*(k-1)**2)/N
        if type in {'newmanreversed', 'newman_r'}:
            phase = phase[::-1]

    # # Explicit construction
    # t = linspace(0, 2*pi, length, endpoint=False)
    # sig = zeros_like(t)
    # for k in arange(N)+N0:
    #     sig += cos(k*t + phase[k-N0])
    # sig *= sqrt(2/N)

    # Inverse FFT construction
    f = zeros(length//2+1, dtype=complex)
    for k in arange(N)+N0:
        f[k] = exp(1j*phase[k-N0])
    sig = irfft(f, length) * length/2 * sqrt(2/N)

    return sig


if __name__ == '__main__':
    from matplotlib import pyplot as plt

    # sampling rate
    fs = 512

    length = fs*1  # seconds

    # impulse, random, rudin, newman
    type = 'rudin'

    # number of tones
    N = 32

    # Starting tone
    N0 = 1

    sig = multitone(type, fs, length, N, N0)

    crest_factor = 20*log10(amax(abs(sig))/rms_flat(sig))
    print('Crest factor: {0:.2f} dB'.format(crest_factor))

    t = linspace(0, 2*pi, length, endpoint=False)
    plt.subplot(2, 1, 1)
    plt.plot(t, sig)
    plt.margins(0.1, 0.1)

    # Use FFT to get the amplitude of the spectrum
    ampl = 1/fs * abs(rfft(sig))

    # FFT frequency bins
    freqs = rfftfreq(fs, 1/fs)

    plt.subplot(2, 1, 2)

    plt.plot(freqs, ampl, '.')
    plt.margins(0.1, 0.1)

Should also try these http://dsp.stackexchange.com/q/18862/29

Summary of several different signals:

Blue is Maximum Length Sequence after going through the analog reconstruction / interpolation process (which ruins the theoretically perfect crest factor):

So for long signals, the sweep-like Newman tone is better than the MLS.