clear all;
close all;
clc;
%
%% Script of the estimation of power spectrums of Jansen' signals
%
%% Description
%
% In that script the power specrtum is estimated for every repetition of
% Jansen' signal and then these spectrums are averaged into one and 
% it is plotted in figure.
%
%% Jansen signals from the sponge experiments with their begginings
%
% sponge_C1_S1_Jansen1_2015.06.16_11.11.27.hdf5 = [16300, 73490, 130800, 187900, 245100];
% sponge_C1_S2_Jansen1_2015.06.16_11.38.46.hdf5 = [21130, 77780, 135100, 193200, 250300];
% sponge_C1_S1_Jansen2_2015.06.16_11.12.52.hdf5 = [27460, 84460, 141300, 198700, 255500];
% sponge_C1_S2_Jansen2_2015.06.16_11.40.12.hdf5 = [36710, 93580, 150500, 209400, 266300];
%
%% Loading of Jansen's signals
%
infoorigsignal = h5info('sponge_C1_S2_Jansen2_2015.06.16_11.40.12.hdf5');
origin = h5read('sponge_C1_S2_Jansen2_2015.06.16_11.40.12.hdf5', '/RawData/Samples');
%
%% Setting of the sampling frequency
%
fs = 4800;
%
%% Plotting of the Jansen's system, 2 inputs and 2 outputs
%
figure
plot(origin(1,:),'r')
hold on
plot(origin(2,:),'b')
hold on
plot(origin(3,:),'g')
hold on
plot(origin(4,:),'y')
%
%% Definition of the vector of begginings of each repetition of Jansen's signal
%
start = [36710, 93580, 150500, 209400, 266300];
% 
%% Definition of the length of each repetition of Jansen's signal
%
difference = 47110;
%
%% For each repetition of Jansen' signal the power spectrum is estimated
%
for l = 1:4;
    for k = 1:length(start)
        ending(k) = start(k) + difference;
        signal{l,k} = double(origin(l,start(k):ending(k)));
        N = length(signal{l,k});
        ws = 2*pi/N;
        wnorm = -pi:ws:pi;
        wnorm = wnorm(1:N);
        w = wnorm*fs;
        win = hamming(N);
        xw = win(:).*signal{l,k}(:);
        Xw{l,k}=fft(xw);
        ampl{l,k}=abs(fftshift(Xw{l,k}));       
    end
end
%
%% These power spectrums from repetitions are everaged, 
% so on power spectrum is obtained for the input 1, 
% for input 2 and for outpu 1 and for output 2
for q = 1:4
  for r = 1:difference;
    for y = 1:5
        average(y,:) = ampl{q,y};
    end
  final{q}(1,r) = mean(average(:,r));
  end
end
