refactored window function and applied window to ACF prior to computi…

…ng secondary spectrum

scintools/dynspec.py

@@ -21,7 +21,7 @@
                          fit_parabola, fit_log_parabola, fitter, \
                          powerspectrum_model
 from scintools.scint_utils import is_valid, svd_model, interp_nan_2d,\
-    centres_to_edges
+    centres_to_edges, get_window
 import scintools.ththmod as thth
 from scipy.interpolate import griddata, interp1d, RectBivariateSpline
 from scipy.signal import convolve2d, medfilt, savgol_filter
@@ -3535,27 +3535,10 @@ def calc_sspec(self, prewhite=False, halve=True, plot=False,
         nf = np.shape(dyn)[0]
         nt = np.shape(dyn)[1]
         dyn = dyn - np.mean(dyn)  # subtract mean
-
+        
         if window is not None:
-            # Window the dynamic spectrum
-            if window.lower() == 'hanning':
-                cw = np.hanning(np.floor(window_frac*nt))
-                sw = np.hanning(np.floor(window_frac*nf))
-            elif window.lower() == 'hamming':
-                cw = np.hamming(np.floor(window_frac*nt))
-                sw = np.hamming(np.floor(window_frac*nf))
-            elif window.lower() == 'blackman':
-                cw = np.blackman(np.floor(window_frac*nt))
-                sw = np.blackman(np.floor(window_frac*nf))
-            elif window.lower() == 'bartlett':
-                cw = np.bartlett(np.floor(window_frac*nt))
-                sw = np.bartlett(np.floor(window_frac*nf))
-            else:
-                print('Window unknown.. Please add it!')
-            chan_window = np.insert(cw, int(np.ceil(len(cw)/2)),
-                                    np.ones([nt-len(cw)]))
-            subint_window = np.insert(sw, int(np.ceil(len(sw)/2)),
-                                      np.ones([nf-len(sw)]))
+            chan_window, subint_window = get_window(nt, nf, window=window,
+                                                    frac=window_frac)
             dyn = np.multiply(chan_window, dyn)
             dyn = np.transpose(np.multiply(subint_window,
                                            np.transpose(dyn)))

scintools/scint_sim.py

@@ -17,6 +17,7 @@
 from scipy.interpolate import griddata
 import scipy.constants as sc
 import matplotlib.pyplot as plt
+from scintools.scint_utils import is_valid, get_window
 
 
 class Simulation():
@@ -470,10 +471,10 @@ def __init__(self, psi=0, phasegrad=0, theta=0, ar=1, alpha=5/3,
         if auto_sampling:
             # calculate to 6 spatial scales along major axis
             self.sp_fac = 6 * ar/spmax
-            # adjust to 81 pixels, doubles at ar=3
-            self.res_fac = (1 + ar/3)*81/nt
-            # triple near core
-            self.core_fac = 3
+            # adjust to 101 pixels, doubles at ar=3
+            self.res_fac = 1 + ar/3
+            # quadruple near core
+            self.core_fac = 4
         else:
             self.sp_fac = spatial_factor
             self.res_fac = resolution_factor
@@ -724,24 +725,36 @@ def plot_acf_efield(self, display=True):
         if display:
             plt.show()
 
-    def calc_sspec(self):
+    def calc_sspec(self, window='hanning', window_frac=1):
         """
         Calculate the secondary spectrum
         """
-        arr = np.fft.fftshift(self.acf)
+        nf, nt = np.shape(self.acf)
+        chan_window, subint_window = get_window(nt, nf, window=window,
+                                                frac=window_frac)
+        arr = np.multiply(chan_window, self.acf)
+        arr = np.transpose(np.multiply(subint_window,
+                                       np.transpose(arr)))
+        arr = np.fft.fftshift(arr)
         arr = np.fft.fft2(arr)
         arr = np.fft.fftshift(arr)
         arr = np.sqrt(np.real(arr * np.conj(arr)))
         self.sspec = 10*np.log10(arr)
 
-    def plot_sspec(self, display=True):
+    def plot_sspec(self, display=True, vmin=None, vmax=None):
         """
         Plots the simulated ACF
         """
         if not hasattr(self, 'sspec'):
             self.calc_sspec()
-
-        plt.pcolormesh(self.tn, self.fn, self.sspec)
+
+        sspec = self.sspec
+        medval = np.median(sspec[is_valid(sspec)*np.array(np.abs(sspec) > 0)])
+        maxval = np.max(sspec[is_valid(sspec)*np.array(np.abs(sspec) > 0)])
+        vmin = medval - 3 if vmin is None else vmin
+        vmax = maxval - 3 if vmax is None else vmax
+
+        plt.pcolormesh(self.tn, self.fn, sspec, vmin=vmin, vmax=vmax)
         plt.colorbar()
         plt.xlabel(r'Delay')
         plt.ylabel(r'Doppler')

scintools/scint_utils.py

@@ -803,6 +803,31 @@ def make_pickle(obj, filepath):
     with open(filepath, 'wb') as f_out:
         for idx in range(0, n_bytes, max_bytes):
             f_out.write(bytes_out[idx:idx+max_bytes])
+    return     
+
+def get_window(nt, nf, window='hanning', frac=0.1):
+    """
+    Returns a window to apply before computing an FFT
+    """
+    if window.lower() == 'hanning':
+        cw = np.hanning(np.floor(frac*nt))
+        sw = np.hanning(np.floor(frac*nf))
+    elif window.lower() == 'hamming':
+        cw = np.hamming(np.floor(frac*nt))
+        sw = np.hamming(np.floor(frac*nf))
+    elif window.lower() == 'blackman':
+        cw = np.blackman(np.floor(frac*nt))
+        sw = np.blackman(np.floor(frac*nf))
+    elif window.lower() == 'bartlett':
+        cw = np.bartlett(np.floor(frac*nt))
+        sw = np.bartlett(np.floor(frac*nf))
+    else:
+        print('Window unknown.. Please add it!')
+    chan_window = np.insert(cw, int(np.ceil(len(cw)/2)),
+                            np.ones([nt-len(cw)]))
+    subint_window = np.insert(sw, int(np.ceil(len(sw)/2)),
+                              np.ones([nf-len(sw)]))
+    return chan_window, subint_window
 
 
 def calculate_curvature_peak_probability(power_data, noise_level, smooth=True,