function out = em27_signoise(sig_dir, noise_dir)

% this function takes ascii spectra and calculates the signal to noise
% ratios (spectra are seperated into directories with signal spectra and
% noise spectra)


% test values
% clear all;
% sig_dir = 'sig/';
% noise_dir = 'noise/';

plotting = 0;

% read in all spectra
list_sig_dir = dir(sig_dir);
sig_spec_filenames = cellstr(char(list_sig_dir.name)); %list of just filenames
sig_spec_filenames = sig_spec_filenames(3:end); %removing the . and ..
list_noise_dir = dir(noise_dir);
noise_spec_filenames = cellstr(char(list_noise_dir.name)); %list of just noise spectra filenames
noise_spec_filenames = noise_spec_filenames(3:end); %removing the . and ..
start_date = noise_spec_filenames{1}(7:14);
end_date = noise_spec_filenames{end}(7:14);

sig_dir_up1_ind = strfind(sig_dir, '/');
if length(sig_dir_up1_ind) > 1;
    sig_dir_up1 = sig_dir(1:sig_dir_up1_ind(end-1));
else
    sig_dir_up1 = []; %empty
end


% % load subfunctions
curr_dir = pwd;
cd('../matlab_subfiles');
    read_asc_spectrum = @read_bin2asc_spectrum;
cd(curr_dir);

if length(sig_spec_filenames) ~= length(noise_spec_filenames);
    error('number of noise spectra is not equal to number of signal spectra')
end

ex_stds = 3; %exclude data outside of this std range
tail_prop = 1-normcdf(ex_stds); %this is proportion of data that should be in tails if data is gaussian
sig_nois = NaN*zeros(length(sig_spec_filenames),2);
for i1=1:length(sig_spec_filenames);
    
    [~, noise_spec_i1]=read_asc_spectrum([noise_dir, noise_spec_filenames{i1}]); %read in noise spectra (don't read in header parameters)
    std_mm = NaN*zeros(20, 3); %should converge within 20 rows, pre-allocating for speed - standard deviation, mean, median
%     out_points=NaN*zeros(20, 1); 
    
    std_mm(1,1) = std(noise_spec_i1.Signal);
    std_mm(1,2) = mean(noise_spec_i1.Signal);
    std_mm(1,3) = median(noise_spec_i1.Signal);
    
    nsignal = noise_spec_i1.Signal;
    out_points = floor(length(nsignal)*tail_prop); %this is the number of points allowed outside of range, can be set manually
    
    if plotting==1; figure(1); clf; hold on;
        plot(noise_spec_i1.Frequency, noise_spec_i1.Signal, '.-')
        nfreq_limits = [min(noise_spec_i1.Frequency), max(noise_spec_i1.Frequency)];
        plot(nfreq_limits, [std_mm(1,3)-std_mm(1,1), std_mm(1,3)-std_mm(1,1)], 'k-');
        plot(nfreq_limits, [std_mm(1,3)+std_mm(1,1)*ex_stds, std_mm(1,3)+std_mm(1,1)*ex_stds], 'k-');
    end
    
    i2=2;
    % calculate noise and remove residual point (not completely blacked out lines)
    while sum(nsignal>(std_mm(i2-1,3)+std_mm(i2-1,1)*ex_stds)) >out_points || sum(nsignal<(std_mm(i2-1,3)-std_mm(i2-1,1)*ex_stds)) >out_points;
        nsignal = nsignal(nsignal<std_mm(i2-1,3)+std_mm(i2-1,1)*ex_stds & nsignal>(std_mm(i2-1,3)-std_mm(i2-1,1)*ex_stds)); % trim some data outside of 3 sigma range
        
        % recalculate means and stds
        std_mm(i2,1) = std(nsignal);
        std_mm(i2,2) = mean(nsignal);
        std_mm(i2,3) = median(nsignal);
        
        if plotting ==1;
            plot(nfreq_limits, [std_mm(i2,3)-std_mm(i2,1)*ex_stds, std_mm(i2,3)-std_mm(i2,1)*ex_stds], 'g-');
            plot(nfreq_limits, [std_mm(i2,3)+std_mm(i2,1)*ex_stds, std_mm(i2,3)+std_mm(i2,1)*ex_stds], 'g-');
        end
        
        i2=i2+1;
    end
    i2=i2-1; %last iteration number didn't count
    % currently removing outliers based on an assumed gaussian
    % distribution. Could also remove based on difference between mean and
    % meadian relative to std.
    % std_mm(:,4) = (std_mm(:,2)-std_mm(:,3))./(std_mm(:,1)); %this typically decreases as the mean approaches the median 
    
    % calculate signal
    [~, signal_spec_i1]=read_asc_spectrum([sig_dir, sig_spec_filenames{i1}]); %read in noise spectra (don't read in header parameters)
    signal_ordered = sort(signal_spec_i1.Signal);
    near_max_sig = signal_ordered(round(length(signal_ordered)*0.99)); %take the 98% highest point, saying top 2% is noise
    
    
    if plotting == 1;
        plot(nfreq_limits, [std_mm(i2,3)-std_mm(i2,1), std_mm(i2,3)-std_mm(i2,1)], 'r-');
        plot(nfreq_limits, [std_mm(i2,3)+std_mm(i2,1), std_mm(i2,3)+std_mm(i2,1)], 'r-');
        figure(3); hist(nsignal, 20)
        title(['number iterations = ', num2str(i2)])
        figure(2); clf; hold on;
        plot(signal_spec_i1.Frequency, signal_spec_i1.Signal, '.-')
        freq_limits = [min(signal_spec_i1.Frequency), max(signal_spec_i1.Frequency)];
        plot(freq_limits, [1 1]*near_max_sig, 'k-')
        pause
    end
    
    
    sig_nois(i1,1) = near_max_sig-std_mm(i2,3); %subtracting noise median from signal
    sig_nois(i1,2) = std_mm(i2,1); %noise
end

sig_nois(:,3) = sig_nois(:,1)./sig_nois(:,2); %signal to noise ratio

% make/write snr file to insert into daily runlogs
snr_rl_filename = [sig_dir_up1, start_date, '_', end_date, '.snrrl']; %for direct use in runlog
snr_filename = [sig_dir_up1, start_date, '_', end_date, '.snr']; %with more decimal places


rl_fid = fopen(snr_rl_filename, 'w');
snr_fid = fopen(snr_filename, 'w');
fprintf(rl_fid,'%s\n',' SNR'); %adds column header
for i1=1:size(sig_nois,1);
    if i1==size(sig_nois,1); %don't print new line at end of file
        fprintf(rl_fid,'%4.f\n',sig_nois(i1,3)); %
        fprintf(snr_fid,'%7.2f',sig_nois(i1,3));
    else
        fprintf(rl_fid,'%4.f\n',sig_nois(i1,3)); %right size for snr spot in runlogs
        fprintf(snr_fid,'%7.2f\r\n',sig_nois(i1,3));
    end
end
fclose('all');

out.snr_filename = snr_filename;
exit
