We're going to try to estimate frequencies in real signals now by using amplitudes of the discrete fourier transform. Before we do that, though, there's one caveat if the frequency that actually exists in our signal does not go through an integer number of periods over the audio that we measured. In that case, the amplitude actually splits over all frequency bins, though it's concentrated on the bins just below and just above the actualy frequency. So peak picking might still be OK
import numpy as np
import matplotlib.pyplot as plt
import librosa
import IPython.display as ipd
def get_dft(x):
"""
Return the non-redundant real/imaginary components
of the DFT, expressed as amplitudes of sines/cosines
involved
Parameters
----------
x: ndarray(N)
A signal
Returns
-------
cos_sum: ndarray(ceil(N/2)), sin_sums(ndarray(ceil(N/2)))
DFT cosine/sine amplitudes
"""
N = len(x)
t = np.linspace(0, 1, N+1)[0:N]
n_freqs = int(np.ceil(N/2))
f = np.fft.fft(x) # O(N log N) time
cos_sums = np.real(f)[0:n_freqs]/(N/2)
sin_sums = -np.imag(f)[0:n_freqs]/(N/2)
return cos_sums, sin_sums
Let's write some code to compute the DFT amplitudes at each frequency index and to plot them againt their corresponding frequencies in hz. We can also estimate the note involved by finding the frequency that has a maximum amplitude and by inverting the note to frequency equation as
Sometimes we get the note wrong by just looking at the peak, but we can at least visually see that the harmonics are spaced in a particular way. Looking for this spacing is in practice how people often visually figure out what the "fundamental frequency" is.
Below is an implementation of this technique that works well on me bowing the A string on my violin
x, sr = librosa.load("A4.wav")
c, s = get_dft(x)
amps = np.sqrt(c**2 + s**2)
ks = np.arange(len(c))
freqs = ks*sr/len(x)
idx = np.argmax(amps)
max_freq = freqs[idx]
# freq = 440*2^(note/12)
# note = log_2(freq/440)*12
note = np.round(np.log2(max_freq/440)*12) % 12
plt.plot(freqs, amps)
plt.scatter([freqs[idx]], [amps[idx]])
plt.title("Max Freq = {:.3f}, Note number = {}".format(max_freq, note))
plt.xlabel("Frequency Hz")
plt.xlim([0, 4000])
ipd.Audio(x, rate=sr)