Cepstral Analysis on NMF Component

Question

I am using Non-negative Matrix Factorization (NMF) on an audio file of a working motor to decompose it into five NMF components. I need to explore the quefrency/cepstral domain of each component, but I am unsure how to proceed.

Since the NMF components are computed after applying the Short-Time Fourier Transform (STFT) to the original signal and are represented as a 2D matrix, I am uncertain about how to transform them into the quefrency domain.

I attempted to generate a cepstrogram for one of my NMF components using the following code:

# STFT Parameters
n_fft = 2048
hop_length = 256
win_length = n_fft
window = 'hann'

# Compute STFT
S = librosa.stft(y, n_fft=n_fft, hop_length=hop_length, win_length=win_length, window=window)
X, X_phase = librosa.magphase(S)  


# Normalization
X_norm = librosa.util.normalize(X)


# NMF Calculation
n_components = 5
model = NMF(n_components=n_components, init='random', beta_loss='kullback-leibler', solver='mu', max_iter=3000, tol=1e-5, l1_ratio=0.8)
W = model.fit_transform(X_norm)
H = model.components_


# Compute the cepstrogram from log magnitude spectrum

c = np.outer(W[:, 4], H[4]) #c is already a absolute value
# Compute log magnitude
log_magnitude = np.log(c + 1e-10)  # Add small value to avoid log(0)

# Compute cepstrum using IFFT
cepstrum = np.fft.ifft(log_magnitude, axis=0).real  # Compute real cepstrum across time frames

# Number of unique cepstral coefficients
NUP = cepstrum.shape[0] // 2  # Keep only the first half

# Extract the meaningful part of the cepstrum
C = cepstrum[:NUP, :]

# Compute quefrency and time vectors
q = np.arange(NUP) / sr  # Quefrency (s)
t = np.arange(C.shape[1]) * (hop_length / sr)  # Time (s)


# ===== APPLY CONDITIONING =====
q_mask = q > 0.0005 # Ignore quefrencies below 0.5 ms
C = C[q_mask, :]  # Filter cepstrum coefficients
q = q[q_mask] * 1000  # Convert quefrency from seconds to milliseconds


plt.figure(figsize=(10, 6))
plt.imshow(C, aspect='auto', origin='lower', extent=[t[0], t[-1], q[0], q[-1]])
plt.xlabel('Time (s)')
plt.ylabel('Quefrency (ms)')
plt.title('Cepstrogram from STFT (Magnitude Spectrum)')
plt.colorbar(label='Amplitude')
plt.show()

I need to determine whether this approach is correct. If it is, how can I extract quefrency domain features such as Linear Frequency Cepstral Coefficients (LFCC), Constant Q Cepstral Coefficients (CQCC), Teager Energy Cepstral Coefficients (TECC), and Gabor Filterbank Cepstral Coefficients (GTFB)?