Google Speech Commands Encoder¶
In [1]:
%reset -f
import time
start = time.time()
from IPython.display import display, HTML
display(HTML('''
<script id="MathJax-script"
async
src="https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-svg.js">
</script>
'''))
import os
import warnings
import numpy as np
import torchaudio
from scipy import signal
from scipy.signal import find_peaks
from torch import Tensor
from torchaudio.datasets import SPEECHCOMMANDS
import torchaudio.transforms as T
from plotly.subplots import make_subplots
import plotly.graph_objects as go
import torch
from pathlib import Path
from noisereduce.torchgate import TorchGate as TG
import plotly
from torchaudio.transforms import Resample
import librosa
from dyned import (
sigma_delta_quantisation,
generate_step_signal,
identify_peaks_and_troughs,
create_spike_plot,
create_rle_spike_plot,
create_lzw_spike_plot,
create_huffman_spike_plot,
create_delta_spike_plot,
create_bwt_spike_plot,
create_lz77_spike_plot,
create_dynedc_v4_spike_plot,
compression_summary,
print_compression_table,
)
plotly.offline.init_notebook_mode(connected=True)
import plotly.io as pio
if os.environ.get("theme", "light") == "dark":
pio.templates.default = "plotly_dark"
else:
pio.templates.default = "plotly"
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
tg = TG(sr=8000, nonstationary=False).to(device)
command_index = 3000
In [2]:
def transform_waveform(waveform: torch.Tensor, waveform_sample_rate: int):
n_fft = 1024 # Size of the FFT typicalvaluesare1024,2048
hop_length = 512 # Hop length typicallynfft//2ornfft//4
n_mels = 64 # Number of mel bands
mel_spectrogram_transform = T.MelSpectrogram(
sample_rate=waveform_sample_rate,
n_fft=n_fft,
hop_length=hop_length,
n_mels=n_mels,
f_min=0,
f_max=waveform_sample_rate // 2,
pad=0,
center=True,
power=2.0
).to(device)
# Apply the transformation
mel_spectrogram = mel_spectrogram_transform(waveform)
# Convert to decibels
mel_spectrogram_db = T.AmplitudeToDB()(mel_spectrogram)
return mel_spectrogram_db
# Load the SPEECHCOMMANDS dataset
class CustomSPEECHCOMMANDS(SPEECHCOMMANDS):
def __init__(self, subset, root="./assets", download=True, transform=None):
super().__init__(root=root, download=download, subset=subset)
self.transform = transform
self.device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
def __getitem__(self, n) -> tuple[Tensor, int, str, str, int, str]:
_waveform, _sample_rate, _label, _speaker_id, _utterance_number = super().__getitem__(n)
_file_path = self._walker[n]
_file_name = Path(_file_path).name
if self.transform is not None:
_waveform = self.transform(_waveform, _sample_rate)
_waveform = _waveform.to(self.device)
return _waveform, _sample_rate, _label, _speaker_id, _utterance_number, _file_name
train_dataset_raw = CustomSPEECHCOMMANDS(subset="training")
test_dataset_raw = CustomSPEECHCOMMANDS(subset="testing")
Cleaning Raw Audio Data¶
In [3]:
def downsample_audio_torch(audio, original_sr, target_sr):
# Ensure audio is a torch tensor
if not isinstance(audio, torch.Tensor):
audio = torch.tensor(audio)
# Get the device of the input audio
device = audio.device
# Move audio to CPU for resampling
audio = audio.cpu()
# Ensure audio is 2D channels,samples
if audio.dim() == 1:
audio = audio.unsqueeze(0)
# Create resampler
resampler = Resample(orig_freq=original_sr, new_freq=target_sr)
# Resample
audio_resampled = resampler(audio)
# Move back to original device
audio_resampled = audio_resampled.to(device)
return audio_resampled
sample_raw = train_dataset_raw[command_index]
audio_data_raw, sample_rate_raw, label_raw, speaker_id_raw, utterance_number_raw, file_name = sample_raw
cleaned_audio_data_raw = tg(audio_data_raw)
print(f"Sample Rate: {sample_rate_raw}")
# Downsample the audio to 8 kHz
target_sr = 8000
audio_data_raw = downsample_audio_torch(audio_data_raw, sample_rate_raw, target_sr)
print(f"Raw Audio Data Shape: {audio_data_raw.shape}")
cleaned_audio_data_raw = tg(audio_data_raw)
print(f"New Sample Rate: {sample_rate_raw}")
print(f"Raw Audio Data Shape: {audio_data_raw.shape}")
print(f"Label: {label_raw}")
print(f"Speaker ID: {speaker_id_raw}")
print(f"Utterance Number: {utterance_number_raw}")
print(f"File Name: {file_name}")
def save_audio(file_path, audio_data, sample_rate):
try:
# Ensure audio_data is a PyTorch tensor
if isinstance(audio_data, np.ndarray):
audio_data = torch.from_numpy(audio_data)
# Ensure audio is 2D (channels, samples)
if audio_data.dim() == 1:
audio_data = audio_data.unsqueeze(0)
# Normalize audio to [-1, 1] range if it's not already
if audio_data.abs().max() > 1:
audio_data = audio_data / audio_data.abs().max()
# Convert to 16-bit PCM
audio_data = (audio_data * 32767).to(torch.int16).cpu()
# Save the audio file
torchaudio.save(file_path, audio_data, sample_rate)
print(f"Audio saved successfully to {file_path}")
except Exception as e:
print(f"Error saving audio file: {str(e)}")
def create_plotly_spectrogram(audio_data, sample_rate, nperseg=256, noverlap=128):
# Compute spectrogram data
f, t, Sxx = signal.spectrogram(audio_data, fs=sample_rate, nperseg=nperseg, noverlap=noverlap)
# Convert to dB
Sxx_db = 10 * np.log10(Sxx) # Adding small value to avoid log(0)
return f, t, Sxx_db
# Plot the raw audio waveform
fig = make_subplots(rows=3, cols=1, vertical_spacing=0.15, subplot_titles=("Raw Audio Waveform", "Mel Spectrogram", "Spectrogram"))
fig.add_trace(go.Scatter(x=np.arange(audio_data_raw.shape[1]), y=audio_data_raw.squeeze().cpu().numpy(), mode="lines"), row=1, col=1)
raw_mel = transform_waveform(audio_data_raw, sample_rate_raw).squeeze().cpu().numpy()
fig.add_trace(go.Heatmap(z=raw_mel,
colorscale='Viridis',
zmin=raw_mel.min(),
zmax=raw_mel.max()), row=2, col=1)
f,t,ssx = create_plotly_spectrogram(audio_data_raw.squeeze().cpu().numpy(), sample_rate_raw)
fig.add_trace(go.Heatmap(z=ssx,
x=t,
y=f,
colorscale='Viridis',
zmin=ssx.min(),
zmax=ssx.max()), row=3, col=1)
fig.update_layout(title=f"Raw Audio Waveform for {label_raw}", xaxis_title="Time", yaxis_title="Amplitude", showlegend=False, height=800, width=1400)
fig.show()
# Update sample rate
sample_rate_raw = target_sr
# Plot the cleaned audio waveform
fig = make_subplots(rows=3, cols=1, vertical_spacing=0.15, subplot_titles=("Cleaned Audio Waveform", "Mel Spectrogram", "Spectrogram"))
fig.add_trace(go.Scatter(x=np.arange(cleaned_audio_data_raw.shape[1]), y=cleaned_audio_data_raw.squeeze().cpu().numpy(), mode="lines"), row=1, col=1)
cleaned_mel = transform_waveform(cleaned_audio_data_raw, sample_rate_raw).squeeze().cpu().numpy()
fig.add_trace(go.Heatmap(z=cleaned_mel,
colorscale='Viridis',
zmin=cleaned_mel.min(),
zmax=cleaned_mel.max()), row=2, col=1)
f,t,ssx = create_plotly_spectrogram(cleaned_audio_data_raw.squeeze().cpu().numpy(), sample_rate_raw)
fig.add_trace(go.Heatmap(z=ssx,
x=t,
y=f,
colorscale='Viridis',
zmin=ssx.min(),
zmax=ssx.max()), row=3, col=1)
fig.update_layout(title=f"Cleaned Audio Waveform for {label_raw} using Spectral Gates", xaxis_title="Time", yaxis_title="Amplitude", showlegend=False, height=800, width=1400)
fig.show()
save_audio("raw_audio.wav", audio_data_raw, sample_rate_raw)
save_audio("cleaned_audio.wav", cleaned_audio_data_raw, sample_rate_raw)
# torchaudio.save("raw_audio.wav", audio_data_raw, sample_rate_raw)
# torchaudio.save("cleaned_audio.wav", cleaned_audio_data_raw, sample_rate_raw)
Sample Rate: 16000 Raw Audio Data Shape: torch.Size([1, 4779]) New Sample Rate: 16000 Raw Audio Data Shape: torch.Size([1, 4779]) Label: bird Speaker ID: 0bfec55f Utterance Number: 0 File Name: 0bfec55f_nohash_0.wav
/tmp/ipykernel_2067593/1997290036.py:74: RuntimeWarning: divide by zero encountered in log10 Sxx_db = 10 * np.log10(Sxx) # Adding small value to avoid log(0)
Audio saved successfully to raw_audio.wav Audio saved successfully to cleaned_audio.wav
Quantisation¶
In [4]:
# Compute FFT of cleaned audio
fft_result = torch.fft.fft(cleaned_audio_data_raw.squeeze())
fft_shifted = torch.fft.fftshift(fft_result)
magnitude_spectrum = torch.abs(fft_shifted)
phase_spectrum = torch.angle(fft_shifted)
magnitude_spectrum = torch.log1p(magnitude_spectrum)
magnitude_spectrum_np = magnitude_spectrum.cpu().numpy()
# Calculate frequency bins
N = len(magnitude_spectrum_np)
freqs = np.fft.fftshift(np.fft.fftfreq(N, d=1/sample_rate_raw))
print(f"Frequency Bins: {freqs.shape}")
# Quantise using DyNED encoder
quantised_magnitude_spectrum, errors = sigma_delta_quantisation(magnitude_spectrum, levels=10)
quantised_step_signal = generate_step_signal(quantised_magnitude_spectrum)
peaks, troughs = identify_peaks_and_troughs(quantised_magnitude_spectrum)
fig = make_subplots(rows=5, cols=1, vertical_spacing=0.1, subplot_titles=("Frequency Spectrum", "Quantised Magnitude Spectrum", "Spikes", "Peaks and Troughs", "Errors"))
fig.add_trace(go.Scatter(x=freqs, y=magnitude_spectrum.cpu(), mode="lines"), row=1, col=1)
fig.add_trace(go.Scatter(x=freqs, y=quantised_magnitude_spectrum, mode="lines"), row=2, col=1)
# Spikes
fig.add_trace(create_spike_plot(freqs, quantised_step_signal), row=3, col=1)
fig.add_trace(go.Scatter(x=freqs[peaks], y=quantised_magnitude_spectrum[peaks], mode='markers', marker=dict(color='green', size=3), name='Peaks'), row=4, col=1)
fig.add_trace(go.Scatter(x=freqs[troughs], y=quantised_magnitude_spectrum[troughs], mode='markers', marker=dict(color='red', size=3), name='Troughs'), row=4, col=1)
fig.add_trace(go.Scatter(x=freqs, y=errors, mode="lines"), row=5, col=1)
fig.update_layout(title="Quantisation of Frequency Spectrum", xaxis_title="Frequency (Hz)", yaxis_title="Magnitude", showlegend=False, height=1080, width=1400)
fig.show()
Frequency Bins: (4608,)
Number of spikes: 2341
Reconstruction¶
In [5]:
from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score
from scipy.interpolate import RBFInterpolator, CubicSpline, KroghInterpolator
def stable_piecewise_interpolation(x, y, segment_size=30, overlap=5):
"""
Perform stable piecewise interpolation with smooth blending.
"""
if len(x) != len(y):
raise ValueError("x and y must have the same length")
result = np.zeros_like(x, dtype=float)
weights = np.zeros_like(x, dtype=float)
for start in range(0, len(x), segment_size - overlap):
end = min(start + segment_size, len(x))
# Create interpolator for this segment
interp = KroghInterpolator(x[start:end], y[start:end])
# Evaluate interpolator on this segment
segment_result = interp(x[start:end])
# Create a weight array for smooth blending
segment_weights = np.ones(end - start)
if start > 0:
segment_weights[:overlap] = np.linspace(0, 1, overlap)
if end < len(x):
segment_weights[-overlap:] = np.linspace(1, 0, overlap)
# Add to result and weights
result[start:end] += segment_result * segment_weights
weights[start:end] += segment_weights
# Normalize by weights
result /= weights
return result
# Assuming quantised_magnitude_spectrum is already defined
quantised_magnitude_spectrum_reshaped = quantised_magnitude_spectrum.reshape(-1, 1)
freqs_reshaped = freqs.reshape(-1, 1) # Ensure this is a 2D array as required by RBFInterpolator
# RBF interpolation
start_time = time.time()
rbfi = RBFInterpolator(freqs_reshaped, quantised_magnitude_spectrum_reshaped, kernel="linear")
freqs_interpolated = np.linspace(freqs.min(), freqs.max(), len(quantised_magnitude_spectrum)) # Smooth interpolation across the frequency range
interpolated_values = rbfi(freqs_interpolated[:, None]).flatten() # Flatten because the output will be 2D otherwise
rbf_time = time.time() - start_time
# Cubic spline interpolation
start_time = time.time()
cube_spline = CubicSpline(freqs, quantised_magnitude_spectrum, bc_type='clamped')
interpolated_values_cubic = cube_spline(freqs_interpolated)
cubic_time = time.time() - start_time
# Radial basis function interpolation
start_time = time.time()
rbf = RBFInterpolator(freqs_reshaped, quantised_magnitude_spectrum_reshaped, kernel='cubic')
interpolated_values_rbf = rbf(freqs_interpolated[:, None]).flatten()
rbf_classic_time = time.time() - start_time
# Krogh interpolation
start_time = time.time()
interpolated_values_krogh = stable_piecewise_interpolation(freqs, quantised_magnitude_spectrum)
krogh_time = time.time() - start_time
# Plot using Plotly
fig = make_subplots(rows=4, cols=1, subplot_titles=("RBF Interpolation", "Cubic Interpolation", "RBF Interpolation (Classic)", "Krogh Interpolation"), vertical_spacing=0.09)
fig.add_trace(go.Scatter(x=freqs, y=magnitude_spectrum.cpu(), mode='lines', name='Quantized Magnitude Spectrum'), row=1, col=1)
fig.add_trace(go.Scatter(x=freqs_interpolated, y=interpolated_values, mode='lines',line = dict(dash='dot'), name='Interpolated Quantized Magnitude Spectrum'), row=1, col=1)
fig.add_trace(go.Scatter(x=freqs, y=magnitude_spectrum.cpu(), mode='lines', name='Quantized Magnitude Spectrum'), row=2, col=1)
fig.add_trace(go.Scatter(x=freqs_interpolated, y=interpolated_values_cubic, mode='lines', line = dict(dash='dot'), name='Interpolated Quantized Magnitude Spectrum'), row=2, col=1)
fig.add_trace(go.Scatter(x=freqs, y=magnitude_spectrum.cpu(), mode='lines', name='Quantized Magnitude Spectrum'), row=3, col=1)
fig.add_trace(go.Scatter(x=freqs_interpolated, y=interpolated_values_rbf, mode='lines', line = dict(dash='dot'), name='Interpolated Quantized Magnitude Spectrum'), row=3, col=1)
fig.add_trace(go.Scatter(x=freqs, y=magnitude_spectrum.cpu(), mode='lines', name='Quantized Magnitude Spectrum'), row=4, col=1)
fig.add_trace(go.Scatter(x=freqs_interpolated, y=interpolated_values_krogh, mode='lines', line = dict(dash='dot'), name='Interpolated Quantized Magnitude Spectrum'), row=4, col=1)
fig.update_layout(title='Interpolated Quantized FFT Magnitude Spectrum', xaxis_title='Frequency (Hz)', yaxis_title='Magnitude (dB)', height=1080, width=1400, showlegend=False)
fig.show()
original_values_at_interpolated_freqs = np.interp(freqs_interpolated, freqs, magnitude_spectrum.cpu())
# Compute interpolation method
# RBF Interpolation Error Calculation
mse_rbfi = mean_squared_error(original_values_at_interpolated_freqs, interpolated_values)
rmse_rbfi = np.sqrt(mse_rbfi)
mae_rbfi = mean_absolute_error(original_values_at_interpolated_freqs, interpolated_values)
r2_rbfi = r2_score(original_values_at_interpolated_freqs, interpolated_values)
# Cubic Spline Interpolation Error Calculation
mse_cubic = mean_squared_error(original_values_at_interpolated_freqs, interpolated_values_cubic)
rmse_cubic = np.sqrt(mse_cubic)
mae_cubic = mean_absolute_error(original_values_at_interpolated_freqs, interpolated_values_cubic)
r2_cubic = r2_score(original_values_at_interpolated_freqs, interpolated_values_cubic)
# Radial Basis Function Interpolation Error Calculation
mse_rbf = mean_squared_error(original_values_at_interpolated_freqs, interpolated_values_rbf)
rmse_rbf = np.sqrt(mse_rbf)
mae_rbf = mean_absolute_error(original_values_at_interpolated_freqs, interpolated_values_rbf)
r2_rbf = r2_score(original_values_at_interpolated_freqs, interpolated_values_rbf)
# Krogh Interpolation Error Calculation
mse_krogh = mean_squared_error(original_values_at_interpolated_freqs, interpolated_values_krogh)
rmse_krogh = np.sqrt(mse_krogh)
mae_krogh = mean_absolute_error(original_values_at_interpolated_freqs, interpolated_values_krogh)
r2_krogh = r2_score(original_values_at_interpolated_freqs, interpolated_values_krogh)
# Print Error Metrics
print(f"RBFInterpolator - RMSE: {rmse_rbfi}, MAE: {mae_rbfi}, R2: {r2_rbfi}, Time: {rbf_time:.4f}s")
print(f"Cubic Spline - RMSE: {rmse_cubic}, MAE: {mae_cubic}, R2: {r2_cubic}, Time: {cubic_time:.4f}s")
print(f"Radial Basis Function - RMSE: {rmse_rbf}, MAE: {mae_rbf}, R2: {r2_rbf}, Time: {rbf_classic_time:.4f}s")
print(f"Krogh Interpolator - RMSE: {rmse_krogh}, MAE: {mae_krogh}, R2: {r2_krogh}, Time: {krogh_time:.4f}s")
fig = make_subplots(rows=1, cols=1)
fig.add_trace(go.Scatter(x=freqs, y=magnitude_spectrum.cpu(), mode='lines', name='Original Data Points'), row=1, col=1)
fig.add_trace(go.Scatter(x=freqs_interpolated, y=interpolated_values, mode='lines', name='RBF Interpolation', line=dict(dash='dot')), row=1, col=1)
fig.add_trace(go.Scatter(x=freqs_interpolated, y=interpolated_values_cubic, name='Cubic Spline', line=dict(dash='dot')), row=1, col=1)
fig.add_trace(go.Scatter(x=freqs_interpolated, y=interpolated_values_rbf, name='RBF', line=dict(dash='dot')), row=1, col=1)
fig.add_trace(go.Scatter(x=freqs_interpolated, y=interpolated_values_krogh, name='Krogh', line=dict(dash='dot')), row=1, col=1)
fig.update_layout(title='Comparison of Interpolation Methods', xaxis_title='Frequency (Hz)', yaxis_title='Magnitude (dB)', height=500, width=1400)
fig.show()
RBFInterpolator - RMSE: 0.18424177718154452, MAE: 0.15118937379884326, R2: 0.9531799743875803, Time: 1.2924s Cubic Spline - RMSE: 0.1842417771814062, MAE: 0.1511893737987078, R2: 0.9531799743876506, Time: 0.0010s Radial Basis Function - RMSE: 0.18424247491317589, MAE: 0.15119319550377452, R2: 0.9531796197679581, Time: 1.3157s Krogh Interpolator - RMSE: 0.18424177025479124, MAE: 0.1511893773920111, R2: 0.9531799779080713, Time: 0.1312s
In [6]:
def reconstruct_speech(interpolated_magnitude, original_phase, scale_factor=1.0):
# Scale the magnitude
scaled_magnitude = interpolated_magnitude * scale_factor
# Combine scaled magnitude with original phase
complex_spectrum = scaled_magnitude * torch.exp(1j * original_phase)
# Perform inverse FFT shift and then inverse FFT
ifft_result = torch.fft.ifft(torch.fft.ifftshift(complex_spectrum))
# Take the real part
reconstructed_audio = torch.real(ifft_result)
return reconstructed_audio
def normalize_amplitude(original, reconstructed):
"""Normalize the amplitude of the reconstructed signal to match the original"""
original_rms = torch.sqrt(torch.mean(original**2))
reconstructed_rms = torch.sqrt(torch.mean(reconstructed**2))
scale_factor = original_rms / reconstructed_rms
return reconstructed * scale_factor
def compare_reconstruction(original, reconstructed, sr, normalize=True):
# Ensure both are on CPU and convert to numpy
original = original.cpu().numpy().squeeze()
reconstructed = reconstructed.cpu().numpy().squeeze()
# Normalize amplitude if requested
if normalize:
reconstructed = normalize_amplitude(torch.from_numpy(original), torch.from_numpy(reconstructed)).numpy()
# Adjust lengths if necessary
min_length = min(len(original), len(reconstructed))
original = original[:min_length]
reconstructed = reconstructed[:min_length]
# Calculate metrics
snr = 10 * np.log10(np.sum(original**2) / np.sum((original - reconstructed)**2))
rmse = np.sqrt(np.mean((original - reconstructed)**2))
# Create time array
time = np.arange(min_length) / sr
# Create plot
fig = make_subplots(rows=2, cols=1, subplot_titles=("Waveforms", "Spectrogram Comparison"))
# Plot waveforms
fig.add_trace(go.Scatter(x=time, y=original, name="Original", line=dict(color="blue")), row=1, col=1)
fig.add_trace(go.Scatter(x=time, y=reconstructed, name="Reconstructed", line=dict(color="red")), row=1, col=1)
# Compute and plot spectrograms
D_original = librosa.stft(original)
D_reconstructed = librosa.stft(reconstructed)
S_db_original = librosa.amplitude_to_db(np.abs(D_original), ref=np.max)
S_db_reconstructed = librosa.amplitude_to_db(np.abs(D_reconstructed), ref=np.max)
fig.add_trace(go.Heatmap(z=S_db_original, colorscale="Viridis", name="Original Spectrogram"), row=2, col=1)
fig.add_trace(go.Heatmap(z=S_db_reconstructed, colorscale="Viridis", name="Reconstructed Spectrogram"), row=2, col=1)
# Update layout
fig.update_layout(height=800, width=1000, title_text=f"Audio Reconstruction Comparison (SNR: {snr:.2f} dB, RMSE: {rmse:.4f})")
fig.update_xaxes(title_text="Time (s)", row=1, col=1)
fig.update_yaxes(title_text="Amplitude", row=1, col=1)
fig.update_xaxes(title_text="Time", row=2, col=1)
fig.update_yaxes(title_text="Frequency", row=2, col=1)
fig.show()
print(f"SNR: {snr:.2f} dB")
print(f"RMSE: {rmse:.4f}")
# Reconstruct the speech signal
reconstructed_audio = reconstruct_speech(torch.tensor(interpolated_values_cubic)*10, phase_spectrum.cpu())
reconstructed_audio = tg(reconstructed_audio.unsqueeze(0).cuda())
compare_reconstruction(cleaned_audio_data_raw, reconstructed_audio, sample_rate_raw, normalize=False)
reconstructed_audio = reconstruct_speech(torch.tensor(interpolated_values_krogh), phase_spectrum.cpu())
reconstructed_audio = tg(reconstructed_audio.unsqueeze(0).cuda())
compare_reconstruction(cleaned_audio_data_raw, reconstructed_audio, sample_rate_raw, normalize=True)
SNR: 4.59 dB RMSE: 0.0915
SNR: 5.49 dB RMSE: 0.0825
Post-Processing¶
In [7]:
def ensure_numpy(audio):
if isinstance(audio, torch.Tensor):
return audio.cpu().numpy()
return np.asarray(audio)
def plot_waveforms(original, processed, sr, frame_length=1024, hop_length=512):
original = ensure_numpy(original).flatten()
processed = ensure_numpy(processed).flatten()
# Find the minimum length
min_length = min(len(original), len(processed))
original = original[:min_length]
processed = processed[:min_length]
# Create subplots
fig = make_subplots(rows=4, cols=1,
subplot_titles=('Full Waveforms', 'First Few Frames', 'Frame Energy', 'Difference (Original - Processed)'))
# Plot full waveforms
time = np.arange(min_length) / sr
fig.add_trace(go.Scatter(x=time, y=original, name='Original', line=dict(color='blue')), row=1, col=1)
fig.add_trace(go.Scatter(x=time, y=processed, name='Processed', line=dict(color='red')), row=1, col=1)
# Plot the first few frames
n_frames = 5
plot_length = min(n_frames * frame_length, min_length)
t = np.arange(plot_length) / sr
fig.add_trace(go.Scatter(x=t, y=original[:plot_length], name='Original', line=dict(color='blue')), row=2, col=1)
fig.add_trace(go.Scatter(x=t, y=processed[:plot_length], name='Processed', line=dict(color='red')), row=2, col=1)
# Plot frame energy
frames_orig = librosa.util.frame(original, frame_length=frame_length, hop_length=hop_length)
frames_proc = librosa.util.frame(processed, frame_length=frame_length, hop_length=hop_length)
energy_orig = np.sum(frames_orig**2, axis=0)
energy_proc = np.sum(frames_proc**2, axis=0)
frame_time = np.arange(len(energy_orig)) * hop_length / sr
fig.add_trace(go.Scatter(x=frame_time, y=energy_orig, name='Original', line=dict(color='blue')), row=3, col=1)
fig.add_trace(go.Scatter(x=frame_time, y=energy_proc, name='Processed', line=dict(color='red')), row=3, col=1)
# Plot difference
difference = original - processed
fig.add_trace(go.Scatter(x=time, y=difference, line=dict(color='green')), row=4, col=1)
# Update layout
fig.update_layout(height=1000, width=1000, title_text="Waveform Comparison")
fig.update_xaxes(title_text="Time (s)", row=1, col=1)
fig.update_xaxes(title_text="Time (s)", row=2, col=1)
fig.update_xaxes(title_text="Time (s)", row=3, col=1)
fig.update_xaxes(title_text="Time (s)", row=4, col=1)
fig.update_yaxes(title_text="Amplitude", row=1, col=1)
fig.update_yaxes(title_text="Amplitude", row=2, col=1)
fig.update_yaxes(title_text="Energy", row=3, col=1)
fig.update_yaxes(title_text="Amplitude", row=4, col=1)
fig.update_yaxes(type="log", row=3, col=1) # Log scale for energy plot
# Show the plot
fig.show()
# Print length information
print(f"Original audio length: {len(original)/sr:.4f} seconds")
print(f"Processed audio length: {len(processed)/sr:.4f} seconds")
if len(original) != len(processed):
print(f"Note: Audio lengths differ by {abs(len(original) - len(processed))/sr:.4f} seconds")
def analyze_audio(original, processed, sr, frame_length=1024, hop_length=512):
original = ensure_numpy(original).flatten()
processed = ensure_numpy(processed).flatten()
# Find the minimum length
min_length = min(len(original), len(processed))
original = original[:min_length]
processed = processed[:min_length]
# Basic stats
print(f"Original audio range: {original.min():.4f} to {original.max():.4f}")
print(f"Processed audio range: {processed.min():.4f} to {processed.max():.4f}")
# Check for silence
silence_threshold = 1e-5
silent_orig = np.sum(np.abs(original) < silence_threshold) / original.size
silent_proc = np.sum(np.abs(processed) < silence_threshold) / processed.size
print(f"Proportion of silence in original: {silent_orig:.2%}")
print(f"Proportion of silence in processed: {silent_proc:.2%}")
# Frame-by-frame analysis
frames_orig = librosa.util.frame(original, frame_length=frame_length, hop_length=hop_length)
frames_proc = librosa.util.frame(processed, frame_length=frame_length, hop_length=hop_length)
for i in range(min(5, frames_orig.shape[1])): # Analyze first 5 frames
energy_orig = np.sum(frames_orig[:, i]**2)
energy_proc = np.sum(frames_proc[:, i]**2)
print(f"Frame {i}: Original energy = {energy_orig:.6f}, Processed energy = {energy_proc:.6f}")
# Usage
plot_waveforms(cleaned_audio_data_raw, reconstructed_audio*10, sample_rate_raw)
analyze_audio(cleaned_audio_data_raw, reconstructed_audio*10, sample_rate_raw)
Original audio length: 0.5760 seconds Processed audio length: 0.5760 seconds Original audio range: -0.9950 to 0.9193 Processed audio range: -1.8276 to 1.4438 Proportion of silence in original: 71.14% Proportion of silence in processed: 55.36% Frame 0: Original energy = 0.000000, Processed energy = 0.005351 Frame 1: Original energy = 0.000000, Processed energy = 0.000000 Frame 2: Original energy = 0.000000, Processed energy = 0.000000 Frame 3: Original energy = 0.000000, Processed energy = 0.000000 Frame 4: Original energy = 0.000000, Processed energy = 0.000000
In [8]:
def ensure_tensor(audio, device='cpu'):
if isinstance(audio, np.ndarray):
return torch.from_numpy(audio).to(device)
return audio.to(device)
def normalize(audio):
"""Normalize audio to range [-1, 1]"""
if isinstance(audio, torch.Tensor):
return audio / torch.max(torch.abs(audio))
return audio / np.max(np.abs(audio))
def peak_match(original, processed):
"""Match the peak amplitude of processed audio to original"""
if isinstance(original, torch.Tensor):
return processed * (torch.max(torch.abs(original)) / torch.max(torch.abs(processed)))
return processed * (np.max(np.abs(original)) / np.max(np.abs(processed)))
def rms_match(original, processed):
"""Match the RMS energy of processed audio to original"""
if isinstance(original, torch.Tensor):
original_rms = torch.sqrt(torch.mean(original**2))
processed_rms = torch.sqrt(torch.mean(processed**2))
else:
original_rms = np.sqrt(np.mean(original**2))
processed_rms = np.sqrt(np.mean(processed**2))
return processed * (original_rms / processed_rms)
def adaptive_boost(original, processed, target_ratio=0.9):
"""Boost processed signal to match a target ratio of original's range"""
if isinstance(original, torch.Tensor):
orig_range = torch.max(original) - torch.min(original)
proc_range = torch.max(processed) - torch.min(processed)
else:
orig_range = np.max(original) - np.min(original)
proc_range = np.max(processed) - np.min(processed)
boost_factor = (orig_range * target_ratio) / proc_range
return processed * boost_factor
def calculate_rms(signal):
"""Calculate Root Mean Square"""
if isinstance(signal, torch.Tensor):
return torch.sqrt(torch.mean(signal**2))
return np.sqrt(np.mean(signal**2))
def compress_dynamic_range(audio, threshold=0.5, ratio=2):
"""Apply dynamic range compression"""
if isinstance(audio, torch.Tensor):
mask = audio.abs() > threshold
compressed = audio.clone()
compressed[mask] = threshold + (compressed[mask] - threshold) / ratio
return compressed
else:
mask = np.abs(audio) > threshold
compressed = audio.copy()
compressed[mask] = threshold + (compressed[mask] - threshold) / ratio
return compressed
def hard_clip(audio, threshold=0.8):
"""Apply hard clipping"""
if isinstance(audio, torch.Tensor):
return torch.clamp(audio, -threshold, threshold)
else:
return np.clip(audio, -threshold, threshold)
def soft_clip(audio, amount=2):
"""Apply soft clipping using tanh function"""
if isinstance(audio, torch.Tensor):
return torch.tanh(amount * audio) / torch.tanh(torch.tensor(amount))
else:
return np.tanh(amount * audio) / np.tanh(amount)
def apply_boost(original, processed, method='normalize'):
original = ensure_tensor(original).flatten()
processed = ensure_tensor(processed).flatten()
if method == 'normalize':
boosted = normalize(processed)
elif method == 'peak_match':
boosted = peak_match(original, processed)
elif method == 'rms_match':
boosted = rms_match(original, processed)
elif method == 'adaptive':
boosted = adaptive_boost(original, processed)
elif method == 'compress':
boosted = compress_dynamic_range(processed)
elif method == 'hard_clip':
boosted = hard_clip(processed)
elif method == 'soft_clip':
boosted = soft_clip(processed)
else:
raise ValueError("Invalid boosting method")
return boosted
def ensure_same_length(original, processed):
"""Adjust processed signal to match the length of the original"""
if len(original) > len(processed):
return original, np.pad(processed, (0, len(original) - len(processed)), 'constant')
elif len(original) < len(processed):
return original, processed[:len(original)]
return original, processed
def calculate_snr(original, processed):
"""Calculate Signal-to-Noise Ratio"""
original, processed = ensure_same_length(original, processed)
noise = original - processed
if isinstance(original, torch.Tensor):
return 10 * torch.log10(torch.sum(original**2) / torch.sum(noise**2))
return 10 * np.log10(np.sum(original**2) / np.sum(noise**2))
def calculate_rmse(original, processed):
"""Calculate Root Mean Square Error"""
original, processed = ensure_same_length(original, processed)
if isinstance(original, torch.Tensor):
return torch.sqrt(torch.mean((original - processed)**2))
return np.sqrt(np.mean((original - processed)**2))
def librosa_effects(audio, sr):
"""Apply various audio effects using librosa"""
effects = {
'Pitch Shift': librosa.effects.pitch_shift(audio, sr=sr, n_steps=2),
'Time Stretch': librosa.effects.time_stretch(audio, rate=1.2),
'Preemphasis': librosa.effects.preemphasis(audio),
'RMS Normalize': librosa.util.normalize(audio, norm=np.inf),
'Harmonic': librosa.effects.harmonic(audio),
'Percussive': librosa.effects.percussive(audio)
}
return effects
def print_stats(name, signal, original):
orig_adjusted, signal_adjusted = ensure_same_length(original, signal)
with warnings.catch_warnings():
warnings.simplefilter("ignore")
snr = calculate_snr(orig_adjusted, signal_adjusted)
rmse = calculate_rmse(orig_adjusted, signal_adjusted)
rms = calculate_rms(signal_adjusted)
snr_str = f"{snr:<10.4f}" if np.isfinite(snr) else "inf "
print(f"{name:<20} {np.min(signal_adjusted):<10.4f} {np.max(signal_adjusted):<10.4f} "
f"{rms:<10.4f} {snr_str} {rmse:<10.4f}")
def compare_all_methods(original, processed, sr):
original = ensure_numpy(original).flatten()
processed = ensure_numpy(processed).flatten()
custom_methods = ['normalize', 'peak_match', 'rms_match', 'adaptive', 'compress', 'hard_clip', 'soft_clip']
custom_boosted = {method: ensure_numpy(apply_boost(torch.from_numpy(original), torch.from_numpy(processed), method)) for method in custom_methods}
# Apply librosa effects
librosa_processed = librosa_effects(processed, sr)
all_methods = {**custom_boosted, **librosa_processed}
# Calculate number of rows
n_rows = (len(all_methods) + 2) // 2 # +1 for original vs processed, rounded up
# Create subplots
fig = make_subplots(rows=n_rows, cols=2,
subplot_titles=['Original vs Processed'] + list(all_methods.keys()),
vertical_spacing=0.05,
horizontal_spacing=0.1)
# Plot original vs processed
fig.add_trace(go.Scatter(y=original, name='Original', line=dict(color='blue')), row=1, col=1)
fig.add_trace(go.Scatter(y=processed, name='Processed (Unboosted)', line=dict(color='red')), row=1, col=1)
# Plot all processed signals
for i, (method, boosted) in enumerate(all_methods.items(), start=2):
row = (i + 1) // 2
col = 2 if i % 2 == 0 else 1
orig_adjusted, boosted_adjusted = ensure_same_length(original, boosted)
fig.add_trace(go.Scatter(y=orig_adjusted, name='Original', line=dict(color='blue'), showlegend=False), row=row, col=col)
fig.add_trace(go.Scatter(y=boosted_adjusted, name=f'Processed ({method})', line=dict(color='green')), row=row, col=col)
# Update layout
fig.update_layout(height=300*n_rows, width=1200,
title_text="Audio Signal Processing Comparison (Boosting)",
showlegend=False)
# Update all subplot axes
for i in range(1, n_rows+1):
for j in range(1, 3):
fig.update_xaxes(title_text="Sample", row=i, col=j)
fig.update_yaxes(title_text="Amplitude", row=i, col=j)
fig.show()
# Print statistics
print("Signal Statistics:")
print(f"{'Method':<20} {'Min':<10} {'Max':<10} {'RMS':<10} {'SNR (dB)':<10} {'RMSE':<10}")
print("-" * 70)
print_stats("Original", original, original)
print_stats("Processed", processed, original)
for method, boosted in all_methods.items():
print_stats(method, boosted, original)
compare_all_methods(cleaned_audio_data_raw, reconstructed_audio*10, sample_rate_raw)
Signal Statistics: Method Min Max RMS SNR (dB) RMSE ---------------------------------------------------------------------- Original -0.9950 0.9193 0.1552 inf 0.0000 Processed -1.8276 1.4438 0.1784 4.5901 0.0915 normalize -1.0000 0.7900 0.0976 5.0113 0.0872 peak_match -0.9950 0.7860 0.0971 4.9912 0.0874 rms_match -1.5894 1.2556 0.1552 5.4864 0.0825 adaptive -0.9625 0.7604 0.0940 4.8573 0.0887 compress -0.6638 0.9719 0.1230 3.4474 0.1043 hard_clip -0.8000 0.8000 0.1605 5.4301 0.0830 soft_clip -1.0359 1.0309 0.2308 1.5604 0.1297 Pitch Shift -1.1226 0.9461 0.1124 -1.8268 0.1915 Time Stretch -0.9322 0.7000 0.0793 -1.0081 0.1743 Preemphasis -2.0658 2.1298 0.2273 -2.4026 0.2046 RMS Normalize -1.0000 0.7900 0.0976 5.0113 0.0872 Harmonic -0.0098 0.0096 0.0004 -0.0000 0.1552 Percussive -1.8276 1.4438 0.1784 4.5905 0.0915
Comparison of Pruning Techniques¶
In [9]:
import torch
import numpy as np
import plotly.graph_objects as go
from plotly.subplots import make_subplots
def print_debug_info(freqs, magnitude_spectrum):
print(f"freqs shape: {freqs.shape}")
print(f"magnitude_spectrum shape: {magnitude_spectrum.shape}")
print(f"freqs dtype: {freqs.dtype}")
print(f"magnitude_spectrum dtype: {magnitude_spectrum.dtype}")
def threshold_pruning(data, threshold):
"""Keep only data points above a certain threshold."""
return data > threshold
def percentile_pruning(data, percentile):
"""Keep only data points above a certain percentile."""
threshold = np.percentile(data.cpu().numpy(), percentile)
return data > threshold
def downsample_pruning(data, factor):
"""Downsample the data by keeping every nth point."""
return torch.arange(len(data)) % factor == 0
def local_maxima_pruning(data, window_size):
"""Keep only local maxima within a sliding window."""
maxima = torch.zeros_like(data, dtype=bool)
for i in range(len(data)):
start = max(0, i - window_size // 2)
end = min(len(data), i + window_size // 2 + 1)
if data[i] == torch.max(data[start:end]):
maxima[i] = True
return maxima
def calculate_metrics(original, pruned, mask):
"""Calculate metrics to quantify data loss."""
retained_percentage = (mask.sum() / len(mask)) * 100
# Reconstruct pruned data to original size for MSE calculation
reconstructed = torch.zeros_like(original)
reconstructed[mask] = pruned
mse = torch.mean((original - reconstructed) ** 2).item()
# PSNR calculation
max_signal = torch.max(original)
psnr = 20 * torch.log10(max_signal / torch.sqrt(torch.tensor(mse))).item()
return retained_percentage, mse, psnr
def demonstrate_pruning(freqs, magnitude_spectrum):
print_debug_info(freqs, magnitude_spectrum)
# Ensure freqs and magnitude_spectrum have the same length
min_length = min(len(freqs), len(magnitude_spectrum))
freqs = freqs[:min_length]
magnitude_spectrum = magnitude_spectrum[:min_length]
fig = make_subplots(rows=5, cols=1, subplot_titles=(
"Original Spectrum", "Threshold Pruning", "90th Percentile Pruning",
"Downsample Pruning (factor=10)", "Local Maxima Pruning (window=20)"),
vertical_spacing=0.1)
freqs_tensor = torch.from_numpy(freqs).to(magnitude_spectrum.device)
# Original data
fig.add_trace(go.Scatter(x=freqs, y=magnitude_spectrum.cpu().numpy(), mode='lines', name='Original'), row=1, col=1)
pruning_methods = [
("Threshold", lambda x: threshold_pruning(x, torch.mean(x) + torch.std(x))),
("Percentile", lambda x: percentile_pruning(x, 90)),
("Downsample", lambda x: downsample_pruning(x, 10)),
("Local Maxima", lambda x: local_maxima_pruning(x, window_size=20))
]
for i, (name, pruning_func) in enumerate(pruning_methods, start=2):
mask = pruning_func(magnitude_spectrum)
pruned_data = magnitude_spectrum[mask]
pruned_freqs = freqs_tensor[mask].cpu().numpy()
fig.add_trace(go.Scatter(x=pruned_freqs, y=pruned_data.cpu().numpy(), mode='lines', name=name), row=i, col=1)
# Calculate and display metrics
retained, mse, psnr = calculate_metrics(magnitude_spectrum, pruned_data, mask)
metrics_text = f"Retained: {retained:.2f}%, MSE: {mse:.4f}, PSNR: {psnr:.2f} dB"
fig.add_annotation(x=0.5, y=1.1, xref=f"x{i}", yref=f"y{i}",
text=metrics_text, showarrow=False)
fig.update_layout(height=1800, width=1000, title_text="Spectrum Pruning Techniques")
for i in range(1, 6):
fig.update_xaxes(title_text="Frequency (Hz)", row=i, col=1)
fig.update_yaxes(title_text="Magnitude", row=i, col=1)
fig.show()
# Example usage
# Assuming magnitude_spectrum is a PyTorch tensor and sample_rate_raw is an integer
freqs = np.fft.fftfreq(len(magnitude_spectrum), d=1/sample_rate_raw)
demonstrate_pruning(freqs, magnitude_spectrum)
freqs shape: (4608,) magnitude_spectrum shape: torch.Size([4608]) freqs dtype: float64 magnitude_spectrum dtype: torch.float32
In [10]:
end = time.time()
print(f"Time taken: {end - start:.4f} seconds")
Time taken: 7.0250 seconds
In [11]:
sample_raw1 = train_dataset_raw[0]
sample_raw2 = train_dataset_raw[30]
sample_raw3 = train_dataset_raw[50]
sample_raw4 = train_dataset_raw[60]
sample_raw5 = train_dataset_raw[90]
sample_raw6 = train_dataset_raw[150]
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline
plt.style.use('default')
# Create figure with subplots
fig, axs = plt.subplots(3, 2, figsize=(8, 12), facecolor='black')
# Samples to plot
samples = [sample_raw1, sample_raw2, sample_raw3, sample_raw4, sample_raw5, sample_raw6]
for idx, (ax, sample) in enumerate(zip(axs.flat, samples)):
# Get audio data from tuple and convert to numpy, removing extra dimension
audio_data = sample[0].cpu().numpy().squeeze()
# Plot the waveform in white
ax.plot(audio_data, color='white', linewidth=0.5)
# Set y-axis limits
ax.set_ylim(-1.1, 1.1)
# Remove axes, labels, and ticks
ax.set_xticks([])
ax.set_yticks([])
for spine in ax.spines.values():
spine.set_visible(False)
# Set subplot background to black
ax.set_facecolor('black')
# Adjust layout and display
plt.tight_layout()
plt.savefig('speech_waveforms.png', dpi=300, bbox_inches='tight', facecolor='black')
plt.show()
Spike Compression¶
In [12]:
# Compression classes are imported from dyned.py
# (RunLengthEncoding, LZWCompressor, HuffmanCompressor, DeltaCompressor,
# BWTransform, LZ77Compressor, DyNEDcCompressor)
# See imports in the first cell.
In [13]:
plt1, cr1, sa1 = create_rle_spike_plot(freqs, quantised_step_signal)
plt2, cr2, sa2 = create_lzw_spike_plot(freqs, quantised_step_signal)
plt3, cr3, sa3 = create_huffman_spike_plot(freqs, quantised_step_signal)
plt4, cr4, sa4 = create_delta_spike_plot(freqs, quantised_step_signal)
plt5, cr5, sa5 = create_bwt_spike_plot(freqs, quantised_step_signal)
plt6, cr6, sa6 = create_lz77_spike_plot(freqs, quantised_step_signal)
plt7, cr7, sa7 = create_dynedc_v4_spike_plot(freqs, quantised_step_signal)
fig = make_subplots(rows=8, cols=1, subplot_titles=(
f"Original Spikes",
f"Run-Length Encoding ({cr1}, {sa1})",
f"LZW Compression ({cr2}, {sa2})",
f"Huffman Compression ({cr3}, {sa3})",
f"Delta Compression ({cr4}, {sa4})",
f"BWT Compression ({cr5}, {sa5})",
f"LZ77 Compression ({cr6}, {sa6})",
f"$\\text{{DyNED}}_{{c}}\\text{{ Compression ({cr7}, {sa7})}}$"
), vertical_spacing=0.09)
fig.add_trace(create_spike_plot(freqs, quantised_step_signal), row=1, col=1)
fig.add_trace(plt1, row=2, col=1)
fig.add_trace(plt2, row=3, col=1)
fig.add_trace(plt3, row=4, col=1)
fig.add_trace(plt4, row=5, col=1)
fig.add_trace(plt5, row=6, col=1)
fig.add_trace(plt6, row=7, col=1)
fig.add_trace(plt7, row=8, col=1)
fig.update_layout(height=1080, width=1400, title="Spike Compression Techniques - <i>bird</i> speech", xaxis_title='Frequency (Hz)', yaxis_title='Magnitude (dB)', showlegend=False)
fig.show()
# Compression summary table with bits-per-spike comparison
summary = compression_summary(quantised_step_signal)
print_compression_table(summary)
Number of spikes: 2341
Run-Length Encoding:
Original data size: 4608 bits
Compressed data size: 11126 bits
Compression ratio: 2.414
Space saving: -141.4%
Number of spikes: 2341
LZW Compression:
Original data size: 4608 bits
Compressed data size: 6680 bits
Compression ratio: 1.450
Space saving: -45.0%
Number of spikes: 2341
Huffman Compression:
Original data size: 4608 bits
Compressed data size: 4543 bits
Compression ratio: 0.986
Space saving: 1.4%
Number of spikes: 2341
Delta Compression:
Original data size: 4608 bits
Compressed data size: 5983 bits
Compression ratio: 1.298
Space saving: -29.8%
Number of spikes: 2341
Burrows-Wheeler Transform:
Original data size: 4608 bits
Compressed data size: 4103 bits
Compression ratio: 0.890
Space saving: 11.0%
Number of spikes: 2341
LZ77:
Original data size: 4608 bits
Compressed data size: 13659 bits
Compression ratio: 2.964
Space saving: -196.4%
Number of spikes: 2341
Mode selected: huffman
DyNEDc Compression:
Original data size: 4608 bits
Compressed data size: 4543 bits
Compression ratio: 0.986
Space saving: 1.4%
Number of spikes: 2341
Spike Train: 4608 bits, 2341 spikes Method Compressed Ratio Bits/Spike Saving -------------------------------------------------------- RLE 11126 2.4145 4.75 -141.4% LZW 6680 1.4497 2.85 -45.0% Huffman 4543 0.9859 1.94 1.4% Delta 5983 1.2984 2.56 -29.8% BWT (skipped) N/A N/A N/A LZ77 13659 2.9642 5.83 -196.4% DyNEDc 4602 0.9987 1.97 0.1% DyNEDc V2 6027 1.3079 2.57 -30.8% DyNEDc V3 4543 0.9859 1.94 1.4% DyNEDc V4 4543 0.9859 1.94 1.4%