format compact
clear fr avg
clf(figure(1))

%% options
tstep = 4*60*1e6;   % [us] time per frame
windowsize = 8;     % [-]  no. of frames for averaging5

%% settings
nplots = 0;         % [-]  number of plots to be shown while running
fmin   = 13/60;     % [Hz] minimal expected breathing frequency
fmax   = 27/60;     % [Hz] maximal expected breathing frequency


%% preparations
% create editable data
tijd2  = tijd;
signal = envelope;

% calculate frames (count & time)
nframes = floor(tijd(end)/tstep);
tb = (0:nframes) * tstep;

% storage vars
fr  = zeros(1,nframes);
avg = zeros(1,nframes);

%% loop
for i=1:nframes
    clear x t
    plotindex = 1;

    % cut previous frame from data
    signal = signal(tijd2>tb(i));
    tijd2 = tijd2(tijd2>tb(i));
    
    % create frame
    x = signal(tijd2<tb(i+1));
    t = tijd2(tijd2<tb(i+1));
    
    % this is where the magic happens
    X = fft(x);

    % calculate frequency variables
    T0=(t(end)-t(1))/1e6;
    f0=1/T0;
    N=max(size(t));
    
    fs=N/T0;
    Ts=1/fs;
    f=-fs/2:f0:fs/2-f0;

    X2=Ts*[X(floor(N/2)+1:N); X(1:floor(N/2))];

%     subplot(nplots,1,plotindex);plotindex=plotindex+1;
%     plot(t/36e8,x*5/1024)
%     xlabel('tijd [h]')
%     ylabel('audio [V]')

%     subplot(nplots,1,plotindex);plotindex=plotindex+1;
%     plot(f,abs(X2))
%     axis([0 50 0 3])
%     xlabel('frequentie [Hz]')
%     ylabel('amplitude')
% end

    % filter frequencies
    X3 = X2(f < fmax);
    f2 = f(f < fmax);
    X3 = X3(f2 > fmin);
    f3 = f2(f2 > fmin);
    
%     subplot(nplots,1,plotindex);plotindex=plotindex+1;
%     plot(f3,abs(X3));
%     xlabel('frequentie [Hz]')
%     ylabel('amplitude')
    
    % assume the highest peak is actual breathing frequency
    fr(i) = f3(X3==max(X3));
    
    % calculate shifting average
    if i < windowsize
        avg(i) = mean(fr(1:i));
    else
        avg(i) = mean(fr(i+1-windowsize:i));
    end
    
%     subplot(nplots,1,plotindex);plotindex=plotindex+1;
%     plot(fr(1:i)*60, 'LineWidth',0.5);hold on
%     plot(avg(1:i)*60,'LineWidth',2);hold off
%     xlabel(['frame [-]'])
%     ylabel('frequentie [/min]')
    
%     pause(0.1)
end

T0 = (tijd(end)-tijd(1))/1e6;
if T0 < 600
    tc = 1e6;
    label = 'tijd [s]';
elseif T0 < 3 * 3600
    tc = 60*1e6;
    label = 'tijd [min]';
else
    tc = 3600*1e6;
    label = 'tijd [h]';
end
tframe = (1:tijd(end)/tstep)*tstep/tc;
clf();
plot(tframe, fr*60, 'LineWidth',0.5);hold on
plot(tframe, avg*60,'LineWidth',2,'Color','r');hold on
xlabel(label)
ylabel('ademhalingsfrequentie [/min]')