I always wanted to work with FFTs and Spectrograms to characterise sounds and their frequencies using python.
Lately I recorded a woodpecker at his work.
To know where to look at I compared the sound with a generated sawtooth from Onlinegenerator.com. I think 17 Hz fits well.
Checking the code with the generated file worked good.
But sadly I cannot find the frequency in the DFT/FFT or spectrogram generated with the code below. No matter which settings I tried, I cannot visualise the frequencies below 20 Hz properly. No major signals cuts on or off. Only a signal at 23 Hz, which seems to be continuous.
from scipy import signal
import scipy.io.wavfile as wav
import matplotlib.pyplot as plt
import numpy as np
# Read WAV File
the_file = 'IMG_0710_short.wav'
samplerate, samples = wav.read(the_file)
#samples_left = samples[:,0]
samples = samples[:,1] # reduce to right channel
timevec = range(len(samples)) # Time vector for plot
timevec = [x / samplerate for x in timevec]
t_max =timevec[-1]
# Lineplot of Signal
dpi = 200
plt.rcParams['figure.dpi']= dpi
plt.plot(timevec,samples)
plt.title(the_file)
plt.ylabel('Amplitude [-]')
plt.xlabel('Time [sec]')
# STFT Settings
nperseg = 0.5 * samplerate # Window Size 1
noverlap = nperseg*0.95 # Overlap
nfft = 1 * samplerate # STFT 2
window = 'hann' # Window Type
timespan = [0, 2] # Calculation Window
fromm = int(len(samples)/t_max*timespan[0])
too = int(len(samples)/t_max*timespan[1])
f, t, Zxx = signal.stft(samples[fromm:too], samplerate, nperseg=nperseg, window=window, noverlap=noverlap, nfft=nfft)
t = t + timespan[0]
cmap=plt.cm.nipy_spectral
vmin = 10
vmax = 18
fig = plt.figure(figsize=(7, 5))
pcm = plt.pcolormesh(t, f, np.log(np.abs(Zxx)), cmap=cmap, vmin=vmin, vmax=vmax)
plt.ylim(0,50)
fig.colorbar(pcm)
plt.title(the_file)
plt.ylabel('Frequency [Hz]')
plt.xlabel('Time [sec]')
# FFT
n = len(samples) # length of the signal
k = np.arange(n)
T = n/samplerate
frq = k/T # two sides frequency range
frq = frq[range(int((n/2)))] # one side frequency range
Y = np.fft.fft(samples, norm='ortho')#/n # fft computing and normalization
Y = Y[range(int(n/2))]
fig, ax = plt.subplots(2, 1)
ax[0].plot(timevec,samples) # plotting the signal
ax[0].set_xlabel('Time')
ax[0].set_ylabel('Amplitude')
ax[1].plot(frq,abs(Y),'r') # plotting the spectrum
ax[1].set_xlabel('Frequency (Hz)')
ax[1].set_ylabel('|Y(freq)|')
ax[1].set_xlim(0,50)
How do I need to set the stft parameters in order to visualise the woodpecker spectrum correctly? Or is there a fundamental problem with stft here?