spectrumgenerator.cpp
12.1 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
// SPDX-License-Identifier: GPL-2.0+
#define _USE_MATH_DEFINES
#include <cmath>
#include <QColor>
#include <QMutex>
#include <QTimer>
#include <fftw3.h>
#include "spectrumgenerator.h"
#include "glscope.h"
#include "settings.h"
#include "utils/printutils.h"
/// \brief Analyzes the data from the dso.
SpectrumGenerator::SpectrumGenerator(const DsoSettingsScope *scope, const DsoSettingsPostProcessing *postprocessing)
: scope(scope), postprocessing(postprocessing) {}
SpectrumGenerator::~SpectrumGenerator() {
if (lastWindowBuffer) fftw_free(lastWindowBuffer);
}
void SpectrumGenerator::process(PPresult *result) {
// Calculate frequencies and spectrums
for (ChannelID channel = 0; channel < result->channelCount(); ++channel) {
DataChannel *const channelData = result->modifyData(channel);
if (channelData->voltage.sample.empty()) {
// Clear unused channels
channelData->spectrum.interval = 0;
channelData->spectrum.sample.clear();
continue;
}
// Calculate new window
size_t sampleCount = channelData->voltage.sample.size();
if (!lastWindowBuffer || lastWindow != postprocessing->spectrumWindow || lastRecordLength != sampleCount) {
if (lastWindowBuffer) fftw_free(lastWindowBuffer);
lastWindowBuffer = fftw_alloc_real(sampleCount);
lastRecordLength = (unsigned)sampleCount;
unsigned int windowEnd = lastRecordLength - 1;
lastWindow = postprocessing->spectrumWindow;
switch (postprocessing->spectrumWindow) {
case Dso::WindowFunction::HAMMING:
for (unsigned int windowPosition = 0; windowPosition < lastRecordLength; ++windowPosition)
*(lastWindowBuffer + windowPosition) = 0.54 - 0.46 * cos(2.0 * M_PI * windowPosition / windowEnd);
break;
case Dso::WindowFunction::HANN:
for (unsigned int windowPosition = 0; windowPosition < lastRecordLength; ++windowPosition)
*(lastWindowBuffer + windowPosition) = 0.5 * (1.0 - cos(2.0 * M_PI * windowPosition / windowEnd));
break;
case Dso::WindowFunction::COSINE:
for (unsigned int windowPosition = 0; windowPosition < lastRecordLength; ++windowPosition)
*(lastWindowBuffer + windowPosition) = sin(M_PI * windowPosition / windowEnd);
break;
case Dso::WindowFunction::LANCZOS:
for (unsigned int windowPosition = 0; windowPosition < lastRecordLength; ++windowPosition) {
double sincParameter = (2.0 * windowPosition / windowEnd - 1.0) * M_PI;
if (sincParameter == 0)
*(lastWindowBuffer + windowPosition) = 1;
else
*(lastWindowBuffer + windowPosition) = sin(sincParameter) / sincParameter;
}
break;
case Dso::WindowFunction::BARTLETT:
for (unsigned int windowPosition = 0; windowPosition < lastRecordLength; ++windowPosition)
*(lastWindowBuffer + windowPosition) =
2.0 / windowEnd * (windowEnd / 2 - std::abs((double)(windowPosition - windowEnd / 2.0)));
break;
case Dso::WindowFunction::TRIANGULAR:
for (unsigned int windowPosition = 0; windowPosition < lastRecordLength; ++windowPosition)
*(lastWindowBuffer + windowPosition) =
2.0 / lastRecordLength *
(lastRecordLength / 2 - std::abs((double)(windowPosition - windowEnd / 2.0)));
break;
case Dso::WindowFunction::GAUSS: {
double sigma = 0.4;
for (unsigned int windowPosition = 0; windowPosition < lastRecordLength; ++windowPosition)
*(lastWindowBuffer + windowPosition) =
exp(-0.5 * pow(((windowPosition - windowEnd / 2) / (sigma * windowEnd / 2)), 2));
} break;
case Dso::WindowFunction::BARTLETTHANN:
for (unsigned int windowPosition = 0; windowPosition < lastRecordLength; ++windowPosition)
*(lastWindowBuffer + windowPosition) = 0.62 -
0.48 * std::abs((double)(windowPosition / windowEnd - 0.5)) -
0.38 * cos(2.0 * M_PI * windowPosition / windowEnd);
break;
case Dso::WindowFunction::BLACKMAN: {
double alpha = 0.16;
for (unsigned int windowPosition = 0; windowPosition < lastRecordLength; ++windowPosition)
*(lastWindowBuffer + windowPosition) = (1 - alpha) / 2 -
0.5 * cos(2.0 * M_PI * windowPosition / windowEnd) +
alpha / 2 * cos(4.0 * M_PI * windowPosition / windowEnd);
} break;
// case Dso::WindowFunction::WINDOW_KAISER:
// TODO WINDOW_KAISER
// double alpha = 3.0;
// for(unsigned int windowPosition = 0; windowPosition <
// lastRecordLength; ++windowPosition)
//*(window + windowPosition) = ;
// break;
case Dso::WindowFunction::NUTTALL:
for (unsigned int windowPosition = 0; windowPosition < lastRecordLength; ++windowPosition)
*(lastWindowBuffer + windowPosition) = 0.355768 -
0.487396 * cos(2 * M_PI * windowPosition / windowEnd) +
0.144232 * cos(4 * M_PI * windowPosition / windowEnd) -
0.012604 * cos(6 * M_PI * windowPosition / windowEnd);
break;
case Dso::WindowFunction::BLACKMANHARRIS:
for (unsigned int windowPosition = 0; windowPosition < lastRecordLength; ++windowPosition)
*(lastWindowBuffer + windowPosition) = 0.35875 -
0.48829 * cos(2 * M_PI * windowPosition / windowEnd) +
0.14128 * cos(4 * M_PI * windowPosition / windowEnd) -
0.01168 * cos(6 * M_PI * windowPosition / windowEnd);
break;
case Dso::WindowFunction::BLACKMANNUTTALL:
for (unsigned int windowPosition = 0; windowPosition < lastRecordLength; ++windowPosition)
*(lastWindowBuffer + windowPosition) = 0.3635819 -
0.4891775 * cos(2 * M_PI * windowPosition / windowEnd) +
0.1365995 * cos(4 * M_PI * windowPosition / windowEnd) -
0.0106411 * cos(6 * M_PI * windowPosition / windowEnd);
break;
case Dso::WindowFunction::FLATTOP:
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);
break;
default: // Dso::WINDOW_RECTANGULAR
for (unsigned int windowPosition = 0; windowPosition < lastRecordLength; ++windowPosition)
*(lastWindowBuffer + windowPosition) = 1.0;
}
}
// Set sampling interval
channelData->spectrum.interval = 1.0 / channelData->voltage.interval / sampleCount;
// Number of real/complex samples
unsigned int dftLength = sampleCount / 2;
// Reallocate memory for samples if the sample count has changed
channelData->spectrum.sample.resize(sampleCount);
// Create sample buffer and apply window
std::unique_ptr<double[]> windowedValues = std::unique_ptr<double[]>(new double[sampleCount]);
for (unsigned int position = 0; position < sampleCount; ++position)
windowedValues[position] = lastWindowBuffer[position] * channelData->voltage.sample[position];
{
// 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 an autocorrelation to get the frequency of the signal
std::unique_ptr<double[]> conjugateComplex = std::move(windowedValues);
// Real values
unsigned int position;
double correctionFactor = 1.0 / dftLength / dftLength;
conjugateComplex[0] = (channelData->spectrum.sample[0] * channelData->spectrum.sample[0]) * correctionFactor;
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
conjugateComplex[dftLength] =
(channelData->spectrum.sample[dftLength] * channelData->spectrum.sample[dftLength]) * correctionFactor;
for (++position; position < sampleCount; ++position) conjugateComplex[position] = 0;
// 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);
fftw_execute(fftPlan);
fftw_destroy_plan(fftPlan);
// Get the frequency from the correlation results
double minimumCorrelation = correlation[0];
double peakCorrelation = 0;
unsigned int peakPosition = 0;
for (unsigned int position = 1; position < sampleCount / 2; ++position) {
if (correlation[position] > peakCorrelation && correlation[position] > minimumCorrelation * 2) {
peakCorrelation = correlation[position];
peakPosition = position;
} else if (correlation[position] < minimumCorrelation)
minimumCorrelation = correlation[position];
}
correlation.reset(nullptr);
// Calculate the frequency in Hz
if (peakPosition)
channelData->frequency = 1.0 / (channelData->voltage.interval * peakPosition);
else
channelData->frequency = 0;
// Finally calculate the real spectrum if we want it
if (scope->spectrum[channel].used) {
// Convert values into dB (Relative to the reference level)
double offset = 60 - postprocessing->spectrumReference - 20 * log10(dftLength);
double offsetLimit = postprocessing->spectrumLimit - postprocessing->spectrumReference;
for (std::vector<double>::iterator spectrumIterator = channelData->spectrum.sample.begin();
spectrumIterator != channelData->spectrum.sample.end(); ++spectrumIterator) {
double value = 20 * log10(fabs(*spectrumIterator)) + offset;
// Check if this value has to be limited
if (offsetLimit > value) value = offsetLimit;
*spectrumIterator = value;
}
}
}
}