%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Fmax  =  100;       % Frequency range for analysis and plots (Hz)
sigmat=0.025; 		% Window size for averaging (sec)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Signal Parameters
RandomSeed=37;		 %RandomSeed>0 resets the random number generator
T = 4.0;             %Length of signal in seconds
dt = 0.001;          %Time step size
upperBandLimit = 2*pi*5; %
lowerBandLimit = 2*pi*0;  %
rms = 0.5;           % RMS signal level
bandwidth = [lowerBandLimit upperBandLimit];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Generate bandlimited white noise
[S,Amps] = genSignal(T,dt,rms,bandwidth,RandomSeed);
Nt = length(S); %
t = dt*[0:Nt-1]; %

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% LIF Neuron Model Parameters
tref = 2;         %
if dt > .0005*tref
    error('ERROR: dt larger than 1/2 Tref'); %
end
tRC  = 20;        %
eps = tRC/tref;   
S0 = 0;           % Firing rate a0 specified at level S0  
a0 = 40;          
S1 = 1.0;         % Firing rate a1 specified at level S1
a1 = 150;

%%% Compute the gain and bias currents
r1 = 1000/(tref*a0);
r2 = 1000/(tref*a1);
f1 = (r1-1)/eps;
f2 = (r2-1)/eps;
gain = (1.0/(exp(f2)-1) - 1.0/(exp(f1)-1))/(S1-S0); %
Sthreshold = S0-1/(gain*(exp(f1)-1));               %
R=1;            %
Vt=1;           %
Jbias = Vt-R*gain*Sthreshold;    %

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
LIFparms = [0.001*tref,0.001*tRC,Vt,Jbias,gain];
tic;
[onOffCnt TonOff] = genOnOFFLIFSpikes(S,dt,LIFparms);  %
runTime = toc;

Non=onOffCnt(1);     %
Noff = onOffCnt(2);  %

fLow = lowerBandLimit/(2*pi);
fHigh = upperBandLimit/(2*pi);
fprintf('Signal RMS=%5.3f,T=%4.1f(sec),dt=%6.4f(sec),BW1=%6.3f (Hz),BW2=%6.3f (Hz)\n',rms,T,dt,fLow,fHigh);
fprintf('tref=%5.1f(ms),tRC=%5.1f(ms),a(0)=%5.1f,a(1)=%5.1f\n',tref,tRC,a0,a1);
fprintf('On Cnt=%d, Off Cnt=%d, runTime=%5.3f\n',Non,Noff,runTime);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% Set up frequency variables for doing the computation
domega = 2*pi/T;  %
omega = [0:domega:2*pi*Fmax];  %
Nomega0 = length(omega);   %
omega = [-fliplr(omega(2:end)),omega];   %
Nt2 = length(omega); %New frequency channel number
dt2 = T/Nt2;         %Time step associated with limited frequency band.
f = omega/(2*pi);   %

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Set up windowing function 
index = find(abs(omega)<3/sigmat); % Restrict the frequency range 
WinOmega = omega(index);
W2 = exp(-WinOmega.^2*sigmat^2); %
W2 = W2/(sum(W2));   %

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Obtain amplitudes of frequency channels for new frequency range.
A = [Amps(end-(Nomega0-2):end),Amps(1:Nomega0)]; % Put zero in the center
PAmps = A.*conj(A); %

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Compute FFT of spike trains 
Ron = sum(exp(-i*TonOff(1,1:Non)'*omega));    %
Roff = sum(exp(-i*TonOff(2,1:Noff)'*omega));  %
R = (Ron-Roff);                               %
PSpikes = R.*conj(R);             % 

%%% Compute the total Spike and Signal powers as a check
Spower  = sum(PSpikes)/Nt2;       % Total Spike Power = (Non+Noff) (For full spectrum!!)
Apower  = (sum(PAmps)/Nt)*(1/Nt); % Total Power/Nt = (rms)^2
fprintf('Non+Noff=%d,Spike Power=%7.3f;Apower=%7.3e,(RMSsig)^2=%7.3e\n',(Non+Noff),Spower,Apower,rms^2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CorPower = A.*conj(R); %

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Average by convolving with the window function W2
WAP    = convn(PAmps,W2,'same');    %
WRP    = convn(PSpikes,W2,'same');  %
WCP    = convn(CorPower,W2,'same'); %

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
OF = WCP./WRP; %
findex = find(real(OF)>.0005);  %
OF2 = OF(findex);   %
f2=f(findex);

figure(1);clf;
plot(f,PAmps,'-','Color',[.8 .8 .8]); %
hold on;
plot(f,((Nt^2*Apower)/(Spower*Nt2))*PSpikes,'k:'); % 

title('???','FontSize',16);
legend('???','???');
xlabel('???');
ylabel('???');

figure(2);clf;
subplot(1,2,1);
plot(f2,real(OF2),'k'); %
xlabel('???');
title('???','FontSize',16);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
FilT = ifft(ifftshift(OF2)); %
index = find(abs(FilT)>0.0005); %
Nf = length(index);
FilT = [FilT(end-Nf:end),FilT(1:Nf)]; %
tf = dt*[-Nf:Nf];
subplot(1,2,2);
plot(tf,real(FilT),'k');
title('???','FontSize',16);
xlabel('???');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Aest = OF.*R; %

DeltaA  = A-Aest; %

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(4);clf;
Sest = (dt/dt2)*real(ifft(ifftshift(Aest))); %
t2 = dt2*[0:length(Sest)-1];
%% Plot results 

plot(t,S,'-','LineWidth',2,'Color',[.8 .8 .8]);
hold on;
plot(t2,Sest,'k:');
y = rms*ones(1,Non);
plot(TonOff(1,1:Non),y,'k.','MarkerSize',1);
y = -rms*ones(2,Noff);
plot(TonOff(2,1:Noff),y,'k.','MarkerSize',1);
title('???','FontSize',16);
legend('???','???','???');
xlabel('???');




