Styjun

I'm processing wav files for amplitude and frequency analysis with FFT, but I am having trouble getting the data out to csv in a time series format.

Using @Beginner's answer heavily from this post: How to convert a .wav file to a spectrogram in python3, I'm able to get the spectrogram output in an image. I'm trying to simplify that somewhat to get to a text output in csv format, but I'm not seeing how to do so. The outcome I'm hoping to achieve would look something like the following:

time_in_ms, amplitude_in_dB, freq_in_kHz

.001, -115, 1

.002, -110, 2

.003, 20, 200

...

19000, 20, 200

For my testing, I have been using http://soundbible.com/2123-40-Smith-Wesson-8x.html, (Notes: I simplified the wav down to a single channel and removed metadata w/ Audacity to get it to work.)

Heavy props to @Beginner for 99.9% of the following, anything nonsensical is surely mine.

import numpy as np
from matplotlib import pyplot as plt
import scipy.io.wavfile as wav
from numpy.lib import stride_tricks

filepath = "40sw3.wav"

""" short time fourier transform of audio signal """
def stft(sig, frameSize, overlapFac=0.5, window=np.hanning):
 win = window(frameSize)
 hopSize = int(frameSize - np.floor(overlapFac * frameSize))

 # zeros at beginning (thus center of 1st window should be for sample nr. 0) 
 samples = np.append(np.zeros(int(np.floor(frameSize/2.0))), sig) 
 # cols for windowing
 cols = np.ceil( (len(samples) - frameSize) / float(hopSize)) + 1
 # zeros at end (thus samples can be fully covered by frames)
 samples = np.append(samples, np.zeros(frameSize))

 frames = stride_tricks.as_strided(samples, shape=(int(cols), frameSize), strides=(samples.strides[0]*hopSize, samples.strides[0])).copy()
 frames *= win

 return np.fft.rfft(frames) 

""" scale frequency axis logarithmically """ 
def logscale_spec(spec, sr=44100, factor=20.):
 timebins, freqbins = np.shape(spec)

 scale = np.linspace(0, 1, freqbins) ** factor
 scale *= (freqbins-1)/max(scale)
 scale = np.unique(np.round(scale))

 # create spectrogram with new freq bins
 newspec = np.complex128(np.zeros([timebins, len(scale)]))
 for i in range(0, len(scale)): 
 if i == len(scale)-1:
 newspec[:,i] = np.sum(spec[:,int(scale[i]):], axis=1)
 else: 
 newspec[:,i] = np.sum(spec[:,int(scale[i]):int(scale[i+1])], axis=1)

 # list center freq of bins
 allfreqs = np.abs(np.fft.fftfreq(freqbins*2, 1./sr)[:freqbins+1])
 freqs = []
 for i in range(0, len(scale)):
 if i == len(scale)-1:
 freqs += [np.mean(allfreqs[int(scale[i]):])]
 else:
 freqs += [np.mean(allfreqs[int(scale[i]):int(scale[i+1])])]

 return newspec, freqs

""" compute spectrogram """
def compute_stft(audiopath, binsize=2**10):
 samplerate, samples = wav.read(audiopath)
 s = stft(samples, binsize)
 sshow, freq = logscale_spec(s, factor=1.0, sr=samplerate)
 ims = 20.*np.log10(np.abs(sshow)/10e-6) # amplitude to decibel
 return ims, samples, samplerate, freq

""" plot spectrogram """
def plot_stft(ims, samples, samplerate, freq, binsize=2**10, plotpath=None, colormap="jet"):
 timebins, freqbins = np.shape(ims)

 plt.figure(figsize=(15, 7.5))
 plt.imshow(np.transpose(ims), origin="lower", aspect="auto", cmap=colormap, interpolation="none")
 plt.colorbar()

 plt.xlabel("time (s)")
 plt.ylabel("frequency (hz)")
 plt.xlim([0, timebins-1])
 plt.ylim([0, freqbins])

 xlocs = np.float32(np.linspace(0, timebins-1, 5))
 plt.xticks(xlocs, ["%.02f" % l for l in ((xlocs*len(samples)/timebins)+(0.5*binsize))/samplerate])
 ylocs = np.int16(np.round(np.linspace(0, freqbins-1, 10)))
 plt.yticks(ylocs, ["%.02f" % freq[i] for i in ylocs])

 if plotpath:
 plt.savefig(plotpath, bbox_inches="tight")
 else:
 plt.show()

 plt.clf()

"" HERE IS WHERE I'm ATTEMPTING TO GET IT OUT TO TXT """
ims, samples, samplerate, freq = compute_stft(filepath)

""" Print lengths """
print('ims len:', len(ims))
print('samples len:', len(samples))
print('samplerate:', samplerate)
print('freq len:', len(freq))

""" Write values to files """
np.savetxt(filepath + '-ims.txt', ims, delimiter=', ', newline='n', header='ims')
np.savetxt(filepath + '-samples.txt', samples, delimiter=', ', newline='n', header='samples')
np.savetxt(filepath + '-frequencies.txt', freq, delimiter=', ', newline='n', header='frequencies')

In terms of values out, the file I'm analyzing is approx 19.1 seconds long and the sample rate is 44100, so I’d expect to have about 842k values for any given variable. But I'm not seeing what I expected. Instead here is what I see:

freqs comes out with just a handful of values, 512 and while they appear to be correct range for expected frequency, they are ordered least to greatest, not in time series like I expected. The 512 values, I assume, is the "fast" in FFT, basically down-sampled...

ims, appears to be amplitude, but values seem too high, although sample size is correct. Should be seeing -50 up to ~240dB.

samples . . . not sure.

In short, can someone advise on how I'd get the FFT out to a text file with time, amp, and freq values for the entire sample set? Is savetxt the correct route, or is there a better way? This code can certainly be used to make a great spectrogram, but how can I just get out the data?

Your output format is too limiting, as the audio spectrum at any interval in time usually contains a range of frequencies. e.g the FFT of a 1024 samples will contain 512 frequency bins for one window of time or time step, each with an amplitude. If you want a time step of one millisecond, then you will have to offset the window of samples you feed each STFT to center the window at that point in your sample vector. Although with an FFT about 23 milliseconds long, that will involve a high overlap of windows. You could use shorter windows, but the time-frequency trade-off will result in proportionately less frequency resolution.