function [R, C, dt_reflection] = menmednormpol(quietwin, IARwins, IARwin, menquants,normalization,polarity,lEfilt,winunique,ptimes,B_output,peakwint,s,FFTscale,noiselevel,Brmsquiet,datestring, stationstring, wintlims, delayguess)

% Function to create windows of spike shapes and quiet period
global fntsize res xsize ysize lmarg bmarg widthv heightv xstr titlestr fstr

global fs

quietwinindex = floor((quietwin(1)/(fs*3600))/IARwin(3))+1;

IARwinsindex = int16(IARwins/(fs*3600)/IARwin(3)+1);

detsubplots = length(IARwinsindex);

dt_reflection = [];

switch detsubplots
    case {1,2}
        R = 2;
        C = 1;
    case {3,4}
        R = 2;
        C = 2;
    case {5,6}
        R = 2;
        C = 3;
    case {7,8,9}
        R = 3;
        C = 3;
    case {10,11,12}
        R = 3;
        C = 4;
    case {13,14,15,16}
        R = 4;
        C = 4;
    case {17,18}
        R = 3;
        C = 6;
    case {19,20}
        R = 4;
        C = 5;
end

if strcmp(normalization,'unnormalized')
    if strcmp(polarity,'normalpolarity')
        if strcmp(lEfilt,'on')
            nn = 9;
            nstring = '';
            pstring = '';
            ptimes(:,9);
        else
            nn = 3;
            nstring = '';
            pstring = '';
        end
    elseif strcmp(polarity,'positivepolarity')
        if strcmp(lEfilt,'on')
            nn = 10;
            nstring = '';
            pstring ='+ve';
            ptimes(:,10)
        else
            nn = 5;
            nstring = '';
            pstring = '+ve';
        end
    end
elseif strcmp(normalization,'normalized')
    if strcmp(polarity,'normalpolarity')
        nn = 4;
        nstring = 'normalized';
        pstring = '';
    elseif strcmp(polarity,'positivepolarity')
        nn = 6;
        nstring = 'normalized';
        pstring = '+ve';
    end
elseif strcmp(normalization,'other')
    if strcmp(polarity,'normalpolarity')
        nn = 7;
        nstring = '';
        pstring = '';
    elseif strcmp(polarity,'positivepolarity')
        nn = 8;
        nstring = '';
        pstring = '+ve';
    end
end
    
    sf1 = figure('Position', [0 0 s(3) s(4)]);%('Position', [0 0 s(3) s(4)]); % New figure
    sf2 = figure('Position', [0 0 s(3) s(4)]);%('Position', [0 0 s(3) s(4)]); % New figure
    switch menquants
        case 'mean'
            % Mean of Quiet peak shape
            peaksxvec = wintlims(1)*peakwint:1/fs:wintlims(2)*peakwint-1/fs;
            if strcmp(winunique,'unique')
                if nn == 9
                    quietpeaks = mean(ptimes{quietwinindex,3}(ptimes{quietwinindex,14},:));
                elseif nn == 10
                    quietpeaks = mean(ptimes{quietwinindex,5}(ptimes{quietwinindex,14},:));
                else
                    quietpeaks = mean(ptimes{quietwinindex,nn}(ptimes{quietwinindex,13},:));
                end
            else
                quietpeaks = mean(cell2mat(ptimes(quietwinindex,nn)));
            end
            % Find reflection peak after superposition peak
            %[qpksabs qlocsabs] = findpeaks(abs(quietpeaks(floor((abs(wintlims(1))*peakwint+delayguess)*fs):end)),'SortStr','descend','NPeaks',3);
            %true_qlocsabs = qlocsabs +  floor((abs(wintlims(1))*peakwint+delayguess)*fs) - 1;
            j=1;
            %dt_reflection(j,:) = (true_qlocsabs - abs(wintlims(1))*peakwint*fs) .* 1/fs;
            %j=j+1;
            % Plot SEA and reflection peak
            figure(sf1);
            ax1(1) = subplot(R,C,1);
            plot(peaksxvec, quietpeaks);
            hold(ax1(1));
            %plot(peaksxvec(true_qlocsabs),quietpeaks(true_qlocsabs),'m*');
            line([peaksxvec(1) peaksxvec(end)],[Brmsquiet(2) Brmsquiet(2)],'Color','k','LineStyle','--');
            line([peaksxvec(1) peaksxvec(end)],[-Brmsquiet(2) -Brmsquiet(2)],'Color','k','LineStyle','--');
            %line([peaksxvec(1) peaksxvec(end)],[noiselevel(quietwinindex) noiselevel(quietwinindex)],'Color','r','LineStyle','--');
            %line([peaksxvec(1) peaksxvec(end)],[-noiselevel(quietwinindex) -noiselevel(quietwinindex)],'Color','k','LineStyle','--');
            %line([0 peakwint],[rms(quietpeaks) rms(quietpeaks)],'Color',[0.5 0 0],'LineStyle','--');
            %line([0 peakwint],[-rms(quietpeaks) -rms(quietpeaks)],'Color',[0.5 0 0],'LineStyle','--');
            %line([0.2*peakwint 0.2*peakwint],[min(quietpeaks(1,:))-1 max(quietpeaks(1,:))+1],'Color','k','LineStyle','--');
            hold off;
            grid on;
            xlabel('Time (s.ss)', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
            ylabel('B (nT)', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axestitle('Superposed peak signatures', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
            title({strcat(datestring, stationstring,'\_','Mean of Quiet (~',num2str(round((quietwin(1)/(fs*3600))/IARwin(3))*IARwin(3)),'-',num2str(round((quietwin(2)/(fs*3600))/IARwin(3))*IARwin(3)),' UTC)'); strcat('peak signatures (',nstring,',',pstring,')')}, 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
            legend('Mean SEA','\sigma_{quiet}');
            % Plot FFT of Quiet peak shape
            [quietfvec, quietfftfreq] = singsidefft(quietpeaks);
            figure(sf2);
            ax2(1)=subplot(R,C,1);
            if strcmp(FFTscale,'log')
                plot(quietfvec, log(quietfftfreq));
                ylabel('log(Magnitude)', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
            else
                plot(quietfvec, quietfftfreq);
                ylabel('Magnitude', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
            end
            grid on;
            xlabel('Frequency (Hz)', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
            title({strcat(datestring, stationstring,'\_','FFT of Mean of Quiet (~',num2str(round((quietwin(1)/(fs*3600))/IARwin(3))*IARwin(3)),'-',num2str(round((quietwin(2)/(fs*3600))/IARwin(3))*IARwin(3)),' UTC)'); strcat('peak signatures (',nstring,',',pstring,')')}, 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
            
            % Mean and FFT of IAR peak shapes
            
            for i = 1:length(IARwinsindex)-1
                if strcmp(winunique,'unique')
                    if nn == 9
                        IARpeaks = mean(ptimes{IARwinsindex(i),3}(ptimes{IARwinsindex(i),14},:));
                    elseif nn == 10
                        IARpeaks = mean(ptimes{IARwinsindex(i),5}(ptimes{IARwinsindex(i),14},:));
                    else
                        IARpeaks = mean(ptimes{IARwinsindex(i),nn}(ptimes{IARwinsindex(i),13},:));
                    end
                else
                    IARpeaks = mean(cell2mat(ptimes(IARwinsindex(i),nn)));
                end
                % Find reflection peak after superposition peak
                [pksabs, locsabs] = findpeaks(abs(IARpeaks(floor((abs(wintlims(1))*peakwint+delayguess)*fs):end)),'SortStr','descend','NPeaks',3);
                true_qlocsabs = locsabs +  floor((abs(wintlims(1))*peakwint+delayguess)*fs) - 1;
                dt_reflection(j,:) = (true_qlocsabs - abs(wintlims(1))*peakwint*fs) .* 1/fs;
                j=j+1;
                % Plot SEA and reflection peak
                figure(sf1);
                ax1(i+1)=subplot(R,C,i+1);
                plot(peaksxvec, IARpeaks);
                hold(ax1(i+1));
                %plot(peaksxvec(true_qlocsabs),IARpeaks(true_qlocsabs),'m*');
                line([peaksxvec(1) peaksxvec(end)],[Brmsquiet(2) Brmsquiet(2)],'Color','k','LineStyle','--');
                line([peaksxvec(1) peaksxvec(end)],[-Brmsquiet(2) -Brmsquiet(2)],'Color','k','LineStyle','--');
                %line([peaksxvec(1) peaksxvec(end)],[noiselevel(IARwinsindex(i)) noiselevel(IARwinsindex(i))],'Color','k','LineStyle','--');
                %line([peaksxvec(1) peaksxvec(end)],[-noiselevel(IARwinsindex(i)) -noiselevel(IARwinsindex(i))],'Color','k','LineStyle','--');
                %line([0 peakwint],[rms(IARpeaks) rms(IARpeaks)],'Color',[0.5 0 0],'LineStyle','--');
                %line([0 peakwint],[-rms(IARpeaks) -rms(IARpeaks)],'Color',[0.5 0 0],'LineStyle','--');
                %line([0.2*peakwint 0.2*peakwint],[min(IARpeaks(1,:))-1 max(IARpeaks(1,:))+1],'Color','k','LineStyle','--');
                hold off;
                grid on;
                xlabel('Time (s.ss)', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
                ylabel('B', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axestitle('Superposed peak signatures', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
                title({strcat(datestring, stationstring,'\_','Mean of IAR (~',num2str(IARwins(i)/(fs*3600)),'-',num2str(IARwins(i+1)/(fs*3600)),' UTC)'); strcat('peak signatures (',nstring,',',pstring,')')}, 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
                legend('Mean SEA','\sigma_{quiet}');
                % Plot FFT of Quiet peak shape
                [IARfvec, IARfftfreq] = singsidefft(IARpeaks);
                figure(sf2);
                ax2(i+1) = subplot(R,C,i+1);
                if strcmp(FFTscale,'log')
                    plot(IARfvec, log(IARfftfreq));
                    ylabel('log(Magnitude)', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
                else
                    plot(IARfvec, IARfftfreq);
                    ylabel('Magnitude', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
                end
                grid on;
                xlabel('Frequency (Hz)', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
                title({strcat(datestring, stationstring,'\_','FFT of Mean of IAR (~',num2str(IARwins(i)/(fs*3600)),'-',num2str(IARwins(i+1)/(fs*3600)),' UTC)'); strcat('peak signatures (',nstring,',',pstring,')')}, 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
            end
        case 'quartiles'
            % Quartiles
            peaksxvec = wintlims(1)*peakwint:1/fs:wintlims(2)*peakwint-1/fs;
            if strcmp(winunique,'unique')
                if nn == 9
                    quietpeaks = quantile(ptimes{quietwinindex,3}(ptimes{quietwinindex,14},:),[.25 .50 .75]);
                elseif nn == 10
                    quietpeaks = quantile(ptimes{quietwinindex,5}(ptimes{quietwinindex,14},:),[.25 .50 .75]);
                else
                    quietpeaks = quantile(ptimes{quietwinindex,nn}(ptimes{quietwinindex,13},:),[.25 .50 .75]);
                end
            else
                quietpeaks = quantile(cell2mat(ptimes(quietwinindex,nn)),[.25 .50 .75]);
            end            
            figure(sf1);
            ax1(1) = subplot(R,C,1);
            if length(peaksxvec) ~= length(quietpeaks)
                plot(quietpeaks);
                skip = 1;
            else
                plot(peaksxvec, quietpeaks);
                skip = 0;
            end
            hold(ax1(1));
            line([peaksxvec(1) peaksxvec(end)],[Brmsquiet(2) Brmsquiet(2)],'Color','k','LineStyle','--');
            line([peaksxvec(1) peaksxvec(end)],[-Brmsquiet(2) -Brmsquiet(2)],'Color','k','LineStyle','--');
            %line([peaksxvec(1) peaksxvec(end)],[noiselevel(quietwinindex) noiselevel(quietwinindex)],'Color','k','LineStyle','--');
            %line([peaksxvec(1) peaksxvec(end)],[-noiselevel(quietwinindex) -noiselevel(quietwinindex)],'Color','k','LineStyle','--');
            if skip == 0
            %line([0 peakwint],[rms(quietpeaks(2,:)) rms(quietpeaks(2,:))],'Color',[0.5 0 0],'LineStyle','--');
            %line([0 peakwint],[-rms(quietpeaks(2,:)) -rms(quietpeaks(2,:))],'Color',[0.5 0 0],'LineStyle','--');
            %line([0.2*peakwint 0.2*peakwint],[min(quietpeaks(1,:))-1 max(quietpeaks(3,:))+1],'Color','k','LineStyle','--');
            end
            hold off;
            grid on;
            xlabel('Time (s.ss)', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
            ylabel('B', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axestitle('Superposed peak signatures', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
            title({strcat(datestring, stationstring,'\_','Quartiles of Quiet (~',num2str(round((quietwin(1)/(fs*3600))/IARwin(3))*IARwin(3)),'-',num2str(round((quietwin(2)/(fs*3600))/IARwin(3))*IARwin(3)),' UTC)'); strcat('peak signatures (',nstring,',',pstring,')')}, 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
            
            % Plot FFT of Quiet peak shape
            [quietfvec, quietfftfreq] = singsidefft(quietpeaks);
            figure(sf2);
            ax2(1)=subplot(R,C,1);
            if strcmp(FFTscale,'log')
                plot(quietfvec, log(quietfftfreq));
                ylabel('log(Magnitude)', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
            else
                plot(quietfvec, quietfftfreq);
                ylabel('Magnitude', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
            end
            grid on;
            xlabel('Frequency (Hz)', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
            title({strcat(datestring, stationstring,'\_','FFT of Quartiles of Quiet (~',num2str(round((quietwin(1)/(fs*3600))/IARwin(3))*IARwin(3)),'-',num2str(round((quietwin(2)/(fs*3600))/IARwin(3))*IARwin(3)),' UTC)'); strcat('peak signatures (',nstring,',',pstring,')')}, 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
            
            % Quartiles and FFT of Quartiles for IAR
            for i = 1:length(IARwinsindex)-1
                if strcmp(winunique,'unique')
                    if nn == 9
                        IARpeaks = quantile(ptimes{IARwinsindex(i),3}(ptimes{IARwinsindex(i),14},:),[.25 .50 .75]);
                    elseif nn == 10
                        IARpeaks = quantile(ptimes{IARwinsindex(i),5}(ptimes{IARwinsindex(i),14},:),[.25 .50 .75]);
                    else
                        IARpeaks = quantile(ptimes{IARwinsindex(i),nn}(ptimes{IARwinsindex(i),13},:),[.25 .50 .75]);
                    end
                else
                    IARpeaks = quantile(cell2mat(ptimes(IARwinsindex(i),nn)),[.25 .50 .75]);
                end         
                figure(sf1);
                ax1(i+1)=subplot(R,C,i+1);
                if length(peaksxvec) ~= length(IARpeaks)
                plot(IARpeaks);
                skip2 = 1;
                else
                plot(peaksxvec, IARpeaks);
                skip2 = 0;
                end
                %plot(peaksxvec, IARpeaks);
                hold(ax1(i+1));
                line([peaksxvec(1) peaksxvec(end)],[Brmsquiet(2) Brmsquiet(2)],'Color','k','LineStyle','--');
                line([peaksxvec(1) peaksxvec(end)],[-Brmsquiet(2) -Brmsquiet(2)],'Color','k','LineStyle','--');
                %line([peaksxvec(1) peaksxvec(end)],[noiselevel(IARwinsindex(i)) noiselevel(IARwinsindex(i))],'Color','k','LineStyle','--');
                %line([peaksxvec(1) peaksxvec(end)],[-noiselevel(IARwinsindex(i)) -noiselevel(IARwinsindex(i))],'Color','k','LineStyle','--');
                if skip2 == 0
                %line([0 peakwint],[rms(IARpeaks(2,:)) rms(IARpeaks(2,:))],'Color',[0.5 0 0],'LineStyle','--');
                %line([0 peakwint],[-rms(IARpeaks(2,:)) -rms(IARpeaks(2,:))],'Color',[0.5 0 0],'LineStyle','--');
                %line([0.2*peakwint 0.2*peakwint],[min(IARpeaks(1,:))-1 max(IARpeaks(3,:))+1],'Color','k','LineStyle','--');
                end
                hold off;
                grid on;
                xlabel('Time (s.ss)', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
                ylabel('B', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axestitle('Superposed peak signatures', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
                title({strcat(datestring, stationstring,'\_','Quartiles of IAR (~',num2str(IARwins(i)/(fs*3600)),'-',num2str(IARwins(i+1)/(fs*3600)),' UTC)'); strcat('peak signatures (',nstring,',',pstring,')')}, 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
                % Plot FFT of IAR peak shape
                [IARfvec, IARfftfreq] = singsidefft(IARpeaks);
                figure(sf2);
                ax2(i+1) = subplot(R,C,i+1);
                if strcmp(FFTscale,'log')
                    plot(IARfvec, log(IARfftfreq));
                    ylabel('log(Magnitude)', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
                else
                    plot(IARfvec, IARfftfreq);
                    ylabel('Magnitude', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
                end
                grid on;
                xlabel('Frequency (Hz)', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
                title({strcat(datestring, stationstring,'\_','FFT of Quartiles of IAR (~',num2str(IARwins(i)/(fs*3600)),'-',num2str(IARwins(i+1)/(fs*3600)),' UTC)'); strcat('peak signatures (',nstring,',',pstring,')')}, 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
            end
        case 'deciles'
            % deciles
            peaksxvec = wintlims(1)*peakwint:1/fs:wintlims(2)*peakwint-1/fs;     
            if strcmp(winunique,'unique')
                if nn == 9
                    quietpeaks = quantile(ptimes{quietwinindex,3}(ptimes{quietwinindex,14},:),[.1 .2 .3 .4 .5 .6 .7 .8 .9]);
                elseif nn == 10
                    quietpeaks = quantile(ptimes{quietwinindex,5}(ptimes{quietwinindex,14},:),[.1 .2 .3 .4 .5 .6 .7 .8 .9]);
                else
                    quietpeaks = quantile(ptimes{quietwinindex,nn}(ptimes{quietwinindex,13},:),[.1 .2 .3 .4 .5 .6 .7 .8 .9]);
                end
            else
                quietpeaks = quantile(cell2mat(ptimes(quietwinindex,nn)),[.1 .2 .3 .4 .5 .6 .7 .8 .9]);
            end 
            figure(sf1);
            ax1(1) = subplot(R,C,1);
            plot(peaksxvec, quietpeaks);
            hold(ax1(1));
            line([peaksxvec(1) peaksxvec(end)],[Brmsquiet(2) Brmsquiet(2)],'Color','k','LineStyle','--');
            line([peaksxvec(1) peaksxvec(end)],[-Brmsquiet(2) -Brmsquiet(2)],'Color','k','LineStyle','--');
            %line([peaksxvec(1) peaksxvec(end)],[noiselevel(quietwinindex) noiselevel(quietwinindex)],'Color','k','LineStyle','--');
            %line([peaksxvec(1) peaksxvec(end)],[-noiselevel(quietwinindex) -noiselevel(quietwinindex)],'Color','k','LineStyle','--');
            %line([0 peakwint],[rms(quietpeaks(5,:)) rms(quietpeaks(5,:))],'Color',[0.5 0 0],'LineStyle','--');
            %line([0 peakwint],[-rms(quietpeaks(5,:)) -rms(quietpeaks(5,:))],'Color',[0.5 0 0],'LineStyle','--');
            %line([0.2*peakwint 0.2*peakwint],[min(quietpeaks(1,:))-1 max(quietpeaks(3,:))+1],'Color','k','LineStyle','--');
            hold off;
            grid on;
            xlabel('Time (s.ss)', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
            ylabel('B', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axestitle('Superposed peak signatures', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
            title({strcat(datestring, stationstring,'\_','Deciles of Quiet (~',num2str(round((quietwin(1)/(fs*3600))/IARwin(3))*IARwin(3)),'-',num2str(round((quietwin(2)/(fs*3600))/IARwin(3))*IARwin(3)),' UTC)'); strcat('peak signatures (',nstring,',',pstring,')')}, 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
            
            % Plot FFT of Quiet peak shape
            [quietfvec, quietfftfreq] = singsidefft(quietpeaks);
            figure(sf2);
            ax2(1)=subplot(R,C,1);
            if strcmp(FFTscale,'log')
                plot(quietfvec, log(quietfftfreq));
                ylabel('log(Magnitude)', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
            else
                plot(quietfvec, quietfftfreq);
                ylabel('Magnitude', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
            end
            grid on;
            xlabel('Frequency (Hz)', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
            title({strcat(datestring, stationstring,'\_','FFT of Deciles of Quiet (~',num2str(round((quietwin(1)/(fs*3600))/IARwin(3))*IARwin(3)),'-',num2str(round((quietwin(2)/(fs*3600))/IARwin(3))*IARwin(3)),' UTC)'); strcat('peak signatures (',nstring,',',pstring,')')}, 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
            
            % Deciles and FFT of Deciles of IAR
            for i = 1:length(IARwinsindex)-1                
                if strcmp(winunique,'unique')
                    if nn == 9
                        IARpeaks = quantile(ptimes{IARwinsindex(i),3}(ptimes{IARwinsindex(i),14},:),[.1 .2 .3 .4 .5 .6 .7 .8 .9]);
                    elseif nn == 10
                        IARpeaks = quantile(ptimes{IARwinsindex(i),5}(ptimes{IARwinsindex(i),14},:),[.1 .2 .3 .4 .5 .6 .7 .8 .9]);
                    else
                        IARpeaks = quantile(ptimes{IARwinsindex(i),nn}(ptimes{IARwinsindex(i),13},:),[.1 .2 .3 .4 .5 .6 .7 .8 .9]);
                    end
                else
                    IARpeaks = quantile(cell2mat(ptimes(IARwinsindex(i),nn)),[.1 .2 .3 .4 .5 .6 .7 .8 .9]);
                end
                figure(sf1);
                ax1(i+1)=subplot(R,C,i+1);
                plot(peaksxvec, IARpeaks);
                hold(ax1(i+1));
                line([peaksxvec(1) peaksxvec(end)],[Brmsquiet(2) Brmsquiet(2)],'Color','k','LineStyle','--');
                line([peaksxvec(1) peaksxvec(end)],[-Brmsquiet(2) -Brmsquiet(2)],'Color','k','LineStyle','--');
                %line([peaksxvec(1) peaksxvec(end)],[noiselevel(IARwinsindex(i)) noiselevel(IARwinsindex(i))],'Color','k','LineStyle','--');
                %line([peaksxvec(1) peaksxvec(end)],[-noiselevel(IARwinsindex(i)) -noiselevel(IARwinsindex(i))],'Color','k','LineStyle','--');
                %line([0 peakwint],[rms(IARpeaks(5,:)) rms(IARpeaks(5,:))],'Color',[0.5 0 0],'LineStyle','--');
                %line([0 peakwint],[-rms(IARpeaks(5,:)) -rms(IARpeaks(5,:))],'Color',[0.5 0 0],'LineStyle','--');
                %line([0.2*peakwint 0.2*peakwint],[min(IARpeaks(1,:))-1 max(IARpeaks(3,:))+1],'Color','k','LineStyle','--');
                hold off;
                grid on;
                xlabel('Time (s.ss)', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
                ylabel('B', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axestitle('Superposed peak signatures', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
                title({strcat(datestring, stationstring,'\_','Deciles of IAR (~',num2str(IARwins(i)/(fs*3600)),'-',num2str(IARwins(i+1)/(fs*3600)),' UTC)'); strcat('peak signatures (',nstring,',',pstring,')')}, 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
                % Plot FFT of IAR peak shape
                [IARfvec, IARfftfreq] = singsidefft(IARpeaks);
                figure(sf2);
                ax2(i+1) = subplot(R,C,i+1);
                if strcmp(FFTscale,'log')
                    plot(IARfvec, log(IARfftfreq));
                    ylabel('log(Magnitude)', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
                else
                    plot(IARfvec, IARfftfreq);
                    ylabel('Magnitude', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
                end
                grid on;
                xlabel('Frequency (Hz)', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
                title({strcat(datestring, stationstring,'\_','FFT of Deciles of IAR (~',num2str(IARwins(i)/(fs*3600)),'-',num2str(IARwins(i+1)/(fs*3600)),' UTC)'); strcat('peak signatures (',nstring,',',pstring,')')}, 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
            end
    end
    
    linkaxes(ax1,'xy');
    linkaxes(ax2,'xy');
    
    
    
    
    % Sums
    % quietpeaks = sum(cell2mat(ptimes(quietwin(1):quietwin(2),3)));
    % IARpeaks1 = sum(cell2mat(ptimes(IARwins(1):IARwins(2),3)));
    % IARpeaks2 = sum(cell2mat(ptimes(IARwins(2):IARwins(3),3)));
    % IARpeaks3 = sum(cell2mat(ptimes(IARwins(3):IARwins(4),3)));
    % IARpeaks4 = sum(cell2mat(ptimes(IARwins(4):IARwins(5),3)));
    
    
    
    % Quartiles of true polarity
    % quietpeaksquart = quantile(cell2mat(ptimes(quietwin(1):quietwin(2),3)),[.25 .50 .75]);
    % IARpeaks1quart = quantile(cell2mat(ptimes(IARwins(1):IARwins(2),3)),[.25 .50 .75]);
    % IARpeaks2quart = quantile(cell2mat(ptimes(IARwins(2):IARwins(3),3)),[.25 .50 .75]);
    % IARpeaks3quart = quantile(cell2mat(ptimes(IARwins(3):IARwins(4),3)),[.25 .50 .75]);
    % IARpeaks4quart = quantile(cell2mat(ptimes(IARwins(4):IARwins(5),3)),[.25 .50 .75]);
    
    % Quartiles of the mean for true polarity
    % quietpeaksquart = quantile(cell2mat(ptimes(quietwin(1):quietwin(2),5)),[.25 .50 .75]);
    % IARpeaks1quart = quantile(cell2mat(ptimes(IARwins(1):IARwins(2),5)),[.25 .50 .75]);
    % IARpeaks2quart = quantile(cell2mat(ptimes(IARwins(2):IARwins(3),5)),[.25 .50 .75]);
    % IARpeaks3quart = quantile(cell2mat(ptimes(IARwins(3):IARwins(4),5)),[.25 .50 .75]);
    % IARpeaks4quart = quantile(cell2mat(ptimes(IARwins(4):IARwins(5),5)),[.25 .50 .75]);
    
    % Quartiles of positive polarity (negative spike windows multiplied by -1)
    % quietpeaksquart = quantile(cell2mat(ptimes(quietwin(1):quietwin(2),4)),[.25 .50 .75]);
    % IARpeaks1quart = quantile(cell2mat(ptimes(IARwins(1):IARwins(2),4)),[.25 .50 .75]);
    % IARpeaks2quart = quantile(cell2mat(ptimes(IARwins(2):IARwins(3),4)),[.25 .50 .75]);
    % IARpeaks3quart = quantile(cell2mat(ptimes(IARwins(3):IARwins(4),4)),[.25 .50 .75]);
    % IARpeaks4quart = quantile(cell2mat(ptimes(IARwins(4):IARwins(5),4)),[.25 .50 .75]);
    
    % Quartiles of the mean of positive polarity (negative spike windows multiplied by -1)
    % quietpeaksquart = quantile(cell2mat(ptimes(quietwin(1):quietwin(2),6)),[.25 .50 .75]);
    % IARpeaks1quart = quantile(cell2mat(ptimes(IARwins(1):IARwins(2),6)),[.25 .50 .75]);
    % IARpeaks2quart = quantile(cell2mat(ptimes(IARwins(2):IARwins(3),6)),[.25 .50 .75]);
    % IARpeaks3quart = quantile(cell2mat(ptimes(IARwins(3):IARwins(4),6)),[.25 .50 .75]);
    % IARpeaks4quart = quantile(cell2mat(ptimes(IARwins(4):IARwins(5),6)),[.25 .50 .75]);
    %
    % Quartiles of normalized positive polarity (negative spike windows multiplied by -1 and normalized to value of peak)
    % quietpeaksquart = quantile(cell2mat(ptimes(quietwin(1):quietwin(2),8)),[.25 .50 .75]);
    % IARpeaks1quart = quantile(cell2mat(ptimes(IARwins(1):IARwins(2),8)),[.25 .50 .75]);
    % IARpeaks2quart = quantile(cell2mat(ptimes(IARwins(2):IARwins(3),8)),[.25 .50 .75]);
    % IARpeaks3quart = quantile(cell2mat(ptimes(IARwins(3):IARwins(4),8)),[.25 .50 .75]);
    % IARpeaks4quart = quantile(cell2mat(ptimes(IARwins(4):IARwins(5),8)),[.25 .50 .75]);
    
    % Quartiles of the mean of normalized positive polarity (negative spike windows multiplied by -1 and normalized to value of peak)
    % quietpeaksquart = quantile(cell2mat(ptimes(quietwin(1):quietwin(2),9)),[.25 .50 .75]);
    % IARpeaks1quart = quantile(cell2mat(ptimes(IARwins(1):IARwins(2),9)),[.25 .50 .75]);
    % IARpeaks2quart = quantile(cell2mat(ptimes(IARwins(2):IARwins(3),9)),[.25 .50 .75]);
    % IARpeaks3quart = quantile(cell2mat(ptimes(IARwins(3):IARwins(4),9)),[.25 .50 .75]);
    % IARpeaks4quart = quantile(cell2mat(ptimes(IARwins(4):IARwins(5),9)),[.25 .50 .75]);
    
    % % Plotting
    % peaksxvec = 0:1/fs:peakwint-1/fs;
    % %allpeaks = sum(cell2mat(ptimes(:,3)));
    % % Plot quartiles
    % lowerlim = min(min([quietpeaksquart IARpeaks1quart IARpeaks2quart IARpeaks3quart IARpeaks4quart],[],2));
    % upperlim = max(max([quietpeaksquart IARpeaks1quart IARpeaks2quart IARpeaks3quart IARpeaks4quart],[],2));
    % f3 = figure('Position', [0 0 s(3) s(4)]);
    % subplot(2,3,1);
    % plot(peaksxvec, quietpeaksquart);
    % ylim([lowerlim upperlim]);
    % xlabel('Time (s.ss)', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
    % ylabel('B', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axestitle('Superposed peak signatures', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
    % title({'Quartiles of the mean of Quiet (~12:00-16:00 UTC)'; 'peak signatures'}, 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
    % subplot(2,3,2);
    % plot(peaksxvec, IARpeaks1quart);
    % ylim([lowerlim upperlim]);
    % xlabel('Time (s.ss)', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
    % ylabel('B', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axestitle('Superposed peak signatures', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
    % title({'Quartiles of the mean of IAR1 (~4:00-5:45 UTC)'; 'peak signatures'}, 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
    % subplot(2,3,3);
    % plot(peaksxvec, IARpeaks2quart);
    % ylim([lowerlim upperlim]);
    % xlabel('Time (s.ss)', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
    % ylabel('B', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axestitle('Superposed peak signatures', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
    % title({'Quartiles of the mean of IAR2 (~5:45-7:30 UTC)'; 'peak signatures'}, 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
    % subplot(2,3,4);
    % plot(peaksxvec, IARpeaks3quart);
    % ylim([lowerlim upperlim]);
    % xlabel('Time (s.ss)', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
    % ylabel('B', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axestitle('Superposed peak signatures', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
    % title({'Quartiles of the mean of IAR3 (~7:30-9:15 UTC)'; 'peak signatures'}, 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
    % subplot(2,3,5);
    % plot(peaksxvec, IARpeaks4quart);
    % ylim([lowerlim upperlim]);
    % xlabel('Time (s.ss)', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
    % ylabel('B', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axestitle('Superposed peak signatures', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
    % title({'Quartiles of the mean of IAR4 (~9:15-11:00 UTC)'; 'peak signatures'}, 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
    
end