Hilbert Huang Spektrum plotten

mit matplotlib, NumPy, pandas, SciPy, SymPy und weiteren mathematischen Programmbibliotheken.
Antworten
Philuup
User
Beiträge: 1
Registriert: Freitag 27. April 2018, 09:02

Freitag 27. April 2018, 09:15

Hallo,

ich benutze das PyEMD package für Python 3.6 und möchte ein 2D Hilbert Spektrum (x: Time, y: Frequency, Colorbar: Amplitude) und ein Hilbert Marginal Spektrum (x: Frequency, y: Amplitude) plotten.

Ich kann durch die Empirische Modenzerlegung (EMD) und die Hilbert-Huang-Transformation (HHT) die momentane Phase, Frequenz und Amplitude der einzelnen intrinsischen Modenfunktionen (IMFs) berechnen, wie im folgenden Code gezeigt:

Code: Alles auswählen

import numpy as np
import pylab as plt
from scipy.signal import hilbert

from PyEMD import EMD


def instant_phase(imfs):
    """Extract analytical signal through Hilbert Transform."""
    analytic_signal = hilbert(imfs)  # Apply Hilbert transform to each row
    phase = np.unwrap(np.angle(analytic_signal))  # Compute angle between img and real
    return phase 

def instant_amplitude(imfs):
    """Extract analytical signal through Hilbert Transform."""
    analytic_signal = hilbert(imfs)  # Apply Hilbert transform to each row
    analytic_amplitude = np.abs(analytic_signal) ## Compute amplitude
    return analytic_amplitude


# Define signal
t = np.linspace(0, 1, 200)
dt = t[1]-t[0]

sin = lambda x, p: np.sin(2*np.pi*x*t+p)
S = 3*sin(18, 0.2)*(t-0.2)**2
S += 5*sin(11, 2.7)
S += 3*sin(14, 1.6)
S += 1*np.sin(4*2*np.pi*(t-0.8)**2)
S += t**2.1 -t

# Compute IMFs with EMD
emd = EMD()
imfs = emd(S, t)

# Extract instantaneous phases, frequencies and amplitude using Hilbert transform 
instant_phases = instant_phase(imfs)
instant_freqs = np.diff(instant_phases)/(2*np.pi*dt)
instant_amplitude = instant_amplitude(imfs)

Jetzt hänge ich aber fest und kann nicht das gesamte Spektrum berechnen bzw plotten, sowie das Marginal Spektrum nach folgenden Formeln:

https://en.wikipedia.org/wiki/Hilbert_spectrum

Hat das jemand schonmal gemacht und kann mir helfen?

Vielen Dank!
Antworten