Commit e900a896 authored by Gregor Moenke's avatar Gregor Moenke

Merge branch 'develop' into 'master'

Fourier estimates

See merge request !37
parents c6863ed0 e759919d
### pyBOAT 0.8.2
- Added time averaging of Wavelet spectra <-> Fourier estimates
- Added Fourier distribution for ensembles
- Reworked FFT visualizations
### pyBOAT 0.8.17
- Fixed crashes during batch analysis with thresholded ridges
......
This source diff could not be displayed because it is too large. You can view the blob instead.
......@@ -3,7 +3,7 @@
import sys,os
import argparse
__version__ = '0.8.17'
__version__ = '0.8.18'
# the object oriented API
from .api import WAnalyzer
......
......@@ -23,7 +23,7 @@ def average_power_distribution(ridge_results, signal_ids = None, exclude_coi = F
about the DataFrame structure
signal_ids : sequence, optional
labels of the analyzed signals, of not given
labels of the analyzed signals, if not given
a numeric sequence of len(ridge_results)
will be used as labels
......@@ -74,8 +74,8 @@ def get_ensemble_dynamics(ridge_results):
'''
Aggregate all the ridge-readouts (period, amplitude and phase)
of a ridge-analyzed signal ensemble and return time-continuous median and
quartiles (as for a time-continuous box plot).
In the case of phase return the 1st order parameter
quartiles (as in a time-continuous box plot).
In the case of phase return the 1st order parameter (length of resultant vector)
as a phase coherence measure over time.
Signals/ridge-readouts of unequal length in time are Ok!
......@@ -91,7 +91,7 @@ def get_ensemble_dynamics(ridge_results):
Returns
--------
A tuple holding 5 data frames, one summary statistic over time
A tuple holding 4 data frames, one summary statistic over time
for period, amplitude, power and phase each.
'''
......@@ -127,6 +127,8 @@ def get_ensemble_dynamics(ridge_results):
return periods_mq1q3, amplitudes_mq1q3, powers_mq1q3, phases_R
if __name__ == '__main__':
'''
......@@ -135,7 +137,7 @@ if __name__ == '__main__':
'''
from pyboat import WAnalyzer, ssg
from pyboat.plotting import plot_ensemble_dynamics, plot_power_distribution
from pyboat.plotting import ensemble_dynamics, power_distribution, Fourier_distribution
# set up analyzing instance
periods = np.linspace(5, 60, 100)
......@@ -144,11 +146,12 @@ if __name__ == '__main__':
# create a bunch of chirp signals
Nsignals = 50 # times 2
T1 = 30 # initial period
Tstart = 30 # initial period
Tend = 50 # period at the end
Nt = 500 # number of samples per signal
signals = [
ssg.create_noisy_chirp(T1 = T1, T2 = T, Nt = Nt, eps = 1)
for T in np.linspace(T1, 50, Nsignals) ]
ssg.create_noisy_chirp(T1 = Tstart, T2 = T, Nt = Nt, eps = 1)
for T in np.linspace(Tstart, Tend, Nsignals) ]
# add the same amount of pure noise
noisy_ones = [
......@@ -159,16 +162,25 @@ if __name__ == '__main__':
# get the the individual ridge readouts
ridge_results = []
for signal in signals:
# store the individual time averaged Wavelet spectra
df_fouriers = pd.DataFrame(columns = np.arange(len(signals)), index = wAn.periods)
for i,signal in enumerate(signals):
wAn.compute_spectrum(signal, sinc_detrend = False, do_plot = False)
rd = wAn.get_maxRidge(smoothing_wsize = 11)
ridge_results.append(rd)
df_fouriers[i] = wAn.get_averaged_spectrum()
# the time-averaged power distribution
powers_series = average_power_distribution( ridge_results, signal_ids = None )
plot_power_distribution(powers_series)
power_distribution(powers_series)
# keeping the pure noise signal out
res = get_ensemble_dynamics(ridge_results[:Nsignals])
plot_ensemble_dynamics(*res)
ensemble_dynamics(*res)
# the Fourier power distribution
Fourier_distribution(df_fouriers)
......@@ -36,7 +36,7 @@ CMAP = "YlGnBu_r"
# CMAP = 'cividis'
# CMAP = 'magma'
# --- max size of signal to plot also the sample points
# --- max size of signal to plot also the sample points as explicit o's
Nmax = 250
# --- define line widths ---
......@@ -119,7 +119,7 @@ def draw_detrended(ax, time_vector, detrended):
return ax2
# --- Fourier Spectrum ------------------------------------------------
# --- Fourier Spectrum --------------------------
def mk_Fourier_ax(fig, time_unit="a.u.", show_periods=False):
......@@ -141,18 +141,12 @@ def mk_Fourier_ax(fig, time_unit="a.u.", show_periods=False):
return ax
def Fourier_spec(ax, fft_freqs, fft_power, show_periods=False):
if show_periods:
# period view, omit the last bin 2/(N*dt)
# skip 0-frequency
if len(fft_freqs) < 1000:
ax.vlines(
1 / fft_freqs[1:],
0,
......@@ -161,7 +155,6 @@ def Fourier_spec(ax, fft_freqs, fft_power, show_periods=False):
alpha=0.8,
color=FOURIER_COLOR,
)
# plot differently for very long signals
else:
ax.plot(
......@@ -172,11 +165,8 @@ def Fourier_spec(ax, fft_freqs, fft_power, show_periods=False):
alpha=0.8,
color=FOURIER_COLOR,
)
else:
if len(fft_freqs) < 1000:
# frequency view
ax.vlines(
fft_freqs,
......@@ -189,7 +179,37 @@ def Fourier_spec(ax, fft_freqs, fft_power, show_periods=False):
else:
ax.plot(fft_freqs, fft_power, "-", lw=1.5, alpha=0.8, color=FOURIER_COLOR)
# --- time averaged Wavelet spectrum -> Fourier estimate
def averaged_Wspec(averaged_Wspec, periods, time_unit = 'a.u', fig = None):
'''
Plots the time averaged Wavelet spectrum, which is a good Fourier estimate.
Parameters
----------
averaged_Wspec : sequence, holding the average power for each period
periods : sequence, the periods of the original Wavelet transform
time_unit : str, the time unit label
fig : matplotlib figure instance, a new figure is created per default
'''
if fig is None:
fig = ppl.figure(figsize = (5, 3.2))
ax = fig.subplots()
ax.set_ylabel("Power (wnp)", fontsize=label_size)
ax.set_xlabel(f"Period ({time_unit})", fontsize=label_size)
ax.plot(periods, averaged_Wspec, lw = SIGNAL_LW, color = FOURIER_COLOR)
ax.fill_between(periods, 0, averaged_Wspec, color = FOURIER_COLOR, alpha = 0.3)
# --- Wavelet spectrum ------
......@@ -422,7 +442,7 @@ def Rice_rule(samples):
Nbins = int(2 * pow(len(samples), 1/3))
return Nbins
def plot_power_distribution(powers, kde = True, fig = None):
def power_distribution(powers, kde = True, fig = None):
'''
Histogram (bin-counts) of Wavelet powers, intended
......@@ -433,6 +453,7 @@ def plot_power_distribution(powers, kde = True, fig = None):
----------
powers : a sequence of floats
kde : bool, Gaussian kde
fig : matplotlib figure instance
'''
......@@ -446,7 +467,7 @@ def plot_power_distribution(powers, kde = True, fig = None):
ax = fig.subplots()
ax.set_ylabel("Counts", fontsize=label_size)
ax.set_xlabel("Average Wavelet Power", fontsize=label_size)
ax.set_xlabel("Average Ridge Power", fontsize=label_size)
ax.tick_params(axis="both", labelsize=tick_label_size)
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
......@@ -494,7 +515,7 @@ def plot_power_distribution(powers, kde = True, fig = None):
fig.tight_layout()
#fig.subplots_adjust(bottom = 0.22, left = 0.15)
def plot_ensemble_dynamics(
def ensemble_dynamics(
periods,
amplitudes,
powers,
......@@ -580,3 +601,53 @@ def plot_ensemble_dynamics(
fig.subplots_adjust(wspace=0.3, left=0.11, top=0.98, right=0.97, bottom=0.14)
#fig.subplots_adjust(bottom = 0.1, left = 0.2, top = 0.95, hspace = 0.1)
def Fourier_distribution(
df_fouriers,
time_unit = 'a.u.',
label = None,
color = FOURIER_COLOR,
fig = None):
'''
Plots the median and quartiles of a distribution of
Fourier power spectra. Typical use case are the
time averaged Wavelet spectra of a population.
Parameters
----------
df_fouriers : DataFrame, columns hold the individual Fourier power spectra,
index holds the periods
time_unit : str, the time unit label
fig : matplotlib figure instance
'''
periods = df_fouriers.index
med = df_fouriers.median(axis = 1)
q1 = df_fouriers.quantile(q = 0.25, axis = 1)
q3 = df_fouriers.quantile(q = 0.75, axis = 1)
if fig is None:
fig = ppl.figure(figsize = (5.6, 4.5))
ax = fig.subplots()
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.grid(axis = 'y')
ax.tick_params(axis="both", labelsize=tick_label_size)
ax.set_xlabel(f'Period ({time_unit})', fontsize = label_size)
ax.set_ylabel('Power (wnp)', fontsize = label_size)
ax.plot(periods, med, lw = 2.5, color = color, label = label)
ax.fill_between(periods, q1, q3, color = color, alpha = 0.3)
if label:
ax.legend()
fig.tight_layout()
......@@ -31,7 +31,7 @@ def create_chirp(T1, T2, Nt):
def create_noisy_chirp(T1, T2, Nt, eps, alpha = 0):
'''
Creates a clean chirp signal,
Creates a noisy chirp signal,
sweeping linearly through the
frequencies:
......
......@@ -3,7 +3,7 @@ import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from PyQt5.QtWidgets import QCheckBox, QTableView, QComboBox, QFileDialog, QAction, QMainWindow, QApplication, QLabel, QLineEdit, QPushButton, QMessageBox, QSizePolicy, QWidget, QVBoxLayout, QHBoxLayout, QDialog, QGroupBox, QFormLayout, QGridLayout, QTabWidget, QTableWidget
from PyQt5.QtWidgets import QCheckBox, QTableView, QComboBox, QFileDialog, QAction, QMainWindow, QApplication, QLabel, QLineEdit, QPushButton, QMessageBox, QSizePolicy, QWidget, QVBoxLayout, QHBoxLayout, QDialog, QGroupBox, QFormLayout, QGridLayout, QTabWidget, QTableWidget, QSpacerItem
from PyQt5.QtGui import QDoubleValidator, QIntValidator, QScreen
from PyQt5.QtCore import Qt, pyqtSignal
......@@ -11,7 +11,7 @@ from matplotlib.backends.backend_qt5agg import NavigationToolbar2QT as Navigatio
from matplotlib.backends.backend_qt5agg import FigureCanvasQTAgg as FigureCanvas
from matplotlib.figure import Figure
from pyboat.ui.util import MessageWindow, posfloatV, posintV
from pyboat.ui.util import MessageWindow, posfloatV, posintV, mkGenericCanvas
from pyboat import core
from pyboat import plotting as pl
......@@ -115,10 +115,9 @@ class WaveletAnalyzer(QWidget):
# no anneal parameters yet
self.anneal_pars = None
#=============Compute Spectrum=============================================
#=============Compute Spectrum================================================
self.modulus, self.wlet = core.compute_spectrum(self.signal, dt, self.periods)
#==========================================================================
#=============================================================================
# Wavelet ridge-readout results
self.ResultWindows = {}
......@@ -150,7 +149,7 @@ class WaveletAnalyzer(QWidget):
# --- Spectrum plotting options ---
spectrum_opt_box = QGroupBox("Plotting Options")
spectrum_opt_box = QGroupBox("Spectrum Plotting Options")
spectrum_opt_box.setSizePolicy(QSizePolicy.Maximum,QSizePolicy.Maximum)
spectrum_opt_layout = QHBoxLayout()
......@@ -162,7 +161,7 @@ class WaveletAnalyzer(QWidget):
self.pmax_edit = QLineEdit()
self.pmax_edit.setToolTip('Rescales the color map of the spectrum')
self.pmax_edit.setMaximumWidth(60)
self.pmax_edit.setMaximumWidth(80)
self.pmax_edit.setValidator(posfloatV)
# retrieve initial power value, axs[1] is the spectrum
......@@ -186,6 +185,22 @@ class WaveletAnalyzer(QWidget):
spectrum_opt_layout.addStretch(0)
spectrum_opt_layout.addWidget(self.cb_coi)
# --- Time average -> asymptotic Fourier ---
time_av_box = QGroupBox("Time Averaging")
estimFourierButton = QPushButton('Estimate Fourier', self)
estimFourierButton.clicked.connect(self.ini_average_spec)
estimFourierButton.setToolTip('Calculates time averaged Wavelet power spectrum')
time_av_layout = QHBoxLayout()
#time_av_layout.addStretch()
time_av_layout.addWidget(estimFourierButton)
#time_av_layout.addStretch()
time_av_box.setLayout(time_av_layout)
# --- Ridge detection and options --
......@@ -210,19 +225,19 @@ class WaveletAnalyzer(QWidget):
power_thresh_edit = QLineEdit()
power_thresh_edit.setToolTip('Sets the minimal power value required to be part of the ridge')
power_thresh_edit.setSizePolicy(QSizePolicy.Fixed, QSizePolicy.Fixed)
power_thresh_edit.setMinimumSize(50,0)
power_thresh_edit.setMinimumSize(60,0)
power_thresh_edit.setValidator(posfloatV)
smooth_label = QLabel("Ridge Smoothing:")
ridge_smooth_edit = QLineEdit()
ridge_smooth_edit.setMinimumSize(50,0)
ridge_smooth_edit.setMinimumSize(60,0)
ridge_smooth_edit.setSizePolicy(QSizePolicy.Fixed, QSizePolicy.Fixed)
ridge_smooth_edit.setValidator( QIntValidator(bottom = 3, top = len(self.signal)) )
# Plot Results
plotResultsButton = QPushButton('Plot Ridge Readout', self)
plotResultsButton.setSizePolicy(QSizePolicy.Fixed, QSizePolicy.Fixed)
# plotResultsButton.setSizePolicy(QSizePolicy.Fixed, QSizePolicy.Fixed)
plotResultsButton.clicked.connect(self.ini_plot_readout)
ridge_opt_layout.addWidget(maxRidgeButton,0,0,1,1)
......@@ -242,13 +257,25 @@ class WaveletAnalyzer(QWidget):
rtool_layout.addWidget(ridge_opt_box)
rtool_layout.addStretch(0)
# put everything together
main_layout = QVBoxLayout()
main_layout.setSpacing(0)
main_layout.addWidget(self.wCanvas)
main_layout.addWidget(ntb)
main_layout.addWidget(spectrum_opt_box)
main_layout.addWidget(ridge_opt_box)
# --- put everything together ---
main_layout = QGridLayout()
#main_layout.setSpacing(0)
main_layout.addWidget(self.wCanvas, 0, 0, 4, 5)
main_layout.addWidget(ntb, 5, 0,1,5)
main_layout.addWidget(spectrum_opt_box, 6, 0, 1, 3)
main_layout.addWidget(time_av_box, 6, 3, 1, 1)
main_layout.addItem(QSpacerItem(2,15), 6, 5)
main_layout.addWidget(ridge_opt_box, 7, 0, 1, 5)
# set stretching (resizing) behavior
main_layout.setRowStretch(0, 1) # plot should stretch
main_layout.setRowMinimumHeight(0, 300) # plot should not get squeezed too much
main_layout.setRowStretch(5, 0) # options shouldn't stretch
main_layout.setRowStretch(6, 0) # options shouldn't stretch
main_layout.setRowStretch(7, 0) # options shouldn't stretch
self.setLayout(main_layout)
# initialize line edits
......@@ -467,6 +494,12 @@ class WaveletAnalyzer(QWidget):
pos_offset = self.w_offset,
DEBUG = self.DEBUG)
self.w_offset += 20
def ini_average_spec(self):
self.avWspecWindow = AveragedWaveletWindow(self.modulus, parent = self)
class mkWaveletCanvas(FigureCanvas):
......@@ -485,83 +518,6 @@ class mkWaveletCanvas(FigureCanvas):
class AnnealConfigWindow(QWidget):
# the signal for the anneal parameters
signal = pyqtSignal('PyQt_PyObject')
def __init__(self, parent, DEBUG):
super().__init__()
# get properly initialized in set_up_anneal
self.parentWaveletWindow = parent
self.DEBUG = DEBUG
def initUI(self, periods):
self.setWindowTitle('Ridge from Simulated Annealing')
self.setGeometry(310,330,350,200)
config_grid = QGridLayout()
ini_per = QLineEdit(str(int(np.mean(periods)))) # start at middle of period interval
ini_T = QLineEdit('10')
Nsteps = QLineEdit('5000')
max_jump = QLineEdit('3')
curve_pen = QLineEdit('0')
# for easy readout and emission
self.edits = {'ini_per' : ini_per, 'ini_T' : ini_T, 'Nsteps' : Nsteps, \
'max_jump' : max_jump, 'curve_pen' : curve_pen}
per_ini_lab = QLabel('Initial period guess')
T_ini_lab = QLabel('Initial temperature')
Nsteps_lab = QLabel('Number of iterations')
max_jump_lab = QLabel('Maximal jumping distance')
curve_pen_lab = QLabel('Curvature cost')
per_ini_lab.setWordWrap(True)
T_ini_lab.setWordWrap(True)
Nsteps_lab.setWordWrap(True)
max_jump_lab.setWordWrap(True)
curve_pen_lab.setWordWrap(True)
OkButton = QPushButton("Run!", self)
OkButton.clicked.connect(self.read_emit_parameters)
# 1st column
config_grid.addWidget( per_ini_lab, 0,0,1,1)
config_grid.addWidget( ini_per, 0,1,1,1)
config_grid.addWidget( T_ini_lab, 1,0,1,1)
config_grid.addWidget( ini_T, 1,1,1,1)
config_grid.addWidget( curve_pen_lab, 2,0,1,1)
config_grid.addWidget( curve_pen, 2,1,1,1)
# 2nd column
config_grid.addWidget( Nsteps_lab, 0,2,1,1)
config_grid.addWidget( Nsteps, 0,3,1,1)
config_grid.addWidget( max_jump_lab, 1,2,1,1)
config_grid.addWidget( max_jump, 1,3,1,1)
config_grid.addWidget( OkButton, 2,3,1,1)
# set main layout
self.setLayout(config_grid)
self.show()
def read_emit_parameters(self):
anneal_pars = {}
for par_key in self.edits:
edit = self.edits[par_key]
anneal_pars[par_key] = float(edit.text())
if self.DEBUG:
print('Anneal pars:', anneal_pars )
self.parentWaveletWindow.do_annealRidge_detection(anneal_pars)
# send to WaveletAnalyzer Window
# self.signal.emit(anneal_pars)
class WaveletReadoutWindow(QWidget):
......@@ -696,4 +652,124 @@ class mkReadoutCanvas(FigureCanvas):
QSizePolicy.Expanding,
QSizePolicy.Expanding)
FigureCanvas.updateGeometry(self)
class AveragedWaveletWindow(QWidget):
def __init__(self, modulus, parent):
super().__init__()
# --- calculate time averaged power spectrum <-> Fourier estimate ---
self.avWspec = np.sum(modulus, axis=1) / modulus.shape[1]
# -------------------------------------------------------------------
# the Wavelet analysis window spawning *this* Widget
self.parentWA = parent
self.initUI()
def initUI(self):
self.setWindowTitle(f'Fourier Spectrum Estimate - {self.parentWA.signal_id}')
self.setGeometry(510,80,550,400)
main_frame = QWidget()
pCanvas = mkGenericCanvas()
pCanvas.setParent(main_frame)
ntb = NavigationToolbar(pCanvas, main_frame)
# plot it
pCanvas.fig.clf()
pl.averaged_Wspec(self.avWspec, self.parentWA.periods, fig = pCanvas.fig)
pCanvas.fig.subplots_adjust(left = 0.15, bottom = 0.17)
main_layout = QGridLayout()
main_layout.addWidget(pCanvas,0,0,9,1)
main_layout.addWidget(ntb,10,0,1,1)
self.setLayout(main_layout)
self.show()
# --- Not used in the public version.. ---
class AnnealConfigWindow(QWidget):
'''
Not used in the public version..
'''
# the signal for the anneal parameters
signal = pyqtSignal('PyQt_PyObject')
def __init__(self, parent, DEBUG):
super().__init__()
# get properly initialized in set_up_anneal
self.parentWaveletWindow = parent
self.DEBUG = DEBUG
def initUI(self, periods):
self.setWindowTitle('Ridge from Simulated Annealing')
self.setGeometry(310,330,350,200)
config_grid = QGridLayout()
ini_per = QLineEdit(str(int(np.mean(periods)))) # start at middle of period interval
ini_T = QLineEdit('10')
Nsteps = QLineEdit('5000')
max_jump = QLineEdit('3')
curve_pen = QLineEdit('0')
# for easy readout and emission
self.edits = {'ini_per' : ini_per, 'ini_T' : ini_T, 'Nsteps' : Nsteps, \
'max_jump' : max_jump, 'curve_pen' : curve_pen}
per_ini_lab = QLabel('Initial period guess')
T_ini_lab = QLabel('Initial temperature')
Nsteps_lab = QLabel('Number of iterations')
max_jump_lab = QLabel('Maximal jumping distance')
curve_pen_lab = QLabel('Curvature cost')
per_ini_lab.setWordWrap(True)
T_ini_lab.setWordWrap(True)
Nsteps_lab.setWordWrap(True)
max_jump_lab.setWordWrap(True)
curve_pen_lab.setWordWrap(True)
OkButton = QPushButton("Run!", self)
OkButton.clicked.connect(self.read_emit_parameters)
# 1st column
config_grid.addWidget( per_ini_lab, 0,0,1,1)
config_grid.addWidget( ini_per, 0,1,1,1)
config_grid.addWidget( T_ini_lab, 1,0,1,1)
config_grid.addWidget( ini_T, 1,1,1,1)
config_grid.addWidget( curve_pen_lab, 2,0,1,1)
config_grid.addWidget( curve_pen, 2,1,1,1)
# 2nd column
config_grid.addWidget( Nsteps_lab, 0,2,1,1)
config_grid.addWidget( Nsteps, 0,3,1,1)
config_grid.addWidget( max_jump_lab, 1,2,1,1)
config_grid.addWidget( max_jump, 1,3,1,1)
config_grid.addWidget( OkButton, 2,3,1,1)
# set main layout
self.setLayout(config_grid)
self.show()
def read_emit_parameters(self):
anneal_pars = {}
for par_key in self.edits:
edit = self.edits[par_key]
anneal_pars[par_key] = float(edit.text())
if self.DEBUG:
print('Anneal pars:', anneal_pars )
self.parentWaveletWindow.do_annealRidge_detection(anneal_pars)
# send to WaveletAnalyzer Window
# self.signal.emit(anneal_pars)