# -*- coding: utf-8 -*- # ------------------------------------------------------------------ # Filename: spectrogram.py # Purpose: Plotting spectrogram of Seismograms. # Author: Christian Sippl, Moritz Beyreuther # Email: sippl@geophysik.uni-muenchen.de # # Copyright (C) 2008-2012 Christian Sippl # Modified by Utpal Kumar, 2021/05 @IESAS # -------------------------------------------------------------------- """ Plotting spectrogram of seismograms. :copyright: The ObsPy Development Team (devs@obspy.org) :license: GNU Lesser General Public License, Version 3 (https://www.gnu.org/copyleft/lesser.html) """ from __future__ import (absolute_import, division, print_function, unicode_literals) from future.builtins import * # NOQA @UnusedWildImport import math import numpy as np from matplotlib import mlab from matplotlib.colors import Normalize from obspy.imaging.cm import obspy_sequential def _nearest_pow_2(x): """ Find power of two nearest to x >>> _nearest_pow_2(3) 2.0 >>> _nearest_pow_2(15) 16.0 :type x: float :param x: Number :rtype: Int :return: Nearest power of 2 to x """ a = math.pow(2, math.ceil(np.log2(x))) b = math.pow(2, math.floor(np.log2(x))) if abs(a - x) < abs(b - x): return a else: return b def compute_spectrogram(data, samp_rate, per_lap=0.9, wlen=None, dbscale=False, mult=8.0): """ Computes and plots spectrogram of the input data. :param data: Input data :type samp_rate: float :param samp_rate: Samplerate in Hz :type per_lap: float :param per_lap: Percentage of overlap of sliding window, ranging from 0 to 1. High overlaps take a long time to compute. :type wlen: int or float :param wlen: Window length for fft in seconds. If this parameter is too small, the calculation will take forever. If None, it defaults to (samp_rate/100.0). :type dbscale: bool :param dbscale: If True 10 * log10 of color values is taken, if False the sqrt is taken. :type mult: float :param mult: Pad zeros to length mult * wlen. This will make the spectrogram smoother. """ # enforce float for samp_rate samp_rate = float(samp_rate) # set wlen from samp_rate if not specified otherwise if not wlen: wlen = samp_rate / 100. npts = len(data) # nfft needs to be an integer, otherwise a deprecation will be raised # XXX add condition for too many windows => calculation takes for ever nfft = int(_nearest_pow_2(wlen * samp_rate)) if nfft > npts: nfft = int(_nearest_pow_2(npts / 8.0)) if mult is not None: mult = int(_nearest_pow_2(mult)) mult = mult * nfft nlap = int(nfft * float(per_lap)) data = data - data.mean() end = npts / samp_rate # Here we call not plt.specgram as this already produces a plot # matplotlib.mlab.specgram should be faster as it computes only the # arrays # XXX mlab.specgram uses fft, would be better and faster use rfft specgram, freq, time = mlab.specgram(data, Fs=samp_rate, NFFT=nfft, pad_to=mult, noverlap=nlap) # db scale and remove zero/offset for amplitude if dbscale: specgram = 10 * np.log10(specgram[1:, :]) else: specgram = np.sqrt(specgram[1:, :]) freq = freq[1:] # calculate half bin width halfbin_time = (time[1] - time[0]) / 2.0 halfbin_freq = (freq[1] - freq[0]) / 2.0 # pcolor expects one bin more at the right end freq = np.concatenate((freq, [freq[-1] + 2 * halfbin_freq])) time = np.concatenate((time, [time[-1] + 2 * halfbin_time])) # center bin time -= halfbin_time freq -= halfbin_freq return specgram, freq, time