generate.py

#!/usr/bin/python

# Copyright (c) 2016, Gordon JC Pearce <gordon@gjcp.net>
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
#
# 1. Redistributions of source code must retain the above copyright notice,
# this list of conditions and the following disclaimer.
#
# 2. Redistributions in binary form must reproduce the above copyright notice,
# this list of conditions and the following disclaimer in the documentation
# and/or other materials provided with the distribution.

# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.

# To use this, just run it with "python generate.py" and play back the
# two generated .wav files - play "naive.wav" first then "antialiased.wav"
# for maximum effect

import numpy as np
import scipy.io.wavfile as wav

s_len = 96000       # how many samples to generate
s_rate = 22050.0    # pick something low enough to show aliasing

f_start = 110       # frequency to start at
f_range = 6000      # uncomfortably close to Nyquist


def polyblep(phase, inc):
    # place a blep at a discontinuity
    if phase < inc:
        # we've just reset
        p = phase / inc
        return p+p - p*p - 1.0
    elif phase > (1.0-inc):
        # we're just about to reset
        p = (phase-1.0) / inc
        return p+p + p*p + 1.0
    else:
        # no need to do anything
        return 0

# define a couple of blank arrays to hold the samples
s_a = np.zeros(s_len)
s_n = np.zeros(s_len)

phase = 0

for i in range(0, s_len):

    ph_inc = (f_start + (f_range*(i/float(s_len)))) / s_rate

    y = (2*phase)-1     # must go from -1 to 1
    s_a[i] = y - polyblep(phase, ph_inc)
    s_n[i] = y

    # phase accumulator
    phase = phase + ph_inc
    if phase > 1.0:
        phase = phase - 1.0

wav.write("antialiased.wav", s_rate, s_a)
wav.write("naive.wav", s_rate, s_n)