From 6f973f1a4fb591c05c986ee7529a4c532364e1c7 Mon Sep 17 00:00:00 2001
From: mhby1g21 <mhby1g21@soton.ac.uk>
Date: Fri, 10 Jan 2025 12:21:01 +0000
Subject: [PATCH] added matlab files for rir analysis from internship repo
---
.../main1_mhby1g21.m | 130 ++++++++++++++++++
.../main2_mhy1g21.m | 91 ++++++++++++
.../main3_mhby1g21.m | 50 +++++++
3 files changed, 271 insertions(+)
create mode 100644 rir analysis objective evaluation/main1_mhby1g21.m
create mode 100644 rir analysis objective evaluation/main2_mhy1g21.m
create mode 100644 rir analysis objective evaluation/main3_mhby1g21.m
diff --git a/rir analysis objective evaluation/main1_mhby1g21.m b/rir analysis objective evaluation/main1_mhby1g21.m
new file mode 100644
index 0000000..175c3d0
--- /dev/null
+++ b/rir analysis objective evaluation/main1_mhby1g21.m
@@ -0,0 +1,130 @@
+% Parameters
+fs = 48000; % Sample rate
+T = 30; % Duration in seconds
+t = (0:1/fs:T-1/fs)';
+
+% Generate an exponential sine sweep
+function sweep = exp_sine_sweep(f0, f1, T, fs)
+ t = (0:1/fs:T-1/fs)';
+ w1 = 2*pi*f0;
+ w2 = 2*pi*f1;
+ K = T*w1/log(w2/w1);
+ L = T/log(w2/w1);
+ sweep = sin(K*(exp(t/L) - 1));
+end
+
+f0 = 20; % Start frequency
+f1 = 20000; % End frequency
+sweep = exp_sine_sweep(f0, f1, T, fs);
+
+% Apply fade-in and fade-out
+fade_length = round(0.1 * fs); % 100 ms fade
+fade = linspace(0, 1, fade_length)';
+sweep(1:fade_length) = sweep(1:fade_length) .* fade;
+sweep(end-fade_length+1:end) = sweep(end-fade_length+1:end) .* flip(fade);
+
+% Generate inverse filter
+inv_filter = flip(sweep) .* exp(-t/T);
+
+% Apply window function to inverse filter
+window = hann(length(inv_filter));
+inv_filter = inv_filter .* window;
+
+% Ensure double precision
+sweep = double(sweep);
+inv_filter = double(inv_filter);
+
+% Normalize both signals
+sweep = sweep / max(abs(sweep));
+inv_filter = inv_filter / max(abs(inv_filter));
+
+% Save sweep as .wav file
+audiowrite('sweep_30sec.wav', sweep, fs);
+disp('Sweep saved as sweep.wav');
+audiowrite('inverse_filter_30sec.wav', inv_filter, fs);
+disp('Inverse filter saved as inverse_filter.wav');
+
+% Use the original sweep as the recorded sweep for perfect impulse demonstration
+recorded_sweep = sweep;
+
+% Convolve recorded sweep with inverse filter to get the RIR
+rir = conv(recorded_sweep, inv_filter, 'same');
+
+% Normalize RIR
+rir = rir / max(abs(rir));
+
+% Plot the sweep and inverse filter
+figure;
+subplot(2,1,1);
+plot(t, sweep);
+xlabel('Time (s)');
+ylabel('Amplitude');
+title('Sweep Signal');
+xlim([0 T]);
+
+subplot(2,1,2);
+plot(t, inv_filter);
+xlabel('Time (s)');
+ylabel('Amplitude');
+title('Inverse Filter');
+xlim([0 T]);
+
+% Plot the RIR
+figure;
+t_rir = (-length(rir)/2:length(rir)/2-1) / fs * 1000; % Time in milliseconds
+plot(t_rir, rir);
+xlabel('Time (ms)');
+ylabel('Amplitude');
+title('Room Impulse Response');
+xlim([-10 10]); % Zoom in to ±10 ms around the center
+
+% Calculate frequency responses
+f = linspace(0, fs/2, length(rir)/2);
+Y_sweep = abs(fft(sweep));
+Y_sweep = Y_sweep(1:length(Y_sweep)/2);
+Y_inv = abs(fft(inv_filter));
+Y_inv = Y_inv(1:length(Y_inv)/2);
+Y_rir = abs(fft(rir));
+Y_rir = Y_rir(1:length(Y_rir)/2);
+
+% Plot the frequency responses
+figure;
+semilogx(f, 20*log10(Y_sweep/max(Y_sweep)), 'b', 'LineWidth', 1.5);
+hold on;
+semilogx(f, 20*log10(Y_inv/max(Y_inv)), 'r', 'LineWidth', 1.5);
+semilogx(f, 20*log10(Y_rir/max(Y_rir)), 'g', 'LineWidth', 1.5);
+xlabel('Frequency (Hz)');
+ylabel('Magnitude (dB)');
+title('Frequency Responses');
+xlim([20 20000]);
+ylim([-60 5]);
+legend('Sweep', 'Inverse Filter', 'RIR');
+grid on;
+
+% Function to calculate RT60 with noise floor cutoff
+function rt60 = calculate_rt60(ir, fs, noise_floor_db)
+ ir_db = 20 * log10(abs(ir) / max(abs(ir)));
+ [~, peak_index] = max(abs(ir));
+ ir_db = ir_db(peak_index:end);
+
+ % Find the point where the decay reaches noise floor
+ decay_end = find(ir_db < noise_floor_db, 1);
+ if isempty(decay_end)
+ decay_end = length(ir_db);
+ end
+
+ % Linear fit to the decay curve
+ t = (0:decay_end-1)' / fs;
+ p = polyfit(t, ir_db(1:decay_end), 1);
+
+ % Calculate RT60
+ rt60 = -60 / p(1);
+end
+
+% Calculate and display RT60
+rt60 = calculate_rt60(rir, fs, -40); % Use -40 dB as noise floor
+disp(['Calculated RT60: ', num2str(rt60), ' seconds']);
+
+% Save RIR as .wav file
+audiowrite('perfect_room_impulse_response_30sec.wav', rir, fs);
+disp('Perfect room Impulse Response saved as perfect_room_impulse_response.wav');
\ No newline at end of file
diff --git a/rir analysis objective evaluation/main2_mhy1g21.m b/rir analysis objective evaluation/main2_mhy1g21.m
new file mode 100644
index 0000000..1bab0b8
--- /dev/null
+++ b/rir analysis objective evaluation/main2_mhy1g21.m
@@ -0,0 +1,91 @@
+% RIR extraction and analysis script
+
+clear all;
+close all;
+
+% Parameters
+fs = 48000; % Sample rate (Hz)
+
+try
+ % Load the original sweep and inverse filter
+ [sweep, fs_sweep] = audioread('sweep.wav');
+ [inv_filter, fs_inv] = audioread('inverse_filter.wav');
+
+ % Ensure the sample rates match
+ assert(fs_sweep == fs && fs_inv == fs, 'Sample rates of sweep and inverse filter must match the specified fs');
+
+ % Read the recorded sweep
+ [recorded_sweep, fs_recorded] = audioread('Recorded/LR_matlabsweep_-20db_attenuation_defaultbaked.wav');
+
+ % Convert recorded sweep to mono if it's stereo
+ if size(recorded_sweep, 2) > 1
+ recorded_sweep = mean(recorded_sweep, 2);
+ end
+
+ % Resample recorded sweep if necessary
+ if fs_recorded ~= fs
+ recorded_sweep = resample(recorded_sweep, fs, fs_recorded);
+ end
+
+ % Print diagnostic information
+ disp(['Length of original sweep: ', num2str(length(sweep))]);
+ disp(['Length of recorded sweep: ', num2str(length(recorded_sweep))]);
+ disp(['Length of inverse filter: ', num2str(length(inv_filter))]);
+
+ % Perform deconvolution in frequency domain
+ N = length(recorded_sweep);
+ F_recorded = fft(recorded_sweep, N);
+ F_inverse = fft(inv_filter, N);
+ F_rir = F_recorded .* F_inverse;
+ rir = real(ifft(F_rir));
+
+ % Normalize RIR
+ rir = rir / max(abs(rir));
+
+ % Plot the full RIR
+ figure;
+ t_rir = (0:length(rir)-1) / fs;
+ plot(t_rir, rir);
+ xlabel('Time (s)');
+ ylabel('Amplitude');
+ title('Full Room Impulse Response');
+
+ % Plot the RIR in dB scale
+ figure;
+ rir_db = 20 * log10(abs(rir) / max(abs(rir)));
+ plot(t_rir, rir_db);
+ xlabel('Time (s)');
+ ylabel('Amplitude (dB)');
+ title('Room Impulse Response (dB scale)');
+ ylim([-60 0]);
+
+ % Plot the frequency response
+ figure;
+ f = (0:N-1) * fs / N;
+ F_rir_mag = abs(F_rir);
+ semilogx(f, 20*log10(F_rir_mag / max(F_rir_mag)));
+ xlabel('Frequency (Hz)');
+ ylabel('Magnitude (dB)');
+ title('Frequency Response of RIR');
+ xlim([20 20000]); % Focus on audible range
+ grid on;
+
+ % Plot the inverse filter
+ figure;
+ t_inv = (0:length(inv_filter)-1) / fs;
+ plot(t_inv, inv_filter);
+ xlabel('Time (s)');
+ ylabel('Amplitude');
+ title('Inverse Filter');
+
+ % Save RIR as .wav file
+ audiowrite('room_impulse_response.wav', rir, fs);
+ disp('Room Impulse Response saved as room_impulse_response.wav');
+
+catch err
+ % Display error message
+ disp('An error occurred:');
+ disp(err.message);
+ % Display the line where the error occurred
+ disp(['Error in line: ', num2str(err.stack(1).line)]);
+end
\ No newline at end of file
diff --git a/rir analysis objective evaluation/main3_mhby1g21.m b/rir analysis objective evaluation/main3_mhby1g21.m
new file mode 100644
index 0000000..1788419
--- /dev/null
+++ b/rir analysis objective evaluation/main3_mhby1g21.m
@@ -0,0 +1,50 @@
+addpath 'IoSR Toolbox' 'octave'
+
+[LS2_sweep, fs2] = audioread("room_impulse_response.wav");
+[sweep, fsst] = audioread("sweep.wav");
+[RT60, DRR, C50, Cfs, EDT] = iosr.acoustics.irStats('room_impulse_response.wav','graph', true, 'spec', 'full', 'y_fit', [-5 -35]);
+
+% Calculating Mean Values
+mean_RT60 = mean(RT60(3:8));
+mean_EDT = mean(EDT(3:8));
+
+% Plot RIR
+t2 = 0:1/fs2:((length(LS2_sweep)-1)/fs2);
+figure;
+plot(t2,LS2_sweep(:,1).'); xlabel("time [s]"); ylabel("Amplitude"); title("RIR from sweep");
+
+% Extract RT60 values for specific frequencies
+freq_indices = find(ismember(Cfs, [500, 1000, 2000, 4000, 8000]));
+RT60_values = RT60(freq_indices);
+frequencies = Cfs(freq_indices);
+
+% Create the graph
+figure;
+plot(1:5, RT60_values, '-o', 'LineWidth', 2, 'MarkerSize', 8, 'MarkerFaceColor', 'auto');
+set(gca, 'XTick', 1:5, 'XTickLabel', {'0.5', '1', '2', '4', '8'});
+xlabel('Frequency (kHz)');
+ylabel('RT60 (s)');
+title('KT - RT60 vs Frequency');
+grid on;
+
+
+% Add value labels on top of each point
+for i = 1:length(RT60_values)
+ text(i, RT60_values(i), sprintf('%.2f', RT60_values(i)), ...
+ 'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'center');
+end
+
+% Display Mean Values
+disp('Mean RT60:');
+disp(mean_RT60);
+
+disp('Mean EDT:');
+disp(mean_EDT);
+
+% Print RT60 and EDT for each octave band contributing to the mean
+disp('RT60 and EDT for each contributing octave band:');
+disp('Frequency (Hz) | RT60 (s) | EDT (s)');
+disp('----------------------------------------');
+for i = 3:8
+ fprintf('%13d | %8.2f | %7.2f\n', Cfs(i), RT60(i), EDT(i));
+end
\ No newline at end of file
--
GitLab