From 9bcef88b8d30d970f3a144b63f2f13832fb718e7 Mon Sep 17 00:00:00 2001 From: oliverhaag Date: Wed, 29 Sep 2010 18:02:20 +0000 Subject: [PATCH] Frequency detection using autocorrelation and spectrums in csv export --- openhantek/ChangeLog | 5 +++++ openhantek/src/dataanalyzer.cpp | 155 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------------------------------------------------------- openhantek/src/exporter.cpp | 12 ++++++++++++ openhantek/src/hantek/control.cpp | 4 ++-- openhantek/src/hantek/device.cpp | 3 ++- openhantek/src/settings.cpp | 4 ++-- 6 files changed, 117 insertions(+), 66 deletions(-) diff --git a/openhantek/ChangeLog b/openhantek/ChangeLog index 48f6918..d088cbf 100644 --- a/openhantek/ChangeLog +++ b/openhantek/ChangeLog @@ -97,3 +97,8 @@ * Some bugfixes for the DSO-5200 * Simplified trigger point calculation and correct handling for fast rate mode * Recalculating pretrigger position after samplerate changes + +2010-09-29 Oliver Haag +* Added spectrums to exported csv files +* Using autocorrelation for frequency detection now +* Some undefined value bugfixes diff --git a/openhantek/src/dataanalyzer.cpp b/openhantek/src/dataanalyzer.cpp index 5e7100c..d794379 100644 --- a/openhantek/src/dataanalyzer.cpp +++ b/openhantek/src/dataanalyzer.cpp @@ -184,55 +184,10 @@ void DataAnalyzer::run() { // Lower priority for spectrum calculation this->setPriority(QThread::LowPriority); - // Calculate peak-to-peak voltage and frequency - for(int channel = 0; channel < this->analyzedData.count(); channel++) { - if(this->settings->scope.voltage[channel].used && this->analyzedData[channel]->samples.voltage.sample) { - bool aboveTrigger = this->analyzedData[channel]->samples.voltage.sample[0] > this->settings->scope.voltage[channel].trigger; - double minimalVoltage, maximalVoltage; - unsigned int firstTrigger = 0, lastTrigger = 0; - int triggerCount = -1; - this->analyzedData[channel]->amplitude = 0; - minimalVoltage = maximalVoltage = this->analyzedData[channel]->samples.voltage.sample[0]; - - for(unsigned int position = 0; position < this->analyzedData[channel]->samples.voltage.count; position++) { - // Check trigger condition - if(aboveTrigger != (this->analyzedData[channel]->samples.voltage.sample[position] > this->settings->scope.voltage[channel].trigger)) { - aboveTrigger = this->analyzedData[channel]->samples.voltage.sample[position] > this->settings->scope.voltage[channel].trigger; - // We measure the time between two trigger slopes - if(aboveTrigger == (this->settings->scope.trigger.slope == Dso::SLOPE_POSITIVE)) { - triggerCount++; - if(triggerCount == 0) - firstTrigger = position; - else - this->analyzedData[channel]->amplitude += maximalVoltage - minimalVoltage; - minimalVoltage = this->analyzedData[channel]->samples.voltage.sample[position]; - maximalVoltage = this->analyzedData[channel]->samples.voltage.sample[position]; - lastTrigger = position; - } - } - else { // Get minimal and maximal voltage - if(this->analyzedData[channel]->samples.voltage.sample[position] < minimalVoltage) - minimalVoltage = this->analyzedData[channel]->samples.voltage.sample[position]; - else if(this->analyzedData[channel]->samples.voltage.sample[position] > maximalVoltage) - maximalVoltage = this->analyzedData[channel]->samples.voltage.sample[position]; - } - } - - // Calculate values - if(triggerCount >= 0) { - this->analyzedData[channel]->amplitude /= triggerCount; - this->analyzedData[channel]->frequency = (double) triggerCount / (lastTrigger - firstTrigger) / this->analyzedData[channel]->samples.voltage.interval; - } - else { - this->analyzedData[channel]->amplitude = maximalVoltage - minimalVoltage; - this->analyzedData[channel]->frequency = 0; - } - } - } - // Calculate spectrums + // Calculate frequencies, peak-to-peak voltages and spectrums for(int channel = 0; channel < this->analyzedData.count(); channel++) { - if(this->settings->scope.spectrum[channel].used && this->analyzedData[channel]->samples.voltage.sample) { + if(this->analyzedData[channel]->samples.voltage.sample) { // Calculate new window if(this->lastWindow != this->settings->scope.spectrumWindow || this->lastBufferSize != this->analyzedData[channel]->samples.voltage.count) { if(this->lastBufferSize != this->analyzedData[channel]->samples.voltage.count) { @@ -325,9 +280,12 @@ void DataAnalyzer::run() { // Set sampling interval this->analyzedData[channel]->samples.spectrum.interval = 1.0 / this->analyzedData[channel]->samples.voltage.interval / this->analyzedData[channel]->samples.voltage.count; + // Number of real/complex samples + unsigned int dftLength = this->analyzedData[channel]->samples.voltage.count / 2; + // Reallocate memory for samples if the sample count has changed - if(this->analyzedData[channel]->samples.spectrum.count != this->analyzedData[channel]->samples.voltage.count / 2) { - this->analyzedData[channel]->samples.spectrum.count = this->analyzedData[channel]->samples.voltage.count / 2; + if(this->analyzedData[channel]->samples.spectrum.count != dftLength) { + this->analyzedData[channel]->samples.spectrum.count = dftLength; if(this->analyzedData[channel]->samples.spectrum.sample) delete[] this->analyzedData[channel]->samples.spectrum.sample; this->analyzedData[channel]->samples.spectrum.sample = new double[this->analyzedData[channel]->samples.voltage.count]; @@ -339,23 +297,98 @@ void DataAnalyzer::run() { windowedValues[position] = this->window[position] * this->analyzedData[channel]->samples.voltage.sample[position]; // Do discrete real to half-complex transformation - // TODO: Reuse plan and use FFTW_MEASURE to get fastest algorithm + /// \todo Check if buffer size is multiple of 2 + /// \todo Reuse plan and use FFTW_MEASURE to get fastest algorithm fftw_plan fftPlan = fftw_plan_r2r_1d(this->analyzedData[channel]->samples.voltage.count, windowedValues, this->analyzedData[channel]->samples.spectrum.sample, FFTW_R2HC, FFTW_ESTIMATE); fftw_execute(fftPlan); fftw_destroy_plan(fftPlan); - // Deallocate sample buffer - delete[] windowedValues; + // Do an autocorrelation to get the frequency of the signal + double *conjugateComplex = windowedValues; // Reuse the windowedValues buffer - // Convert values into dB (Relative to the reference level) - double offset = 60 - this->settings->scope.spectrumReference - 20 * log10(sqrt(this->analyzedData[channel]->samples.spectrum.count)); - double offsetLimit = this->settings->scope.spectrumLimit - this->settings->scope.spectrumReference; - for(unsigned int position = 0; position < this->analyzedData[channel]->samples.spectrum.count; position++) { - this->analyzedData[channel]->samples.spectrum.sample[position] = 20 * log10(fabs(this->analyzedData[channel]->samples.spectrum.sample[position])) + offset; - - // Check if this value has to be limited - if(offsetLimit > this->analyzedData[channel]->samples.spectrum.sample[position]) - this->analyzedData[channel]->samples.spectrum.sample[position] = offsetLimit; + // Real values + unsigned int position; + double correctionFactor = 1.0 / dftLength / dftLength; + conjugateComplex[0] = (this->analyzedData[channel]->samples.spectrum.sample[0] * this->analyzedData[channel]->samples.spectrum.sample[0]) * correctionFactor; + for(position = 1; position < dftLength; position++) + conjugateComplex[position] = (this->analyzedData[channel]->samples.spectrum.sample[position] * this->analyzedData[channel]->samples.spectrum.sample[position] + this->analyzedData[channel]->samples.spectrum.sample[this->analyzedData[channel]->samples.voltage.count - position] * this->analyzedData[channel]->samples.spectrum.sample[this->analyzedData[channel]->samples.voltage.count - position]) * correctionFactor; + // Complex values, all zero for autocorrelation + conjugateComplex[dftLength] = (this->analyzedData[channel]->samples.spectrum.sample[dftLength] * this->analyzedData[channel]->samples.spectrum.sample[dftLength]) * correctionFactor; + for(position++; position < this->analyzedData[channel]->samples.voltage.count; position++) + conjugateComplex[position] = 0; + + // Do half-complex to real inverse transformation + double *correlation = new double[this->analyzedData[channel]->samples.voltage.count]; + fftPlan = fftw_plan_r2r_1d(this->analyzedData[channel]->samples.voltage.count, conjugateComplex, correlation, FFTW_HC2R, FFTW_ESTIMATE); + fftw_execute(fftPlan); + fftw_destroy_plan(fftPlan); + delete[] conjugateComplex; + + // Calculate peak-to-peak voltage + double minimalVoltage, maximalVoltage; + minimalVoltage = maximalVoltage = this->analyzedData[channel]->samples.voltage.sample[0]; + + for(unsigned int position = 1; position < this->analyzedData[channel]->samples.voltage.count; position++) { + if(this->analyzedData[channel]->samples.voltage.sample[position] < minimalVoltage) + minimalVoltage = this->analyzedData[channel]->samples.voltage.sample[position]; + else if(this->analyzedData[channel]->samples.voltage.sample[position] > maximalVoltage) + maximalVoltage = this->analyzedData[channel]->samples.voltage.sample[position]; + } + + this->analyzedData[channel]->amplitude = maximalVoltage - minimalVoltage; + + // Get the frequency from the correlation results + double correlationLimit = pow(sqrt(maximalVoltage - minimalVoltage) / 2, 4); + bool newPeak = false; // Ignore correlation without offset (position = 0) + double bestPeak = 0, lastPeak = 0; + unsigned int bestPeakPosition = 0, currentPeakPosition = 0; + + for(unsigned int position = 1; position < this->analyzedData[channel]->samples.voltage.count; position++) { + if(correlation[position] < correlationLimit) { + // Check if there was a good peak before + if(currentPeakPosition) { + // Is this really a better correlation and not just a secondary peak of the first one? + if(lastPeak > bestPeak * 1.2) { + bestPeak = lastPeak; + bestPeakPosition = currentPeakPosition; + } + currentPeakPosition = 0; + } + newPeak = true; + } + else if((currentPeakPosition || newPeak) && correlation[position] > lastPeak) { + // We want this peak, store it + lastPeak = correlation[position]; + currentPeakPosition = position; + newPeak = false; + } + } + delete[] correlation; + + // Check if there's a possible peak available that wasn't finished + if(currentPeakPosition && currentPeakPosition < this->analyzedData[channel]->samples.voltage.count - 1 && lastPeak > bestPeak * 1.2) { + bestPeak = lastPeak; + bestPeakPosition = currentPeakPosition; + } + + // Calculate the frequency in Hz + if(bestPeakPosition) + this->analyzedData[channel]->frequency = 1.0 / (this->analyzedData[channel]->samples.voltage.interval * bestPeakPosition); + else + this->analyzedData[channel]->frequency = 0; + + // Finally calculate the real spectrum if we want it + if(this->settings->scope.spectrum[channel].used) { + // Convert values into dB (Relative to the reference level) + double offset = 60 - this->settings->scope.spectrumReference - 20 * log10(dftLength); + double offsetLimit = this->settings->scope.spectrumLimit - this->settings->scope.spectrumReference; + for(unsigned int position = 0; position < this->analyzedData[channel]->samples.spectrum.count; position++) { + this->analyzedData[channel]->samples.spectrum.sample[position] = 20 * log10(fabs(this->analyzedData[channel]->samples.spectrum.sample[position])) + offset; + + // Check if this value has to be limited + if(offsetLimit > this->analyzedData[channel]->samples.spectrum.sample[position]) + this->analyzedData[channel]->samples.spectrum.sample[position] = offsetLimit; + } } } else if(this->analyzedData[channel]->samples.spectrum.sample) { diff --git a/openhantek/src/exporter.cpp b/openhantek/src/exporter.cpp index 5307b61..dabbdd3 100644 --- a/openhantek/src/exporter.cpp +++ b/openhantek/src/exporter.cpp @@ -384,6 +384,18 @@ bool Exporter::doExport() { // Finally a newline csvStream << '\n'; } + + if(this->settings->scope.spectrum[channel].used) { + // Start with channel name and the sample interval + csvStream << "\"" << this->settings->scope.spectrum[channel].name << "\"," << this->dataAnalyzer->data(channel)->samples.spectrum.interval; + + // And now all magnitudes in dB + for(unsigned int position = 0; position < this->dataAnalyzer->data(channel)->samples.spectrum.count; position++) + csvStream << "," << this->dataAnalyzer->data(channel)->samples.spectrum.sample[position]; + + // Finally a newline + csvStream << '\n'; + } } csvFile.close(); diff --git a/openhantek/src/hantek/control.cpp b/openhantek/src/hantek/control.cpp index 9b68574..9e5ce71 100644 --- a/openhantek/src/hantek/control.cpp +++ b/openhantek/src/hantek/control.cpp @@ -745,7 +745,7 @@ namespace Hantek { scalerId = gainId % 3; this->sampleRange[channel] = 0xff; break; - case 1: + default: /// \todo Use calibration data to get the DSO-5200 sample ranges if(gainId == GAIN_10MV) { scalerId = 1; @@ -909,7 +909,7 @@ namespace Hantek { maximum = 0xfd; break; - case 1: + default: // The range is the same as used for the offsets for 10 bit models minimum = ((unsigned short int) *((unsigned char *) &(this->channelLevels[channel][this->gain[channel]][OFFSET_START])) << 8) + *((unsigned char *) &(this->channelLevels[channel][this->gain[channel]][OFFSET_START]) + 1); maximum = ((unsigned short int) *((unsigned char *) &(this->channelLevels[channel][this->gain[channel]][OFFSET_START])) << 8) + *((unsigned char *) &(this->channelLevels[channel][this->gain[channel]][OFFSET_END]) + 1); diff --git a/openhantek/src/hantek/device.cpp b/openhantek/src/hantek/device.cpp index 006ffb7..e9aaf86 100644 --- a/openhantek/src/hantek/device.cpp +++ b/openhantek/src/hantek/device.cpp @@ -165,10 +165,11 @@ namespace Hantek { libusb_close(this->handle); ssize_t deviceCount = libusb_get_device_list(this->context, &deviceList); - if (deviceCount < 0) + if(deviceCount < 0) return tr("Failed to get device list: %3").arg(Helper::libUsbErrorString(errorCode)); // Iterate through all usb devices + this->model = MODEL_UNKNOWN; for(ssize_t deviceIterator = 0; deviceIterator < deviceCount; deviceIterator++) { device = deviceList[deviceIterator]; // Get device descriptor diff --git a/openhantek/src/settings.cpp b/openhantek/src/settings.cpp index d669d32..41364d3 100644 --- a/openhantek/src/settings.cpp +++ b/openhantek/src/settings.cpp @@ -56,8 +56,8 @@ DsoSettings::DsoSettings() { this->scope.trigger.special = false; // General this->scope.physicalChannels = 0; - this->scope.spectrumLimit = 0.0; - this->scope.spectrumReference = 20.0; + this->scope.spectrumLimit = -20.0; + this->scope.spectrumReference = 0.0; this->scope.spectrumWindow = Dso::WINDOW_HANN; -- libgit2 0.21.4