def matrix_smooth(P):
    '''
    Purpose : smooth a spectrogram
    
    Parameters
    ----------
    hd : numpy array (2D)
        - spectrogram data
        
    Returns
    -------
    S2 : numpy array (2D)
        - smoothed spectrogram data
        
    Diary
    ---------------------------
    2019-09-09 (v01) 
        - Fuh borrowed this script from Dr. Tzu-Hao (Harry) Lin, written in Matlab
    2021-01-01 (v02) 
        - Fuh provided annotations and saved this script as a fucntion under my shared "functions" subfolder, in Matlab
    2021-10-08 (v03)
        - Fuh rewrote this script in Python
        
    ToDo
    ----
    None
                
    '''

    import numpy as np
    
    # initialize output array (2D)
    S2 = np.full(P.shape, np.nan)
    
    # spectrogram smoothing
    smoothing_factor = np.array([[1.0, 2.0, 1.0], [2.0, 4.0, 2.0], [1.0, 2.0, 1.0]]) # smoothing matrix
    sum_sm = np.sum(smoothing_factor)       # sum of all smoothing factors
    
    # create a smoothing array for each time bin (i.e., column)
    (a, b) = P.shape
    P_tmp = np.full((a+2, b+2), 0.0)
    P_tmp[1:-1, 1:-1] = P_tmp[1:-1, 1:-1] + P
    smoothing_tmp = np.tile(smoothing_factor, (a+2, 1))
    
    # smooth input array, one column at a time
    for m in range(1, b+1):
        tmp1 = P_tmp[:,m-1:m+2]        
        tmp2 = np.concatenate((P_tmp[1:, m-1:m+2], np.full((1,3), 0.0)), axis=0)        
        tmp3 = np.concatenate((P_tmp[2:, m-1:m+2], np.full((2,3), 0.0)), axis=0)
        P_tmp_m = np.concatenate((tmp1, tmp2, tmp3), axis=1).T.reshape((3,-1), order='F').T
        
        tmp4 = np.sum(P_tmp_m * smoothing_tmp, axis=1)  
        P_tmp_m = np.sum(tmp4.reshape((3,-1), order='F').T, axis=1) / sum_sm        
        S2[:,m-1] = P_tmp_m[:-2]
        
    return S2
