Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- import numpy as np
- from matplotlib import pyplot as plt
- import scipy.io.wavfile as wav
- from numpy.lib import stride_tricks
- colormaps = "Accent, Accent_r, Blues, Blues_r, BrBG, BrBG_r, BuGn, BuGn_r, BuPu, BuPu_r, CMRmap, CMRmap_r, Dark2, Dark2_r, " \
- "GnBu, GnBu_r, Greens, Greens_r, Greys, Greys_r, OrRd, OrRd_r, Oranges, Oranges_r, PRGn, PRGn_r, Paired, Paired_r, " \
- "Pastel1, Pastel1_r, Pastel2, Pastel2_r, PiYG, PiYG_r, PuBu, PuBuGn, PuBuGn_r, PuBu_r, PuOr, PuOr_r, PuRd, PuRd_r, " \
- "Purples, Purples_r, RdBu, RdBu_r, RdGy, RdGy_r, RdPu, RdPu_r, RdYlBu, RdYlBu_r, RdYlGn, RdYlGn_r, Reds, Reds_r, " \
- "Set1, Set1_r, Set2, Set2_r, Set3, Set3_r, Spectral, Spectral_r, Vega10, Vega10_r, Vega20, Vega20_r, Vega20b, " \
- "Vega20b_r, Vega20c, Vega20c_r, Wistia, Wistia_r, YlGn, YlGnBu, YlGnBu_r, YlGn_r, YlOrBr, YlOrBr_r, YlOrRd, " \
- "YlOrRd_r, afmhot, afmhot_r, autumn, autumn_r, binary, binary_r, bone, bone_r, brg, brg_r, bwr, bwr_r, cool, " \
- "cool_r, coolwarm, coolwarm_r, copper, copper_r, cubehelix, cubehelix_r, flag, flag_r, gist_earth, gist_earth_r, " \
- "gist_gray, gist_gray_r, gist_heat, gist_heat_r, gist_ncar, gist_ncar_r, gist_rainbow, gist_rainbow_r, gist_stern, " \
- "gist_stern_r, gist_yarg, gist_yarg_r, gnuplot, gnuplot2, gnuplot2_r, gnuplot_r, gray, gray_r, hot, hot_r, hsv, " \
- "hsv_r, inferno, inferno_r, jet, jet_r, magma, magma_r, nipy_spectral, nipy_spectral_r, ocean, ocean_r, pink, " \
- "pink_r, plasma, plasma_r, prism, prism_r, rainbow, rainbow_r, seismic, seismic_r, spectral, spectral_r, spring, " \
- "spring_r, summer, summer_r, tab10, tab10_r, tab20, tab20_r, tab20b, tab20b_r, tab20c, tab20c_r, terrain, terrain_r, " \
- "viridis, viridis_r, winter, winter_r"
- """ 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
- """ plot spectrogram"""
- def plotstft(audiopath, binsize=2**10, plotpath=None, colormap="jet"):
- samplerate, samples = wav.read(audiopath)
- s = stft(samples, binsize)
- sshow, freq = logscale_spec(s, factor=2.0, sr=samplerate)
- ims = 20.*np.log10(np.abs(sshow)/10e-6) # amplitude to decibel
- timebins, freqbins = np.shape(ims)
- print("timebins: ", timebins)
- print("freqbins: ", freqbins)
- 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()
- return ims
- ims = plotstft('/home/rendier/Ptolemy/Pharos/audio/demo3.wav') # You will need to provide your own wave file here
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement