Source code for spectrum.spectrogram

from spectrum import Periodogram, pmtm
import numpy as np

__all__ = ["Spectrogram"]

[docs]class Spectrogram(object): """Simple example of spectrogram .. plot:: from spectrum import Spectrogram, dolphin_filename, readwav data, samplerate = readwav(dolphin_filename) p = Spectrogram(data, ws=128, W=4096, sampling=samplerate) p.periodogram() p.plot() .. warning:: this is a prototype and need careful checking about x/y axis """ def __init__(self, signal, ws=128, W=4096, sampling=1, channel=1): if len(signal.shape) == 1: self.signal = signal else: self.signal = signal[:,channel-1] self.W = W = ws self._start_y = 10 self.sampling = sampling self.duration = len(self.signal) / float(self.sampling)
[docs] def plot(self, filename=None, vmin=None, vmax=None, cmap='jet_r'): import pylab pylab.clf() pylab.imshow(-np.log10(self.results[self._start_y:,:]), origin="lower", aspect="auto", cmap=cmap, vmin=vmin, vmax=vmax) pylab.colorbar() # Fix xticks XMAX = float(self.results.shape[1]) # The max integer on xaxis xpos = list(range(0, int(XMAX), int(XMAX/5))) xx = [int(this*100)/100 for this in np.array(xpos) / XMAX * self.duration] pylab.xticks(xpos, xx, fontsize=16) # Fix yticks YMAX = float(self.results.shape[0]) # The max integer on xaxis ypos = list(range(0, int(YMAX), int(YMAX/5))) yy = [int(this) for this in np.array(ypos) / YMAX * self.sampling] pylab.yticks(ypos, yy, fontsize=16) #pylab.yticks([1000,2000,3000,4000], [5500,11000,16500,22000], fontsize=16) #pylab.title("%s echoes" % filename.replace(".png", ""), fontsize=25) pylab.xlabel("Time (seconds)", fontsize=25) pylab.ylabel("Frequence (Hz)", fontsize=25) pylab.tight_layout() if filename: pylab.savefig(filename)
[docs] def periodogram(self): W = self.W ws = N = int(len(self.signal)/ws) self.results = np.zeros((W*2+1, N-8)) print("Duration: %s" % self.duration) print("W: %s" % W) print("ws: %s" % ws) print("Computing %s TFs" % N) for i in range(N-8): data = self.signal[i*ws:i*ws+W] p = Periodogram(data, sampling=self.sampling, NFFT=W*4) p() self.results[:,i] = p.psd print("done")
[docs] def pmtm(self): W = self.W ws = N = int(len(self.signal) / ws) self.results = np.zeros((W+1, N-8)) for i in range(N-8): data = self.signal[i*ws:i*ws+W] a = pmtm(data, 4, NFFT=W*4, show=False) Sk = np.mean(abs(a[0].transpose())**2 * a[1], axis=1) self.results[:, i] = Sk[0:self.W+1] print(i, N) print("done")