diff --git a/Source/Processors/DataThreads/RHD2000Thread.cpp b/Source/Processors/DataThreads/RHD2000Thread.cpp index c7564316bc3fad95c4614cf64a1b92f7c0c3810d..7bc9ce5e0a5e0f27b24938ed02ff489d24991d43 100644 --- a/Source/Processors/DataThreads/RHD2000Thread.cpp +++ b/Source/Processors/DataThreads/RHD2000Thread.cpp @@ -124,11 +124,6 @@ RHD2000Thread::RHD2000Thread(SourceNode* sn) : DataThread(sn), // automatically find connected headstages scanPorts(); // things would appear to run more smoothly if this were done after the editor has been created - // to perform electrode impedance measurements at very low frequencies. - const int maxNumBlocks = 120; - int numStreams = 8; - allocateDoubleArray3D(amplifierPreFilter, numStreams, 32, SAMPLES_PER_DATA_BLOCK * maxNumBlocks); - // probably better to do this with a thread, but a timer works for now: // startTimer(10); // initialize the board in the background dacStream = new int[8]; @@ -1601,499 +1596,6 @@ bool RHD2000Thread::updateBuffer() } - -/***********************************/ -/* Below is code for impedance measurements */ - - - - -// Update electrode impedance measurement frequency, after checking that -// requested test frequency lies within acceptable ranges based on the -// amplifier bandwidth and the sampling rate. See impedancefreqdialog.cpp -// for more information. -float RHD2000Thread::updateImpedanceFrequency(float desiredImpedanceFreq, bool& impedanceFreqValid) -{ - int impedancePeriod; - double lowerBandwidthLimit, upperBandwidthLimit; - float actualImpedanceFreq; - - upperBandwidthLimit = actualUpperBandwidth / 1.5; - lowerBandwidthLimit = actualLowerBandwidth * 1.5; - if (dspEnabled) - { - if (actualDspCutoffFreq > actualLowerBandwidth) - { - lowerBandwidthLimit = actualDspCutoffFreq * 1.5; - } - } - - if (desiredImpedanceFreq > 0.0) - { - impedancePeriod = (boardSampleRate / desiredImpedanceFreq); - if (impedancePeriod >= 4 && impedancePeriod <= 1024 && - desiredImpedanceFreq >= lowerBandwidthLimit && - desiredImpedanceFreq <= upperBandwidthLimit) - { - actualImpedanceFreq = boardSampleRate / impedancePeriod; - impedanceFreqValid = true; - } - else - { - actualImpedanceFreq = 0.0; - impedanceFreqValid = false; - } - } - else - { - actualImpedanceFreq = 0.0; - impedanceFreqValid = false; - } - - return actualImpedanceFreq; -} - - -// Reads numBlocks blocks of raw USB data stored in a queue of Rhd2000DataBlock -// objects, loads this data into this SignalProcessor object, scaling the raw -// data to generate waveforms with units of volts or microvolts. -int RHD2000Thread::loadAmplifierData(queue<Rhd2000DataBlock>& dataQueue, - int numBlocks, int numDataStreams) -{ - - int block, t, channel, stream, i, j; - int indexAmp = 0; - int indexAux = 0; - int indexSupply = 0; - int indexAdc = 0; - int indexDig = 0; - int numWordsWritten = 0; - - int bufferIndex; - int16 tempQint16; - uint16 tempQuint16; - int32 tempQint32; - - bool triggerFound = false; - const double AnalogTriggerThreshold = 1.65; - - - for (block = 0; block < numBlocks; ++block) - { - - // Load and scale RHD2000 amplifier waveforms - // (sampled at amplifier sampling rate) - for (t = 0; t < SAMPLES_PER_DATA_BLOCK; ++t) - { - for (channel = 0; channel < 32; ++channel) - { - for (stream = 0; stream < numDataStreams; ++stream) - { - // Amplifier waveform units = microvolts - amplifierPreFilter[stream][channel][indexAmp] = 0.195 * - (dataQueue.front().amplifierData[stream][channel][t] - 32768); - } - } - ++indexAmp; - } - // We are done with this Rhd2000DataBlock object; remove it from dataQueue - dataQueue.pop(); - } - - return 0; -} - -#define PI 3.14159265359 -#define TWO_PI 6.28318530718 -#define DEGREES_TO_RADIANS 0.0174532925199 -#define RADIANS_TO_DEGREES 57.2957795132 - -// Return the magnitude and phase (in degrees) of a selected frequency component (in Hz) -// for a selected amplifier channel on the selected USB data stream. -void RHD2000Thread::measureComplexAmplitude(std::vector<std::vector<std::vector<double>>>& measuredMagnitude, - std::vector<std::vector<std::vector<double>>>& measuredPhase, - int capIndex, int stream, int chipChannel, int numBlocks, - double sampleRate, double frequency, int numPeriods) -{ - int period = (sampleRate / frequency); - int startIndex = 0; - int endIndex = startIndex + numPeriods * period - 1; - - // Move the measurement window to the end of the waveform to ignore start-up transient. - while (endIndex < SAMPLES_PER_DATA_BLOCK * numBlocks - period) - { - startIndex += period; - endIndex += period; - } - - double iComponent, qComponent; - - // Measure real (iComponent) and imaginary (qComponent) amplitude of frequency component. - amplitudeOfFreqComponent(iComponent, qComponent, amplifierPreFilter[stream][chipChannel], - startIndex, endIndex, sampleRate, frequency); - // Calculate magnitude and phase from real (I) and imaginary (Q) components. - measuredMagnitude[stream][chipChannel][capIndex] = - sqrt(iComponent * iComponent + qComponent * qComponent); - measuredPhase[stream][chipChannel][capIndex] = - RADIANS_TO_DEGREES *atan2(qComponent, iComponent); -} - -// Returns the real and imaginary amplitudes of a selected frequency component in the vector -// data, between a start index and end index. -void RHD2000Thread::amplitudeOfFreqComponent(double& realComponent, double& imagComponent, - const std::vector<double>& data, int startIndex, - int endIndex, double sampleRate, double frequency) -{ - int length = endIndex - startIndex + 1; - const double k = TWO_PI * frequency / sampleRate; // precalculate for speed - - // Perform correlation with sine and cosine waveforms. - double meanI = 0.0; - double meanQ = 0.0; - for (int t = startIndex; t <= endIndex; ++t) - { - meanI += data.at(t) * cos(k * t); - meanQ += data.at(t) * -1.0 * sin(k * t); - } - meanI /= (double) length; - meanQ /= (double) length; - - realComponent = 2.0 * meanI; - imagComponent = 2.0 * meanQ; -} - - - -// Given a measured complex impedance that is the result of an electrode impedance in parallel -// with a parasitic capacitance (i.e., due to the amplifier input capacitance and other -// capacitances associated with the chip bondpads), this function factors out the effect of the -// parasitic capacitance to return the acutal electrode impedance. -void RHD2000Thread::factorOutParallelCapacitance(double& impedanceMagnitude, double& impedancePhase, - double frequency, double parasiticCapacitance) -{ - // First, convert from polar coordinates to rectangular coordinates. - double measuredR = impedanceMagnitude * cos(DEGREES_TO_RADIANS * impedancePhase); - double measuredX = impedanceMagnitude * sin(DEGREES_TO_RADIANS * impedancePhase); - - double capTerm = TWO_PI * frequency * parasiticCapacitance; - double xTerm = capTerm * (measuredR * measuredR + measuredX * measuredX); - double denominator = capTerm * xTerm + 2 * capTerm * measuredX + 1; - double trueR = measuredR / denominator; - double trueX = (measuredX + xTerm) / denominator; - - // Now, convert from rectangular coordinates back to polar coordinates. - impedanceMagnitude = sqrt(trueR * trueR + trueX * trueX); - impedancePhase = RADIANS_TO_DEGREES * atan2(trueX, trueR); -} - -// This is a purely empirical function to correct observed errors in the real component -// of measured electrode impedances at sampling rates below 15 kS/s. At low sampling rates, -// it is difficult to approximate a smooth sine wave with the on-chip voltage DAC and 10 kHz -// 2-pole lowpass filter. This function attempts to somewhat correct for this, but a better -// solution is to always run impedance measurements at 20 kS/s, where they seem to be most -// accurate. -void RHD2000Thread::empiricalResistanceCorrection(double& impedanceMagnitude, double& impedancePhase, - double boardSampleRate) -{ - // First, convert from polar coordinates to rectangular coordinates. - double impedanceR = impedanceMagnitude * cos(DEGREES_TO_RADIANS * impedancePhase); - double impedanceX = impedanceMagnitude * sin(DEGREES_TO_RADIANS * impedancePhase); - - // Emprically derived correction factor (i.e., no physical basis for this equation). - impedanceR /= 10.0 * exp(-boardSampleRate / 2500.0) * cos(TWO_PI * boardSampleRate / 15000.0) + 1.0; - - // Now, convert from rectangular coordinates back to polar coordinates. - impedanceMagnitude = sqrt(impedanceR * impedanceR + impedanceX * impedanceX); - impedancePhase = RADIANS_TO_DEGREES * atan2(impedanceX, impedanceR); -} - -void RHD2000Thread::runImpedanceTest(Array<int>& streams, Array<int>& channels, Array<float>& magnitudes, Array<float>& phases) -{ - int commandSequenceLength, stream, channel, capRange; - double cSeries; - vector<int> commandList; - int triggerIndex; // dummy reference variable; not used - queue<Rhd2000DataBlock> bufferQueue; // dummy reference variable; not used - int numdataStreams = evalBoard->getNumEnabledDataStreams(); - - bool rhd2164ChipPresent = false; - Array<int> enabledStreams; - for (stream = 0; stream < MAX_NUM_DATA_STREAMS; ++stream) - { - - if (evalBoard->isStreamEnabled(stream)) - { - enabledStreams.add(stream); - } - - if (chipId[stream] == CHIP_ID_RHD2164_B) - { - rhd2164ChipPresent = true; - } - } - - bool validImpedanceFreq; - float actualImpedanceFreq = updateImpedanceFrequency(1000.0, validImpedanceFreq); - if (!validImpedanceFreq) - { - return ; - } - // Create a command list for the AuxCmd1 slot. - commandSequenceLength = chipRegisters.createCommandListZcheckDac(commandList, actualImpedanceFreq, 128.0); - evalBoard->uploadCommandList(commandList, Rhd2000EvalBoard::AuxCmd1, 1); - evalBoard->selectAuxCommandLength(Rhd2000EvalBoard::AuxCmd1, - 0, commandSequenceLength - 1); - if (fastTTLSettleEnabled) - { - evalBoard->enableExternalFastSettle(false); - } - - evalBoard->selectAuxCommandBank(Rhd2000EvalBoard::PortA, - Rhd2000EvalBoard::AuxCmd1, 1); - evalBoard->selectAuxCommandBank(Rhd2000EvalBoard::PortB, - Rhd2000EvalBoard::AuxCmd1, 1); - evalBoard->selectAuxCommandBank(Rhd2000EvalBoard::PortC, - Rhd2000EvalBoard::AuxCmd1, 1); - evalBoard->selectAuxCommandBank(Rhd2000EvalBoard::PortD, - Rhd2000EvalBoard::AuxCmd1, 1); - - // Select number of periods to measure impedance over - int numPeriods = (0.020 * actualImpedanceFreq); // Test each channel for at least 20 msec... - if (numPeriods < 5) numPeriods = 5; // ...but always measure across no fewer than 5 complete periods - double period = boardSampleRate / actualImpedanceFreq; - int numBlocks = ceil((numPeriods + 2.0) * period / 60.0); // + 2 periods to give time to settle initially - if (numBlocks < 2) numBlocks = 2; // need first block for command to switch channels to take effect. - - actualDspCutoffFreq = chipRegisters.setDspCutoffFreq(desiredDspCutoffFreq); - actualLowerBandwidth = chipRegisters.setLowerBandwidth(desiredLowerBandwidth); - actualUpperBandwidth = chipRegisters.setUpperBandwidth(desiredUpperBandwidth); - chipRegisters.enableDsp(dspEnabled); - chipRegisters.enableZcheck(true); - commandSequenceLength = chipRegisters.createCommandListRegisterConfig(commandList, false); - // Upload version with no ADC calibration to AuxCmd3 RAM Bank 1. - evalBoard->uploadCommandList(commandList, Rhd2000EvalBoard::AuxCmd3, 3); - evalBoard->selectAuxCommandLength(Rhd2000EvalBoard::AuxCmd3, 0, commandSequenceLength - 1); - evalBoard->selectAuxCommandBank(Rhd2000EvalBoard::PortA, Rhd2000EvalBoard::AuxCmd3, 3); - evalBoard->selectAuxCommandBank(Rhd2000EvalBoard::PortB, Rhd2000EvalBoard::AuxCmd3, 3); - evalBoard->selectAuxCommandBank(Rhd2000EvalBoard::PortC, Rhd2000EvalBoard::AuxCmd3, 3); - evalBoard->selectAuxCommandBank(Rhd2000EvalBoard::PortD, Rhd2000EvalBoard::AuxCmd3, 3); - - evalBoard->setContinuousRunMode(false); - evalBoard->setMaxTimeStep(SAMPLES_PER_DATA_BLOCK * numBlocks); - - // Create matrices of doubles of size (numStreams x 32 x 3) to store complex amplitudes - // of all amplifier channels (32 on each data stream) at three different Cseries values. - std::vector<std::vector<std::vector<double>>> measuredMagnitude; - std::vector<std::vector<std::vector<double>>> measuredPhase; - - measuredMagnitude.resize(evalBoard->getNumEnabledDataStreams()); - measuredPhase.resize(evalBoard->getNumEnabledDataStreams()); - for (int i = 0; i < evalBoard->getNumEnabledDataStreams(); ++i) - { - measuredMagnitude[i].resize(32); - measuredPhase[i].resize(32); - for (int j = 0; j < 32; ++j) - { - measuredMagnitude[i][j].resize(3); - measuredPhase[i][j].resize(3); - } - } - - - - double distance, minDistance, current, Cseries; - double impedanceMagnitude, impedancePhase; - - const double bestAmplitude = 250.0; // we favor voltage readings that are closest to 250 uV: not too large, - // and not too small. - const double dacVoltageAmplitude = 128 * (1.225 / 256); // this assumes the DAC amplitude was set to 128 - const double parasiticCapacitance = 14.0e-12; // 14 pF: an estimate of on-chip parasitic capacitance, - // including 10 pF of amplifier input capacitance. - double relativeFreq = actualImpedanceFreq / boardSampleRate; - - int bestAmplitudeIndex; - - // We execute three complete electrode impedance measurements: one each with - // Cseries set to 0.1 pF, 1 pF, and 10 pF. Then we select the best measurement - // for each channel so that we achieve a wide impedance measurement range. - for (capRange = 0; capRange < 3; ++ capRange) - { - - - switch (capRange) - { - case 0: - chipRegisters.setZcheckScale(Rhd2000Registers::ZcheckCs100fF); - cSeries = 0.1e-12; - cout << "setting capacitance to 0.1pF" << endl; - break; - case 1: - chipRegisters.setZcheckScale(Rhd2000Registers::ZcheckCs1pF); - cSeries = 1.0e-12; - cout << "setting capacitance to 1pF" << endl; - break; - case 2: - chipRegisters.setZcheckScale(Rhd2000Registers::ZcheckCs10pF); - cSeries = 10.0e-12; - cout << "setting capacitance to 10pF" << endl; - break; - } - - // Check all 32 channels across all active data streams. - for (channel = 0; channel < 32; ++channel) - { - cout << "running impedance on channel " << channel << endl; - - chipRegisters.setZcheckChannel(channel); - commandSequenceLength = - chipRegisters.createCommandListRegisterConfig(commandList, false); - // Upload version with no ADC calibration to AuxCmd3 RAM Bank 1. - evalBoard->uploadCommandList(commandList, Rhd2000EvalBoard::AuxCmd3, 3); - - evalBoard->run(); - while (evalBoard->isRunning()) - { - - } - queue<Rhd2000DataBlock> dataQueue; - evalBoard->readDataBlocks(numBlocks, dataQueue); - loadAmplifierData(dataQueue, numBlocks, numdataStreams); - for (stream = 0; stream < numdataStreams; ++stream) - { - if (chipId[stream] != CHIP_ID_RHD2164_B) - { - measureComplexAmplitude(measuredMagnitude, measuredPhase, - capRange, stream, channel, numBlocks, boardSampleRate, - actualImpedanceFreq, numPeriods); - } - } - - // If an RHD2164 chip is plugged in, we have to set the Zcheck select register to channels 32-63 - // and repeat the previous steps. - if (rhd2164ChipPresent) - { - chipRegisters.setZcheckChannel(channel + 32); // address channels 32-63 - commandSequenceLength = - chipRegisters.createCommandListRegisterConfig(commandList, false); - // Upload version with no ADC calibration to AuxCmd3 RAM Bank 1. - evalBoard->uploadCommandList(commandList, Rhd2000EvalBoard::AuxCmd3, 3); - - evalBoard->run(); - while (evalBoard->isRunning()) - { - - } - evalBoard->readDataBlocks(numBlocks, dataQueue); - loadAmplifierData(dataQueue, numBlocks, numdataStreams); - - for (stream = 0; stream < evalBoard->getNumEnabledDataStreams(); ++stream) - { - if (chipId[stream] == CHIP_ID_RHD2164_B) - { - measureComplexAmplitude(measuredMagnitude, measuredPhase, - capRange, stream, channel, numBlocks, boardSampleRate, - actualImpedanceFreq, numPeriods); - } - } - } - } - } - - streams.clear(); - channels.clear(); - magnitudes.clear(); - phases.clear(); - - for (stream = 0; stream < evalBoard->getNumEnabledDataStreams(); ++stream) - { - for (channel = 0; channel < 32; ++channel) - { - if (1) - { - minDistance = 9.9e99; // ridiculously large number - for (capRange = 0; capRange < 3; ++capRange) - { - // Find the measured amplitude that is closest to bestAmplitude on a logarithmic scale - distance = abs(log(measuredMagnitude[stream][channel][capRange] / bestAmplitude)); - if (distance < minDistance) - { - bestAmplitudeIndex = capRange; - minDistance = distance; - } - } - switch (bestAmplitudeIndex) - { - case 0: - Cseries = 0.1e-12; - break; - case 1: - Cseries = 1.0e-12; - break; - case 2: - Cseries = 10.0e-12; - break; - } - - // Calculate current amplitude produced by on-chip voltage DAC - current = TWO_PI * actualImpedanceFreq * dacVoltageAmplitude * Cseries; - - // Calculate impedance magnitude from calculated current and measured voltage. - impedanceMagnitude = 1.0e-6 * (measuredMagnitude[stream][channel][bestAmplitudeIndex] / current) * - (18.0 * relativeFreq * relativeFreq + 1.0); - - // Calculate impedance phase, with small correction factor accounting for the - // 3-command SPI pipeline delay. - impedancePhase = measuredPhase[stream][channel][bestAmplitudeIndex] + (360.0 * (3.0 / period)); - - // Factor out on-chip parasitic capacitance from impedance measurement. - factorOutParallelCapacitance(impedanceMagnitude, impedancePhase, actualImpedanceFreq, - parasiticCapacitance); - - // Perform empirical resistance correction to improve accuarcy at sample rates below - // 15 kS/s. - empiricalResistanceCorrection(impedanceMagnitude, impedancePhase, - boardSampleRate); - - streams.add(enabledStreams[stream]); - channels.add(channel); - magnitudes.add(impedanceMagnitude); - phases.add(impedancePhase); - - if (impedanceMagnitude > 1000000) - cout << "stream " << stream << " channel " << 1+channel << " magnitude: " << String(impedanceMagnitude/1e6,2) << " mOhm , phase : " <<impedancePhase << endl; - else - cout << "stream " << stream << " channel " << 1+channel << " magnitude: " << String(impedanceMagnitude/1e3,2) << " kOhm , phase : " <<impedancePhase << endl; - - } - } - } - - evalBoard->setContinuousRunMode(false); - evalBoard->setMaxTimeStep(0); - evalBoard->flush(); - - // Switch back to flatline - evalBoard->selectAuxCommandBank(Rhd2000EvalBoard::PortA, Rhd2000EvalBoard::AuxCmd1, 0); - evalBoard->selectAuxCommandBank(Rhd2000EvalBoard::PortB, Rhd2000EvalBoard::AuxCmd1, 0); - evalBoard->selectAuxCommandBank(Rhd2000EvalBoard::PortC, Rhd2000EvalBoard::AuxCmd1, 0); - evalBoard->selectAuxCommandBank(Rhd2000EvalBoard::PortD, Rhd2000EvalBoard::AuxCmd1, 0); - evalBoard->selectAuxCommandLength(Rhd2000EvalBoard::AuxCmd1, 0, 1); - - evalBoard->selectAuxCommandBank(Rhd2000EvalBoard::PortA, Rhd2000EvalBoard::AuxCmd3, - fastSettleEnabled ? 2 : 1); - evalBoard->selectAuxCommandBank(Rhd2000EvalBoard::PortB, Rhd2000EvalBoard::AuxCmd3, - fastSettleEnabled ? 2 : 1); - evalBoard->selectAuxCommandBank(Rhd2000EvalBoard::PortC, Rhd2000EvalBoard::AuxCmd3, - fastSettleEnabled ? 2 : 1); - evalBoard->selectAuxCommandBank(Rhd2000EvalBoard::PortD, Rhd2000EvalBoard::AuxCmd3, - fastSettleEnabled ? 2 : 1); - - if (fastTTLSettleEnabled) - { - evalBoard->enableExternalFastSettle(true); - } -} - int RHD2000Thread::getChannelFromHeadstage(int hs, int ch) { int channelCount = 0; @@ -2237,4 +1739,504 @@ Rhd2000EvalBoard::BoardDataSource RHDHeadstage::getDataStream(int index) bool RHDHeadstage::isPlugged() { return (numStreams > 0); +} + +/***********************************/ +/* Below is code for impedance measurements */ + +RHDImpedanceMeasure::RHDImpedanceMeasure(RHD2000Thread* b, ImpedanceData* d) : Thread(""), data(d), board(b) +{ + // to perform electrode impedance measurements at very low frequencies. + const int maxNumBlocks = 120; + int numStreams = 8; + allocateDoubleArray3D(amplifierPreFilter, numStreams, 32, SAMPLES_PER_DATA_BLOCK * maxNumBlocks); +} + + + + +// Update electrode impedance measurement frequency, after checking that +// requested test frequency lies within acceptable ranges based on the +// amplifier bandwidth and the sampling rate. See impedancefreqdialog.cpp +// for more information. +float RHDImpedanceMeasure::updateImpedanceFrequency(float desiredImpedanceFreq, bool& impedanceFreqValid) +{ + int impedancePeriod; + double lowerBandwidthLimit, upperBandwidthLimit; + float actualImpedanceFreq; + + upperBandwidthLimit = board->actualUpperBandwidth / 1.5; + lowerBandwidthLimit = board->actualLowerBandwidth * 1.5; + if (board->dspEnabled) + { + if (board->actualDspCutoffFreq > board->actualLowerBandwidth) + { + lowerBandwidthLimit = board->actualDspCutoffFreq * 1.5; + } + } + + if (desiredImpedanceFreq > 0.0) + { + impedancePeriod = (board->boardSampleRate / desiredImpedanceFreq); + if (impedancePeriod >= 4 && impedancePeriod <= 1024 && + desiredImpedanceFreq >= lowerBandwidthLimit && + desiredImpedanceFreq <= upperBandwidthLimit) + { + actualImpedanceFreq = board->boardSampleRate / impedancePeriod; + impedanceFreqValid = true; + } + else + { + actualImpedanceFreq = 0.0; + impedanceFreqValid = false; + } + } + else + { + actualImpedanceFreq = 0.0; + impedanceFreqValid = false; + } + + return actualImpedanceFreq; +} + + +// Reads numBlocks blocks of raw USB data stored in a queue of Rhd2000DataBlock +// objects, loads this data into this SignalProcessor object, scaling the raw +// data to generate waveforms with units of volts or microvolts. +int RHDImpedanceMeasure::loadAmplifierData(queue<Rhd2000DataBlock>& dataQueue, + int numBlocks, int numDataStreams) +{ + + int block, t, channel, stream, i, j; + int indexAmp = 0; + int indexAux = 0; + int indexSupply = 0; + int indexAdc = 0; + int indexDig = 0; + int numWordsWritten = 0; + + int bufferIndex; + int16 tempQint16; + uint16 tempQuint16; + int32 tempQint32; + + bool triggerFound = false; + const double AnalogTriggerThreshold = 1.65; + + + for (block = 0; block < numBlocks; ++block) + { + + // Load and scale RHD2000 amplifier waveforms + // (sampled at amplifier sampling rate) + for (t = 0; t < SAMPLES_PER_DATA_BLOCK; ++t) + { + for (channel = 0; channel < 32; ++channel) + { + for (stream = 0; stream < numDataStreams; ++stream) + { + // Amplifier waveform units = microvolts + amplifierPreFilter[stream][channel][indexAmp] = 0.195 * + (dataQueue.front().amplifierData[stream][channel][t] - 32768); + } + } + ++indexAmp; + } + // We are done with this Rhd2000DataBlock object; remove it from dataQueue + dataQueue.pop(); + } + + return 0; +} + +#define PI 3.14159265359 +#define TWO_PI 6.28318530718 +#define DEGREES_TO_RADIANS 0.0174532925199 +#define RADIANS_TO_DEGREES 57.2957795132 + +// Return the magnitude and phase (in degrees) of a selected frequency component (in Hz) +// for a selected amplifier channel on the selected USB data stream. +void RHDImpedanceMeasure::measureComplexAmplitude(std::vector<std::vector<std::vector<double>>>& measuredMagnitude, + std::vector<std::vector<std::vector<double>>>& measuredPhase, + int capIndex, int stream, int chipChannel, int numBlocks, + double sampleRate, double frequency, int numPeriods) +{ + int period = (sampleRate / frequency); + int startIndex = 0; + int endIndex = startIndex + numPeriods * period - 1; + + // Move the measurement window to the end of the waveform to ignore start-up transient. + while (endIndex < SAMPLES_PER_DATA_BLOCK * numBlocks - period) + { + startIndex += period; + endIndex += period; + } + + double iComponent, qComponent; + + // Measure real (iComponent) and imaginary (qComponent) amplitude of frequency component. + amplitudeOfFreqComponent(iComponent, qComponent, amplifierPreFilter[stream][chipChannel], + startIndex, endIndex, sampleRate, frequency); + // Calculate magnitude and phase from real (I) and imaginary (Q) components. + measuredMagnitude[stream][chipChannel][capIndex] = + sqrt(iComponent * iComponent + qComponent * qComponent); + measuredPhase[stream][chipChannel][capIndex] = + RADIANS_TO_DEGREES *atan2(qComponent, iComponent); +} + +// Returns the real and imaginary amplitudes of a selected frequency component in the vector +// data, between a start index and end index. +void RHDImpedanceMeasure::amplitudeOfFreqComponent(double& realComponent, double& imagComponent, + const std::vector<double>& data, int startIndex, + int endIndex, double sampleRate, double frequency) +{ + int length = endIndex - startIndex + 1; + const double k = TWO_PI * frequency / sampleRate; // precalculate for speed + + // Perform correlation with sine and cosine waveforms. + double meanI = 0.0; + double meanQ = 0.0; + for (int t = startIndex; t <= endIndex; ++t) + { + meanI += data.at(t) * cos(k * t); + meanQ += data.at(t) * -1.0 * sin(k * t); + } + meanI /= (double)length; + meanQ /= (double)length; + + realComponent = 2.0 * meanI; + imagComponent = 2.0 * meanQ; +} + + + +// Given a measured complex impedance that is the result of an electrode impedance in parallel +// with a parasitic capacitance (i.e., due to the amplifier input capacitance and other +// capacitances associated with the chip bondpads), this function factors out the effect of the +// parasitic capacitance to return the acutal electrode impedance. +void RHDImpedanceMeasure::factorOutParallelCapacitance(double& impedanceMagnitude, double& impedancePhase, + double frequency, double parasiticCapacitance) +{ + // First, convert from polar coordinates to rectangular coordinates. + double measuredR = impedanceMagnitude * cos(DEGREES_TO_RADIANS * impedancePhase); + double measuredX = impedanceMagnitude * sin(DEGREES_TO_RADIANS * impedancePhase); + + double capTerm = TWO_PI * frequency * parasiticCapacitance; + double xTerm = capTerm * (measuredR * measuredR + measuredX * measuredX); + double denominator = capTerm * xTerm + 2 * capTerm * measuredX + 1; + double trueR = measuredR / denominator; + double trueX = (measuredX + xTerm) / denominator; + + // Now, convert from rectangular coordinates back to polar coordinates. + impedanceMagnitude = sqrt(trueR * trueR + trueX * trueX); + impedancePhase = RADIANS_TO_DEGREES * atan2(trueX, trueR); +} + +// This is a purely empirical function to correct observed errors in the real component +// of measured electrode impedances at sampling rates below 15 kS/s. At low sampling rates, +// it is difficult to approximate a smooth sine wave with the on-chip voltage DAC and 10 kHz +// 2-pole lowpass filter. This function attempts to somewhat correct for this, but a better +// solution is to always run impedance measurements at 20 kS/s, where they seem to be most +// accurate. +void RHDImpedanceMeasure::empiricalResistanceCorrection(double& impedanceMagnitude, double& impedancePhase, + double boardSampleRate) +{ + // First, convert from polar coordinates to rectangular coordinates. + double impedanceR = impedanceMagnitude * cos(DEGREES_TO_RADIANS * impedancePhase); + double impedanceX = impedanceMagnitude * sin(DEGREES_TO_RADIANS * impedancePhase); + + // Emprically derived correction factor (i.e., no physical basis for this equation). + impedanceR /= 10.0 * exp(-boardSampleRate / 2500.0) * cos(TWO_PI * boardSampleRate / 15000.0) + 1.0; + + // Now, convert from rectangular coordinates back to polar coordinates. + impedanceMagnitude = sqrt(impedanceR * impedanceR + impedanceX * impedanceX); + impedancePhase = RADIANS_TO_DEGREES * atan2(impedanceX, impedanceR); +} + +void RHDImpedanceMeasure::run() +{ + int commandSequenceLength, stream, channel, capRange; + double cSeries; + vector<int> commandList; + int triggerIndex; // dummy reference variable; not used + queue<Rhd2000DataBlock> bufferQueue; // dummy reference variable; not used + int numdataStreams = board->evalBoard->getNumEnabledDataStreams(); + + bool rhd2164ChipPresent = false; + Array<int> enabledStreams; + for (stream = 0; stream < MAX_NUM_DATA_STREAMS; ++stream) + { + + if (board->evalBoard->isStreamEnabled(stream)) + { + enabledStreams.add(stream); + } + + if (board->chipId[stream] == CHIP_ID_RHD2164_B) + { + rhd2164ChipPresent = true; + } + } + + bool validImpedanceFreq; + float actualImpedanceFreq = updateImpedanceFrequency(1000.0, validImpedanceFreq); + if (!validImpedanceFreq) + { + return; + } + // Create a command list for the AuxCmd1 slot. + commandSequenceLength = board->chipRegisters.createCommandListZcheckDac(commandList, actualImpedanceFreq, 128.0); + board->evalBoard->uploadCommandList(commandList, Rhd2000EvalBoard::AuxCmd1, 1); + board->evalBoard->selectAuxCommandLength(Rhd2000EvalBoard::AuxCmd1, + 0, commandSequenceLength - 1); + if (board->fastTTLSettleEnabled) + { + board->evalBoard->enableExternalFastSettle(false); + } + + board->evalBoard->selectAuxCommandBank(Rhd2000EvalBoard::PortA, + Rhd2000EvalBoard::AuxCmd1, 1); + board->evalBoard->selectAuxCommandBank(Rhd2000EvalBoard::PortB, + Rhd2000EvalBoard::AuxCmd1, 1); + board->evalBoard->selectAuxCommandBank(Rhd2000EvalBoard::PortC, + Rhd2000EvalBoard::AuxCmd1, 1); + board->evalBoard->selectAuxCommandBank(Rhd2000EvalBoard::PortD, + Rhd2000EvalBoard::AuxCmd1, 1); + + // Select number of periods to measure impedance over + int numPeriods = (0.020 * actualImpedanceFreq); // Test each channel for at least 20 msec... + if (numPeriods < 5) numPeriods = 5; // ...but always measure across no fewer than 5 complete periods + double period = board->boardSampleRate / actualImpedanceFreq; + int numBlocks = ceil((numPeriods + 2.0) * period / 60.0); // + 2 periods to give time to settle initially + if (numBlocks < 2) numBlocks = 2; // need first block for command to switch channels to take effect. + + board->actualDspCutoffFreq = board->chipRegisters.setDspCutoffFreq(board->desiredDspCutoffFreq); + board->actualLowerBandwidth = board->chipRegisters.setLowerBandwidth(board->desiredLowerBandwidth); + board->actualUpperBandwidth = board->chipRegisters.setUpperBandwidth(board->desiredUpperBandwidth); + board->chipRegisters.enableDsp(board->dspEnabled); + board->chipRegisters.enableZcheck(true); + commandSequenceLength = board->chipRegisters.createCommandListRegisterConfig(commandList, false); + // Upload version with no ADC calibration to AuxCmd3 RAM Bank 1. + board->evalBoard->uploadCommandList(commandList, Rhd2000EvalBoard::AuxCmd3, 3); + board->evalBoard->selectAuxCommandLength(Rhd2000EvalBoard::AuxCmd3, 0, commandSequenceLength - 1); + board->evalBoard->selectAuxCommandBank(Rhd2000EvalBoard::PortA, Rhd2000EvalBoard::AuxCmd3, 3); + board->evalBoard->selectAuxCommandBank(Rhd2000EvalBoard::PortB, Rhd2000EvalBoard::AuxCmd3, 3); + board->evalBoard->selectAuxCommandBank(Rhd2000EvalBoard::PortC, Rhd2000EvalBoard::AuxCmd3, 3); + board->evalBoard->selectAuxCommandBank(Rhd2000EvalBoard::PortD, Rhd2000EvalBoard::AuxCmd3, 3); + + board->evalBoard->setContinuousRunMode(false); + board->evalBoard->setMaxTimeStep(SAMPLES_PER_DATA_BLOCK * numBlocks); + + // Create matrices of doubles of size (numStreams x 32 x 3) to store complex amplitudes + // of all amplifier channels (32 on each data stream) at three different Cseries values. + std::vector<std::vector<std::vector<double>>> measuredMagnitude; + std::vector<std::vector<std::vector<double>>> measuredPhase; + + measuredMagnitude.resize(board->evalBoard->getNumEnabledDataStreams()); + measuredPhase.resize(board->evalBoard->getNumEnabledDataStreams()); + for (int i = 0; i < board->evalBoard->getNumEnabledDataStreams(); ++i) + { + measuredMagnitude[i].resize(32); + measuredPhase[i].resize(32); + for (int j = 0; j < 32; ++j) + { + measuredMagnitude[i][j].resize(3); + measuredPhase[i][j].resize(3); + } + } + + + + double distance, minDistance, current, Cseries; + double impedanceMagnitude, impedancePhase; + + const double bestAmplitude = 250.0; // we favor voltage readings that are closest to 250 uV: not too large, + // and not too small. + const double dacVoltageAmplitude = 128 * (1.225 / 256); // this assumes the DAC amplitude was set to 128 + const double parasiticCapacitance = 14.0e-12; // 14 pF: an estimate of on-chip parasitic capacitance, + // including 10 pF of amplifier input capacitance. + double relativeFreq = actualImpedanceFreq / board->boardSampleRate; + + int bestAmplitudeIndex; + + // We execute three complete electrode impedance measurements: one each with + // Cseries set to 0.1 pF, 1 pF, and 10 pF. Then we select the best measurement + // for each channel so that we achieve a wide impedance measurement range. + for (capRange = 0; capRange < 3; ++capRange) + { + + + switch (capRange) + { + case 0: + board->chipRegisters.setZcheckScale(Rhd2000Registers::ZcheckCs100fF); + cSeries = 0.1e-12; + cout << "setting capacitance to 0.1pF" << endl; + break; + case 1: + board->chipRegisters.setZcheckScale(Rhd2000Registers::ZcheckCs1pF); + cSeries = 1.0e-12; + cout << "setting capacitance to 1pF" << endl; + break; + case 2: + board->chipRegisters.setZcheckScale(Rhd2000Registers::ZcheckCs10pF); + cSeries = 10.0e-12; + cout << "setting capacitance to 10pF" << endl; + break; + } + + // Check all 32 channels across all active data streams. + for (channel = 0; channel < 32; ++channel) + { + cout << "running impedance on channel " << channel << endl; + + board->chipRegisters.setZcheckChannel(channel); + commandSequenceLength = + board->chipRegisters.createCommandListRegisterConfig(commandList, false); + // Upload version with no ADC calibration to AuxCmd3 RAM Bank 1. + board->evalBoard->uploadCommandList(commandList, Rhd2000EvalBoard::AuxCmd3, 3); + + board->evalBoard->run(); + while (board->evalBoard->isRunning()) + { + + } + queue<Rhd2000DataBlock> dataQueue; + board->evalBoard->readDataBlocks(numBlocks, dataQueue); + loadAmplifierData(dataQueue, numBlocks, numdataStreams); + for (stream = 0; stream < numdataStreams; ++stream) + { + if (board->chipId[stream] != CHIP_ID_RHD2164_B) + { + measureComplexAmplitude(measuredMagnitude, measuredPhase, + capRange, stream, channel, numBlocks, board->boardSampleRate, + actualImpedanceFreq, numPeriods); + } + } + + // If an RHD2164 chip is plugged in, we have to set the Zcheck select register to channels 32-63 + // and repeat the previous steps. + if (rhd2164ChipPresent) + { + board->chipRegisters.setZcheckChannel(channel + 32); // address channels 32-63 + commandSequenceLength = + board->chipRegisters.createCommandListRegisterConfig(commandList, false); + // Upload version with no ADC calibration to AuxCmd3 RAM Bank 1. + board->evalBoard->uploadCommandList(commandList, Rhd2000EvalBoard::AuxCmd3, 3); + + board->evalBoard->run(); + while (board->evalBoard->isRunning()) + { + + } + board->evalBoard->readDataBlocks(numBlocks, dataQueue); + loadAmplifierData(dataQueue, numBlocks, numdataStreams); + + for (stream = 0; stream < board->evalBoard->getNumEnabledDataStreams(); ++stream) + { + if (board->chipId[stream] == CHIP_ID_RHD2164_B) + { + measureComplexAmplitude(measuredMagnitude, measuredPhase, + capRange, stream, channel, numBlocks, board->boardSampleRate, + actualImpedanceFreq, numPeriods); + } + } + } + } + } + + data->streams.clear(); + data->channels.clear(); + data->magnitudes.clear(); + data->phases.clear(); + + for (stream = 0; stream < board->evalBoard->getNumEnabledDataStreams(); ++stream) + { + for (channel = 0; channel < 32; ++channel) + { + if (1) + { + minDistance = 9.9e99; // ridiculously large number + for (capRange = 0; capRange < 3; ++capRange) + { + // Find the measured amplitude that is closest to bestAmplitude on a logarithmic scale + distance = abs(log(measuredMagnitude[stream][channel][capRange] / bestAmplitude)); + if (distance < minDistance) + { + bestAmplitudeIndex = capRange; + minDistance = distance; + } + } + switch (bestAmplitudeIndex) + { + case 0: + Cseries = 0.1e-12; + break; + case 1: + Cseries = 1.0e-12; + break; + case 2: + Cseries = 10.0e-12; + break; + } + + // Calculate current amplitude produced by on-chip voltage DAC + current = TWO_PI * actualImpedanceFreq * dacVoltageAmplitude * Cseries; + + // Calculate impedance magnitude from calculated current and measured voltage. + impedanceMagnitude = 1.0e-6 * (measuredMagnitude[stream][channel][bestAmplitudeIndex] / current) * + (18.0 * relativeFreq * relativeFreq + 1.0); + + // Calculate impedance phase, with small correction factor accounting for the + // 3-command SPI pipeline delay. + impedancePhase = measuredPhase[stream][channel][bestAmplitudeIndex] + (360.0 * (3.0 / period)); + + // Factor out on-chip parasitic capacitance from impedance measurement. + factorOutParallelCapacitance(impedanceMagnitude, impedancePhase, actualImpedanceFreq, + parasiticCapacitance); + + // Perform empirical resistance correction to improve accuarcy at sample rates below + // 15 kS/s. + empiricalResistanceCorrection(impedanceMagnitude, impedancePhase, + board->boardSampleRate); + + data->streams.add(enabledStreams[stream]); + data->channels.add(channel); + data->magnitudes.add(impedanceMagnitude); + data->phases.add(impedancePhase); + + if (impedanceMagnitude > 1000000) + cout << "stream " << stream << " channel " << 1 + channel << " magnitude: " << String(impedanceMagnitude / 1e6, 2) << " mOhm , phase : " << impedancePhase << endl; + else + cout << "stream " << stream << " channel " << 1 + channel << " magnitude: " << String(impedanceMagnitude / 1e3, 2) << " kOhm , phase : " << impedancePhase << endl; + + } + } + } + + board->evalBoard->setContinuousRunMode(false); + board->evalBoard->setMaxTimeStep(0); + board->evalBoard->flush(); + + // Switch back to flatline + board->evalBoard->selectAuxCommandBank(Rhd2000EvalBoard::PortA, Rhd2000EvalBoard::AuxCmd1, 0); + board->evalBoard->selectAuxCommandBank(Rhd2000EvalBoard::PortB, Rhd2000EvalBoard::AuxCmd1, 0); + board->evalBoard->selectAuxCommandBank(Rhd2000EvalBoard::PortC, Rhd2000EvalBoard::AuxCmd1, 0); + board->evalBoard->selectAuxCommandBank(Rhd2000EvalBoard::PortD, Rhd2000EvalBoard::AuxCmd1, 0); + board->evalBoard->selectAuxCommandLength(Rhd2000EvalBoard::AuxCmd1, 0, 1); + + board->evalBoard->selectAuxCommandBank(Rhd2000EvalBoard::PortA, Rhd2000EvalBoard::AuxCmd3, + board->fastSettleEnabled ? 2 : 1); + board->evalBoard->selectAuxCommandBank(Rhd2000EvalBoard::PortB, Rhd2000EvalBoard::AuxCmd3, + board->fastSettleEnabled ? 2 : 1); + board->evalBoard->selectAuxCommandBank(Rhd2000EvalBoard::PortC, Rhd2000EvalBoard::AuxCmd3, + board->fastSettleEnabled ? 2 : 1); + board->evalBoard->selectAuxCommandBank(Rhd2000EvalBoard::PortD, Rhd2000EvalBoard::AuxCmd3, + board->fastSettleEnabled ? 2 : 1); + + if (board->fastTTLSettleEnabled) + { + board->evalBoard->enableExternalFastSettle(true); + } } \ No newline at end of file diff --git a/Source/Processors/DataThreads/RHD2000Thread.h b/Source/Processors/DataThreads/RHD2000Thread.h index 8ca1a68b55e408ad996ad7b87f4434ccce32a9a5..3cba97cd8e94ec4a2c066628e119526a8df74c99 100644 --- a/Source/Processors/DataThreads/RHD2000Thread.h +++ b/Source/Processors/DataThreads/RHD2000Thread.h @@ -44,6 +44,16 @@ class SourceNode; class RHDHeadstage; +class RHDImpedanceMeasure; + +struct ImpedanceData +{ + Array<int> streams; + Array<int> channels; + Array<float> magnitudes; + Array<float> phases; +}; + /** Communicates with the RHD2000 Evaluation Board from Intan Technologies @@ -53,8 +63,8 @@ class RHDHeadstage; */ class RHD2000Thread : public DataThread, public Timer - { + friend class RHDImpedanceMeasure; public: RHD2000Thread(SourceNode* sn); ~RHD2000Thread(); @@ -83,22 +93,11 @@ public: void setDSPOffset(bool state); int setNoiseSlicerLevel(int level); - void runImpedanceTest(Array<int>& stream, Array<int>& channel, Array<float>& magnitude, Array<float>& phase); void setFastTTLSettle(bool state, int channel); void setTTLoutputMode(bool state); void setDAChpf(float cutoff, bool enabled); void scanPorts(); - float updateImpedanceFrequency(float desiredImpedanceFreq, bool& impedanceFreqValid); - int loadAmplifierData(queue<Rhd2000DataBlock>& dataQueue, - int numBlocks, int numDataStreams); - void measureComplexAmplitude(std::vector<std::vector<std::vector<double>>>& measuredMagnitude, - std::vector<std::vector<std::vector<double>>>& measuredPhase, - int capIndex, int stream, int chipChannel, int numBlocks, - double sampleRate, double frequency, int numPeriods); - void amplitudeOfFreqComponent(double& realComponent, double& imagComponent, - const std::vector<double>& data, int startIndex, - int endIndex, double sampleRate, double frequency); int getNumEventChannels(); void enableAdcs(bool); @@ -136,12 +135,7 @@ private: Rhd2000Registers chipRegisters; Rhd2000DataBlock* dataBlock; - std::vector<std::vector<std::vector<double>>> amplifierPreFilter; - void factorOutParallelCapacitance(double& impedanceMagnitude, double& impedancePhase, - double frequency, double parasiticCapacitance); - void empiricalResistanceCorrection(double& impedanceMagnitude, double& impedancePhase, - double boardSampleRate); - int numChannels; + int numChannels; bool deviceFound; float thisSample[256]; @@ -227,4 +221,31 @@ private: JUCE_DECLARE_NON_COPYABLE_WITH_LEAK_DETECTOR(RHDHeadstage); }; +class RHDImpedanceMeasure : public Thread +{ +public: + RHDImpedanceMeasure(RHD2000Thread* b, ImpedanceData* d); + void run(); +private: + void measureComplexAmplitude(std::vector<std::vector<std::vector<double>>>& measuredMagnitude, + std::vector<std::vector<std::vector<double>>>& measuredPhase, + int capIndex, int stream, int chipChannel, int numBlocks, + double sampleRate, double frequency, int numPeriods); + void amplitudeOfFreqComponent(double& realComponent, double& imagComponent, + const std::vector<double>& data, int startIndex, + int endIndex, double sampleRate, double frequency); + float updateImpedanceFrequency(float desiredImpedanceFreq, bool& impedanceFreqValid); + void factorOutParallelCapacitance(double& impedanceMagnitude, double& impedancePhase, + double frequency, double parasiticCapacitance); + void empiricalResistanceCorrection(double& impedanceMagnitude, double& impedancePhase, + double boardSampleRate); + int loadAmplifierData(queue<Rhd2000DataBlock>& dataQueue, + int numBlocks, int numDataStreams); + + std::vector<std::vector<std::vector<double>>> amplifierPreFilter; + + ImpedanceData* data; + RHD2000Thread* board; +}; + #endif // __RHD2000THREAD_H_2C4CBD67__