% ### EXturbine.m ###    9/13/13
% Code developed by Chris Ni to explore 'overlap' projection from two
% nested oscillators (http://www.yorku.ca/cberge/images/Ni2013CUPCsmall.pdf)

clear
% -------
% User parameters
tmax= 20;   % time          [s]
SR= 80;    % sample rate   [Hz]
f1= 0.2;       % sinusoid freq [Hz]
f2= 0.3;   % sinusoid freq [Hz]
p1= 0;    % phase angle   [rad]
p2= 0;    % phase angle   [rad]
D1= 0.5;    % angular width [rad]
D2= 0.7;    % angular width [rad]
refP= pi;   % reference phase for determining I
SNR= 0.0;  % signal-to-noise factor
trials= 90; % # of trials to repeat noise runs
maxNoise= 0.1; % [0.04]
linewidth = 3;               % adjustable plotting parameter
% -------
% allocate more data points for autocorrelation
tmax= 2*tmax;
N= tmax*SR;
t= linspace(0,tmax,N);  % digitize time
storeI= [];
indx=1;
for mm=1:trials    
    % create 'frozen' noise
    noise1= SNR*randn(numel(t),1)'; %position noise
    noise2= SNR*randn(numel(t),1)';
    % create time waveforms for each phase oscillator
    % [kludge? not sure if two variables are needed for each here]
    % blade 1
    randW= (1+((mm-1)/trials)*maxNoise*randn(1,1)); %up to 1+1*max
    w1(mm)= randW*2*pi*f1;
    a1b = mod((w1(mm))*t + p1 + D1/2,2*pi)+ noise1; % leading edge
    a1e = mod(a1b-D1+p1,2*pi)+ noise1;                % trailng edge
    % blade 2
    a2b = mod((2*pi*f2)*t + p2 + D2/2,2*pi)+ noise2; % leading edge
    a2e = mod(a2b-D2+p2,2*pi)+ noise2;                % trailng edge
    % determine 'intensity' output: I=0 if POV is not blocked, I=1 if it is
    % blocked
    for n=1:N
        if a1b(n) >= refP && a1e(n) <= refP || a2b(n) >= refP && a2e(n) <= refP
            I(n) = 0;   % blocked
            %tP(indx,:)= [mm t(n) 1] ;
            %indx=indx+1;
        else
            I(n) = 1;   % not blocked
        end
    end
    storeI(mm,:)= I;
end

% +++
% develop raster plot to visualize
A = storeI;                  % a convenient synonym
% invert color scheme (switch 0s and 1s) to make it easier...
A = abs(A-1);
[nrows,ncols] = size(A);
x = 0:(ncols-1);
figure(1); clf;
for m=1:nrows
  y = A(m,:);
  dy = [0,diff(y)];
  ups = find(dy>0);
  downs = find(dy<0);
  % deal with the ends...
  if (ups(1)>downs(1))
    ups = [1,ups];
  end
  if (ups(end)>downs(end))
    downs = [downs,ncols];
  end
  if (numel(ups) ~= numel(downs))
    warning('rethink this!');
  end
  
  for n=1:numel(ups)
    x1 = x(ups(n));
    x2 = x(downs(n));
    plot([x1,x2]/SR,[m,m],'k-','LineWidth',linewidth); hold on
  end
end
xlim([x(1),x(end)/SR]);
ylim([1,nrows]);
title('Noise Representation')
xlabel('Time [s]')
ylabel('Trial #')

% +++
figure(2); clf;
plot(linspace(1,trials,trials),w1,'ko-')
hold on; grid on;
xlabel('Trial #')
ylabel('w1')
plot(linspace(1,trials,trials),w1(1)+linspace(0,1,trials)*maxNoise,'r-')
plot(linspace(1,trials,trials),w1(1)-linspace(0,1,trials)*maxNoise,'r-')
legend('w1 variation','STD bound','Location','NorthWest')


return
% -----------------------------------
% take half of total segment to (auto)correlate with full segment
% [so to maintain a baseline amplitude]
for n=1:1+N/2
    s(n) = I(n);
end

% visualize
fHand= figure(1); clf;
% graph model
subplot(311);
grid on;hold on;
plot(t(1:N/2),a1b(1:N/2),'.-','Color',[0 0 1])
plot(t(1:N/2),a1e(1:N/2),'.-','Color',[0 0 0.5])
plot(t(1:N/2),a2b(1:N/2),'.-','Color',[1 0 0])
plot(t(1:N/2),a2e(1:N/2),'.-','Color',[0.5 0 0])
plot(t(1:N/2),refP,'g-')
xlabel('Time [s]')
ylabel('Phase [rad]')
axis([0 tmax/2 0 7])
% intensity output
subplot(312); plot(t(1:N/2),I(1:N/2),'rd-')
grid on;hold on;
title('Model Output')
xlabel('Time [s]')
ylabel('Intensity')

% ----------------
% determine if periodicity exist based upon autocorrelation
% slide function segment
for n=1:N/2 
    % do the convolution
    c(n)= sum(s.*I(n:n+numel(s)-1));
    % normalize re initial meas.
    if n==1
        ref= c(n);
        c(n)= c(n)/c(1);
    else
        c(n)= c(n)/ref;
    end
    
end

% intensity output autocorrelation
figure(fHand)
subplot(313);
plot(t(1:N/2),c,'b.-')
grid on;hold on;
title('Autocorrelation')
xlabel('time [s]')
ylabel('ACF')

% kludge way to visualize periodicity
indxP= find(c>=0.95);
plot(t(indxP),c(indxP),'rd','LineWidth',3)