import numpy as np

def multipleFilterAnalysis(data,alfa,cfreq,dt,ndec):
    """
    Perform a multiple filter analysis of data.
    data:  Array of detrended and demeaned data whose length is power of 2 
    alfa:  Width parameter of Gaussian bandpass filter
    cfreq: Array of center frequencies of Gaussian filter
    dt: sampling interval of data
    ndec: decimation factor for instantaneous amplitude output
    """
    npts = len(data)
    nd = int(pow(2, np.ceil(np.log(npts)/np.log(2))))   # find next higher power of 2 of npts
    ftd = np.fft.rfft(data,nd)                   # Fourier transform of entire data set (pos. freq.)
                                                 # data are padded with zeros since npts <= nd
    freq = np.fft.rfftfreq(nd,dt)                # Fourier frequencies (positive frequencies)
    nt = int(np.ceil(npts/ndec))
    mfa = np.zeros((len(cfreq),nt))              # numpy array for MFA result
    for jf,cf in enumerate(cfreq):
        hg = np.exp(-alfa*((freq-cf)/cf)**2)     # Gaussian filter (use f instead of omega here)
        fk = hg*ftd                              # multiply FT of data with filter
        qk = np.complex(0,1)*fk                  # FT of Hilbert transform
        ftk = np.fft.irfft(fk)                   # filtered data
        qtk = np.fft.irfft(qk)                   # Hilbert transform of filtered data
        at = np.sqrt(ftk**2+qtk**2)              # instantaneous amplitude
        mfa[jf,:] = at[0:npts:ndec]              # store decimated original result
    return mfa
#-----------------------------------------------------------------------
# normalize multiple filter result either along time or frequency axis
def normalizeMFT(mfa,mode,exp):
    """
    Normalize the result of the mutiple filtering operation.
    mfa: array with instantaneous amplitudes versus frequency (mfa(f,t))
    mode: normalization mode: if 'time', normalize along time axis
          else normalize along frequency axis
    exp:  exponent for modifying inst amp using a power less than 1
    """
    nf,nt = mfa.shape
    if mode == 'time':
        for jf in range(0,nf):
            mfamax = np.amax(mfa[jf,:])
            mfa[jf,:] = np.power(mfa[jf,:]/mfamax+1.e-10,exp)          
    else:
        for jt in range(0,nt):
            mfamax = np.amax(mfa[:,jt])
            mfa[:,jt] = np.power(mfa[:,jt]/mfamax+1.e-10,exp)
    return mfa