%% DATA IMPORT 
mic_raw = importdata('mic_raw.mat');

fs_mic = 2500000; %MHz Microphone Sampling Frequency


%% ROI Extraction

start_time = 0.5; % Choose your ROI lower limit
end_time = 1.5;  % Choose your ROI upper limit
start_index = floor(start_time * fs_mic);
end_index = floor(end_time * fs_mic);
numFiles_mic = size(mic_raw,2); 

num_samples_roi = end_index - start_index + 1;
mic_roi = zeros(num_samples_roi, numFiles_mic);  

time_roi = time_mic(start_index:end_index);

for i = 1:numFiles_mic
    mic_roi(:,i) = mic_raw(start_index:end_index, i);
end


%% PSD COMPUTATION
% psd_cell = cell(1, 3);  
f_kHz = [];
window = 1024; % Window length for STFT
noverlap = window; % Overlap 
nfft = window; % Number of FFT points

for i = 1:3
      [psd_cell{i}, f] = pwelch(mic_raw(:,i), window, noverlap, nfft, fs_mic);
    if isempty(f_kHz)
        f_kHz = f / 1000;  
    end
end
clear window noverlap nfft fs

%%%%% PSD PLOT %%%%%
figure;
hold on;

for i = 1:numFiles_mic
    plot(f_kHz, 10*log10(psd_cell{i}), 'LineWidth', 2,color="k");
end

legend('P_L = 1250 W', 'P_L = 1350 W', 'P_L = 2000 W')
xlabel('Frequency (kHz)', 'FontSize', 34, 'FontName', 'Times New Roman');  
ylabel('Power/Frequency (dB/Hz)', 'FontSize', 34, 'FontName', 'Times New Roman');
set(gca, 'FontSize', 44, 'fontname', 'Times New Roman', 'XMinorTick', 'on', 'YMinorTick', 'on');

grid off;
box on;
hold off;


%% Band-stop filter

f_min = 8000 ; % Hz, lower limit band-stop frequency
f_max = 26000 ; % Hz, upper limit band-stop frequency

for i=1:numFiles_mic
    mic_roi_stop(:,i) = bandstop(mic_roi(:,i),[fmin fmax],fs_mic,'Steepness',0.99,'StopbandAttenuation',40);
end
clear i 

%% PSD Quotient 

weld_1 = mic_roi(:,2); % Select the two weld tests you want to compare
Weld_2 = mic_roi(:,3);

psd_quotient(weld_1,weld_2)
grid minor

%% Sound power calculation

for i = 1:numFiles_mic
    rms_value = rms(mic_roi(:,i)); % RMS
    sound_power(i) = rms_value^2; % Sound power
end