<reponame>fjeng/aeptools
def estimatef0(hd='C:', logfilepath='C:/data/GroupData_NMF_feasibility_python/', logfilename='NMF_feasibility_log.json', file='', win=50.0, gap=0.5, freqResolution=1.0, upperFreq=1000.0, low=92.4, high=149.8, se_absFlag=False, r_method=1, f0ContourPath=''):
    '''
    Purpose : estimatef0 (spectral maximum method) for FFR recordings
    
    parameters
    ----------
    hd : string
        - hard drive
    logiflepath : string
        - file path for log book
    logfilename : string
        - file name for log book
    file : a core.frame.Pandas object 
        - a file to be processed
    win : float64
        - window size in ms for calculating spectrogram 
    gap : float64
        - step size in ms for calculating spectrogram
    freqResolution : float64
        - frequency resolution in Hz for calculating spectrogram 
    upperFreq : float64
        - the highest frequency to take out from the spectrogram output [default 1000]
    low : float
        - the lower boundary when searching candidates for f0 location
    high : float
        - the higher boundary when searching candidates for f0 location
    se_absFlag : boolean
        - flag indicating whether to take an absolute value of se, or not (True yes, False no) 
    r_method : integer
        - method to calculate Tracking Accuracy (0 ==> linear regression, r between [0 1]; 1 ==> correlcation coefficient, r between [-1 +1])
    f0ContourPath : string
        - file path where f0 contours are saved
    
    Returns
    -------
    None

    Diary
    -----
    2021-09-09 (v01) Fuh
        - rewrote _03estimatef0 script in Python
    2021-10-15 (v02)
        - rewrote this script as a function
        
    ToDo
    ----
    None
        
    '''
    
    # clear console, clear variables
    from IPython import get_ipython         # The 'IPython' is a Python built-in module, for interactive computing in different languages
    get_ipython().magic('clear')            # clear: clear the console 
    get_ipython().magic('reset -f')         # -f : force reset without asking for confirmation (i.e., a 'hard' reset of the kernel)
    
    # close all figures
    import matplotlib.pyplot as plt         
    plt.close('all')
    
    # import modules
    from pathlib import Path                    # pathlib is a built-in library in since Python 3.4. The 'Path' class is a new, correct, and easy way to work with file paths and file names.
    from scipy import signal                    # for resample_poly(), firwin(), lfilter(), etc.
    from scipy.io import wavfile                # for wavfile.read()
    import sys                                  # import the Python built-in module 'sys' for sys.exit(), sys.platform
    # import pandas as pd                         # import Pandas for DataFrame
    import numpy as np
    import json
    
    # # initialize parameters
    # hd = 'E:'                               # hard drive (Windows)
    # if sys.platform == 'darwin':            # macOS
    #     # hd = '/Volumes/AEPLAB'            # Fuh Mac computer, VPN S drive
    #     # hd = '/Users/groverlab'           # Fuh Mac computer, local HD
    #     hd = '/Users/groverlab/OneDrive - Ohio University'   # Fuh Mac computer, OneDrive - Ohio University
    #     # hd = '/Volumes/FUHBACKUP'         # Fuh Mac computer, USB thumb drive
    
    # # --- read in par ---
    # logfilepath = hd + '/data/GroupData_NMF_feasibility_python/'   
    # logfilename = 'NMF_feasibility_log.json'
    # logfilepathname = logfilepath + logfilename             
    # par = pd.read_json(logfilepathname, orient='table')
    
    # hide the toolbar when showing new figures
    plt.rcParams['toolbar'] = 'None'       # This will hide toolbar for all figures, unless it is turned back on. ['toolbar2' (default), 'None', 'toolmanager']
    
    # define a dictionary for converting numbers to color maps (for ploting spectrograms)
    dict_colorMap = {
        0: 'gray_r',
        1: 'jet',
        }
    
    # define a dictionary for converting stimToken to a tuple (stimpath, stimname)
    dict_stimToken = {
        'i2_150ms': ('/Stim2/RecordedToneStimuli/i/origcycN_150ms/', 'i2_150ms.wav'),   
        'i2_250ms': ('/Stim2/RecordedToneStimuli/i/origcycN/', 'i2_250ms.wav'),   
        'i2_CycAmpN': ('/Stim2/RecordedToneStimuli/i/cycAmpN/', 'i2_CycAmpN.wav'),   
        }
    
    # stores a list of stimTokens that have been processed in this subroutine
    # TODO: return the accumulated stimTokens list back to the calling script, so that the same stimToken will not be processed twice.
    stimTokens = []
    
    # data flags
    outputStimf0Flag = True
    outputRespf0Flag = True
    outputf0RegressionFlag = True
    
    # figure flags
    outputFigureStimSpectrogramFlag = False
    outputFigureRespSpectrogramFlag = False
    outputFigureStimf0Flag = False
    outputFigureRespf0Flag = False
    outputFigureRegressionFlag = False
    outputFiguref0ComparisonFlag = True
    setFigureLocationFlag = True        
    space = 10                                  # space between adjacent figures in pixels 
    dx = 640                                    # width of a figure in pixels
    dy = 480                                    # height of a figure in pixels
    colorMap_stimulus = 1                       # define color map for stimulus spectrogram
    colorMap_recording = 1                      # define color map for recording spectrogram
    
    # # spectrogram parameters
    # win = 50.0             # window size in ms for calculating spectrogram 
    # gap = 0.5              # step size in ms for calculating spectrogram [changed from 1 to 0.5 ms for 170 ms di1234 tokens]
    # centerFlag = True      # time array is based on the center point of each time bin, or not [default = 1]
    # freqResolution = 1.0   # have a frequency resolution of 0.1 Hz here to increase accuracy in estimating f0 --- it will take a long time to run through, but the data will be worth it.
    #                        # 20000/(2^12) # frequency resolution in Hz for narrow-band spectrogram 
    #                        # old default freqResolution = fs / (2^12) = 20000 / (2^12) = 4.88 Hz
    #                        # new default freqResolution = fs / fs = 1 Hz
    # upperFreq = 1000.0     # frequency range to take out from the spectrogram output [default 1000, i.e., take out 1-1000 Hz data out of a total of a 10,000 Hz data range
    
    backend = plt.get_backend()  # get graphics backend
    
    polarityNames = ['']                  # ['', '_a', '_b']   # a list of polarity names to be included for analysis
    npolarities = len(polarityNames)      # number of polarity conditions
    
    # low = 92.4             # 102.4 from "f0_contours_i1234_150ms_origCycN_spec_150ms.txt"
    # high = 149.8           # 139.8
    
    # se_absFlag = False     # flag to take an absolute value of se (1 yes, 0 no)
    
    # r_method = 1           # method to calculate Tracking Accuracy (0 ==> linear regression, r between [0 1] 1 ==> correlcation coefficient, r between [-1 +1])
    
    # f0ContourPath = hd + '/Stim2/RecordedToneStimuli/i/f0contours_150ms/'
            
    # # clear all local variables, except a few, before a new iteration starts. 
    # keepList = ['get_ipython', 'plt', 'Path', 'signal', 'wavfile', 'sys', 'pd', 'np', 'json', 
    #             'hd', 'par', 'file', 'space', 'dx', 'dy', 'win', 'gap', 'freqResolution', 'upperFreq', 'backend', 'polarityNames', 'npolarities',
    #             'low', 'high', 'r_method', 'f0ContourPath']  
    # keepPrefix = ('keep', 'var', 'log', 'dict', 'colorMap', 'stim')     # keep keep*, var*, log*, dict*, etc
    # keepSuffix = ('Flag')  # keep *Flag, etc
    # varListLocal = list(locals().keys())    # Get local variables ==> dictionary keys ==> convert dict into a list
    # for var in varListLocal:
    #     if (var not in keepList) and (not var.startswith(keepPrefix)) and (not var.endswith(keepSuffix)) and (var[0] != '_'):   # if (1) var is not in the keep list, (2) var does not start with a certain prefix, (3) var does not end with a certain suffix, and (4) var is not a system 'hidden' variable with a prefix '_', ==> delete it.
    #         exec(f'del {var}')        

    # define filepath, filename, and related variables
    filepath = hd + '/data/' + file.directory + '/'
    filepre = file.subjectCode + 'ses' + str(file.sessionNumber) + '_'  
    filenamestr = filepre + file.fileNumberString
    filepathstr = filepath + filepre + file.fileNumberString 
    
    # read in header
    if file.labviewIndex == 1:  # whether the recording is made from Labview or not (1 yes, 0 no). When no, it is recorded by using Neuroscan.
        header = np.load(filepathstr + '_header.npy', allow_pickle=True).item()
        # jsonFile = open(filepathstr + '_header.json', 'r')
        # header = json.load(jsonFile)
        # jsonFile.close()
        
        softwareversion = header['softwareversion']
        if softwareversion >= 2.17:        
            nchans = int(header['AI_nchannels'])
            if nchans <= 0:                 # nchans must be >= 1
                nchans = 1                  
            fs = header['AI_samplingRate']
            si = header['AO_silentInterval'] / 1000.0
            isi = header['AI_nsamplesReadPerChannel'] / fs 
            stimPolarity = int(header['stim_polarity'])
            if stimPolarity == 2:        # alternating polarity
                stimdur = (isi / 2) - si
            else:                        # original or inverted polarity
                stimdur = isi - si
            dBSPL = header['AO0_dBSPL'] 
            tailFlag = int(header['tail_flag'])
            stimPathName = header['AO0_stimPath']
            stimdur_ms = stimdur * 1000.0              # stim dur in ms
    
        else:
            sys.exit('not recognizable AEP software version...')            
    else:    
        fs = 20000.0
        stimdur = 0.150            
        stimdur_ms = stimdur * 1000.0
        si = file.si / 1000.0      
        isi = stimdur + si
        stimPolarity = 2
        
    # -------------- initialize parameters -------------------------------
    # nrep = file.nRepetitions  # number of repetitions
    # overlap = win - gap     
    # window = np.hanning(np.floor(win * fs / 1000.0))      # window is a hanning window of 50 ms
    # noverlap = int(np.floor(overlap * fs / 1000.0))       # overlaps noverlap samples of each segment, 
    # nfft = int(np.floor(win * fs / 1000.0))               # default 1000
    # pad_to = int(np.floor((1.0 / freqResolution) * fs))   # convert freqResolution to number of data points [default nfft = 2^12]. 
    # #         # zero padding each segment to have nfft number of data points 
    # #         # when nfft = fs, frequency resolution = 1 Hz

    # take only 1 out of every 'interval' data points along the frequency axis, so that frequency resolution is 10 Hz on the spectrogram plot (for Python and Veusz).
    freq_stepsize = 10.0                                      # step size on the frequency axis [default freq step size = 10 Hz]
    interval = int(np.floor(freq_stepsize / freqResolution))  # make frequency interval to be 10 Hz, for Phtyon and Veusz plotting purposes

    # ---- read stimulus token, first occurrence only ------------
    stimToken = file.stimToken              # stimulus token entered on MATLAB logbook, column 8
    if dict_stimToken.get(stimToken) is not None:
        stimpath, stimname = dict_stimToken[stimToken]
    else:
        sys.exit('not recognizable stimulus token...')
    
    stimfs, stim = wavfile.read(hd + stimpath + stimname)

    # make stimulus token to have the same sampling rate as recording fs [default = 20000]
    up = 100                                            # set 'up' sampling at 100 (artibutrarily)
    down = int(up * (stimfs / fs))                 # compute 'down' sampling, to get desired new fs (relative to the old fs)
    stim = signal.resample_poly(stim, up, down)         # resample y to have a total of npts data points
    if len(stim) > int(np.floor(fs * stimdur)):
        stim = stim[0 : int(np.floor(fs * stimdur))]     
        
    # normalize stim amplitudes (time waveform) 
    stim = stim / np.amax(stim)

    # -------------Read in stimulus and recording spectrograms ----------
    for  chan in range(nchans):          
        for polarity in range(npolarities):    
            
            # close all figures
            plt.close('all')   
            
            # append chanstr and polaritystr to filepathstr and filenamestr
            chanstr = '' if nchans==1 else '_chan' + str(chan)
            polaritystr = polarityNames[polarity]                  
            filepathstr = filepath + filepre + file.fileNumberString + chanstr + polaritystr  
            filenamestr = filepre + file.fileNumberString + chanstr + polaritystr
            filename = filenamestr + '.npy'

            annotationStimulus = stimToken
            annotationRecording = filename + ', ' + file.experimentString
                       
            print('Reading file: ' + filepath + filename, flush=True)
    
            # --- read in stimulus spectrogram data ---
            T_ms = np.load(logfilepath + file.stimToken + '_stim_spectrogram_time_ms.npy')    
            F = np.load(logfilepath + file.stimToken + '_stim_spectrogram_freq.npy') 
            nV_S = np.load(logfilepath + file.stimToken + '_stim_spectrogram_nV.npy')     
            nV_S_veusz = np.load(logfilepath + file.stimToken + '_stim_spectrogram_nV_veusz.npy')  
    
            # compute (left, right, bottom, top) boundaries for imshow(), for stimulus and response
            left = win/2.0 - gap/2.0              
            right = (stimdur_ms-win/2.0) + gap/2.0
            bottom = -(freqResolution/2.0) * interval
            top = F[-1] + (freqResolution/2.0)*interval
            
            # define x tick labels
            if int(stimdur*1000) == 150:
                xticks = (25,50,75,100,125)
            elif int(stimdur*1000) == 250:
                xticks = (25,75,125,175,225)
            else:
                sys.exit('Unrecognizable stimulus duration')
                      
            # --- plot stimulus spectrogram ---
            fig101 = plt.figure(101)
            plt.imshow(nV_S_veusz, origin='lower', extent=(left, right, bottom, top), aspect='auto', cmap=dict_colorMap.get(colorMap_stimulus)) # , vmin=0, vmax=8000)   # set zero origin at lower left corner, define the extent of axis boundaries (so that the x and y tick numbers will be correct), set aspect ratio to 'auto' (so that the spectrogram image will fill in the entire space), set color map to 'jet' for color and 'gray_r' for gray scale, set vmin and vmax for the lower and upper limits along the z axis (i.e., amplitude axis)
            plt.title(annotationStimulus)     
            plt.xlabel('Time (ms)')
            plt.ylabel('Frequency (Hz)')
            plt.xticks(ticks=xticks)
            # plt.colorbar()               
            if setFigureLocationFlag:
                if 'Qt' in backend:          # Qt (e.g., Qt5Agg)
                    mgr = plt.get_current_fig_manager()              # get a manager for current figure
                    mgr.window.setGeometry(space, space, dx, dy)     # set geometry of current figure (x, y, width, height) in pixels
                else:
                    sys.exit('Unrecognizable graphical backend')
            plt.show(block=False)    
            
            # --- read in response spectrogram data ---
            t_ms = np.load(filepathstr + '_spectrogram_time_ms.npy')    
            f = np.load(filepathstr + '_spectrogram_freq.npy')    
            nV = np.load(filepathstr + '_spectrogram_nV.npy')   
            nV_veusz = np.load(filepathstr + '_spectrogram_nV_veusz.npy')  
    
            # --- plot response spectrogram ---
            fig102 = plt.figure(102)
            plt.imshow(nV_veusz, origin='lower', extent=(left, right, bottom, top), aspect='auto', cmap=dict_colorMap.get(colorMap_recording)) # , vmin=0, vmax=8000)   # set zero origin at lower left corner, define the extent of axis boundaries (so that the x and y tick numbers will be correct), set aspect ratio to 'auto' (so that the spectrogram image will fill in the entire space), set color map to 'jet' for color and 'gray_r' for gray scale, set vmin and vmax for the lower and upper limits along the z axis (i.e., amplitude axis)
            plt.title(annotationRecording)
            plt.xlabel('Time (ms)')
            plt.ylabel('Frequency (Hz)')
            plt.xticks(ticks=xticks)
            # plt.colorbar()               
            if setFigureLocationFlag:
                if 'Qt' in backend:          # Qt (e.g., Qt5Agg)
                    mgr = plt.get_current_fig_manager()                     # get a manager for current figure
                    mgr.window.setGeometry(space+dx+space, space, dx, dy)   # set geometry of current figure (x, y, width, height) in pixels
                else:
                    sys.exit('Unrecognizable graphical backend')
            plt.show(block=False)
            
            ###############################
            # --- Find spectral maximum ---
            ###############################
                
            #----------------- extract f0s for stimulus ----------------
            if stimToken not in stimTokens:           # compute stimulus f0s, only upon the first encounter of a stimToken
                stimTokens.append(stimToken)
                highIndices = np.argwhere(F <= high)  # find the upper frequency index of the searching range
                highIndex = highIndices[-1, 0] + 1
                lowIndices = np.argwhere(F >= low)    # find the lower frequency index of the searching range
                lowIndex = lowIndices[0, 0]
                
                F_range = F[lowIndex:highIndex]
                nV_S_range = nV_S[lowIndex:highIndex, :]
                
                # in the range of T_ms, find the maximum y and its index
                f0amps = np.full(T_ms.shape, None)
                indices = np.full(T_ms.shape, None)
                for i in range(len(T_ms)):
                    f0amps[i] = np.amax(nV_S_range[:, i])
                    indices[i] = np.argmax(nV_S_range[:, i])
                    
                # use the index to find the corresponding frequency
                f0s = np.full(T_ms.shape, None)
                for j in range(len(T_ms)):
                    f0s[j] = F_range[indices[j]]
                f0s = f0s.astype(float)            # force f0s to float64
                    
                # output f0s
                np.save(f0ContourPath + 'f0_' + stimToken + '.npy', f0s)  
            else:
                f0s = np.load(f0ContourPath + 'f0_' + stimToken + '.npy', allow_pickle=True)  
                                    
            # --- plot stimulus f0s ---
            fig103 = plt.figure(103)
            plt.plot(T_ms, f0s, 'k')         # {b: blue, k: black, r: red, etc}, {o: circle, +: plus, x: cross, etc}
            plt.title('F0s of the stimulus, ' + stimToken)
            plt.xlabel('Time (ms)')
            plt.ylabel('Frequency (Hz)')
            if setFigureLocationFlag:
                if 'Qt' in backend:          # Qt (e.g., Qt5Agg)
                    mgr = plt.get_current_fig_manager()              # get a manager for current figure
                    mgr.window.setGeometry(space, int(space+dy*6/10+space), dx, dy)     # set geometry of current figure (x, y, width, height) in pixels
                elif 'wx' in backend:        
                    print('wx backend')   
                else:
                    sys.exit('Unrecognizable graphical backend')
            plt.show(block=False)
    
            #----------------- extract f0s for a recording ----------------
            highIndices = np.argwhere(f <= high)  # find the upper frequency index of the searching range
            highIndex = highIndices[-1, 0] + 1
            lowIndices = np.argwhere(f >= low)    # find the lower frequency index of the searching range
            lowIndex = lowIndices[0, 0]
            
            f_range = f[lowIndex:highIndex]
            nV_range = nV[lowIndex:highIndex, :]
            
            # in the range of T_ms, find the maximum y and its index
            f0amp = np.full(t_ms.shape, None)
            index = np.full(t_ms.shape, None)
            for i in range(len(t_ms)):
                f0amp[i] = np.amax(nV_range[:, i])
                index[i] = np.argmax(nV_range[:, i])
                
            # use the index to find the corresponding frequency
            f0 = np.full(t_ms.shape, None)
            for j in range(len(t_ms)):
                f0[j] = f_range[index[j]]
            f0 = f0.astype(float)             # force f0 to float64
            
            # calculate Spectral Amplitude (SA)
            sa = np.mean(abs(f0amp))        
    
            # --- plot response f0s ---
            fig104 = plt.figure(104)
            plt.plot(t_ms, f0, 'r')
            plt.title('F0s of a recording, ' + annotationRecording)
            plt.xlabel('Time (ms)')
            plt.ylabel('Frequency (Hz)')
            if setFigureLocationFlag:
                if 'Qt' in backend:          # Qt (e.g., Qt5Agg)
                    mgr = plt.get_current_fig_manager()              # get a manager for current figure
                    mgr.window.setGeometry(space+dx+space, int(space+dy*6/10+space), dx, dy)     # set geometry of current figure (x, y, width, height) in pixels
                elif 'wx' in backend:        
                    print('wx backend')   
                else:
                    sys.exit('Unrecognizable graphical backend')
            plt.show(block=False)
                
            # ---- comparison of f0 contours ----
            fig105 = plt.figure(105)
            plt.plot(T_ms, f0s, 'k', label='Stimulus')
            plt.plot(t_ms, f0, 'r', label='Response')
            plt.title('F0 contours, ' + annotationRecording)
            plt.xlabel('Time (ms)')
            plt.ylabel('Frequency (Hz)')
            plt.legend()
            if setFigureLocationFlag:
                if 'Qt' in backend:          # Qt (e.g., Qt5Agg)
                    mgr = plt.get_current_fig_manager()              # get a manager for current figure
                    mgr.window.setGeometry(space*2, space*8, dx, dy)     # set geometry of current figure (x, y, width, height) in pixels
                elif 'wx' in backend:        
                    print('wx backend')   
                else:
                    sys.exit('Unrecognizable graphical backend')
            plt.show(block=False)
    
            # --- compute Frequency Error (FE) ---
            X = np.vstack([f0s, np.ones(len(f0s))]).transpose()   # append ones column to x               
            (b, _, _, _) = np.linalg.lstsq(X, f0, rcond=None)   # perform linear regression       
            yestimate = b[0] * f0s + b[1]    # compute y estimate      
            fe = np.mean(abs(yestimate - f0))     # compute frequency error (fe)
            
            # --- compute Tracking Accuracy (TA) ---
            r = np.corrcoef(f0s, f0)[0, 1]      # correlation coefficient (TA) between [-1.0 1.0]
            if r_method == 0:
                r = abs(r)
            R2 = r ** 2                      # R2 = r^2 = coefficient of determination
    
            fig106 = plt.figure(106)
            plt.plot(f0s, f0, 'r', label='F0s')
            plt.plot(f0s, yestimate, 'b', label='Fitted line')
            plt.title('Linear regression, ' + annotationRecording)
            plt.xlabel('F0 of stimulus (Hz)')
            plt.ylabel('F0 of response (Hz)')
            plt.legend()
            if setFigureLocationFlag:
                if 'Qt' in backend:          # Qt (e.g., Qt5Agg)
                    mgr = plt.get_current_fig_manager()              # get a manager for current figure
                    mgr.window.setGeometry(space*4 + dx, space*8, dx, dy)     # set geometry of current figure (x, y, width, height) in pixels
                elif 'wx' in backend:        
                    print('wx backend')   
                else:
                    sys.exit('Unrecognizable graphical backend')
            plt.show(block=False)
            
            # --- compute Slope Error (SE) ---
            X = np.vstack([t_ms, np.ones(len(t_ms))]).transpose()     # append ones column to x               
            (b_s, _, _, _) = np.linalg.lstsq(X, f0s, rcond=None)      # stimulus f0 regression
            (b_ffr, _, _, _) = np.linalg.lstsq(X, f0, rcond=None)     # FFR f0 regression
            se = (b_ffr[0] - b_s[0]) * 1000.0                         # slope error in Hz/sec
            if se_absFlag:
                se = abs(se)
    
            # --- compute RMS amplitude ---
            resp = np.load(filepathstr + '_averaged_' + str(int(stimdur_ms)) + 'ms_recording.npy', allow_pickle=True)    # read time waveform
            resp = resp.astype(float)           # force to float64 precision
            rms = np.sqrt(np.mean(resp**2))     # compute RMS amplitude
    
            # display results on console 
            print('Frequency Error (Hz) = %3.8f' % fe)
            print('Slope Error (Hz/s) = %3.8f' % se)
            print('Tracking Accuracy = %3.8f' % r)
            print('Spectral Amplitude (nV) = %3.8f' % sa)
            print('RMS Amplitude (nV) = %3.8f' % rms)
            
            # --- store results in data (i.e., a dictionary)---
            data = {
                'filepathstr': filepathstr,
                'stimToken': stimToken,
                'fe': fe,
                'se': se,
                'ta': r,
                'sa': sa,
                'rms': rms,
                'b_s_slope': b_s[0] * 1000.0,      # slope of stimulus f0s in Hz/second
                'b_s_intercept': b_s[1],    
                'b_ffr_slope': b_ffr[0] * 1000.0,  # slope of FFR f0s in Hz/second
                'b_ffr_intercept': b_ffr[1],  
                'b0': b[0],                   
                'b1': b[1],    
                'R2': R2,                          # coefficient of determination   
                }
        
            # --- Output data ---
            if outputStimf0Flag and file[0] == 0:  # stimulus f0s
                np.save(logfilepath + stimToken + '_stimf0.npy', f0s)
            
            if outputRespf0Flag:                   # response f0s and f0 amplitudes
                np.save(filepathstr + '_respf0.npy', f0)
                np.save(filepathstr + '_f0amp.npy', f0amp)
            
            if outputf0RegressionFlag:             # regression data
                # save regression data (dictionary) as .json
                jsonString = json.dumps(data)                                  # convert dict to a json string
                jsonFile = open(filepathstr + '_f0_LinearRegress.json', 'w')   # open a json file
                jsonFile.write(jsonString)                                     # write string to file
                jsonFile.close()                                               # close file
    
                # # read regression data (dictionary) back from .json          # debug only
                # fileObject = open(filepathstr + '_f0_LinearRegress.json', 'r') # debug only
                # jsonContent = fileObject.read()                                        # debug only
                # data_back = json.loads(jsonContent)                                    # debug only

                # save regression data (dictionary) as .npy
                np.save(filepathstr + '_f0_LinearRegress.npy', data)  
                
                # # read regression data (dictionary) back from .npy      # debug only
                # data_back = np.load(filepathstr + '_f0_LinearRegress.npy', allow_pickle=True).item()  # debug only
                
            # --- output figures ---
            logfigurepath = Path(logfilepath + 'images/')   # create a Path object        
            if not logfigurepath.is_dir():                  # Check and see if this Path object points to an existing directory. If not, make a new directory.
                logfigurepath.mkdir()
            logfigurepath = str(logfigurepath) + '/'        # convert back to a regular string
            
            if outputFigureStimSpectrogramFlag:
                fig101.savefig(logfigurepath + file.stimToken + '_spectrogram.jpg', dpi=100)  
                
            if outputFigureStimf0Flag:
                fig103.savefig(logfigurepath + file.stimToken + '_f0contour.jpg', dpi=100)  
    
            if outputFigureRespSpectrogramFlag:
                fig102.savefig(logfigurepath + filenamestr + '_spectrogram.jpg', dpi=100)  
    
            if outputFigureRespf0Flag:
                fig104.savefig(logfigurepath + filenamestr + '_f0contour.jpg', dpi=100)  
    
            if outputFiguref0ComparisonFlag:
                fig105.savefig(logfigurepath + filenamestr + '_f0_Comparison.jpg', dpi=100)  
    
            if outputFigureRegressionFlag:
                fig106.savefig(logfigurepath + filenamestr + '_f0_Regression.jpg', dpi=100)  
            
            plt.pause(1)  # pause for a second, allowing the user time to view the results. Most importantly, this pause() allows/forces the figures to be displayed, before being closed at the beginning of the next iteration.

            if npolarities > 1: print('Done with polarity ' + polarityNames[polarity], flush=True)
        if nchans > 1: print('Done with chan ' + str(chan), flush=True)
    print('Done with ' + filename, flush=True)












def estimatef0(hd='C:', logfilepath='C:/data/GroupData_NMF_feasibility_python/', logfilename='NMF_feasibility_log.json', file='', win=50.0, gap=0.5, freqResolution=1.0, upperFreq=1000.0, low=92.4, high=149.8, se_absFlag=False, r_method=1, f0ContourPath=''):
    '''
    Purpose : estimatef0 (spectral maximum method) for FFR recordings
    
    parameters
    ----------
    hd : string
        - hard drive
    logiflepath : string
        - file path for log book
    logfilename : string
        - file name for log book
    file : a core.frame.Pandas object 
        - a file to be processed
    win : float64
        - window size in ms for calculating spectrogram 
    gap : float64
        - step size in ms for calculating spectrogram
    freqResolution : float64
        - frequency resolution in Hz for calculating spectrogram 
    upperFreq : float64
        - the highest frequency to take out from the spectrogram output [default 1000]
    low : float
        - the lower boundary when searching candidates for f0 location
    high : float
        - the higher boundary when searching candidates for f0 location
    se_absFlag : boolean
        - flag indicating whether to take an absolute value of se, or not (True yes, False no) 
    r_method : integer
        - method to calculate Tracking Accuracy (0 ==> linear regression, r between [0 1]; 1 ==> correlcation coefficient, r between [-1 +1])
    f0ContourPath : string
        - file path where f0 contours are saved
    
    Returns
    -------
    None

    Diary
    -----
    2021-09-09 (v01) Fuh
        - rewrote _03estimatef0 script in Python
    2021-10-15 (v02)
        - rewrote this script as a function
        
    ToDo
    ----
    None
        
    '''
    
    # clear console, clear variables
    from IPython import get_ipython         # The 'IPython' is a Python built-in module, for interactive computing in different languages
    get_ipython().magic('clear')            # clear: clear the console 
    get_ipython().magic('reset -f')         # -f : force reset without asking for confirmation (i.e., a 'hard' reset of the kernel)
    
    # close all figures
    import matplotlib.pyplot as plt         
    plt.close('all')
    
    # import modules
    from pathlib import Path                    # pathlib is a built-in library in since Python 3.4. The 'Path' class is a new, correct, and easy way to work with file paths and file names.
    from scipy import signal                    # for resample_poly(), firwin(), lfilter(), etc.
    from scipy.io import wavfile                # for wavfile.read()
    import sys                                  # import the Python built-in module 'sys' for sys.exit(), sys.platform
    # import pandas as pd                         # import Pandas for DataFrame
    import numpy as np
    import json
    
    # # initialize parameters
    # hd = 'E:'                               # hard drive (Windows)
    # if sys.platform == 'darwin':            # macOS
    #     # hd = '/Volumes/AEPLAB'            # Fuh Mac computer, VPN S drive
    #     # hd = '/Users/groverlab'           # Fuh Mac computer, local HD
    #     hd = '/Users/groverlab/OneDrive - Ohio University'   # Fuh Mac computer, OneDrive - Ohio University
    #     # hd = '/Volumes/FUHBACKUP'         # Fuh Mac computer, USB thumb drive
    
    # # --- read in par ---
    # logfilepath = hd + '/data/GroupData_NMF_feasibility_python/'   
    # logfilename = 'NMF_feasibility_log.json'
    # logfilepathname = logfilepath + logfilename             
    # par = pd.read_json(logfilepathname, orient='table')
    
    # hide the toolbar when showing new figures
    plt.rcParams['toolbar'] = 'None'       # This will hide toolbar for all figures, unless it is turned back on. ['toolbar2' (default), 'None', 'toolmanager']
    
    # define a dictionary for converting numbers to color maps (for ploting spectrograms)
    dict_colorMap = {
        0: 'gray_r',
        1: 'jet',
        }
    
    # define a dictionary for converting stimToken to a tuple (stimpath, stimname)
    dict_stimToken = {
        'i2_150ms': ('/Stim2/RecordedToneStimuli/i/origcycN_150ms/', 'i2_150ms.wav'),   
        'i2_250ms': ('/Stim2/RecordedToneStimuli/i/origcycN/', 'i2_250ms.wav'),   
        'i2_CycAmpN': ('/Stim2/RecordedToneStimuli/i/cycAmpN/', 'i2_CycAmpN.wav'),   
        }
    
    # stores a list of stimTokens that have been processed in this subroutine
    # TODO: return the accumulated stimTokens list back to the calling script, so that the same stimToken will not be processed twice.
    stimTokens = []
    
    # data flags
    outputStimf0Flag = True
    outputRespf0Flag = True
    outputf0RegressionFlag = True
    
    # figure flags
    outputFigureStimSpectrogramFlag = False
    outputFigureRespSpectrogramFlag = False
    outputFigureStimf0Flag = False
    outputFigureRespf0Flag = False
    outputFigureRegressionFlag = False
    outputFiguref0ComparisonFlag = True
    setFigureLocationFlag = True        
    space = 10                                  # space between adjacent figures in pixels 
    dx = 640                                    # width of a figure in pixels
    dy = 480                                    # height of a figure in pixels
    colorMap_stimulus = 1                       # define color map for stimulus spectrogram
    colorMap_recording = 1                      # define color map for recording spectrogram
    
    # # spectrogram parameters
    # win = 50.0             # window size in ms for calculating spectrogram 
    # gap = 0.5              # step size in ms for calculating spectrogram [changed from 1 to 0.5 ms for 170 ms di1234 tokens]
    # centerFlag = True      # time array is based on the center point of each time bin, or not [default = 1]
    # freqResolution = 1.0   # have a frequency resolution of 0.1 Hz here to increase accuracy in estimating f0 --- it will take a long time to run through, but the data will be worth it.
    #                        # 20000/(2^12) # frequency resolution in Hz for narrow-band spectrogram 
    #                        # old default freqResolution = fs / (2^12) = 20000 / (2^12) = 4.88 Hz
    #                        # new default freqResolution = fs / fs = 1 Hz
    # upperFreq = 1000.0     # frequency range to take out from the spectrogram output [default 1000, i.e., take out 1-1000 Hz data out of a total of a 10,000 Hz data range
    
    backend = plt.get_backend()  # get graphics backend
    
    polarityNames = ['']                  # ['', '_a', '_b']   # a list of polarity names to be included for analysis
    npolarities = len(polarityNames)      # number of polarity conditions
    
    # low = 92.4             # 102.4 from "f0_contours_i1234_150ms_origCycN_spec_150ms.txt"
    # high = 149.8           # 139.8
    
    # se_absFlag = False     # flag to take an absolute value of se (1 yes, 0 no)
    
    # r_method = 1           # method to calculate Tracking Accuracy (0 ==> linear regression, r between [0 1] 1 ==> correlcation coefficient, r between [-1 +1])
    
    # f0ContourPath = hd + '/Stim2/RecordedToneStimuli/i/f0contours_150ms/'
            
    # # clear all local variables, except a few, before a new iteration starts. 
    # keepList = ['get_ipython', 'plt', 'Path', 'signal', 'wavfile', 'sys', 'pd', 'np', 'json', 
    #             'hd', 'par', 'file', 'space', 'dx', 'dy', 'win', 'gap', 'freqResolution', 'upperFreq', 'backend', 'polarityNames', 'npolarities',
    #             'low', 'high', 'r_method', 'f0ContourPath']  
    # keepPrefix = ('keep', 'var', 'log', 'dict', 'colorMap', 'stim')     # keep keep*, var*, log*, dict*, etc
    # keepSuffix = ('Flag')  # keep *Flag, etc
    # varListLocal = list(locals().keys())    # Get local variables ==> dictionary keys ==> convert dict into a list
    # for var in varListLocal:
    #     if (var not in keepList) and (not var.startswith(keepPrefix)) and (not var.endswith(keepSuffix)) and (var[0] != '_'):   # if (1) var is not in the keep list, (2) var does not start with a certain prefix, (3) var does not end with a certain suffix, and (4) var is not a system 'hidden' variable with a prefix '_', ==> delete it.
    #         exec(f'del {var}')        

    # define filepath, filename, and related variables
    filepath = hd + '/data/' + file.directory + '/'
    filepre = file.subjectCode + 'ses' + str(file.sessionNumber) + '_'  
    filenamestr = filepre + file.fileNumberString
    filepathstr = filepath + filepre + file.fileNumberString 
    
    # read in header
    if file.labviewIndex == 1:  # whether the recording is made from Labview or not (1 yes, 0 no). When no, it is recorded by using Neuroscan.
        header = np.load(filepathstr + '_header.npy', allow_pickle=True).item()
        # jsonFile = open(filepathstr + '_header.json', 'r')
        # header = json.load(jsonFile)
        # jsonFile.close()
        
        softwareversion = header['softwareversion']
        if softwareversion >= 2.17:        
            nchans = int(header['AI_nchannels'])
            if nchans <= 0:                 # nchans must be >= 1
                nchans = 1                  
            fs = header['AI_samplingRate']
            si = header['AO_silentInterval'] / 1000.0
            isi = header['AI_nsamplesReadPerChannel'] / fs 
            stimPolarity = int(header['stim_polarity'])
            if stimPolarity == 2:        # alternating polarity
                stimdur = (isi / 2) - si
            else:                        # original or inverted polarity
                stimdur = isi - si
            dBSPL = header['AO0_dBSPL'] 
            tailFlag = int(header['tail_flag'])
            stimPathName = header['AO0_stimPath']
            stimdur_ms = stimdur * 1000.0              # stim dur in ms
    
        else:
            sys.exit('not recognizable AEP software version...')            
    else:    
        fs = 20000.0
        stimdur = 0.150            
        stimdur_ms = stimdur * 1000.0
        si = file.si / 1000.0      
        isi = stimdur + si
        stimPolarity = 2
        
    # -------------- initialize parameters -------------------------------
    # nrep = file.nRepetitions  # number of repetitions
    # overlap = win - gap     
    # window = np.hanning(np.floor(win * fs / 1000.0))      # window is a hanning window of 50 ms
    # noverlap = int(np.floor(overlap * fs / 1000.0))       # overlaps noverlap samples of each segment, 
    # nfft = int(np.floor(win * fs / 1000.0))               # default 1000
    # pad_to = int(np.floor((1.0 / freqResolution) * fs))   # convert freqResolution to number of data points [default nfft = 2^12]. 
    # #         # zero padding each segment to have nfft number of data points 
    # #         # when nfft = fs, frequency resolution = 1 Hz

    # take only 1 out of every 'interval' data points along the frequency axis, so that frequency resolution is 10 Hz on the spectrogram plot (for Python and Veusz).
    freq_stepsize = 10.0                                      # step size on the frequency axis [default freq step size = 10 Hz]
    interval = int(np.floor(freq_stepsize / freqResolution))  # make frequency interval to be 10 Hz, for Phtyon and Veusz plotting purposes

    # ---- read stimulus token, first occurrence only ------------
    stimToken = file.stimToken              # stimulus token entered on MATLAB logbook, column 8
    if dict_stimToken.get(stimToken) is not None:
        stimpath, stimname = dict_stimToken[stimToken]
    else:
        sys.exit('not recognizable stimulus token...')
    
    stimfs, stim = wavfile.read(hd + stimpath + stimname)

    # make stimulus token to have the same sampling rate as recording fs [default = 20000]
    up = 100                                            # set 'up' sampling at 100 (artibutrarily)
    down = int(up * (stimfs / fs))                 # compute 'down' sampling, to get desired new fs (relative to the old fs)
    stim = signal.resample_poly(stim, up, down)         # resample y to have a total of npts data points
    if len(stim) > int(np.floor(fs * stimdur)):
        stim = stim[0 : int(np.floor(fs * stimdur))]     
        
    # normalize stim amplitudes (time waveform) 
    stim = stim / np.amax(stim)

    # -------------Read in stimulus and recording spectrograms ----------
    for  chan in range(nchans):          
        for polarity in range(npolarities):    
            
            # close all figures
            plt.close('all')   
            
            # append chanstr and polaritystr to filepathstr and filenamestr
            chanstr = '' if nchans==1 else '_chan' + str(chan)
            polaritystr = polarityNames[polarity]                  
            filepathstr = filepath + filepre + file.fileNumberString + chanstr + polaritystr  
            filenamestr = filepre + file.fileNumberString + chanstr + polaritystr
            filename = filenamestr + '.npy'

            annotationStimulus = stimToken
            annotationRecording = filename + ', ' + file.experimentString
                       
            print('Reading file: ' + filepath + filename, flush=True)
    
            # --- read in stimulus spectrogram data ---
            T_ms = np.load(logfilepath + file.stimToken + '_stim_spectrogram_time_ms.npy')    
            F = np.load(logfilepath + file.stimToken + '_stim_spectrogram_freq.npy') 
            nV_S = np.load(logfilepath + file.stimToken + '_stim_spectrogram_nV.npy')     
            nV_S_veusz = np.load(logfilepath + file.stimToken + '_stim_spectrogram_nV_veusz.npy')  
    
            # compute (left, right, bottom, top) boundaries for imshow(), for stimulus and response
            left = win/2.0 - gap/2.0              
            right = (stimdur_ms-win/2.0) + gap/2.0
            bottom = -(freqResolution/2.0) * interval
            top = F[-1] + (freqResolution/2.0)*interval
            
            # define x tick labels
            if int(stimdur*1000) == 150:
                xticks = (25,50,75,100,125)
            elif int(stimdur*1000) == 250:
                xticks = (25,75,125,175,225)
            else:
                sys.exit('Unrecognizable stimulus duration')
                      
            # --- plot stimulus spectrogram ---
            fig101 = plt.figure(101)
            plt.imshow(nV_S_veusz, origin='lower', extent=(left, right, bottom, top), aspect='auto', cmap=dict_colorMap.get(colorMap_stimulus)) # , vmin=0, vmax=8000)   # set zero origin at lower left corner, define the extent of axis boundaries (so that the x and y tick numbers will be correct), set aspect ratio to 'auto' (so that the spectrogram image will fill in the entire space), set color map to 'jet' for color and 'gray_r' for gray scale, set vmin and vmax for the lower and upper limits along the z axis (i.e., amplitude axis)
            plt.title(annotationStimulus)     
            plt.xlabel('Time (ms)')
            plt.ylabel('Frequency (Hz)')
            plt.xticks(ticks=xticks)
            # plt.colorbar()               
            if setFigureLocationFlag:
                if 'Qt' in backend:          # Qt (e.g., Qt5Agg)
                    mgr = plt.get_current_fig_manager()              # get a manager for current figure
                    mgr.window.setGeometry(space, space, dx, dy)     # set geometry of current figure (x, y, width, height) in pixels
                else:
                    sys.exit('Unrecognizable graphical backend')
            plt.show(block=False)    
            
            # --- read in response spectrogram data ---
            t_ms = np.load(filepathstr + '_spectrogram_time_ms.npy')    
            f = np.load(filepathstr + '_spectrogram_freq.npy')    
            nV = np.load(filepathstr + '_spectrogram_nV.npy')   
            nV_veusz = np.load(filepathstr + '_spectrogram_nV_veusz.npy')  
    
            # --- plot response spectrogram ---
            fig102 = plt.figure(102)
            plt.imshow(nV_veusz, origin='lower', extent=(left, right, bottom, top), aspect='auto', cmap=dict_colorMap.get(colorMap_recording)) # , vmin=0, vmax=8000)   # set zero origin at lower left corner, define the extent of axis boundaries (so that the x and y tick numbers will be correct), set aspect ratio to 'auto' (so that the spectrogram image will fill in the entire space), set color map to 'jet' for color and 'gray_r' for gray scale, set vmin and vmax for the lower and upper limits along the z axis (i.e., amplitude axis)
            plt.title(annotationRecording)
            plt.xlabel('Time (ms)')
            plt.ylabel('Frequency (Hz)')
            plt.xticks(ticks=xticks)
            # plt.colorbar()               
            if setFigureLocationFlag:
                if 'Qt' in backend:          # Qt (e.g., Qt5Agg)
                    mgr = plt.get_current_fig_manager()                     # get a manager for current figure
                    mgr.window.setGeometry(space+dx+space, space, dx, dy)   # set geometry of current figure (x, y, width, height) in pixels
                else:
                    sys.exit('Unrecognizable graphical backend')
            plt.show(block=False)
            
            ###############################
            # --- Find spectral maximum ---
            ###############################
                
            #----------------- extract f0s for stimulus ----------------
            if stimToken not in stimTokens:           # compute stimulus f0s, only upon the first encounter of a stimToken
                stimTokens.append(stimToken)
                highIndices = np.argwhere(F <= high)  # find the upper frequency index of the searching range
                highIndex = highIndices[-1, 0] + 1
                lowIndices = np.argwhere(F >= low)    # find the lower frequency index of the searching range
                lowIndex = lowIndices[0, 0]
                
                F_range = F[lowIndex:highIndex]
                nV_S_range = nV_S[lowIndex:highIndex, :]
                
                # in the range of T_ms, find the maximum y and its index
                f0amps = np.full(T_ms.shape, None)
                indices = np.full(T_ms.shape, None)
                for i in range(len(T_ms)):
                    f0amps[i] = np.amax(nV_S_range[:, i])
                    indices[i] = np.argmax(nV_S_range[:, i])
                    
                # use the index to find the corresponding frequency
                f0s = np.full(T_ms.shape, None)
                for j in range(len(T_ms)):
                    f0s[j] = F_range[indices[j]]
                f0s = f0s.astype(float)            # force f0s to float64
                    
                # output f0s
                np.save(f0ContourPath + 'f0_' + stimToken + '.npy', f0s)  
            else:
                f0s = np.load(f0ContourPath + 'f0_' + stimToken + '.npy', allow_pickle=True)  
                                    
            # --- plot stimulus f0s ---
            fig103 = plt.figure(103)
            plt.plot(T_ms, f0s, 'k')         # {b: blue, k: black, r: red, etc}, {o: circle, +: plus, x: cross, etc}
            plt.title('F0s of the stimulus, ' + stimToken)
            plt.xlabel('Time (ms)')
            plt.ylabel('Frequency (Hz)')
            if setFigureLocationFlag:
                if 'Qt' in backend:          # Qt (e.g., Qt5Agg)
                    mgr = plt.get_current_fig_manager()              # get a manager for current figure
                    mgr.window.setGeometry(space, int(space+dy*6/10+space), dx, dy)     # set geometry of current figure (x, y, width, height) in pixels
                elif 'wx' in backend:        
                    print('wx backend')   
                else:
                    sys.exit('Unrecognizable graphical backend')
            plt.show(block=False)
    
            #----------------- extract f0s for a recording ----------------
            highIndices = np.argwhere(f <= high)  # find the upper frequency index of the searching range
            highIndex = highIndices[-1, 0] + 1
            lowIndices = np.argwhere(f >= low)    # find the lower frequency index of the searching range
            lowIndex = lowIndices[0, 0]
            
            f_range = f[lowIndex:highIndex]
            nV_range = nV[lowIndex:highIndex, :]
            
            # in the range of T_ms, find the maximum y and its index
            f0amp = np.full(t_ms.shape, None)
            index = np.full(t_ms.shape, None)
            for i in range(len(t_ms)):
                f0amp[i] = np.amax(nV_range[:, i])
                index[i] = np.argmax(nV_range[:, i])
                
            # use the index to find the corresponding frequency
            f0 = np.full(t_ms.shape, None)
            for j in range(len(t_ms)):
                f0[j] = f_range[index[j]]
            f0 = f0.astype(float)             # force f0 to float64
            
            # calculate Spectral Amplitude (SA)
            sa = np.mean(abs(f0amp))        
    
            # --- plot response f0s ---
            fig104 = plt.figure(104)
            plt.plot(t_ms, f0, 'r')
            plt.title('F0s of a recording, ' + annotationRecording)
            plt.xlabel('Time (ms)')
            plt.ylabel('Frequency (Hz)')
            if setFigureLocationFlag:
                if 'Qt' in backend:          # Qt (e.g., Qt5Agg)
                    mgr = plt.get_current_fig_manager()              # get a manager for current figure
                    mgr.window.setGeometry(space+dx+space, int(space+dy*6/10+space), dx, dy)     # set geometry of current figure (x, y, width, height) in pixels
                elif 'wx' in backend:        
                    print('wx backend')   
                else:
                    sys.exit('Unrecognizable graphical backend')
            plt.show(block=False)
                
            # ---- comparison of f0 contours ----
            fig105 = plt.figure(105)
            plt.plot(T_ms, f0s, 'k', label='Stimulus')
            plt.plot(t_ms, f0, 'r', label='Response')
            plt.title('F0 contours, ' + annotationRecording)
            plt.xlabel('Time (ms)')
            plt.ylabel('Frequency (Hz)')
            plt.legend()
            if setFigureLocationFlag:
                if 'Qt' in backend:          # Qt (e.g., Qt5Agg)
                    mgr = plt.get_current_fig_manager()              # get a manager for current figure
                    mgr.window.setGeometry(space*2, space*8, dx, dy)     # set geometry of current figure (x, y, width, height) in pixels
                elif 'wx' in backend:        
                    print('wx backend')   
                else:
                    sys.exit('Unrecognizable graphical backend')
            plt.show(block=False)
    
            # --- compute Frequency Error (FE) ---
            X = np.vstack([f0s, np.ones(len(f0s))]).transpose()   # append ones column to x               
            (b, _, _, _) = np.linalg.lstsq(X, f0, rcond=None)   # perform linear regression       
            yestimate = b[0] * f0s + b[1]    # compute y estimate      
            fe = np.mean(abs(yestimate - f0))     # compute frequency error (fe)
            
            # --- compute Tracking Accuracy (TA) ---
            r = np.corrcoef(f0s, f0)[0, 1]      # correlation coefficient (TA) between [-1.0 1.0]
            if r_method == 0:
                r = abs(r)
            R2 = r ** 2                      # R2 = r^2 = coefficient of determination
    
            fig106 = plt.figure(106)
            plt.plot(f0s, f0, 'r', label='F0s')
            plt.plot(f0s, yestimate, 'b', label='Fitted line')
            plt.title('Linear regression, ' + annotationRecording)
            plt.xlabel('F0 of stimulus (Hz)')
            plt.ylabel('F0 of response (Hz)')
            plt.legend()
            if setFigureLocationFlag:
                if 'Qt' in backend:          # Qt (e.g., Qt5Agg)
                    mgr = plt.get_current_fig_manager()              # get a manager for current figure
                    mgr.window.setGeometry(space*4 + dx, space*8, dx, dy)     # set geometry of current figure (x, y, width, height) in pixels
                elif 'wx' in backend:        
                    print('wx backend')   
                else:
                    sys.exit('Unrecognizable graphical backend')
            plt.show(block=False)
            
            # --- compute Slope Error (SE) ---
            X = np.vstack([t_ms, np.ones(len(t_ms))]).transpose()     # append ones column to x               
            (b_s, _, _, _) = np.linalg.lstsq(X, f0s, rcond=None)      # stimulus f0 regression
            (b_ffr, _, _, _) = np.linalg.lstsq(X, f0, rcond=None)     # FFR f0 regression
            se = (b_ffr[0] - b_s[0]) * 1000.0                         # slope error in Hz/sec
            if se_absFlag:
                se = abs(se)
    
            # --- compute RMS amplitude ---
            resp = np.load(filepathstr + '_averaged_' + str(int(stimdur_ms)) + 'ms_recording.npy', allow_pickle=True)    # read time waveform
            resp = resp.astype(float)           # force to float64 precision
            rms = np.sqrt(np.mean(resp**2))     # compute RMS amplitude
    
            # display results on console 
            print('Frequency Error (Hz) = %3.8f' % fe)
            print('Slope Error (Hz/s) = %3.8f' % se)
            print('Tracking Accuracy = %3.8f' % r)
            print('Spectral Amplitude (nV) = %3.8f' % sa)
            print('RMS Amplitude (nV) = %3.8f' % rms)
            
            # --- store results in data (i.e., a dictionary)---
            data = {
                'filepathstr': filepathstr,
                'stimToken': stimToken,
                'fe': fe,
                'se': se,
                'ta': r,
                'sa': sa,
                'rms': rms,
                'b_s_slope': b_s[0] * 1000.0,      # slope of stimulus f0s in Hz/second
                'b_s_intercept': b_s[1],    
                'b_ffr_slope': b_ffr[0] * 1000.0,  # slope of FFR f0s in Hz/second
                'b_ffr_intercept': b_ffr[1],  
                'b0': b[0],                   
                'b1': b[1],    
                'R2': R2,                          # coefficient of determination   
                }
        
            # --- Output data ---
            if outputStimf0Flag and file[0] == 0:  # stimulus f0s
                np.save(logfilepath + stimToken + '_stimf0.npy', f0s)
            
            if outputRespf0Flag:                   # response f0s and f0 amplitudes
                np.save(filepathstr + '_respf0.npy', f0)
                np.save(filepathstr + '_f0amp.npy', f0amp)
            
            if outputf0RegressionFlag:             # regression data
                # save regression data (dictionary) as .json
                jsonString = json.dumps(data)                                  # convert dict to a json string
                jsonFile = open(filepathstr + '_f0_LinearRegress.json', 'w')   # open a json file
                jsonFile.write(jsonString)                                     # write string to file
                jsonFile.close()                                               # close file
    
                # # read regression data (dictionary) back from .json          # debug only
                # fileObject = open(filepathstr + '_f0_LinearRegress.json', 'r') # debug only
                # jsonContent = fileObject.read()                                        # debug only
                # data_back = json.loads(jsonContent)                                    # debug only

                # save regression data (dictionary) as .npy
                np.save(filepathstr + '_f0_LinearRegress.npy', data)  
                
                # # read regression data (dictionary) back from .npy      # debug only
                # data_back = np.load(filepathstr + '_f0_LinearRegress.npy', allow_pickle=True).item()  # debug only
                
            # --- output figures ---
            logfigurepath = Path(logfilepath + 'images/')   # create a Path object        
            if not logfigurepath.is_dir():                  # Check and see if this Path object points to an existing directory. If not, make a new directory.
                logfigurepath.mkdir()
            logfigurepath = str(logfigurepath) + '/'        # convert back to a regular string
            
            if outputFigureStimSpectrogramFlag:
                fig101.savefig(logfigurepath + file.stimToken + '_spectrogram.jpg', dpi=100)  
                
            if outputFigureStimf0Flag:
                fig103.savefig(logfigurepath + file.stimToken + '_f0contour.jpg', dpi=100)  
    
            if outputFigureRespSpectrogramFlag:
                fig102.savefig(logfigurepath + filenamestr + '_spectrogram.jpg', dpi=100)  
    
            if outputFigureRespf0Flag:
                fig104.savefig(logfigurepath + filenamestr + '_f0contour.jpg', dpi=100)  
    
            if outputFiguref0ComparisonFlag:
                fig105.savefig(logfigurepath + filenamestr + '_f0_Comparison.jpg', dpi=100)  
    
            if outputFigureRegressionFlag:
                fig106.savefig(logfigurepath + filenamestr + '_f0_Regression.jpg', dpi=100)  
            
            plt.pause(1)  # pause for a second, allowing the user time to view the results. Most importantly, this pause() allows/forces the figures to be displayed, before being closed at the beginning of the next iteration.

            if npolarities > 1: print('Done with polarity ' + polarityNames[polarity], flush=True)
        if nchans > 1: print('Done with chan ' + str(chan), flush=True)
    print('Done with ' + filename, flush=True)
Purpose : estimatef0 (spectral maximum method) for FFR recordings          parameters     ----------     hd : string         - hard drive     logiflepath : string         - file path for log book     logfilename : string         - file name for log book     file : a core.frame.Pandas object          - a file to be processed     win : float64         - window size in ms for calculating spectrogram      gap : float64         - step size in ms for calculating spectrogram     freqResolution : float64         - frequency resolution in Hz for calculating spectrogram      upperFreq : float64         - the highest frequency to take out from the spectrogram output [default 1000]     low : float         - the lower boundary when searching candidates for f0 location     high : float         - the higher boundary when searching candidates for f0 location     se_absFlag : boolean         - flag indicating whether to take an absolute value of se, or not (True yes, False no)      r_method : integer         - method to calculate Tracking Accuracy (0 ==> linear regression, r between [0 1]; 1 ==> correlcation coefficient, r between [-1 +1])     f0ContourPath : string         - file path where f0 contours are saved          Returns     -------     None      Diary     -----     2021-09-09 (v01) Fuh         - rewrote _03estimatef0 script in Python     2021-10-15 (v02)         - rewrote this script as a function              ToDo     ----     None # clear console, clear variables # The 'IPython' is a Python built-in module, for interactive computing in different languages # clear: clear the console # -f : force reset without asking for confirmation (i.e., a 'hard' reset of the kernel) # close all figures # import modules # pathlib is a built-in library in since Python 3.4. The 'Path' class is a new, correct, and easy way to work with file paths and file names. # for resample_poly(), firwin(), lfilter(), etc. # for wavfile.read() # import the Python built-in module 'sys' for sys.exit(), sys.platform # import pandas as pd                         # import Pandas for DataFrame # # initialize parameters # hd = 'E:'                               # hard drive (Windows) # if sys.platform == 'darwin':            # macOS #     # hd = '/Volumes/AEPLAB'            # Fuh Mac computer, VPN S drive #     # hd = '/Users/groverlab'           # Fuh Mac computer, local HD #     hd = '/Users/groverlab/OneDrive - Ohio University'   # Fuh Mac computer, OneDrive - Ohio University #     # hd = '/Volumes/FUHBACKUP'         # Fuh Mac computer, USB thumb drive # # --- read in par --- # logfilepath = hd + '/data/GroupData_NMF_feasibility_python/' # logfilename = 'NMF_feasibility_log.json' # logfilepathname = logfilepath + logfilename # par = pd.read_json(logfilepathname, orient='table') # hide the toolbar when showing new figures # This will hide toolbar for all figures, unless it is turned back on. ['toolbar2' (default), 'None', 'toolmanager'] # define a dictionary for converting numbers to color maps (for ploting spectrograms) # define a dictionary for converting stimToken to a tuple (stimpath, stimname) # stores a list of stimTokens that have been processed in this subroutine # TODO: return the accumulated stimTokens list back to the calling script, so that the same stimToken will not be processed twice. # data flags # figure flags # space between adjacent figures in pixels # width of a figure in pixels # height of a figure in pixels # define color map for stimulus spectrogram # define color map for recording spectrogram # # spectrogram parameters # win = 50.0             # window size in ms for calculating spectrogram # gap = 0.5              # step size in ms for calculating spectrogram [changed from 1 to 0.5 ms for 170 ms di1234 tokens] # centerFlag = True      # time array is based on the center point of each time bin, or not [default = 1] # freqResolution = 1.0   # have a frequency resolution of 0.1 Hz here to increase accuracy in estimating f0 --- it will take a long time to run through, but the data will be worth it. #                        # 20000/(2^12) # frequency resolution in Hz for narrow-band spectrogram #                        # old default freqResolution = fs / (2^12) = 20000 / (2^12) = 4.88 Hz #                        # new default freqResolution = fs / fs = 1 Hz # upperFreq = 1000.0     # frequency range to take out from the spectrogram output [default 1000, i.e., take out 1-1000 Hz data out of a total of a 10,000 Hz data range # get graphics backend # ['', '_a', '_b']   # a list of polarity names to be included for analysis # number of polarity conditions # low = 92.4             # 102.4 from "f0_contours_i1234_150ms_origCycN_spec_150ms.txt" # high = 149.8           # 139.8 # se_absFlag = False     # flag to take an absolute value of se (1 yes, 0 no) # r_method = 1           # method to calculate Tracking Accuracy (0 ==> linear regression, r between [0 1] 1 ==> correlcation coefficient, r between [-1 +1]) # f0ContourPath = hd + '/Stim2/RecordedToneStimuli/i/f0contours_150ms/' # # clear all local variables, except a few, before a new iteration starts. # keepList = ['get_ipython', 'plt', 'Path', 'signal', 'wavfile', 'sys', 'pd', 'np', 'json', #             'hd', 'par', 'file', 'space', 'dx', 'dy', 'win', 'gap', 'freqResolution', 'upperFreq', 'backend', 'polarityNames', 'npolarities', #             'low', 'high', 'r_method', 'f0ContourPath'] # keepPrefix = ('keep', 'var', 'log', 'dict', 'colorMap', 'stim')     # keep keep*, var*, log*, dict*, etc # keepSuffix = ('Flag')  # keep *Flag, etc # varListLocal = list(locals().keys())    # Get local variables ==> dictionary keys ==> convert dict into a list # for var in varListLocal: #     if (var not in keepList) and (not var.startswith(keepPrefix)) and (not var.endswith(keepSuffix)) and (var[0] != '_'):   # if (1) var is not in the keep list, (2) var does not start with a certain prefix, (3) var does not end with a certain suffix, and (4) var is not a system 'hidden' variable with a prefix '_', ==> delete it. #         exec(f'del {var}') # define filepath, filename, and related variables # read in header # whether the recording is made from Labview or not (1 yes, 0 no). When no, it is recorded by using Neuroscan. # jsonFile = open(filepathstr + '_header.json', 'r') # header = json.load(jsonFile) # jsonFile.close() # nchans must be >= 1 # alternating polarity # original or inverted polarity # stim dur in ms # -------------- initialize parameters ------------------------------- # nrep = file.nRepetitions  # number of repetitions # overlap = win - gap # window = np.hanning(np.floor(win * fs / 1000.0))      # window is a hanning window of 50 ms # noverlap = int(np.floor(overlap * fs / 1000.0))       # overlaps noverlap samples of each segment, # nfft = int(np.floor(win * fs / 1000.0))               # default 1000 # pad_to = int(np.floor((1.0 / freqResolution) * fs))   # convert freqResolution to number of data points [default nfft = 2^12]. # #         # zero padding each segment to have nfft number of data points # #         # when nfft = fs, frequency resolution = 1 Hz # take only 1 out of every 'interval' data points along the frequency axis, so that frequency resolution is 10 Hz on the spectrogram plot (for Python and Veusz). # step size on the frequency axis [default freq step size = 10 Hz] # make frequency interval to be 10 Hz, for Phtyon and Veusz plotting purposes # ---- read stimulus token, first occurrence only ------------ # stimulus token entered on MATLAB logbook, column 8 # make stimulus token to have the same sampling rate as recording fs [default = 20000] # set 'up' sampling at 100 (artibutrarily) # compute 'down' sampling, to get desired new fs (relative to the old fs) # resample y to have a total of npts data points # normalize stim amplitudes (time waveform) # -------------Read in stimulus and recording spectrograms ---------- # close all figures # append chanstr and polaritystr to filepathstr and filenamestr # --- read in stimulus spectrogram data --- # compute (left, right, bottom, top) boundaries for imshow(), for stimulus and response # define x tick labels # --- plot stimulus spectrogram --- # , vmin=0, vmax=8000)   # set zero origin at lower left corner, define the extent of axis boundaries (so that the x and y tick numbers will be correct), set aspect ratio to 'auto' (so that the spectrogram image will fill in the entire space), set color map to 'jet' for color and 'gray_r' for gray scale, set vmin and vmax for the lower and upper limits along the z axis (i.e., amplitude axis) # plt.colorbar() # Qt (e.g., Qt5Agg) # get a manager for current figure # set geometry of current figure (x, y, width, height) in pixels # --- read in response spectrogram data --- # --- plot response spectrogram --- # , vmin=0, vmax=8000)   # set zero origin at lower left corner, define the extent of axis boundaries (so that the x and y tick numbers will be correct), set aspect ratio to 'auto' (so that the spectrogram image will fill in the entire space), set color map to 'jet' for color and 'gray_r' for gray scale, set vmin and vmax for the lower and upper limits along the z axis (i.e., amplitude axis) # plt.colorbar() # Qt (e.g., Qt5Agg) # get a manager for current figure # set geometry of current figure (x, y, width, height) in pixels ############################### # --- Find spectral maximum --- ############################### #----------------- extract f0s for stimulus ---------------- # compute stimulus f0s, only upon the first encounter of a stimToken # find the upper frequency index of the searching range # find the lower frequency index of the searching range # in the range of T_ms, find the maximum y and its index # use the index to find the corresponding frequency # force f0s to float64 # output f0s # --- plot stimulus f0s --- # {b: blue, k: black, r: red, etc}, {o: circle, +: plus, x: cross, etc} # Qt (e.g., Qt5Agg) # get a manager for current figure # set geometry of current figure (x, y, width, height) in pixels #----------------- extract f0s for a recording ---------------- # find the upper frequency index of the searching range # find the lower frequency index of the searching range # in the range of T_ms, find the maximum y and its index # use the index to find the corresponding frequency # force f0 to float64 # calculate Spectral Amplitude (SA) # --- plot response f0s --- # Qt (e.g., Qt5Agg) # get a manager for current figure # set geometry of current figure (x, y, width, height) in pixels # ---- comparison of f0 contours ---- # Qt (e.g., Qt5Agg) # get a manager for current figure # set geometry of current figure (x, y, width, height) in pixels # --- compute Frequency Error (FE) --- # append ones column to x # perform linear regression # compute y estimate # compute frequency error (fe) # --- compute Tracking Accuracy (TA) --- # correlation coefficient (TA) between [-1.0 1.0] # R2 = r^2 = coefficient of determination # Qt (e.g., Qt5Agg) # get a manager for current figure # set geometry of current figure (x, y, width, height) in pixels # --- compute Slope Error (SE) --- # append ones column to x # stimulus f0 regression # FFR f0 regression # slope error in Hz/sec # --- compute RMS amplitude --- # read time waveform # force to float64 precision # compute RMS amplitude # display results on console # --- store results in data (i.e., a dictionary)--- # slope of stimulus f0s in Hz/second # slope of FFR f0s in Hz/second # coefficient of determination # --- Output data --- # stimulus f0s # response f0s and f0 amplitudes # regression data # save regression data (dictionary) as .json # convert dict to a json string # open a json file # write string to file # close file # # read regression data (dictionary) back from .json          # debug only # fileObject = open(filepathstr + '_f0_LinearRegress.json', 'r') # debug only # jsonContent = fileObject.read()                                        # debug only # data_back = json.loads(jsonContent)                                    # debug only # save regression data (dictionary) as .npy # # read regression data (dictionary) back from .npy      # debug only # data_back = np.load(filepathstr + '_f0_LinearRegress.npy', allow_pickle=True).item()  # debug only # --- output figures --- # create a Path object # Check and see if this Path object points to an existing directory. If not, make a new directory. # convert back to a regular string # pause for a second, allowing the user time to view the results. Most importantly, this pause() allows/forces the figures to be displayed, before being closed at the beginning of the next iteration.