plot power spectrum of a signal with wavelet transform

946 Views Asked by At

Using this dataset https://philharmonia.co.uk/resources/sound-samples/ I'm trying to plot the power spectrum of note played by a specific instrument.

I'm using librosa to load the audio file and get some information with this code

import librosa
import numpy as np
import matplotlib.pyplot as plt
import pywt

y, sr = librosa.load(file_path)
duration = librosa.get_duration(y=y, sr=sr)
delta_t = duration / len(y)
t0=0
time = np.arange(0, len(y)) * delta_t + t0

I'm also following this https://ataspinar.com/2018/12/21/a-guide-for-using-the-wavelet-transform-in-machine-learning/ guide to plot the power spectrum and I'm using pywavelet library.

The problem that I have with this code is RuntimeWarning: divide by zero encountered in log2 and the plot is not shown.

def plot_wavelet(time, signal, scales, 
                 waveletname = 'cmor', 
                 cmap = plt.cm.seismic, 
                 title = 'Wavelet Transform (Power Spectrum) of signal', 
                 ylabel = 'Period (years)', 
                 xlabel = 'Time'):
    
    dt = time[1] - time[0]
    print("dt ", dt)
    [coefficients, frequencies] = pywt.cwt(signal, scales, waveletname, dt)
    power = (abs(coefficients)) ** 2
    period = 1. / frequencies
    levels = [0.0625, 0.125, 0.25, 0.5, 1, 2, 4, 8]
    contourlevels = np.log2(levels)
    
    fig, ax = plt.subplots(figsize=(15, 10))
    im = ax.contourf(time, np.log2(period), np.log2(power), contourlevels, extend='both',cmap=cmap)
    

    ax.set_title(title, fontsize=20)
    ax.set_ylabel(ylabel, fontsize=18)
    ax.set_xlabel(xlabel, fontsize=18)
    
    yticks = 2**np.arange(np.ceil(np.log2(period.min())), np.ceil(np.log2(period.max())))
    ax.set_yticks(np.log2(yticks))
    ax.set_yticklabels(yticks)
    ax.invert_yaxis()
    ylim = ax.get_ylim()
    ax.set_ylim(ylim[0], -1)
    
    cbar_ax = fig.add_axes([0.95, 0.5, 0.03, 0.25])
    fig.colorbar(im, cax=cbar_ax, orientation="vertical")
    plt.show()

scales = np.arange(1, 128)
plot_wavelet(time=time, signal=y, scales=scales, waveletname='gaus5')

Noting that some values in the power array are at -inf. How can I solve?

0

There are 0 best solutions below