Skip to content

Commit

Permalink
Correct spectrum display: use magnitude instead of real part. Limit m…
Browse files Browse the repository at this point in the history 
…ax frequencybase according to samplerate

Signed-off-by: Martin <Ho-Ro@users.noreply.github.com>
  • Loading branch information
@Ho-Ro
Ho-Ro committed May 24, 2019
1 parent d78d3e2 commit 229f6e9
Show file tree
Hide file tree
Showing 2 changed files with 34 additions and 17 deletions.
5 changes: 5 additions & 0 deletions openhantek/src/docks/HorizontalDock.cpp
View file
Expand Up @@ -14,6 +14,7 @@
#include "HorizontalDock.h"
#include "dockwindows.h"

#include "viewconstants.h"
#include "scopesettings.h"
#include "sispinbox.h"
#include "utils/printutils.h"
Expand Down Expand Up @@ -124,6 +125,10 @@ double HorizontalDock::setSamplerate(double samplerate) {
//printf( "HD::setSamplerate( %g )\n", samplerate );
QSignalBlocker blocker(samplerateSiSpinBox);
samplerateSiSpinBox->setValue(samplerate);
double maxFreqBase = samplerate / DIVS_TIME / 2;
frequencybaseSiSpinBox->setMaximum( maxFreqBase );
if (frequencybaseSiSpinBox->value() > maxFreqBase )
setFrequencybase( maxFreqBase );
return samplerateSiSpinBox->value();
}

Expand Down
46 changes: 29 additions & 17 deletions openhantek/src/post/spectrumgenerator.cpp
View file
Expand Up @@ -124,12 +124,12 @@ void SpectrumGenerator::process(PPresult *result) {
0.1365995 * cos(4 * M_PI * windowPosition / windowEnd) -
0.0106411 * cos(6 * M_PI * windowPosition / windowEnd);
break;
case Dso::WindowFunction::FLATTOP:
case Dso::WindowFunction::FLATTOP: // wikipedia.de
for (unsigned int windowPosition = 0; windowPosition < lastRecordLength; ++windowPosition)
*(lastWindowBuffer + windowPosition) = 1.0 - 1.93 * cos(2 * M_PI * windowPosition / windowEnd) +
1.29 * cos(4 * M_PI * windowPosition / windowEnd) -
0.388 * cos(6 * M_PI * windowPosition / windowEnd) +
0.032 * cos(8 * M_PI * windowPosition / windowEnd);
0.028 * cos(8 * M_PI * windowPosition / windowEnd);
break;
default: // Dso::WINDOW_RECTANGULAR
for (unsigned int windowPosition = 0; windowPosition < lastRecordLength; ++windowPosition)
Expand Down Expand Up @@ -168,15 +168,15 @@ void SpectrumGenerator::process(PPresult *result) {
channelData->ac = sqrt( ac2 ); // rms of AC component
channelData->rms = sqrt( dc * dc + ac2 ); // total rms = U eff

{
// Do discrete real to half-complex transformation
/// \todo Check if record length is multiple of 2
/// \todo Reuse plan and use FFTW_MEASURE to get fastest algorithm
fftw_plan fftPlan = fftw_plan_r2r_1d(sampleCount, windowedValues.get(),
&channelData->spectrum.sample.front(), FFTW_R2HC, FFTW_ESTIMATE);
fftw_execute(fftPlan);
fftw_destroy_plan(fftPlan);
}
// Do discrete real to half-complex transformation
/// \todo Check if record length is multiple of 2
/// \todo Reuse plan and use FFTW_MEASURE to get fastest algorithm

fftw_plan fftPlan;
fftPlan = fftw_plan_r2r_1d(sampleCount, windowedValues.get(),
&channelData->spectrum.sample.front(), FFTW_R2HC, FFTW_ESTIMATE);
fftw_execute(fftPlan);
fftw_destroy_plan(fftPlan);

// Do an autocorrelation to get the frequency of the signal
// HORO:
Expand All @@ -190,22 +190,34 @@ void SpectrumGenerator::process(PPresult *result) {
// Real values
unsigned int position;
double correctionFactor = 1.0 / dftLength / dftLength;
// correct the (half-)compled values in spectrum (1st part real forward), (2nd part imag backwards) -> magnitude
std::vector<double>::iterator fwdIterator = channelData->spectrum.sample.begin();
std::vector<double>::iterator backIterator = channelData->spectrum.sample.end();
conjugateComplex[0] = (channelData->spectrum.sample[0] * channelData->spectrum.sample[0]) * correctionFactor;
for (position = 1; position < dftLength; ++position)
fwdIterator++;
*backIterator-- = 0; // clear mirrored 2nd half of spectrum
for (position = 1; position < dftLength; ++position) {
conjugateComplex[position] =
(channelData->spectrum.sample[position] * channelData->spectrum.sample[position] +
channelData->spectrum.sample[sampleCount - position] *
channelData->spectrum.sample[sampleCount - position]) *
correctionFactor;
// Complex values, all zero for autocorrelation
// convert complex to magnitude
*fwdIterator++ = sqrt( *fwdIterator * *fwdIterator + *backIterator * *backIterator );
*backIterator-- = 0; // clear mirrored 2nd half of spectrum
}
conjugateComplex[dftLength] =
(channelData->spectrum.sample[dftLength] * channelData->spectrum.sample[dftLength]) * correctionFactor;
for (++position; position < sampleCount; ++position) conjugateComplex[position] = 0;
// Complex values, all zero for autocorrelation
for (++position; position < sampleCount; ++position) {
conjugateComplex[position] = 0;
}
// clear mirrored 2nd half of spectrum
channelData->spectrum.sample.resize( dftLength );

// Do half-complex to real inverse transformation
std::unique_ptr<double[]> correlation = std::unique_ptr<double[]>(new double[sampleCount]);
fftw_plan fftPlan =
fftw_plan_r2r_1d(sampleCount, conjugateComplex.get(), correlation.get(), FFTW_HC2R, FFTW_ESTIMATE);
fftPlan = fftw_plan_r2r_1d(sampleCount, conjugateComplex.get(), correlation.get(), FFTW_HC2R, FFTW_ESTIMATE);
fftw_execute(fftPlan);
fftw_destroy_plan(fftPlan);

Expand All @@ -215,7 +227,7 @@ void SpectrumGenerator::process(PPresult *result) {
unsigned int peakPosition = 0;

for (unsigned int position = 1; position < sampleCount / 2; ++position) {
if ( correlation[position] > 1.5 * peakCorrelation && correlation[position] > minimumCorrelation ) {
if ( correlation[position] > peakCorrelation && correlation[position] > minimumCorrelation ) {
peakCorrelation = correlation[position];
peakPosition = position;
} else if (correlation[position] < minimumCorrelation)
Expand Down

0 comments on commit 229f6e9

Please sign in to comment.