function [h, hl, freq_windowed] = chs4p3_sub(mode, data, timeseq, rms_time, fig, f, p1)

% subroutine called by CHS4P3, returns two handles
% Plot Graphs to accompany notes on initial Sounder tests

L = length(timeseq);
fs = 500;  % sample freq
X = linspace(0,fs,L+1);
X = X(1:L);

if mode==4, 
    L2 = 276;
    window = [ zeros((L-L2)/2,1); blackman(L2); zeros((L-L2)/2,1) ];
    % this gives some strange answers - needs investigating
    % The idea was that the window would simply be a 3kHz(ish) filter on 
    % each point, but it doesnt seem to work loke that,
    % there are strange dips in the response that shouldnt be there 
else
    window = cosinewindow(L, 'sym', [0.35875 0.48829 0.14128 0.01174]);
    %window = hamming(L);
end
coherent_gain = sum(window)/L;   % we scale the freq spectrum by this
NEB = L * sum(window.^2) ./ (sum(window)).^2

freqseq = abs( fft(timeseq) ) / (L/2) ;
time_windowed = timeseq .* window;
freq_windowed = abs(fft(time_windowed)/(L/2));
freq_windowed = freq_windowed / coherent_gain;

%B = 2048; % width of bar in points
% for latest graph, set B = 256, or 977Hz
B = 64;
K = L/B;  % number of bars in full spectrum

freq_energy = freq_windowed;
meansq = sum (reshape(freq_energy.^2,B,K)) / 2;
% now we scale this by the NEB because otherwise we over-estimate
meansq = meansq/NEB;
sum_meansq = sum(meansq)/2  % allows for reflection in fs/2

signal_meansq = rms_time.^2
ratio_of_above = sum_meansq/signal_meansq
meansq = meansq/sum_meansq;    % remove ';' to print values of histogram bars
%floor(meansq(1:K/2)*100 + 0.5)
widthX = fs/K;  % the width of the bars in kHz
barX = linspace(fs/K/2, (fs - fs/K)/2, K/2);


% get this from calling routine: f = [195 201]                                     
% work out signal between these two points
s = floor(f * L/fs);
range = s(2) + 1 - s(1);
meansq_over_specified_range = sum (freq_energy([s(1):s(2)]).^2) / 2 /NEB /sum_meansq

% ============== now plot frequency data

figure(fig);

h(1) = SUBPLOT(2,1,1);
hl(1) = semilogy(X,freq_windowed, p1, 'linewidth',1.5 );
xlabel('frequency (kHz)','fontsize',12); ylabel('received signal (mV)','fontsize',12)
title(['data: ', data],'fontsize',12)

h(2) = SUBPLOT(2,1,2);

if mode==1,
    bar(barX,meansq(1:K/2), 1, 'c')
    xlabel('frequency (kHz)','fontsize',12); ylabel('proportion of energy','fontsize',12)
elseif mode==2,
    hl(2) = semilogy(X,freqseq, p1, 'linewidth',1.5 );
    xlabel('frequency (kHz)','fontsize',12); ylabel('received signal (mV)','fontsize',12)
else
    hl(2) = semilogy(X,freq_windowed, p1, 'linewidth',1.5 );
    xlabel('frequency (kHz)','fontsize',12); ylabel('received signal (mV)','fontsize',12)
end

set(h, 'LineWidth',1.5,'FontSize',12,'TickDir','out')
set(h, 'PlotBoxAspectRatio', [3 1 1])
%set(h,'YScale', 'lin')
