Отношение сигнал/шум SRN низкое (-1, -2) в аудиофайлах, использующих пользовательский stft/istft и полный класс PLSI, SR ⇐ Python

[Anonymous]

** Исследование решения PLSI по использованию STFT и ISTFT для измерения соотношения сигнал-шум в аудиофайлах (3) для обеспечения точности **
** Постановка проблемы **
Используйте специально закодированное вероятностное латентно-семантическое индексирование (PLSI) для шумоподавления речи. Преобразуйте два обучающих сигнала trs.wav и trn.wav с помощью пользовательской настройки STFT. Создайте систему шумоподавления речи, используя алгоритм PLSI
(не используя NMF). Восстановите речевой сигнал во временной области, применив обратное STFT из tes.wav
** Неправильные результаты для PLSI и плохое отношение сигнал-шум (-2) **
Проблема этого PLSI заключается в том, что сигнал/шум в аудиофайлах равен -2 (низкая точность сигнала восстановления), и необходимо найти решение, позволяющее восстановить этот сигнал/шум в аудиофайлах в положительном (от +2 до +20).
** Были созданы PLSI, STFT и ISTST **
Закодированы были пользовательские PLSI, STFT, ISTAF, функции hann_window, update_W_rule, update_H_rule, Spectral_subtraction, Recovery_signal_with_spectral_subtraction и другие, которые включены в решение кода. Также тестировалось с librosa.stft, librosa.istft.
** Ожидаемый результат заключался в восстановлении сигнала и положительном соотношении сигнал/шум (SRN) (от +1 до + 25), но результаты по этому коду привели только к отрицательному соотношению сигнал-шум (SRN) (-2, -4, что является плохим восстановлением сигнала. Этот вопрос заключается в том, чтобы выяснить, почему соотношение сигнал-шум (SRN) отрицательное (-2). , -4) и как я могу восстановить сигнал и получить положительное отношение сигнал/шум (SRN) (от +1 до +25).
** Вот код для поиска любые обновления о логических ошибках **
** Базовый рабочий процесс для аудиофайлов TRN и TRS **
trs = base_path + '/trs.wav'
trn = base_path + '/trn.wav'
tex = base_path + '/tex.wav'
tes = base_path + '/tes.wav'

audio_trs, s_rate_trs = librosa.load(trs, sr=None)
audio_trn, s_rate_trn = librosa.load(trn, sr=None)
audio_tes, s_rate_tes = librosa.load(tes, sr=None)
audio_tex, s_rate_tex = librosa.load(tex, sr=None)

N = 1024
hop_length = 128
num_components = 30
num_iterations = 20
num_topics = 5
num_features=1000
n_fft=2048

#audio_trs
stft_s_hat = librosa.stft(audio_trs, n_fft=2048, hop_length=512)
magnitude_trs = np.abs(stft_s_hat)
phase_trs = np.exp(1j * np.angle(stft_s_hat))

plsi_model = PLSI(num_topics=30, num_iterations=50) # original low SRN had num_topics=30, iterations=500
W, H = plsi_model.fit_transform(magnitude_trs) # original low SRN had np.abs(stft_s_hat)
noise_indices = [0, 1] # This needs to be determined somehow
H_cleaned = plsi_model.remove_noise_topics(H, noise_indices)
cleaned_spectrum = plsi_model.reconstruct_spectrogram(W, H_cleaned)

cleaned_stft_full = cleaned_spectrum * np.exp(1j * np.angle(stft_s_hat))
recovered_signal = librosa.istft(cleaned_stft_full, hop_length=512)

snr = compute_signal_to_noise_rate(audio_tes, recovered_signal)
print(f" signal to noise rate (SNR): {snr}")
result: signal to noise rate (SNR): -2.3489089966432846

#audio_trn
stft_s_hat = librosa.stft(audio_trn, n_fft=2048, hop_length=512)
magnitude_trn = np.abs(stft_s_hat)
phase_trs = np.exp(1j * np.angle(stft_s_hat))

plsi_model = PLSI(num_topics=30, num_iterations=50) # original low SRN had num_topics=30, iterations=500
W, H = plsi_model.fit_transform(magnitude_trn) # original low SRN had np.abs(stft_s_hat)noise_indices = [0, 1] # H_cleaned = plsi_model.remove_noise_topics(H, noise_indices)
cleaned_spectrum = plsi_model.reconstruct_spectrogram(W, H_cleaned)spectrum)
cleaned_stft_full = cleaned_spectrum * np.exp(1j * np.angle(stft_s_hat))
recovered_signal = librosa.istft(cleaned_stft_full, hop_length=512)

snr = compute_signal_to_noise_rate(audio_tes, recovered_signal)
print(f" signal to noise rate (SNR): {snr}")
signal to noise rate (SNR): -2.0276653852458826

** Окончательное восстановление сигнала осуществляется с помощью этого кода **
display(Audio(audio_trn, rate = s_rate_trn))
display(Audio(audio_trs, rate = s_rate_trs))
display(Audio(recovered_signal, rate = s_rate_trn))

** Базовый класс PLSI и вспомогательные функции **
class PLSI:
def __init__(self, num_topics, num_iterations, lambda_reg=0.01, psi_reg=0.01):
self.num_topics = num_topics
self.num_iterations = num_iterations
self.W = None
self.H = None
self.topic_prob = None
self.doc_topic_prob = None
self.topic_word_prob = None
self.lambda_reg = lambda_reg
self.psi_reg = psi_reg

def fit(self, X):
if X.size == 0 or X.shape[0] == 0 or X.shape[1] == 0:
raise ValueError("Input data X is empty or has invalid dimensions.")

# Initialize W and H matrices with non-negative values
self.W = np.random.rand(X.shape[0], self.num_topics)
self.H = np.random.rand(self.num_topics, X.shape[1])

# Normalize W and H safely
max_W = np.max(self.W, axis=0, keepdims=True)
max_H = np.max(self.H, axis=0, keepdims=True)
self.W = np.divide(self.W, max_W, where=(max_W != 0))
self.H = np.divide(self.H, max_H, where=(max_H != 0))

# Iteratively update W and H using the update rules
for iteration in range(self.num_iterations):
self.W = update_W_rule(X, self.W, self.H)
self.H = update_H_rule(X, self.W, self.H)
return self.W, self.H

def fit_emm(self, data, lambda_reg=0.01, psi_reg=0.01):
num_documents, num_words = data.shape

self.doc_topic_prob = np.random.rand(num_documents, self.num_topics)
self.doc_topic_prob /= self.doc_topic_prob.sum(axis=1, keepdims=True)

self.topic_word_prob = np.random.rand(self.num_topics, num_words)
self.topic_word_prob /= self.topic_word_prob.sum(axis=1, keepdims=True)

self.topic_prob = np.random.rand(num_documents, num_words, self.num_topics)
self.topic_prob /= self.topic_prob.sum(axis=2, keepdims=True)

self.W = np.random.rand(num_documents, self.num_topics)
self.H = np.random.rand(self.num_topics, num_words)

self.W /= np.max(self.W, axis=0)
self.H /= np.max(self.H, axis=0)

# Expectation-Maximization algorithm
for iteration in range(self.num_iterations):
#print(f'Iteration {iteration + 1}/{self.num_iterations}')

# Expectation step
for e in range(num_documents):
for n in range(num_words):
# Calculate the probability of word 'v' at position 'n' in document 'd'
prob_word_topic = np.dot(self.doc_topic_prob[e, :], self.topic_word_prob[:, n])
prob_topic_word = self.topic_prob[e, n, :] * prob_word_topic
prob_topic_word /= np.sum(prob_topic_word)
self.topic_prob[e, n, :] = prob_topic_word

# M-step with regularization
for k in range(self.num_topics):
for w in range(num_words):
self.topic_word_prob[k, w] = np.sum(data[:, w] * self.topic_prob[:, w, k])
# Regularization for topic-word distribution
self.topic_word_prob[k, w] -= psi_reg * (1 - np.sum(self.topic_word_prob[k, :]))
self.topic_word_prob[k, :] /= np.sum(self.topic_word_prob[k, :])

for d in range(num_documents):
for k in range(self.num_topics):
self.doc_topic_prob[d, k] = np.sum(data[d, :] * self.topic_prob[d, :, k])
# Regularization for document-topic distribution
self.doc_topic_prob[d, k] -= lambda_reg * (1 - np.sum(self.doc_topic_prob[d, :]))
self.doc_topic_prob[d, :] /= np.sum(self.doc_topic_prob[d, :])

self.W = self.doc_topic_prob
self.H = self.topic_word_prob.T # Transposed to align dimensions correctly

def transform(self, X):
# Assuming W is already trained and we want to get H for new data
new_H = np.random.rand(self.num_topics, X.shape[1])
new_H /= np.max(new_H, axis=0)
for iteration in range(self.num_iterations):
new_H = update_H_rule(X, self.W, new_H)
return new_H

def fit_transform(self, X):
self.fit(X)
return self.W, self.H

def transform_to_matrix(texts, vocabulary):
# Initialize the document-term matrix with zeros
matrix = np.zeros((len(texts), len(vocabulary)))
for row, text in enumerate(texts):
tokens = tokenize(text)
for token in tokens:
if token in vocabulary:
col = vocabulary[token]
matrix[row, col] += 1
return matrix

def remove_noise_topics(self, H, noise_indices):
""" Reduce or remove noise-related topics before reconstructing the signal. """
H_cleaned = np.copy(H)
H_cleaned[noise_indices, :] = 0 # Set noise topic rows to zero
return H_cleaned

def reconstruct_spectrogram(self, W, H):
#Reconstruct the spectrogram from W and H matrices, excluding noise
return np.dot(W, H)

def tokenize(text):
return text.lower().split()

def build_vocabulary(texts):
vocabulary = {}
for text in texts:
tokens = tokenize(text)
for token in tokens:
if token not in vocabulary:
vocabulary[token] = len(vocabulary)
return vocabulary

def transform_to_matrix(texts, vocabulary):
matrix = np.zeros((len(texts), len(vocabulary)))
for row, text in enumerate(texts):
tokens = tokenize(text)
for token in tokens:
if token in vocabulary:
col = vocabulary[token]
matrix[row, col] += 1
return matrix

def compute_count_vectorizer(texts):
vocabulary = build_vocabulary(texts)
matrix = transform_to_matrix(texts, vocabulary)
return matrix, sorted(vocabulary, key=vocabulary.get)

def check_non_negativity(X):
return np.all(X >= 0)

def hann_window(N):
return 0.5 - 0.5 * np.cos(2 * np.pi * np.arange(N) / N)

def update_W_rule(X, W, H):
product = np.dot(W, H)
product[product == 0] = 1e-10
adjustment = X / product
numerator = np.dot(adjustment, H.T)
denominator = np.dot(np.ones((X.shape[0], 1)), np.sum(H, axis=1).reshape(1, -1))
denominator[denominator == 0] = 1e-10
W_new = W * (numerator / denominator)
return W_new

def update_H_rule(X, W, H):
product = np.dot(W, H)
product[product == 0] = 1e-10
adjustment = X / product
numerator = np.dot(W.T, adjustment)
denominator = np.dot(W.T, np.ones(X.shape))
denominator[denominator == 0] = 1e-10
H_new = H * (numerator / denominator)
return H_new

Сигнальные функции
def signal_reverse(X, N):
X_out = X[:(X.shape[0] - int(N/2)),0]
for col in range(1,X.shape[1]):
X_col = X[(X.shape[0] - int(N/2)):,col-1] + X[0:int(N/2),col]
X_out = np.concatenate((X_out,X_col))

X_out = X_out.reshape(-1)

return X_out

def signal_inverse(X, N):
return np.array([[np.exp(1j*(2*np.pi*f*n/N)) for f in range(0,N)] for n in range(0,N)])/N

** Функция SRN **
def compute_signal_to_noise_rate(signal, estimate):
min_length = min(len(signal), len(estimate))
signal = signal[:min_length]
estimate = estimate[:min_length]
signal_power = np.sum(signal**2)
noise_power = np.sum((signal - estimate)**2)
if noise_power == 0:
return float('inf') # Prevent division by zero
snr = 10 * np.log10(signal_power / noise_power)
return snr

Мои собственные STFT и ISTFT, но тестируются с помощью librose
def istft(stft_matrix, n_fft=2048, hop_length=512):
num_freqs, num_frames = stft_matrix.shape
expected_signal_len = n_fft + hop_length * (num_frames - 1)
signal = np.zeros(expected_signal_len)
window = np.hanning(n_fft) # Window function to reduce artifacts

for i in range(num_frames):
start = i * hop_length
end = start + n_fft
ifft_segment = np.fft.irfft(stft_matrix[:, i], n=n_fft)
windowed_ifft = ifft_segment * window # Apply window function
signal[start:end] += windowed_ifft
return signal

def stft(signal, n_fft=2048, hop_length=512):
signal = np.asarray(signal, dtype=float)
S = len(signal)
num_frames = max(1, 1 + (S - n_fft) // hop_length)
if S < n_fft:
signal = np.pad(signal, (0, n_fft - S), mode='constant')

stft_matrix = np.zeros((n_fft // 2 + 1, num_frames), dtype=np.complex64)
window = np.hanning(n_fft)

for i in range(num_frames):
start = i * hop_length
end = start + n_fft
signal_segment = signal[start:end] if end

Подробнее здесь: https://stackoverflow.com/questions/783 ... stft-istft