Spatial Audio Equalization with (Active) Noise Filtering
Spatial Audio Equalization with Active Noise Cancellation (3D Audio System)
Milestone 1: Load and Prepare Audio
% Load an audio signal (e.g., a mono audio file)
filename = 'audio/Cmaj7.mp3';
b = stereo2mono(filename); % check if it's mono or not
if b == "False"
[audioSignal, audioFs] = audioread([filename(1 : end - 4) , '_mono.wav']);
else
[audioSignal, audioFs] = audioread(filename);
end
sound(audioSignal, audioFs)
Milestone 2: 3D Spatial Audio using HRTF
% Load HRIR data for HRTF (choose a random subject for binaural listening)
subjects = {'003', '008', '009', '010', '011', '012', '015', '017'};
subj = subjects{randi(length(subjects))};
subject_folder = strcat('standard_hrir_database/subject_', subj, '/');
load([subject_folder 'hrir_final.mat']); % Load hrir_l, hrir_r
% Define azimuth and elevation for 3D spatial sound
desired_azimuth = 200; % 30 degrees
desired_elevation = 0; % Horizontal direction (no elevation)
% Find closest azimuth and elevation in HRIR dataset
azim_v = -80:5:80;
elev_v = -45:5:90;
[~, az_idx] = min(abs(azim_v - desired_azimuth));
[~, el_idx] = min(abs(elev_v - desired_elevation));
% Extract HRIR for the chosen azimuth and elevation
hrir_left = squeeze(hrir_l(el_idx, az_idx, :));
hrir_right = squeeze(hrir_r(el_idx, az_idx, :));
% Convolve audio signal with HRIR for spatialization
y_left = conv(audioSignal, hrir_left);
y_right = conv(audioSignal, hrir_right);
% Combine left and right channels
spatialSignal = [y_left, y_right];
% Play the spatialized audio
sound(spatialSignal, audioFs);
audiowrite('spatialized_audio.wav', spatialSignal, audioFs); % Save audio
Milestone 3: Dynamic Audio Equalization using STFT
% Define parameters for STFT
windowSize = 1024;
overlap = windowSize / 2;
[S, F, T] = stft(audioSignal, audioFs, 'Window', hanning(windowSize), 'OverlapLength', overlap);
% Frequency bands: Bass (< 300Hz), Mid (300Hz-5kHz), Treble (> 5kHz)
BassAmp = 6.62; % Boost for bass
MidAmp = 2.0; % Boost for mids
TrebleAmp = 1.5; % Slight boost for treble
% Boost bass frequencies below 300 Hz
bassIndices = find(F < 300);
S(bassIndices, :) = BassAmp * S(bassIndices, :);
% Modify mid and treble frequencies
midIndices = find(F >= 300 & F < 5000);
S(midIndices, :) = MidAmp * S(midIndices, :);
trebleIndices = find(F >= 5000);
S(trebleIndices, :) = TrebleAmp * S(trebleIndices, :);
% Reconstruct signal from modified STFT
equalizedSignal = istft(S, audioFs, 'Window', hanning(windowSize), 'OverlapLength', overlap);
equalizedSignal = real(equalizedSignal) / max(abs(equalizedSignal)); % Normalize
% Play the equalized signal
sound(equalizedSignal, audioFs);
% Measure Low-frequency boost and SNR improvement
P_before = bandpower(audioSignal, audioFs, [0 300]); % Power before boost
P_after = bandpower(equalizedSignal, audioFs, [0 300]); % Power after boost
low_freq_boost = 10 * log10(P_after / P_before); % Boost in dB
% Simulate noise and calculate SNR
noise = 0.05 * randn(size(audioSignal));
eqnoise = 0.05 * randn(size(equalizedSignal));
snr_before = snr(audioSignal + noise, noise);
snr_after = snr(equalizedSignal + eqnoise, eqnoise);
fprintf('Low-frequency boost: %.2f dB\n', low_freq_boost);
fprintf('SNR improvement: %.2f dB\n', snr_after - snr_before);
Milestone 4: Active Noise Cancellation using NLMS Algorithm
% Generate noise and noisy signal
t = 0:1/audioFs:(length(equalizedSignal)-1)/audioFs; % Time vector
clean_signal = equalizedSignal; % Use equalized signal as clean signal
noise_signal = 0.1 * randn(size(equalizedSignal)); % Generate white noise
noisy_signal = clean_signal + noise_signal; % Create noisy signal
sound(noisy_signal, audioFs) % play noisy audio
% pause(size(noisy_signal)/audioFs);
% Define variables
mu = 1e-4; % Learning rate
filter_order = 64; % Length of the adaptive filter
epsilon = 1e-6; % Regularization term to prevent division by zero
% Initialize NLMS variables
weights = zeros(1, filter_order); % Initialize filter coefficients
error_signal = zeros(size(noisy_signal)); % Pre-allocate error signal
% Apply NLMS for noise cancellation
for i = filter_order:length(noisy_signal)
% Extract input vector (noise segment) of size filter_order
x = noisy_signal(i:-1:i-filter_order+1);
% Ensure x is a row vector for element-wise multiplication
x = x(:)'; % Force it to be a row vector
% Calculate filter output (dot product of weights and input vector)
y = dot(weights, x); % Alternatively, y = weights * x';
% Compute the error (difference between clean and filtered output)
error_signal(i) = clean_signal(i) - y;
% Update filter weights using the NLMS rule
norm_factor = x * x' + epsilon; % Normalize by input power
weights = weights + (2 * mu * error_signal(i) * x) / norm_factor;
end
% Play the noise-canceled signal
sound(error_signal, audioFs);
% Measure noise reduction in dB
P_noise_before = bandpower(noise_signal); % Power of noise before
P_noise_after = bandpower(error_signal); % Power of noise after
noise_reduction_dB = 10 * log10(P_noise_before / P_noise_after); % dB reduction
% Display noise reduction
fprintf('Noise reduction: %.2f dB\n', noise_reduction_dB);
% Visualizing the original audio signal
figure;
time_audio = (0:length(audioSignal)-1) / audioFs;
subplot(3,1,1);
plot(time_audio, audioSignal);
title('Original Audio Signal (Time Domain)');
xlabel('Time (s)');
ylabel('Amplitude');
grid on;
% Visualizing the HRIRs for left and right ears
time_hrir = (0:length(hrir_left)-1) / audioFs; % Time vector for HRIRs
subplot(3,1,2);
plot(time_hrir, hrir_left, 'b', time_hrir, hrir_right, 'r');
title('HRIR for Left (blue) and Right (red) Ears');
xlabel('Time (s)');
ylabel('Amplitude');
legend('HRIR Left', 'HRIR Right');
grid on;
% Visualizing the convolved audio signals (left and right ear signals)
time_convolved = (0:length(y_left)-1) / audioFs;
subplot(3,1,3);
plot(time_convolved, y_left, 'b', time_convolved, y_right, 'r');
title('Convolved Audio Signal (Spatialized Left and Right Ears)');
xlabel('Time (s)');
ylabel('Amplitude');
legend('Left Ear', 'Right Ear');
grid on;
% Compute Interaural Time Difference (ITD) and Interaural Level Difference (ILD)
% ITD: Time delay between left and right signals
[corr, lag] = xcorr(y_left, y_right); % Cross-correlation to find ITD
[~, max_idx] = max(abs(corr)); % Index of maximum correlation
ITD = lag(max_idx) / audioFs; % Convert lag to seconds
% ILD: Difference in intensity between left and right signals
ILD = 20 * log10(rms(y_right) / rms(y_left)); % ILD in decibels (dB)
% Display ITD and ILD
fprintf('Interaural Time Difference (ITD): %.4f seconds\n', ITD);
fprintf('Interaural Level Difference (ILD): %.4f dB\n', ILD);
% Additional plot to show ITD visually (cross-correlation)
figure;
plot(lag / audioFs, corr);
title('Cross-correlation Between Left and Right Ear Signals (for ITD)');
xlabel('Lag (seconds)');
ylabel('Correlation');
grid on;
