%
disp(' ')
disp(' SPL.m   ver 1.5  October 19, 2012 ')
disp(' by Tom Irvine  Email: tomirvine@aol.com ')
disp(' ')
disp(' This program calculates the sound pressure level ')
disp(' of a time history ')
disp(' ')
disp(' The time history must be in a two-column matrix format: ')
disp(' Time(sec)  & amplitude(psi) ')
disp(' ')
%
close all;
%
clear n;
clear Y;
clear amp;
clear tim;
clear N;
clear df;
clear dt;
clear Mag;
clear full;
clear freq;
clear complex_FFT;
clear THM;
clear mu;
clear mean;
clear length;
%
fig_num=1;
%
[t,f,dt,sr,tmx,tmi,~,ncontinue]=enter_time_history();
THM=[t f];
%
disp(' Select input data pressure unit. ');
disp(' 1=psi   2=Pascal   3=micro Pascal  ');
%   
ipress=input(' ');   
%
if(ipress==1)
       p_unit=sprintf('psi');
	   ref = 2.9e-09;
end
if(ipress==2)
       p_unit=sprintf('Pa');    
	   ref = 20.e-06;
end
if(ipress==3)
       p_unit=sprintf('micro Pa');        
	   ref = 20.;
end
%
x_label=sprintf('Time(sec)');
y_label=sprintf('Pressure(%s)',p_unit);
t_string=sprintf('Time History');
[fig_num]=plot_TH(fig_num,x_label,y_label,t_string,THM);
%
disp(' ');
st=input(' Enter starting time (sec) ');
%
disp(' ');
te=input(' Enter ending time (sec) ');
%
j=1;
jfirst=1;
jlast=max(size(THM));
for i=1:max(size(THM))
    if(THM(i,1)<st)
        jfirst=i;
    end
    if(THM(i,1)>te)
        jlast=i;
        break;
    end
end
%
tim=double(THM(jfirst:jlast,1));
amp=double(THM(jfirst:jlast,2));    
%
n = max(size(amp));
%
N=2^(floor(log(n)/log(2)));
%
out4 = sprintf(' time history length = %d ',n);
disp(out4)
disp(' ');
%
mu=mean(amp);
amp=amp-mu;
%
if(ncontinue==1)
%
    [dt,df,mmm,NW,io]=FFT_advise(tim,amp);
%
    [mk,freq,time_a,dt,NW]=FFT_time_freq_set(mmm,NW,dt,df,tmi,io);
%
    [store,store_p,freq_p,max_a,max_f]=FFT_core_seg(NW,mmm,mk,freq,amp,io);                               
%
    store=store';
%
    sz=size(store);
    imax=sz(1);
    jmax=sz(2);
    full=zeros(imax,1);
    for i=1:imax
 %
        ms=0;  
        for j=1:jmax
            ms=ms+0.5*store(i,j)^2;
        end
 %
        full(i)=sqrt(ms/jmax);   % rms
    end
%
    full=sqrt(2)*full;  % peak
%
    [fl,fc,fu,imax]=one_third_octave_frequencies();
%
    [band_rms]=convert_FFT_to_one_third(freq,fl,fu,full);                   
%
    [splevel,oaspl]=convert_one_third_octave_mag_to_dB(band_rms,ref);
%
    [pf,pspl,oaspl]=trim_acoustic_SPL(fc,splevel,ref);  
%
    imax=length(pf);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
    disp(' ');
    disp(' Apply Weighting Network? ');
    disp('   1=No (physics & engineering) ');
    disp('   2=Yes (human hearing) ');
%
    iwn=input(' ');
%
    if(iwn==2)
%     
        [pf,pspl,oaspl,scale_label]=acoustic_weighting(pf,pspl,imax,ref);        
%
    end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%       
%
    out1=sprintf('\n OASPL = %8.4g dB  ref 20 micro Pa \n',oaspl);
    disp(out1);
%
    if(iwn==1)
        p_rms = ref*(10^(oaspl/20.));
        out2=sprintf('  pressure in air = %8.4g %s rms  \n',p_rms,p_unit);
        disp(out2);
    end
%
    figure(fig_num);
    fig_num=fig_num+1;
    plot(pf,pspl);
    ylabel('SPL (dB)');
    xlabel('Center Frequency (Hz)');
%
    if(iwn==1)
        out7=sprintf('One-Third Octave SPL   OASPL=%7.4g dB  Ref= 20 micro Pa',oaspl);
    else
        out7=sprintf('One-Third Octave SPL %s-scale  OASPL=%7.4g dB  Ref= 20 micro Pa',scale_label,oaspl);        
    end
%    
    title(out7);
%
    grid on;
    set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','lin');
%
    disp(' Write the SPL data to a external file? ')
    choice = input(' 1=yes 2=no ');
%
    if(choice == 1)
        disp(' ');   
        disp(' Enter the output filename ');
        filename1 = input(' ','s');
        fid1 = fopen(filename1,'w');  
        for k=1:max(size(pf))
            fprintf(fid1,'%11.4e \t %11.4e \n',pf(k), pspl(k));    
        end
        fclose(fid1);
    end    
%
    pf=fix_size(pf);
    pspl=fix_size(pspl);
%
    disp(' ');
    disp(' Matlab filename:  sound_pressure_level');
    sound_pressure_level=[pf pspl];
%
end