Python/Musical intervals (numpy matplotlib)

This code creates three graphs of musical intervals that deviate from just intonation by a user-selected frequency. The graphs exhibit something analogous to an acoustical beat, except that these beats involve only phase shifts because each wave is a pure sinusoidal. To ensure that all three timescales are identical, it is important that the first graph plotted have the lowest beat frequency (here it is the perfect fifth.)

The code requires numpy and matplotlib. The modules scipy and sound file are called, but may be removed since they are not used in this version of the code. A sample output is shown below:

Plots

edit
Beat pattern for a tritone that is off by 4.32 cents.
Beat pattern for a perfect fifth that is off by 4.32 cents.
Beat pattern for a minor sixth that is off by 4.32 cents.

Code

edit

import numpy as np #------------------------numerical operations on large arrays
import matplotlib.pyplot as plt #-----------creates  and displays plots
import scipy, sys,os, time, math #----------scipy not used here
from scipy import signal #------------------signal not used here
import soundfile as sf #--------------------soundfile not used here
start_time = time.time()#-------------------lets me see run time for code
#FORMATTING CONSTANTS ----------------------these determine plot's appearance
figHigh,figWide=3,50#-----------------------figure Height and Width
Amp=1; samplerate=44100#--------------------standard CD sampling (not used)
dt=1/samplerate#----------------------------dt to samplerate conversion
slw=1#-------------------------------------slw signal's line width
barx=1#------------------------------------barx quasi-period marker width
barb=1.4#-----------------------------------barb beat marker width
fs=10#--------------------------------------fontsize, txp (latter not used?)
show_images=False#--------------------------show_images: option to show images
graph_dots=False#---------------------------graph_dots plots dots (not lines)
beat_dots=False#----------------------------beat_dots marks beats with dots
blu,pur="#0072B2","#CC79A7"#----------------colorblind friendly colors
ver,sky,yel,bgr="#D55E00","#56b4e9","#F0E442","#009E73"
MAR=2#--------------------------------------MAR-normalized margin
#ESSENTIAL CONSTANTS------------------------these determine what gets plotted
f_00=100#-----------------------------------Longest quasiperiod
tmar=MAR/f_00#------------------------------margin at start and stop
Df_q=.5#------------------------------------Df_p change p-tone pitch
Df_p=0#-------------------------------------THIS SHOULD ALWAYS EQUAL ZERO
BTS=1#--------------------------------------BTS=beats_to_show

for count in range(3):#---------------------counts the three intervals plotted 
    if count==0:#-------------------------- 3/2 fifth
        interval="fifth 3-2"#--------------- p, q, f_0 defined
        p , q , f_0  = 3 , 2 , f_00#-------- 1/f_0 = quasiperiod
        q_0=q
    if count==1:#------------------------------ 7/5 tritone
        interval="tritone 7-5"
        p  , q  = 7 , 5 
        f_0=f_00*q_0/q
    if count==2:#------------------------------ 8/5 minor sixth
        interval="min_6th 8-5"
        p  , q = 8,  5 
        f_0=f_00*q_0/q
    print("count =%i: %s"%(count,interval))
    T_0=1/f_0#---------------------------------- T_0 quasiperiod                                   
    f_p= p*f_0#--------------------------------- Defines p-wave pitch
    f_q= q*f_0
    f_b=abs(p*Df_q - q*Df_p)#----------------- f_b=2*beat freq
    T_b=1/f_b#---------------------------------- Defines beat period
    #Next, I verify that the first graph requires the longest time period
    if count==0: T_M=T_b#----------------------- Check for margin error
    if T_b>T_M:sys.exit("Margin error: i=%i"%i)#-
    cents=abs( 1200*( np.log2(1+Df_p/f_p)
                -np.log2(1+Df_q/f_q) )  )#------ deviation from just (cents)
    tmax=BTS*T_b+2*tmar#------------------------ tmax=time plotted
    twopi=2*np.pi# (twopi=2*pi)----------------- define twopi=2*pi
    om_p,om_q=(f_p+Df_p)*twopi,(f_q+Df_q)*twopi#-define omega_p omega_q   
    datapoints=int(round(tmax/dt))#------------- datapoints
    print("%i datapoints"%(datapoints))#-------- not standard for CD players
    t=np.linspace(-tmar,#------------------------create numpy array for time
                  BTS*T_b+tmar,#-----------------Define timespan with margins
                  num=datapoints)   
    yp=Amp*np.cos( om_p*(t) )#------------------ yp (numpy array)
    yq=Amp*np.cos( om_q*(t))#------------------- y = yp + yq
    y=yp+yq#------------------------------------ (all are numpy arrays)
    plt.figure(figsize=(figWide,figHigh))#------ plt=variable name of plot
    i=0# --------------------------------------- initialize beat line location
    while i*T_0 < tmax:#------------------------ Plot  quasiperiod markers
        plt.axvline(x=1/f_0+i/f_0, color=pur,linewidth=barx)
        i+=1
    if graph_dots:#------------------------------ Plot signal (dots or line)
        plt.scatter(t,y,s=dotsize, color="k")#--- I used dots to verify lines
    else:
        plt.plot(t,y, color='k', linewidth=slw)#--Lines easier to see
    plt.xlim(-tmar, BTS*T_M+tmar)#----------------Sets limits on axis     
    plt.axhline(y=0, color='k', linewidth=slw/4)# Plot axis (y=0)
    a0=Amp #------------------------------------- Plot y_p and y_q signals
    a1=Amp;  offset=a0+2*a1 #-------------------- Offset signals downward
    arr=np.empty(datapoints); arr.fill(offset)#-- Numpy array speeds things up
    plt.plot(t,a1*yp-arr, color=ver, linewidth=slw)#Plot p-wave
    plt.plot(t,a1*yq-arr, color=bgr, linewidth=slw)#Plot q-wave
    i=0#----------------------------------------- Plot beats T_b=1/f_b
    while i*T_b<tmax:# Plots vertical markers for the beats
        plt.axvline(x=i*T_b, color=sky, linewidth=barb)
        i+=1
    i=0#----------------------------------------- Label quasiperiods T_0
    while i*T_0<tmax:
        s=int(i)#---------------------------------- s=string representing i
        plt.annotate(s, xy=(i/f_0,2.25*Amp),#------ Places s on t axis at i*T_0  
            xycoords='data', fontsize=fs, horizontalalignment='center',
            color=blu, verticalalignment='bottom')
        i+=1
    titlestring="%s %5.1f to %5.1f"%(interval,f_p+Df_p,f_q)
    titlestring+=" (%5.2f cents off)"%(cents)#-------Title for graph
    path2image=titlestring+".svg"#-------------------Names svg file
    path2image=os.path.join("images", path2image)#---Location for image
    plt.xlabel(titlestring,fontsize=fs,loc="left")#---TITLE
    plt.savefig(path2image,dpi=1600,bbox_inches="tight")
    if show_images: plt.show()#----------------------Show image
    print("T_b=%f for %s"%(T_b,interval))
    print("---------------- %s seconds" % (time.time() - start_time))

soundfile and playsound (pip installed packages)

edit

I discoverd that soundfile cannot create large OGG files, probably due their compression. So instead I made WAV files. And I installed playsound so I could hear the interval each time I ran the code. The latest version of the code currently resides at: