Added Fourier distribution ouputs

CHANGELOG.md

+### pyBOAT 0.8.18
+
+- Added time averaging of Wavelet spectra <-> Fourier estimate
+- Reworked FFT visualizations
+
 ### pyBOAT 0.8.17
 
 - Fixed crashes during batch analysis with thresholded ridges

pyboat/__init__.py

@@ -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

pyboat/ensemble_measures.py

@@ -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)

pyboat/plotting.py

@@ -183,6 +183,22 @@ def Fourier_spec(ax, fft_freqs, fft_power, show_periods=False):
 
 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))
 
@@ -426,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
@@ -451,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)
@@ -499,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,
@@ -585,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()

pyboat/ssg.py

@@ -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:
 

pyboat/ui/batch_process.py

@@ -88,9 +88,13 @@ class BatchProcessWindow(QWidget):
         self.cb_power_dis.setToolTip('Show time-averaged wavelet power of the ensemble')
         self.cb_plot_ens_dynamics = QCheckBox('Ensemble Dynamics')
         self.cb_plot_ens_dynamics.setToolTip('Show period, amplitude and phase distribution over time')
+        self.cb_plot_Fourier_dis = QCheckBox('Fourier Spectra Distribution')
+        self.cb_plot_Fourier_dis.setToolTip('Distribution of the time averaged Wavelet spectra')
+        
         lo = QGridLayout()
         lo.addWidget(self.cb_power_dis,0,0)
         lo.addWidget(self.cb_plot_ens_dynamics,1,0)
+        lo.addWidget(self.cb_plot_Fourier_dis,2,0)        
         plotting_options.setLayout(lo)
 
         # -- Save Out Results --
@@ -111,6 +115,9 @@ class BatchProcessWindow(QWidget):
         self.cb_sorted_powers.setToolTip("Saves the time-averaged powers in descending order")
         self.cb_save_ensemble_dynamics = QCheckBox('Ensemble Dynamics')
         self.cb_save_ensemble_dynamics.setToolTip("Saves each period, amplitude and phase summary statistics to a csv table")
+
+        self.cb_save_Fourier_dis = QCheckBox('Fourier Distribution')
+        self.cb_save_Fourier_dis.setToolTip("Saves median and quartiles of the Fourier power spectra")
         
         home = expanduser("~")
         OutPath_label = QLabel('Export to:')
@@ -138,6 +145,7 @@ class BatchProcessWindow(QWidget):
         #lo.addWidget(line1, 3,0)        
         lo.addWidget(self.cb_sorted_powers,4,0)
         lo.addWidget(self.cb_save_ensemble_dynamics,5,0)
+        lo.addWidget(self.cb_save_Fourier_dis,6,0)        
         #lo.addWidget(line2, 6,0)        
         lo.addWidget(PathButton,7,0)                
         lo.addWidget(self.OutPath_edit,8,0)
@@ -198,9 +206,8 @@ class BatchProcessWindow(QWidget):
             if OutPath is None:
                 return
 
-        # is a dictionary holding the ridge-data
-        # for each signal and the signal_id as key
-        ridge_results = self.do_the_loop()
+        # TODO: parallelize
+        ridge_results, df_fouriers = self.do_the_loop()
 
         # check for empty ridge_results
         if not ridge_results:
@@ -209,17 +216,14 @@ class BatchProcessWindow(QWidget):
 
         
         
-        # compute the time-averaged powers
+        # --- compute the time-averaged powers ---
+        
         if self.cb_power_dis.isChecked() or self.cb_sorted_powers.isChecked():
             
             powers_series = em.average_power_distribution(ridge_results.values(),
                                                           ridge_results.keys(),
                                                           exclude_coi = True)
 
-            # sort by power, descending
-            powers_series.sort_values(
-                ascending = False,
-                inplace = True)
 
         if self.cb_power_dis.isChecked():            
             # plot the distribution
@@ -228,9 +232,11 @@ class BatchProcessWindow(QWidget):
 
         # save out the sorted average powers
         if self.export_options.isChecked() and self.cb_sorted_powers.isChecked():
-            fname = f'{OutPath}/average_powers_{dataset_name}.csv'
+            fname = f'{OutPath}/average-ridge-powers_{dataset_name}.csv'
             powers_series.to_csv(fname, sep = ',', index = True, header = False)
-        # compute summary statistics over time
+            
+        # --- compute summary statistics over time ---
+        
         if self.cb_plot_ens_dynamics.isChecked() or self.cb_save_ensemble_dynamics.isChecked():
             # res is a tuple of  DataFrames, one each for
             # periods, amplitude, power and phase
@@ -254,6 +260,27 @@ class BatchProcessWindow(QWidget):
                 df.index.name = 'time'
                 df.to_csv(fname, sep = ',', float_format = '%.3f')
 
+        # --- Fourier Distribution Outputs ---
+        
+        if self.cb_plot_Fourier_dis.isChecked():
+
+            self.fdw = FourierDistributionWindow(df_fouriers,
+                                                 self.parentDV.time_unit,
+                                                 dataset_name)
+            
+
+        if self.export_options.isChecked() and self.cb_save_Fourier_dis.isChecked():
+
+            fname = f'{OutPath}/fourier-distribution_{dataset_name}.csv'
+
+            # save out median and quartiles of Fourier powers
+            df_fdis = pd.DataFrame(index = df_fouriers.index)
+            df_fdis['Median'] = df_fouriers.median(axis = 1)
+            df_fdis['Q1'] = df_fouriers.quantile(q = 0.25, axis = 1)
+            df_fdis['Q3'] = df_fouriers.quantile(q = 0.75, axis = 1)
+            
+            df_fdis.to_csv(fname, sep = ',', float_format = '%.3f')
+            
         if self.debug:
             print(list(ridge_results.items())[:2])
 
@@ -378,8 +405,12 @@ class BatchProcessWindow(QWidget):
         # retrieve batch settings
         power_thresh = self.get_thresh()
         rsmooth = self.get_ridge_smooth()
-        
+
+        # results get stored here
         ridge_results = {}
+        df_fouriers = pd.DataFrame(index = periods)
+        df_fouriers.index.name = 'period'
+        
         for i, signal_id in enumerate(self.parentDV.df):
 
             # log to terminal
@@ -423,11 +454,16 @@ class BatchProcessWindow(QWidget):
                 tvec,
                 power_thresh,
                 smoothing_wsize = rsmooth)
-            
+
+            # from ridge thresholding..
             if ridge_data.empty:
                 EmptyRidge += 1
             else:
                 ridge_results[signal_id] = ridge_data                            
+
+            # time average the spectrum, all have shape len(periods)!
+            averaged_Wspec = np.mean(modulus, axis = 1)
+            df_fouriers[signal_id] = averaged_Wspec
             
             # -- Save out individual results --
             
@@ -470,7 +506,7 @@ class BatchProcessWindow(QWidget):
                 f'{EmptyRidge} ridge readouts entirely below threshold..',
                 'Discarded ridges')
                 
-        return ridge_results
+        return ridge_results, df_fouriers
 
 
 class PowerDistributionWindow(QWidget):
@@ -486,7 +522,7 @@ class PowerDistributionWindow(QWidget):
     def initUI(self, dataset_name):
 
         self.setWindowTitle(f'Ridge Power Distribution - {dataset_name}')
-        self.setGeometry(510,80,550,400)
+        self.setGeometry(410,220,550,400)
 
         main_frame = QWidget()
         pCanvas = mkGenericCanvas()        
@@ -495,7 +531,7 @@ class PowerDistributionWindow(QWidget):
 
         # plot it
         pCanvas.fig.clf()
-        pl.plot_power_distribution(self.powers, fig = pCanvas.fig)
+        pl.power_distribution(self.powers, fig = pCanvas.fig)
         pCanvas.fig.subplots_adjust(left = 0.15, bottom = 0.17)
         
         main_layout = QGridLayout()
@@ -519,7 +555,7 @@ class EnsembleDynamicsWindow(QWidget):
     def initUI(self, dataset_name):
 
         self.setWindowTitle(f'Ensemble Dynamics - {dataset_name}')
-        self.setGeometry(510,80,700,480)
+        self.setGeometry(210,80,700,480)
 
         main_frame = QWidget()
         Canvas = mkGenericCanvas()        
@@ -527,7 +563,7 @@ class EnsembleDynamicsWindow(QWidget):
         ntb = NavigationToolbar(Canvas, main_frame)        
 
         Canvas.fig.clf()
-        pl.plot_ensemble_dynamics(*self.results,
+        pl.ensemble_dynamics(*self.results,
                                   dt = self.dt,
                                   time_unit = self.time_unit,
                                   fig = Canvas.fig)
@@ -540,6 +576,39 @@ class EnsembleDynamicsWindow(QWidget):
         self.setLayout(main_layout)
         self.show()
 
+class FourierDistributionWindow(QWidget):
+    def __init__(self, df_fouriers, time_unit, dataset_name = ''):
+        super().__init__()
+
+        self.time_unit = time_unit
+        self.df_fouriers = df_fouriers
+        
+        self.initUI(dataset_name)
+
+    def initUI(self, dataset_name):
+
+        self.setWindowTitle(f'Fourier Power Distribution - {dataset_name}')
+        self.setGeometry(510,330,550,400)
+
+        main_frame = QWidget()
+        Canvas = mkGenericCanvas()        
+        Canvas.setParent(main_frame)
+        ntb = NavigationToolbar(Canvas, main_frame)        
+
+        Canvas.fig.clf()
+        pl.Fourier_distribution(self.df_fouriers,
+                                  time_unit = self.time_unit,
+                                  fig = Canvas.fig)
+        
+        Canvas.fig.subplots_adjust(wspace = 0.3, left = 0.1, top = 0.98,
+                        right = 0.95, bottom = 0.15)
+        main_layout = QGridLayout()
+        main_layout.addWidget(Canvas,0,0,9,1)
+        main_layout.addWidget(ntb,10,0,1,1)
+        
+        self.setLayout(main_layout)
+        self.show()
+        
 
         
 

pyboat/ui/start_menu.py

@@ -50,7 +50,6 @@ class MainWindow(QMainWindow):
 
         self.nViewers = 0
         self.DataViewers = {}  # no Viewers yet
-        self.detr = {}
         self.initUI()
 
     def initUI(self):
@@ -74,7 +73,7 @@ class MainWindow(QMainWindow):
         go_to_doc = QAction("&Documentation..", self)
         go_to_doc.triggered.connect(self.open_doc_link)
 
-        # online help in low left corner
+        # online help in lower left corner
         self.statusBar()
 
         # the menu bar
@@ -111,6 +110,7 @@ class MainWindow(QMainWindow):
 
         synsigButton = QPushButton("Start Generator", self)
         synsigButton.setStyleSheet("background-color: orange")
+        synsigButton.setStatusTip("Start the synthetic signal generator")        
         synsigButton.clicked.connect(self.init_synsig_generator)
 
         # quitButton = QPushButton("Quit", self)