MuseOSC_EEG.cpp

/*
 * MuseOSC_EEG.cpp
 *
 *  Created on: 7.7.2015
 *      Author: Tomas
 */


#include "MuseOSC_EEG.h"
#include <math.h>
#include <unistd.h>
#include <chrono>

#ifdef __linux__
#include <pthread.h>
#include <sys/resource.h>
#endif

#ifdef WINOS
#include <windows.h>
#endif

#include <dinrhiw.h>

#include "spectral_entropy.h"
#include "spectral_analysis.h"
#include "spectrum_phase_difference.h"

#include "oscpkt.hh"
#include "udp.hh"

using namespace oscpkt;
using namespace std::chrono;

namespace whiteice {
namespace resonanz {

MuseOSC_EEG::MuseOSC_EEG(const unsigned int portNum,
			 const float analyze_sampling_hz,
			 const bool useICA) : 
  port(portNum)
{
  worker_thread = nullptr;
  analyze_thread = nullptr;
  
  running = false;
  hasConnection = false;
  quality = 0.0f;
  
  value.resize(this->getNumberOfSignals());
  for(auto& v : value) v = 0.0f;

  latest_sample_seen_t = 0LL;

  if(analyze_sampling_hz >= 5.0f)
    this->default_analyze_sampling_hz = analyze_sampling_hz;

  extra_value.resize(default_analyze_sampling_hz);
  for(auto& v : extra_value) v = 0.0f;
  
  ica_extra_value.resize(default_analyze_sampling_hz);
  for(auto& v : ica_extra_value) v = 0.0f;
  
  this->calculate_ica_components = useICA;
  
  try{
    std::unique_lock<std::mutex> lock(connection_mutex);
    running = true;
    
    worker_thread = new std::thread(&MuseOSC_EEG::muse_loop, this);
    
    analyze_thread = new std::thread(&MuseOSC_EEG::sample_and_analyze_eeg_loop, this);

    // waits for up to 5 seconds for connection to be established
    // or silently fails (device has no connection when accessed)

    std::chrono::system_clock::time_point now =
      std::chrono::system_clock::now();
    auto end_time = now + std::chrono::seconds(5);

    while(hasConnection == true || std::chrono::system_clock::now() >= end_time)
      connection_cond.wait_until(lock, end_time);
  }
  catch(std::exception){
    running = false;
    throw std::runtime_error("MuseOSC_EEG: couldn't create worker thread.");
  }
}


MuseOSC_EEG::~MuseOSC_EEG()
{
  running = false;
  
  if(worker_thread != nullptr){
    worker_thread->join();
    delete worker_thread;
  }

  if(analyze_thread != nullptr){
    analyze_thread->join();
    delete analyze_thread;
  }
  
  worker_thread = nullptr;
  analyze_thread = nullptr;
}

/*
 * Returns unique DataSource name
 */
std::string MuseOSC_EEG::getDataSourceName() const
{
  return "Interaxon Muse 2/S (Mind Monitor App)";
}

/**
 * Returns true if connection and data collection to device is currently working.
 */
bool MuseOSC_EEG::connectionOk() const
{
  if(hasConnection == false)
    return false; // there is no connection
  
  long long ms_since_epoch = (long long)duration_cast<milliseconds>(system_clock::now().time_since_epoch()).count();
  
  if((ms_since_epoch - latest_sample_seen_t) > 2000)
    return false; // latest sample is more than 2000ms second old => bad connection/data
  else
    return true; // sample within 2000ms => good
}

/**
 * returns current value
 */
bool MuseOSC_EEG::data(std::vector<float>& x) const
{
  std::lock_guard<std::mutex> lock(data_mutex);
  
  if(this->connectionOk() == false || value.size() != this->getNumberOfSignals()){

    if(x.size() != this->getNumberOfSignals()){  
      x.resize(this->getNumberOfSignals());
      
      for(unsigned int i=0;i<x.size();i++)
	x[i] = 0.0f;
    }

    return false;
  }
  
  x = value;
  
  return true;
}

bool MuseOSC_EEG::getSignalNames(std::vector<std::string>& names) const
{
  names.resize(6+3+3);

  // delta, theta, alpha, beta, gamma
  
  names[0] = "Muse: Delta";
  names[1] = "Muse: Theta";
  names[2] = "Muse: Alpha";
  names[3] = "Muse: Beta";
  names[4] = "Muse: Gamma";
  names[5] = "Muse: Total Power";
  
  names[6] = "Muse: Spectral Entropy";
  names[7] = "Muse: Phase Difference";
  names[8] = "Muse: Phase Slope Index";

  names[9]  = "Muse: Spectral Entropy (ICA signals)";
  names[10] = "Muse: Phase Difference (ICA signals)";
  names[11] = "Muse: Phase Slope Index (ICA signals)";
  
  return true;
}

unsigned int MuseOSC_EEG::getNumberOfSignals() const
{
  return 6+3+3;/*+1 no phase difference for now*/ // entropy is calculated signal
}

bool MuseOSC_EEG::extradata(std::vector<float>& x) const
{
  std::lock_guard<std::mutex> lock(extra_data_mutex);
  
  if(this->connectionOk() == false || extra_value.size() != this->getExtraNumberOfSignals()){

    if(x.size() != this->getExtraNumberOfSignals()){
      x.resize(this->getExtraNumberOfSignals());
      
      for(auto& xi : x) xi = 0.0;
    }

    return false;
  }

  x.resize(this->getExtraNumberOfSignals());

  for(auto& xi : x) xi = 0.0;

  for(unsigned int i=0;i<extra_value.size()&&i<x.size();i++)
    x[i] = extra_value[i];

  for(unsigned int i=0;i<ica_extra_value.size()&&(i+extra_value.size())<x.size();i++)
    x[i+extra_value.size()] = ica_extra_value[i];
  
  return true;
}

unsigned int MuseOSC_EEG::getExtraNumberOfSignals() const
{
  // in case of 100 Hz sampling frequency => two 100 FFT power values
  return 2*(this->default_analyze_sampling_hz);
}

void MuseOSC_EEG::muse_loop() // worker thread loop
{
  // sets MuseOSC internal thread high priority thread
  // so that connection doesn't timeout
  if(calculate_ica_components == true){
    
    if(getpriority(PRIO_PGRP, 0) > 0) setpriority(PRIO_PGRP, 0, 0);
    if(getpriority(PRIO_PGRP, 0) > -19) setpriority(PRIO_PGRP, 0, -19);
    
    sched_param sch_params;
    int policy = SCHED_OTHER;
    
    pthread_getschedparam(pthread_self(),
			  &policy, &sch_params);
    
    policy = SCHED_OTHER; // DEADLINE; // FIFO, OTHER
    sch_params.sched_priority = 48; // sched_get_priority_min(policy);
	
    if(pthread_setschedparam(pthread_self(),
			     policy, &sch_params) != 0)
      logging.error("MuseOSC_EEG::muse_loop(): setting REAL-TIME priority FAILED");
    else
      logging.info("MuseOSC_EEG::muse_loop(): setting REAL-TIME priority DONE");
    
#ifdef WINOS
    logging.info("windows os: MuseOSC_EEG::muse_loop(): setting thread high priority");
    SetPriorityClass(GetCurrentProcess(), ABOVE_NORMAL_PRIORITY_CLASS);
    SetPriorityClass(GetCurrentProcess(), HIGH_PRIORITY_CLASS);
    SetPriorityClass(GetCurrentProcess(), REALTIME_PRIORITY_CLASS);
#endif
  }
  
  hasConnection = false;
  
  UdpSocket sock;
  
  while(running){
    sock.bindTo(port);
    if(sock.isOk()) break;
    sock.close();
    sleep(1);
  }
  
  PacketReader pr;
  PacketWriter pw;
  
  std::vector<int> connectionQuality;
  
  float delta = 0.0f, theta = 0.0f, alpha = 0.0f, beta = 0.0f, gamma = 0.0f;
  bool hasNewData = false;
  
  while(running){
    
    if(hasConnection && hasNewData){ // updates data
      float q = 0.0f;
      for(auto qi : connectionQuality)
	if(qi > 0) q++;
      q = q / connectionQuality.size();
      
      quality = q;
      
      std::vector<float> v;
      
      
      // converts absolute power (logarithmic bels) to [0,1] value by saturating
      // values using tanh(t) this limits effective range of the values to [-0.1, 1.2]
      // TODO: calculate statistics of delta, theta, alpha, beta, gamma to optimally saturate values..
      
      auto t = delta; t = (1 + tanh(2*(t - 0.6)))/2.0;
      v.push_back(t);
      
      t = theta; t = (1 + tanh(2*(t - 0.6)))/2.0;
      v.push_back(t);
      
      t = alpha; t = (1 + tanh(2*(t - 0.6)))/2.0;
      v.push_back(t);
      
      t = beta; t = (1 + tanh(2*(t - 0.6)))/2.0;
      v.push_back(t);
      
      t = gamma; t = (1 + tanh(2*(t - 0.6)))/2.0;
      v.push_back(t);
      
      // calculates total power in decibels [sums power terms together]
      float total = 0.0f;

      {
	if(delta/10.0f < 1e6f) total += pow(10.0f, delta/10.0f);
	else total += total += pow(10.0f, 1e6f/10.0f);

	if(theta/10.0f < 1e6f) total += pow(10.0f, theta/10.0f);
	else total += total += pow(10.0f, 1e6f/10.0f);
	
	if(alpha/10.0f < 1e6f) total += pow(10.0f, alpha/10.0f);
	else total += total += pow(10.0f, 1e6f/10.0f);

	if(beta/10.0f < 1e6f) total += pow(10.0f, beta/10.0f);
	else total += total += pow(10.0f, 1e6f/10.0f);
	
	if(gamma/10.0f < 1e6f) total += pow(10.0f, gamma/10.0f);
	else total += total += pow(10.0f, 1e6f/10.0f);

	if(isnan(total) != 0 || isinf(total) != 0){
	  total = 0.50f; // bad data: uses dummy value
	}
	else if(total < 0.0f){ // cannot calculate negative logarithm.. => dummy value
	  total = 0.50f;
	}
	else if(total > 1e30f){
	  total = 10.0f * +30.0f;
	}
	else if(total < 1e-30f){
	  total = 10.0f * -30.0f;
	}
	else{
	  total = 10.0f * log10(total);
	}
      }

      t = total; t = (1 + tanh(2*(t - 7.0)))/2.0;
      v.push_back(t);

      {
	std::lock_guard<std::mutex> lock(extra_values_mutex);
	
	// adds spectral entropy
	v.push_back(global_spectral_entropy);
	
	// adds phase difference
	v.push_back(global_phase_difference);
	
	// adds phase slope index
	v.push_back(global_phase_slope_index);
	
	// adds ICA spectral entropy
	v.push_back(global_ica_spectral_entropy);
	
	// adds ICA phase difference
	v.push_back(global_ica_phase_difference);
	
	// adds ICA phase slope index
	v.push_back(global_ica_phase_slope_index);
      }
      
      // gets current time
      auto ms_since_epoch = duration_cast<milliseconds>(system_clock::now().time_since_epoch()).count();
      
      std::lock_guard<std::mutex> lock(data_mutex);
      value = v;
      latest_sample_seen_t = (long long)ms_since_epoch;
      
      // printf("EEQ: D:%.2f T:%.2f A:%.2f B:%.2f G:%.2f [QUALITY: %.2f]\n", delta, theta, alpha, beta, gamma, q);
      // printf("EEQ POWER: %.2f [QUALITY %.2f]\n", log(exp(delta)+exp(theta)+exp(alpha)+exp(beta)+exp(gamma)), q);
      // fflush(stdout);
      
      hasNewData = false; // this data point has been processed
    }
    
    if(sock.receiveNextPacket(10)){
      
      pr.init(sock.packetData(), sock.packetSize());
      Message* msg = NULL;
      
      while(pr.isOk() && ((msg = pr.popMessage()) != 0)){
	Message::ArgReader r = msg->match("/muse/elements/is_good");
	
	if(r.isOk() == true){ // matched
	  // there are 4 ints telling connection quality
	  std::vector<int> quality;
	  
	  while(r.nbArgRemaining()){
	    if(r.isInt32()){
	      int32_t i;
	      r = r.popInt32(i);
	      quality.push_back((int)i);
	    }
	    else{
	      r = r.pop();
	    }
	  }
	  
	  if(quality.size() > 0){
	    connectionQuality = quality;
	  }
	  
	  bool connection = false;

	  if(quality.size() == 4){
	    unsigned int counter = 0;
	    
	    for(auto q : quality)
	      if(q > 0) counter++;

	    if(counter >= 4) connection = true; // MUST HAVE FULL CONNECTION TO ALL EEG CHANNELS 
	  }

	  {
	    std::lock_guard<std::mutex> lock(connection_mutex);
	    hasConnection = connection;
	    connection_cond.notify_all();
	  }
	}
	else{

	  r = msg->match("/muse/elements/horseshoe");

	  if(r.isOk()){

	    std::vector<float> quality;

	    while(r.nbArgRemaining()){
	      if(r.isFloat()){
		float f = 0.0f;
		r = r.popFloat(f);
		quality.push_back(f);
	      }
	      else{
		r = r.pop();
	      }
	    }

	    if(quality.size() > 0){
	      
	      std::vector<int> q;
	      for(auto qi : quality)
		if(qi == 1.00f) q.push_back(1);
		else q.push_back(0);
	      
	      connectionQuality = q;
	    }

	    
	    bool connection = false;

	    if(quality.size() == 4){
	      unsigned int counter = 0;
	      
	      for(auto q : quality)
		if(q > 0) counter++;
	      
	      if(counter >= 4) connection = true; // MUST HAVE FULL CONNECTION TO ALL EEG CHANNELS 
	    }

	    
	    {
	      std::lock_guard<std::mutex> lock(connection_mutex);
	      hasConnection = connection;
	      connection_cond.notify_all();
	    }
	    
	  }
	  
	}

	

	// gets relative frequency bands powers..
	std::vector<float> f(connectionQuality.size());
	
	r = msg->match("/muse/elements/delta_absolute");
	if(r.isOk() && r.nbArgRemaining() == connectionQuality.size()){

	  for(unsigned int i=0;i<f.size();i++)
	    r = r.popFloat(f[i]);

	  if(r.isOkNoMoreArgs()){
	    std::vector<float> samples;
	    for(unsigned int i=0;i<f.size();i++)
	      if(connectionQuality[i]) samples.push_back(f[i]);
	    
	    float mean = 0.0f;
	    for(auto& si :samples) mean += pow(10.0f, si/10.0f);
	    if(samples.size()) mean /= samples.size();

	    if(mean > 0.0f){
	      delta = 10.0f * log10(mean);
	    }
	    else{
	      delta = 10.0f * (-100.0f);
	    }
	    
	    hasNewData = true;
	  }
	  
	}
	
	r = msg->match("/muse/elements/theta_absolute");
	if(r.isOk() && r.nbArgRemaining() == connectionQuality.size()){
	  
	  for(unsigned int i=0;i<f.size();i++)
	    r = r.popFloat(f[i]);

	  if(r.isOkNoMoreArgs()){
	    std::vector<float> samples;
	    for(unsigned int i=0;i<f.size();i++)
	      if(connectionQuality[i]) samples.push_back(f[i]);
	    
	    float mean = 0.0f;
	    for(auto& si :samples) mean += pow(10.0f, si/10.0f);
	    if(samples.size()) mean /= samples.size();

	    if(mean > 0.0f){
	      theta = 10.0f * log10(mean);
	    }
	    else{
	      theta = 10.0f * (-100.0f);
	    }
	    
	    hasNewData = true;
	  }
	}
	
	r = msg->match("/muse/elements/alpha_absolute");
	if(r.isOk() && r.nbArgRemaining() == connectionQuality.size()){
	  
	  for(unsigned int i=0;i<f.size();i++)
	    r = r.popFloat(f[i]);

	  if(r.isOkNoMoreArgs()){
	    std::vector<float> samples;
	    for(unsigned int i=0;i<f.size();i++)
	      if(connectionQuality[i]) samples.push_back(f[i]);
	    
	    float mean = 0.0f;
	    for(auto& si :samples) mean += pow(10.0f, si/10.0f);
	    if(samples.size()) mean /= samples.size();

	    if(mean > 0.0f){
	      alpha = 10.0f * log10(mean);
	    }
	    else{
	      alpha = 10.0f * (-100.0f);
	    }
	    
	    hasNewData = true;
	  }
	}
	
	r = msg->match("/muse/elements/beta_absolute");
	if(r.isOk() && r.nbArgRemaining() == connectionQuality.size()){

	  for(unsigned int i=0;i<f.size();i++)
	    r = r.popFloat(f[i]);
	  
	  if(r.isOkNoMoreArgs()){
	    std::vector<float> samples;
	    for(unsigned int i=0;i<f.size();i++)
	      if(connectionQuality[i]) samples.push_back(f[i]);
	    
	    float mean = 0.0f;
	    for(auto& si :samples) mean += pow(10.0f, si/10.0f);
	    if(samples.size()) mean /= samples.size();

	    if(mean > 0.0f){
	      beta = 10.0f * log10(mean);
	    }
	    else{
	      beta = 10.0f * (-100.0f);
	    }
	    
	    hasNewData = true;
	  }
	}
	
	r = msg->match("/muse/elements/gamma_absolute");
	if(r.isOk() && r.nbArgRemaining() == connectionQuality.size()){

	  for(unsigned int i=0;i<f.size();i++)
	    r = r.popFloat(f[i]);
	  
	  if(r.isOkNoMoreArgs()){
	    std::vector<float> samples;
	    for(unsigned int i=0;i<f.size();i++)
	      if(connectionQuality[i]) samples.push_back(f[i]);
	    
	    float mean = 0.0f;
	    for(auto& si :samples) mean += pow(10.0f, si/10.0f);
	    if(samples.size()) mean /= samples.size();

	    if(mean > 0.0f){
	      gamma = 10.0f * log10(mean);
	    }
	    else{
	      gamma = 10.0f * (-100.0f);
	    }
	    
	    hasNewData = true;
	  }
	}

	r = msg->match("/muse/eeg"); // Mind Monitor App EEG-values 

	if(r.isOk() && r.nbArgRemaining() >= 4){ // typicalle there should be 6 values.. 

	  // this assumes first 4 values are the REAL EEG values in MindMonitor App
	  
	  std::vector<float> f(4);

	  if(f.size() > connectionQuality.size()) f.resize(connectionQuality.size()); 

	  for(unsigned int i=0;i<f.size();i++){
	    r = r.popFloat(f[i]);
	    if(!connectionQuality[i]) f[i] = 0.0; // sets bad data to static zero
	  }

	  {
	    std::lock_guard<std::mutex> lock(eeg_mutex);

	    if(f.size() >= 1) eegTP9  = f[0];
	    if(f.size() >= 2) eegAF7  = f[1];
	    if(f.size() >= 3) eegAF8  = f[2];
	    if(f.size() >= 4) eegTP10 = f[3];
	  }
	}

      }
    }
  
    if(sock.isOk() == false){
      // tries to reconnect the socket to port
      sock.close();
      sleep(1);
      sock.bindTo(port);
    }
    
  }
  
  sock.close();
}

void MuseOSC_EEG::sample_and_analyze_eeg_loop() // worker thread loop
{
  // sets MuseOSC internal thread high priority thread
  // so that calculations are in sync
  if(calculate_ica_components == true){
    if(getpriority(PRIO_PGRP, 0) > 0) setpriority(PRIO_PGRP, 0, 0);
    if(getpriority(PRIO_PGRP, 0) > -19) setpriority(PRIO_PGRP, 0, -19);
    
    sched_param sch_params;
    int policy = SCHED_OTHER;
    
    pthread_getschedparam(pthread_self(),
			  &policy, &sch_params);
    
    policy = SCHED_OTHER; // DEADLINE; // FIFO, OTHER
    sch_params.sched_priority = 40; // sched_get_priority_min(policy);
	
    if(pthread_setschedparam(pthread_self(),
			     policy, &sch_params) != 0)
      logging.error("MuseOSC_EEG::sample_and_analyze_eeg_loop(): setting REAL-TIME priority FAILED");
    else
      logging.info("MuseOSC_EEG::sample_and_analyze_eeg_loop(): setting REAL-TIME priority DONE");
    
#ifdef WINOS
    logging.info("windows os: MuseOSC_EEG::sample_and_analyze_eeg_loop(): setting thread high priority");
    SetPriorityClass(GetCurrentProcess(), ABOVE_NORMAL_PRIORITY_CLASS);
    SetPriorityClass(GetCurrentProcess(), HIGH_PRIORITY_CLASS);
    SetPriorityClass(GetCurrentProcess(), REALTIME_PRIORITY_CLASS);
#endif
  }  

  
  float spectral_entropy  = 0.50f;
  float phase_difference  = 1.0f;
  float phase_slope_index = 0.0f;

  float ica_spectral_entropy = 0.50f;
  float ica_phase_difference = 1.0f;
  float ica_phase_slope_index = 0.0f;
  
  std::vector< std::vector<float> > EEG;
  std::vector< std::vector<double> > EEGd;

  std::vector< std::vector<float> > IC; // ICA signals of the EEG signals
  std::vector< std::vector<double> > ICd;

  EEG.resize(4); // 4 EEG channels in device
  EEGd.resize(4); // 4 EEG channels in device
  IC.resize(4);
  ICd.resize(4);

  // was: 30 Hz sampling frequency!! (NOTE: was 30 Hz, check it works..)
  // [100 Hz is TOO SLOW WITH ICA signals, try: 50 Hz]
  // (100 Hz should be enough but not good for normal laptops)
  float SAMPLING_HZ = 50.0;

  if(this->default_analyze_sampling_hz >= 5.0f)
    SAMPLING_HZ = this->default_analyze_sampling_hz;
  
  const unsigned int EEG_LENGTH = (unsigned int)SAMPLING_HZ;

  {
    char buf[80];

    sprintf(buf, "MuseOSC_EEG::sample_and_analyze_eeg_loop(): sampling frequency: %f Hz",
	    SAMPLING_HZ);
    
    logging.info(buf);
  }
  
  for(unsigned int i=0;i<IC.size();i++){
    IC[i].resize(EEG_LENGTH);
    ICd[i].resize(EEG_LENGTH);
  }

  // timing
  auto start = std::chrono::high_resolution_clock::now();

  while(running){

    {
      std::lock_guard<std::mutex> lock(extra_values_mutex);

      global_spectral_entropy  = spectral_entropy;
      global_phase_difference  = phase_difference;
      global_phase_slope_index = phase_slope_index;

      global_ica_spectral_entropy  = ica_spectral_entropy;
      global_ica_phase_difference  = ica_phase_difference;
      global_ica_phase_slope_index = ica_phase_slope_index;
    }
    
    auto elapsed = std::chrono::high_resolution_clock::now() - start;
    
    const long long microseconds_elapsed =
      std::chrono::duration_cast<std::chrono::microseconds>(elapsed).count();

    const long long microseconds_sleep_time = 1000000/SAMPLING_HZ;
    
    if(microseconds_elapsed > microseconds_sleep_time){
      if(calculate_ica_components == true)
	whiteice::logging.warn("MuseOSC_EEG::sample_and_analyze_eeg_loop(): Warning sampling out of sync!");
    }
    else{
      std::this_thread::sleep_for
	(std::chrono::microseconds(microseconds_sleep_time - microseconds_elapsed));
      
      // sleep 1000000us/100 Hz = 10ms between updates => 100 Hz polling/sampling interval
      // std::this_thread::sleep_for(std::chrono::milliseconds(1000/SAMPLING_HZ));
    }
    
    start = std::chrono::high_resolution_clock::now();
    
    
    {
      {
	std::lock_guard<std::mutex> lock(eeg_mutex);

	if(isnan(eegTP9) == 0 && isinf(eegTP9) == 0){
	  EEG[0].push_back(eegTP9);
	  EEGd[0].push_back(eegTP9);
	}
	else{
	  //printf("WARN: Bad EEG data: %f\n", eegTP9);

	  int end_eeg = EEG[0].size();
	  int start_eeg = EEG[0].size() - (int)SAMPLING_HZ;

	  float max = 0.0f;

	  if(start_eeg < 0) start_eeg = 0;

	  for(int i=start_eeg;i<end_eeg;i++){
	    if(fabs(EEG[0][i]) > max) max = fabs(EEG[0][i]);
	  }

	  float value = rng.uniform().c[0]*max;

	  EEG[0].push_back(value);
	  EEGd[0].push_back(value);
	}

	if(isnan(eegAF7) == 0 && isinf(eegAF7) == 0){
	  EEG[1].push_back(eegAF7);
	  EEGd[1].push_back(eegAF7);
	}
	else{
	  //printf("WARN: Bad EEG data: %f\n", eegAF7);

	  int end_eeg = EEG[1].size();
	  int start_eeg = EEG[1].size() - (int)SAMPLING_HZ;

	  float max = 0.0f;

	  if(start_eeg < 0) start_eeg = 0;

	  for(int i=start_eeg;i<end_eeg;i++){
	    if(fabs(EEG[1][i]) > max) max = fabs(EEG[1][i]);
	  }

	  float value = rng.uniform().c[0]*max;

	  EEG[1].push_back(value);
	  EEGd[1].push_back(value);
	}

	if(isnan(eegAF8) == 0 && isinf(eegAF8) == 0){
	  EEG[2].push_back(eegAF8);
	  EEGd[2].push_back(eegAF8);
	}
	else{
	  //printf("WARN: Bad EEG data: %f\n", eegAF8);
	  
	  int end_eeg = EEG[2].size();
	  int start_eeg = EEG[2].size() - (int)SAMPLING_HZ;

	  float max = 0.0f;

	  if(start_eeg < 0) start_eeg = 0;

	  for(int i=start_eeg;i<end_eeg;i++){
	    if(fabs(EEG[2][i]) > max) max = fabs(EEG[2][i]);
	  }

	  float value = rng.uniform().c[0]*max;

	  EEG[2].push_back(value);
	  EEGd[2].push_back(value);
	}

	if(isnan(eegTP10) == 0 && isinf(eegTP10) == 0){
	  EEG[3].push_back(eegTP10);
	  EEGd[3].push_back(eegTP10);
	}
	else{
	  //printf("WARN: Bad EEG data: %f\n", eegTP10);

	  int end_eeg = EEG[3].size();
	  int start_eeg = EEG[3].size() - (int)SAMPLING_HZ;

	  float max = 0.0f;

	  if(start_eeg < 0) start_eeg = 0;

	  for(int i=start_eeg;i<end_eeg;i++){
	    if(fabs(EEG[3][i]) > max) max = fabs(EEG[3][i]);
	  }

	  float value = rng.uniform().c[0]*max;

	  EEG[3].push_back(value);
	  EEGd[3].push_back(value);
	}
	
      }
    }

    {
      while(EEG[0].size() > EEG_LENGTH){
	for(unsigned int i=0;i<EEG.size();i++){
	  EEG[i].erase(EEG[i].begin()); // removes first element [SLOW]
	}
      }

      while(EEGd[0].size() > EEG_LENGTH){
	for(unsigned int i=0;i<EEGd.size();i++){
	  EEGd[i].erase(EEGd[i].begin()); // removes first element [SLOW]
	}
      }
      
      if(EEG[0].size() >= EEG_LENGTH){
	// calculate spectral entropy from FFT of 1-sec EEG
	
	std::vector<float> spectral_entropys;
	
	for(unsigned int a=0;a<EEG.size();a++){
	  spectral_entropys.push_back(spectral_entropy_eeg(EEG[a]));
	}

	spectral_entropy = 0.0f;
	
	for(const auto& se : spectral_entropys){
	  spectral_entropy += se; 
	}
	
	if(spectral_entropys.size() > 0)
	  spectral_entropy /= spectral_entropys.size();
	else
	  spectral_entropy = 0.50f;

	// calculate extra signals from 1-sec EEG (power spectrum)
	{
	  std::vector<float> power;
	  power.resize(EEG_LENGTH);
	  for(auto& p : power) p = 0.0f;

	  for(unsigned int i=0;i<EEGd.size();i++){
	    std::vector<double> P;
	    double hz_per_sample = 0.0;

	    P.resize(EEG_LENGTH);

	    for(auto& pi : P) pi = 0.0;

	    if(power_spectral_analysis(EEGd[i], SAMPLING_HZ, P, hz_per_sample)){
	      for(unsigned int n=0;n<P.size();n++){
		power[n] += sqrt(fabs(P[n]))/EEG.size();
	      }
	    }
	  }

	  float psum = 0.0f;
	  for(const auto& p : power) psum += p;
	  if(psum) for(auto& p : power) p /= psum;

	  // updates extra values (power spectrum)
	  {
	    std::lock_guard<std::mutex> lock(extra_data_mutex);
	    
	    extra_value = power;
	  }
	  
	}
	    
	    
	// calculate phase difference from FFT of 100-point EEG
	
	std::vector<float> phase_diffs;
	
	for(unsigned int a=0;a<EEG.size();a++){
	  for(unsigned int b=0;b<EEG.size();b++){
	    float pd = 0.0f;
	    
	    if(a == b) continue;
	    
	    if(spectral_phase_difference(EEG[a], EEG[b], pd))
	      phase_diffs.push_back(pd);
	  }
	}
	    
	phase_difference = 0.0f;
	
	for(unsigned int i=0;i<phase_diffs.size();i++){
	  phase_difference +=
	    (fabs(phase_diffs[i])/180.0f); // [0,1] range value
	}
	
	if(phase_diffs.size() > 0){
	  phase_difference /= phase_diffs.size();
	}
	else{
	  phase_difference = 1.00f;
	}
	
	
	// calculate phase slope difference from FFT of 64-point EEG
	
	std::vector<float> psi;
	
	for(unsigned int a=0;a<EEG.size();a++){
	  for(unsigned int b=0;b<EEG.size();b++){
	    float pd = 0.0f;
	    
	    if(a == b) continue;
	    
	    if(spectral_phase_slope_index(EEG[a], EEG[b], SAMPLING_HZ, pd))
	      psi.push_back(pd);
	  }
	}
	
	phase_slope_index = 0.0f;
	
	// average absolute PSI [high is better => high intelligence????]
	for(unsigned int i=0;i<psi.size();i++){
	  phase_slope_index += (fabs(psi[i])); // psi range???
	}
	
	if(psi.size() > 0){
	  phase_slope_index /= psi.size();
	}
	else{
	  phase_slope_index = 0.00f; // psi default value ???
	}
	
	// Calculates ICA signals and from ICA signals mean key figures
	if(calculate_ica_components){
	  whiteice::math::matrix<> PCA;
	  whiteice::math::vertex<> m;
	  whiteice::math::blas_real<float> origvar, reducedvar;

	  std::vector< whiteice::math::vertex<> > xdata;
	  whiteice::math::vertex<> x;

	  xdata.resize(EEG[0].size());
	  x.resize(EEG.size());

	  for(unsigned int i=0;i<EEG[0].size();i++){
	    for(unsigned int j=0;j<EEG.size();j++){
	      x[j] = EEG[j][i];
	    }
	    
	    xdata[i] = x;
	  }

	  // This calculates linear ICA from the measurements.
	  // FIXME: should calculate convolutional ICA instead with different time-delays to measurement points..
	  
	  if(whiteice::math::pca(xdata, EEG.size(), PCA, m, origvar, reducedvar, true, true) == false){
	    continue; // stop computations..
	  }
	  else{
	    for(unsigned int i=0;i<xdata.size();i++){
	      xdata[i] = PCA*(xdata[i] - m);
	    }
	  }

	  whiteice::math::matrix<> ICA;
	  
	  whiteice::math::blas_real<float> ica_tolerance = 0.01f; // converges more quickly..
	  const unsigned int ICA_MAXITERS = 100; // stops always quickly [adjust to proper value]

	  // FastICA algorithm is still slow for real-time use..
	  if(whiteice::math::ica(xdata, ICA, false, ica_tolerance, ICA_MAXITERS) == false){
	    continue; // stop computations..
	  }
	  else{
	    for(unsigned int i=0;i<xdata.size();i++){
	      xdata[i] = ICA*xdata[i];
	    }
	  }

	  // now converts format back to regular floats for computing key figures..

	  for(unsigned int i=0;i<xdata.size();i++){
	    for(unsigned int j=0;j<EEG.size();j++){
	      whiteice::math::convert(IC[j][i], xdata[i][j]);
	      whiteice::math::convert(ICd[j][i], xdata[i][j]);
	    }
	  }

	  std::vector<float> spectral_entropys;
		
	  for(unsigned int a=0;a<EEG.size();a++){
	    spectral_entropys.push_back(spectral_entropy_eeg(IC[a]));
	  }

	  ica_spectral_entropy = 0.0f;
	
	  for(const auto& se : spectral_entropys){
	    ica_spectral_entropy += se; 
	  }
	
	  if(spectral_entropys.size() > 0)
	    ica_spectral_entropy /= spectral_entropys.size();
	  else
	    ica_spectral_entropy = 0.50f;

	  // calculate extra signals from 1-sec EEG (power spectrum)
	  {
	    std::vector<float> power;
	    power.resize(EEG_LENGTH);
	    for(auto& p : power) p = 0.0f;

	    for(unsigned int i=0;i<ICd.size();i++){
	      std::vector<double> P;
	      double hz_per_sample = 0.0;

	      P.resize(EEG_LENGTH);
	      
	      for(auto& pi : P) pi = 0.0;
	      
	      if(power_spectral_analysis(ICd[i], SAMPLING_HZ, P, hz_per_sample)){
		for(unsigned int n=0;n<P.size();n++){
		  power[n] += sqrt(fabs(P[n]))/EEG.size();
		}
	      }
	    }
	    
	    float psum = 0.0f;
	    for(const auto& p : power) psum += p;
	    if(psum) for(auto& p : power) p /= psum;

	    // updates extra values (power spectrum)
	    {
	      std::lock_guard<std::mutex> lock(extra_data_mutex);
	      
	      ica_extra_value = power; // NOTE: mean power spectrum of ICA signals is not currently used..
	    }
	    
	  }
	    
	    
	  // calculate phase difference from FFT of 100-point ICA signals
	  
	  std::vector<float> phase_diffs;
	  
	  for(unsigned int a=0;a<IC.size();a++){
	    for(unsigned int b=0;b<IC.size();b++){
	      float pd = 0.0f;
	      
	      if(a == b) continue;
	      
	      if(spectral_phase_difference(IC[a], IC[b], pd))
		phase_diffs.push_back(pd);
	    }
	  }
	    
	  ica_phase_difference = 0.0f;
	
	  for(unsigned int i=0;i<phase_diffs.size();i++){
	    ica_phase_difference +=
	      (fabs(phase_diffs[i])/180.0f); // [0,1] range value
	  }
	  
	  if(phase_diffs.size() > 0){
	    ica_phase_difference /= phase_diffs.size();
	  }
	  else{
	    ica_phase_difference = 1.00f;
	  }
	  
	  
	  // calculate phase slope difference from FFT of 64-point EEG
	
	  std::vector<float> psi;
	  
	  for(unsigned int a=0;a<IC.size();a++){
	    for(unsigned int b=0;b<IC.size();b++){
	      float pd = 0.0f;
	      
	      if(a == b) continue;
	      
	      if(spectral_phase_slope_index(IC[a], IC[b], SAMPLING_HZ, pd))
		psi.push_back(pd);
	    }
	  }
	  
	  ica_phase_slope_index = 0.0f;
	
	  // average absolute PSI [high is better => high intelligence????]
	  for(unsigned int i=0;i<psi.size();i++){
	    ica_phase_slope_index += (fabs(psi[i])); // psi range???
	  }
	  
	  if(psi.size() > 0){
	    ica_phase_slope_index /= psi.size();
	  }
	  else{
	    ica_phase_slope_index = 0.00f; // psi default value ???
	  }
	  
	}
	
      }
      else{
	spectral_entropy = 0.50f; // mean spectral entropy value (50%)
	phase_difference = 1.00f; // maximum phase difference (100%)
	phase_slope_index = 0.0f; // psi default value ????

	ica_spectral_entropy = 0.50f; // mean spectral entropy value (50%)
	ica_phase_difference = 1.00f; // maximum phase difference (100%)
	ica_phase_slope_index = 0.0f; // psi default value ????
      }
      
    }
    
  }
  
}

  
} /* namespace resonanz */
} /* namespace whiteice */