In [1]:

    

%pylab inline
pylab.rcParams['figure.figsize'] = (18, 6)

import numpy as np

sample_rate = 48000
tau = 2.0 * np.pi


    

    




Populating the interactive namespace from numpy and matplotlib

In [2]:

    

def plot_spectrum(signal, sample_rate):
    # http://samcarcagno.altervista.org/blog/basic-sound-processing-python/
    size = len(signal)
    fft_res = np.fft.fft(signal)
    fft_unique_points = int(ceil((size+1)/2.0))
    p = abs(fft_res[0:fft_unique_points])
    p = p / float(size)
    p = p ** 2
    if size % 2 == 1:
        p[1:len(p)] = p[1:len(p)] * 2
    else:
        p[1:len(p) - 1] = p[1:len(p) - 1] * 2
    freq_axis_array = arange(0, fft_unique_points, 1.0) * (sample_rate / size)
    plot(freq_axis_array/1000, 10*log10(p))
    #xscale('log')
    #xlim(xmin=20)
    xlabel('Frequency (kHz)')
    ylabel('Power (dB)')


    

In [3]:

    

def gen_additive_saw(osc_freq, sample_rate, nsamples):
    nyquist = sample_rate / 2.0
    harmonics = (osc_freq * x for x in range(1, int(math.ceil(nyquist / osc_freq))))
    
    end_phases = [((tau * f) / sample_rate) * nsamples
                  for f in harmonics]
    
    return np.sum([np.sin(np.linspace(0, end_phase, nsamples)) / (i + 1.0)
                   for (i, end_phase) in enumerate(end_phases)], axis=0) / 2.0


    

In [4]:

    

plot(gen_additive_saw(220, sample_rate, 800))


    

    
Out[4]:





[<matplotlib.lines.Line2D at 0x10b1cc7d0>]






    

In [5]:

    

plot_spectrum(gen_additive_saw(220, sample_rate, 48000), sample_rate)


    

In [6]:

    

def gen_naive_saw(osc_freq, sample_rate, nsamples):
    end_phase = (osc_freq / float(sample_rate)) * nsamples
    return 2.0 * (0.5 - np.mod(np.linspace(0, end_phase, nsamples), 1.0))


    

In [7]:

    

plot(gen_naive_saw(220, sample_rate, 800))


    

    
Out[7]:





[<matplotlib.lines.Line2D at 0x10bb7b610>]






    

In [8]:

    

plot_spectrum(gen_naive_saw(220, sample_rate, 48000), sample_rate)


    

In [30]:

    

def polyblep_2sample(phase, phase_step):
    if phase < phase_step:
        t = phase / phase_step
        blep = (t*2) - (t**2) - 1.0
        print("blep!\nphase: %s\nphase-step: %s\nblep: %s" % (phase, phase_step, blep))
        return blep
    
    if phase > (1.0 - phase_step):
        t = (phase - 1.0) / phase_step
        blep = (t**2) + (t*2) + 1.0
        print("blep!\nphase: %s\nphase-step: %s\nblep: %s" % (phase, phase_step, blep))
        return blep
    
    return 0.0

def gen_polyblep_saw(osc_freq, sample_rate, nsamples):
    end_phase = (osc_freq / float(sample_rate)) * nsamples
    
    phases = np.mod(np.linspace(0, end_phase, nsamples), 1.0)
    
    step = phases[1] - phases[0]
    print("phases[0]: %s\nphases[1]: %s" % (phases[0], phases[1]))
    print("step: %s" % step)
    polyblep = [polyblep_2sample(x, step) for x in phases]
    
    return (1.0 - (2.0 * phases)) + polyblep


    

In [32]:

    

plot(gen_polyblep_saw(440, sample_rate, 250))


    

    




phases[0]: 0.0
phases[1]: 0.00920348058902
step: 0.00920348058902
blep!
phase: 0.0
phase-step: 0.00920348058902
blep: -1.0
blep!
phase: 0.993975903614
phase-step: 0.00920348058902
blep: 0.119338842975
blep!
phase: 0.00317938420348
phase-step: 0.00920348058902
blep: -0.428429752066
blep!
phase: 0.997155287818
phase-step: 0.00920348058902
blep: 0.477355371901
blep!
phase: 0.00635876840696
phase-step: 0.00920348058902
blep: -0.0955371900826






    
Out[32]:





[<matplotlib.lines.Line2D at 0x10bd1df50>]






    

In [11]:

    

plot_spectrum(gen_polyblep_saw(880, sample_rate, 48000), sample_rate)


    

In [12]:

    

import scipy.signal

plot_spectrum(
    scipy.signal.resample_poly(gen_polyblep_saw(220, sample_rate * 2, 96000), 1, 2),
    sample_rate)


    

In [13]:

    

plot_spectrum(gen_additive_saw(220, sample_rate, 48000), sample_rate)


    

In [14]:

    

plot_spectrum(gen_additive_saw(220, sample_rate, 48000), sample_rate)
plot_spectrum(
    scipy.signal.resample_poly(gen_polyblep_saw(220, sample_rate * 2, 96000), 1, 2),
    sample_rate)


    

In [15]:

    

from IPython.display import Audio

freq = 3520
seconds = 2


    

In [16]:

    

Audio(gen_additive_saw(freq, sample_rate, sample_rate * seconds), rate=sample_rate)


    

    
Out[16]:





                
              

In [17]:

    

Audio(gen_naive_saw(freq, sample_rate, sample_rate * seconds), rate=sample_rate)


    

    
Out[17]:





                
              

In [18]:

    

Audio(gen_polyblep_saw(freq, sample_rate, sample_rate * seconds), rate=sample_rate)


    

    
Out[18]:





                
              

In [19]:

    

Audio(
    scipy.signal.resample_poly(gen_polyblep_saw(freq, sample_rate * 2, sample_rate * 2 * seconds), 1, 2),
    rate=sample_rate)


    

    
Out[19]:





                
              

In [ ]: