Fourier analysis of piano hammer shapes and positions where they hit the string

pianohammers.py

# Fourier analysis of piano hammer shapes on string.
# Hammers of different shapes hit piano string at different points,
# what is the frequency spectrum?

import numpy as np
import matplotlib.pyplot as plt

# line + scatter plot of multiple data sets a
def linescatter(a, titles=["", "", ""], labels=[], styles=[], fpath="", cm=plt.cm.rainbow):
    f = plt.figure(figsize=(10, 10))
    plt.rcParams.update({'font.size': 22})
    colors = cm(np.linspace(0, 1, len(a)))
    plt.rcParams['axes.prop_cycle'] = plt.cycler(color=colors)

    i = 0
    for ai in a:
        alabel = str(i+1)
        amarker = ""
        alinestyle = "-"
        al = 1.
        
        if (len(labels) > 0):
            alabel = labels[i]
            
        if (len(styles) > 0):
            if (styles[i] == 1):
                amarker = ""
                alinestyle = "-"
                al = 0.5
            if (styles[i] == 2):
                amarker = "o"
                alinestyle = ""
            if (styles[i] == 3):
                amarker = "o"
                alinestyle = ""
                al = 0.5
            
        if (len(a) > 1):
            plt.plot(ai[:, 0], ai[:, 1], label=alabel, marker=amarker, linestyle=alinestyle, alpha=al)
        else:
            plt.plot(ai[:, 0], ai[:, 1], marker=amarker, linestyle=alinestyle)
        i += 1
    
    plt.title(titles[0])
    plt.xlabel(titles[1])
    plt.ylabel(titles[2])
    
    plt.grid(True)
    plt.legend(fontsize='xx-small')

    if (len(fpath) > 0):
        f.savefig(fpath, bbox_inches='tight') # save to file

    plt.show()
    plt.close()

# generate hammer shapes
def hammers(n, l, w, hx):
    x = np.linspace(-l/2., l/2., n)
    wi = abs(x-hx) <= w
    h1 = np.zeros(n)
    h2 = np.zeros(n)
    
    h1[wi] = 1. # flat hammer
    h2[wi] = np.sqrt(1. - ((x[wi]-hx)/w)**2) # circular hammer
    h3 = np.exp(- 4 * np.log(2.) * (x-hx)**2 / w**2) # Gaussian hammer

    h1[0] = 0.
    h1[-1] = 0.
    h2[0] = 0.
    h2[-1] = 0.
    h3[0] = 0.
    h3[-1] = 0.
    
    linescatter([np.array([x, h1]).T, np.array([x, h2]).T, np.array([x, h3]).T], titles=["Hammer shapes", "x", "y"])

    return h1, h2, h3

# calculate spectrum for hammers
def spectrum(n, l, h):
    bw = np.zeros(n)
    bn = np.zeros(n)

    for ni in range(n): # calculate Fourier coefficients
        fi = 0.
        for nxi in range(n):
            x = (nxi / n - 0.5) * l
            fi += h[nxi] * np.cos((ni+1) * np.pi * x / l) * l / n
            
        bw[ni] = ni / n
        bn[ni] = 2. / l * fi

    return np.array([bw, bn]).T

# calculate spectra for all hammers with hx offset
def spectra(n, l, hx):
    h1, h2, h3 = hammers(n, l, w, hx)

    f1 = spectrum(n, l, h1)
    f2 = spectrum(n, l, h2)
    f3 = spectrum(n, l, h3)

    return f1, f2, f3

n = 1000 # number of points
l = 1. # length of string
w = l / 50. # half width of hammer

f1, f2, f3 = spectra(n, l, 0.) # centered
fh1, fh2, fh3 = spectra(n, l, l / 4. - w) # halfway
fe1, fe2, fe3 = spectra(n, l, l / 2. - w) # end

linescatter([np.abs(f1)], titles=["Spectrum", "$\omega$", "$|A|$"])
linescatter([np.abs(f2)], titles=["Spectrum", "$\omega$", "$|A|$"])
linescatter([np.abs(f3)], titles=["Spectrum", "$\omega$", "$|A|$"])

f1[:, 1] = np.abs(f1[:, 1])**2
fh1[:, 1] = np.abs(fh1[:, 1])**2
fe1[:, 1] = np.abs(fe1[:, 1])**2
linescatter([f1[0:200], fh1[0:200], fe1[0:200]], titles=["Spectrum", "$\omega$", "$|A|^2$"], labels=["Centered", "Half", "End"], styles=[1, 1, 1])

print("h1 |A|^2:", f1[0, 1], fh1[0, 1], fe1[0, 1])

f2[:, 1] = np.abs(f2[:, 1])**2
fh2[:, 1] = np.abs(fh2[:, 1])**2
fe2[:, 1] = np.abs(fe2[:, 1])**2
linescatter([f2[0:200], fh2[0:200], fe2[0:200]], titles=["Spectrum", "$\omega$", "$|A|^2$"], labels=["Centered", "Half", "End"], styles=[1, 1, 1])

print("h2 |A|^2:", f2[0, 1], fh2[0, 1], fe2[0, 1])

f3[:, 1] = np.abs(f3[:, 1])**2
fh3[:, 1] = np.abs(fh3[:, 1])**2
fe3[:, 1] = np.abs(fe3[:, 1])**2
linescatter([f3[0:200], fh3[0:200], fe3[0:200]], titles=["Spectrum", "$\omega$", "$|A|^2$"], labels=["Centered", "Half", "End"], styles=[1, 1, 1])

print("h3 |A|^2:", f3[0, 1], fh3[0, 1], fe3[0, 1])