%%%%%MATLAB CODE BY BHUPATINDRA MALLA TO PROCESS NEAR-FIELD PRESSURE DATA%%%%
%%%%%%%%%%%%%%%%%%%%%%%Updated 06/25/2022%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc; clear all; close all

%%%READ DATA IN .mat FORMAT
%%% Load Pressure Data
root='C:\Users\Bhupa\Desktop\'
fdir=[root 'Mat p\']

file_detect = dir(fdir);
files=file_detect(3:23)% Since first two items are delimiters

na=31 %Total number of mics in the axial direction
nr=length(files) %Total number of mic array positions in radial direction

N=4096 % Block Size. Change value according to the frequency resolution desired.

mean_Xf=zeros(N,na,nr);
NB=zeros(N,na,nr);
prms=zeros(1,na,nr);

tic
for k = 1:length(files)%Total number of array positions
    
    load([fdir files(k).name])% Opening one array position at a time.
  
         for i = 1:na %Number of mics in one array position
            pt(:,1) = p(:,i);%%Just 1 mic at a time

        %%%Set up SIGNAL PARAMETERS for performing FFT
            fs = 204800; %Acquisition Frequency
            h = 1/fs; %Time resolution
            T = h*N; %Time window
            df = 1/T; %Frequency resolution
            fh = fs/2; %Nyquist frequency
            t = 0:h:T-h ;  %discrete time points
            f = 0:df:fs-df; f=f'; %discrete frequency points
         %%%Set up FILTER PARAMETERS 
                     Wl = 350; %Hz
                     Wh = 100000; %Hz
                    [b,a] = butter(5,[Wl Wh]./fh,'bandpass');
                     p(:,i) = filtfilt(b,a,pt(:,1)); %Apply Butterworth Filter
            %        p(:,i) = pt(:,i); %FOR NO FILTER
        %     prms(:,i) = norm(p(:,i))/sqrt(length(p(:,i)));

            w=rectwin(N); %sets rectangular window
            %w=hanning(N); %sets hanning window
         blks = length(pt)/N;% number of blocks               
                Xf = zeros(N,blks);
                pt_finish = 0;
                for n = 1:blks
                    pt_start = pt_finish + 1;
                    pt_finish = n*N;
                    Xf(:,n)  = 2*abs(fft(p(pt_start:pt_finish,i).*w(:,1)))/(N);       
                end        
                mean_Xf(:,i,k)=mean(Xf(:,:),2);
                NB(:,i,k) = 20*log10((mean_Xf(:,i,k)/(2^0.5))./20e-6);%%%Averaged Narrowband SPL in dB
                prms(1,i,k) = norm(p(:,i))/sqrt(length(p(:,i)));%For OASPL
                clear Xf 

         end
    clear p pt
end

 NB=NB(1:size(NB,1)/2,:,:);
 NB=permute(NB,[2 3 1]);
 f = f(1:size(f,1)/2); 
 OASPL_dB = 20.*log10(prms./20e-6);

toc
%%%%SAVE Processed DATA
dir_to_save=[fdir 'Processed\']
save([dir_to_save 'Processed_BS4096.mat'],'f','NB','OASPL_dB');

beep

