In a previous article about the Fourier Series calculation we have illustrated how different numbers of harmonics approximate artificial periodic functions. The present post applies the results to the analysis of instrument sounds, namely sounds of a saxophone.

When playing a stationary tone on a saxophone, we hear a constant sound. Hence, we can assume its waveform is periodic, since we could start to listen to the tone at any time and would still hear the same tone. So, the waveform needs to repeat itself over and over again. In this case, it should be possible to expand the waveform into sines and cosines of harmonic frequencies and reconstruct the original signal from them.

We want to verify this is assumption with this post. Let us start with functions to calculate the Fourier series, `fourierSeries`

and for reconstructing a signal from its Fourier series coefficients, `reconstruct`

.

```
def fourierSeries(period, N):
"""Calculate the Fourier series coefficients up to the Nth harmonic"""
result = []
T = len(period)
t = np.arange(T)
for n in range(N+1):
an = 2/T*(period * np.cos(2*np.pi*n*t/T)).sum()
bn = 2/T*(period * np.sin(2*np.pi*n*t/T)).sum()
result.append((an, bn))
return np.array(result)
```

```
def reconstruct(P, anbn):
"""Sum up sines and cosines according to the coefficients to
produce a reconstruction of the original waveform"""
result = 0
t = np.arange(P)
for n, (a, b) in enumerate(anbn):
if n == 0:
a = a/2
result = result + a*np.cos(2*np.pi*n*t/P) + b * np.sin(2*np.pi*n*t/P)
return result
```

First of all, let's load three sounds of saxophone sounds from the internet:

```
def loadSound(url):
R = requests.get(url)
rate, data = wavfile.read(BytesIO(R.content))
data = data * 1.0 / (abs(data).max())
if len(data.shape) > 1: # for stereo data, use only first channel
data = data[:,0]
return data, rate
# the URLs of the wav files
urls = ["http://www.colorado.edu/physics/phys4830/phys4830_fa01/sounds/nash1.wav",
"http://cd.textfiles.com/sbsw/INSTRMNT/TENORSA.WAV",
"http://www.gentrythomas.com/downloads/GASPAR/Programs/Fruity%20Loop%20Packs%20and%20Drum%20Kits/extra%20Snyth/C3sax.wav"]
sounds = []
for url in urls:
sound, rate = loadSound(url)
sounds.append((sound, rate))
```

```
# Utility function two display two audios side by side in the notebook
def audioSideBySide(name1, audio1, name2, audio2):
text = '
```%s %s %s %s

' % (name1, name2, audio1._repr_html_(), audio2._repr_html_())
display(HTML(text))

Let's listen to the downloaded sounds of the played saxophone:

```
for n, s in enumerate(sounds):
display(HTML("Saxophone sound %d" % (n+1) + Audio(data=s[0], rate=s[1])._repr_html_()))
```

```
t = np.arange(len(sounds[0][0])) / sounds[0][1]
plt.subplot(211)
plt.plot(t, sounds[0][0])
plt.subplot(212)
plt.plot(t, sounds[0][0])
```

```
def extractPeriod(data, rate, t_start, t_end):
t = np.arange(0,len(data))/rate
plt.plot(t, data)
duration = t_end - t_start
plt.xlabel('$t$'); plt.ylabel('$x(t)$');
sample_start = int(t_start * rate)
sample_end = int(t_end*rate)
period = data[sample_start:sample_end]
audioSideBySide("Original", Audio(data=data,rate=rate),
"Extracted period", Audio(np.tile(period, int(1/duration)), rate=rate))
return period, rate
periods = []
rates = []
# The manually found start and end time of one period of the signal
periodBounds = [(0.8005,0.8034),
(0.8044,0.8075),
(0.20555,0.2132)]
for n, (S, P) in enumerate(zip(sax, periodBounds)):
plt.subplot(len(sax), 1, n+1)
period, rate = extractPeriod(S[0], S[1], P[0], P[1])
periods.append(period); rates.append(rate)
```