added matlab files for rir analysis from internship repo

rir analysis objective evaluation/main1_mhby1g21.m

+% 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

rir analysis objective evaluation/main2_mhy1g21.m

+% 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

rir analysis objective evaluation/main3_mhby1g21.m

+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