1 /***************************************************************************
2 psd.cpp: Power Spectra for KST
3 -------------------
4 begin : Fri Feb 10 2002
5 copyright : (C) 2002 by C. Barth Netterfield
6 email : netterfield@astro.utoronto.ca
7 ***************************************************************************/
8
9 /***************************************************************************
10 * *
11 * This program is free software; you can redistribute it and/or modify *
12 * it under the terms of the GNU General Public License as published by *
13 * the Free Software Foundation; either version 2 of the License, or *
14 * (at your option) any later version. *
15 * *
16 ***************************************************************************/
17
18 /** A class for handling power spectra for kst
19 *@author C. Barth Netterfield
20 */
21
22 #include "psd.h"
23
24 #include <assert.h>
25 #include <math.h>
26
27 #include <QXmlStreamWriter>
28
29
30 #include <qdebug.h>
31
32 #include "dialoglauncher.h"
33 #include "datacollection.h"
34 #include "debug.h"
35 #include "psdcalculator.h"
36 #include "objectstore.h"
37
38 #include "dataobjectscriptinterface.h"
39
40 using namespace std;
41
42 extern "C" void rdft(int n, int isgn, double *a);
43
44 namespace Kst {
45
46 const QString PSD::staticTypeString = "Power Spectrum";
47 const QString PSD::staticTypeTag = "powerspectrum";
48
49 static const QLatin1String& INVECTOR = QLatin1String("I");
50 static const QLatin1String& SVECTOR = QLatin1String("S");
51 static const QLatin1String& FVECTOR = QLatin1String("F");
52
53 #define KSTPSDMAXLEN 27
54
PSD(ObjectStore * store)55 PSD::PSD(ObjectStore *store)
56 : DataObject(store) {
57 _changed = true;
58 _typeString = staticTypeString;
59 _type = "PowerSpectrum";
60 _initializeShortName();
61
62 Q_ASSERT(store);
63 VectorPtr ov = store->createObject<Vector>();
64 ov->setProvider(this);
65 ov->setSlaveName("f");
66 ov->resize(2);
67 _fVector = _outputVectors.insert(FVECTOR, ov).value();
68
69 ov = store->createObject<Vector>();
70 ov->setProvider(this);
71 ov->setSlaveName("psd");
72 ov->resize(2);
73 _sVector = _outputVectors.insert(SVECTOR, ov).value();
74 }
75
_initializeShortName()76 void PSD::_initializeShortName() {
77 _shortName = 'S'+QString::number(_psdnum);
78 if (_psdnum>max_psdnum)
79 max_psdnum = _psdnum;
80 _psdnum++;
81 }
82
83
createScriptInterface()84 ScriptInterface* PSD::createScriptInterface() {
85 return new SpectrumSI(this);
86 }
87
88
change(VectorPtr in_V,double in_freq,bool in_average,int in_averageLen,bool in_apodize,bool in_removeMean,const QString & in_VUnits,const QString & in_RUnits,ApodizeFunction in_apodizeFxn,double in_gaussianSigma,PSDType in_output)89 void PSD::change(VectorPtr in_V,
90 double in_freq, bool in_average, int in_averageLen, bool in_apodize,
91 bool in_removeMean, const QString& in_VUnits, const QString& in_RUnits,
92 ApodizeFunction in_apodizeFxn, double in_gaussianSigma, PSDType in_output) {
93
94 if (in_V) {
95 _inputVectors[INVECTOR] = in_V;
96 }
97 _Frequency = in_freq;
98 _Average = in_average;
99 _Apodize = in_apodize;
100 _apodizeFxn = in_apodizeFxn;
101 _gaussianSigma = in_gaussianSigma;
102 _prevOutput = PSDUndefined;
103 _RemoveMean = in_removeMean;
104 _vectorUnits = in_VUnits;
105 _rateUnits = in_RUnits;
106 _Output = in_output;
107 _averageLength = in_averageLen;
108
109 _last_n_subsets = 0;
110 _last_n_new = 0;
111 _last_n_new = 0;
112
113 _PSDLength = 1;
114
115 _fVector->resize(_PSDLength);
116 _sVector->resize(_PSDLength);
117
118 _changed = true;
119 updateVectorLabels();
120 }
121
122
~PSD()123 PSD::~PSD() {
124 _sVector = 0L;
125 _fVector = 0L;
126 }
127
128
curveHints() const129 const CurveHintList *PSD::curveHints() const {
130 _curveHints->clear();
131 _curveHints->append(new CurveHint(tr("PSD Curve"), _fVector->shortName(),
132 _sVector->shortName()));
133 return _curveHints;
134 }
135
136
internalUpdate()137 void PSD::internalUpdate() {
138 writeLockInputsAndOutputs();
139
140 VectorPtr iv = _inputVectors[INVECTOR];
141
142 const int v_len = iv->length();
143
144 _last_n_new += iv->numNew();
145 assert(_last_n_new >= 0);
146
147 int n_subsets = (v_len)/_PSDLength;
148
149 // determine if the PSD needs to be updated.
150 // if not using averaging, then we need at least _PSDLength/16 new data points.
151 // if averaging, then we want enough new data for a complete subset.
152 // ... unless we are counting from end at fixed length (scrolling data).
153 bool scrolling_data = (_last_n == iv->length());
154 if ( (!_changed) && ((_last_n_new < _PSDLength/16) ||
155 (_Average && scrolling_data && (_last_n_new < _PSDLength/16)) ||
156 (_Average && !scrolling_data && (n_subsets - _last_n_subsets < 1))) &&
157 iv->length() != iv->numNew()) {
158 unlockInputsAndOutputs();
159 return;
160 }
161
162 _changed = false;
163
164 _adjustLengths();
165
166 double *psd = _sVector->raw_V_ptr();
167 double *f = _fVector->raw_V_ptr();
168
169 int i_samp;
170 for (i_samp = 0; i_samp < _PSDLength; ++i_samp) {
171 f[i_samp] = i_samp * 0.5 * _Frequency / (_PSDLength - 1);
172 }
173 //f[0] = -1E-280; // really 0 (this shouldn't be needed...)
174
175 _psdCalculator.calculatePowerSpectrum(iv->noNanValue(), v_len, psd, _PSDLength, _RemoveMean, _Average, _averageLength, _Apodize, _apodizeFxn, _gaussianSigma, _Output, _Frequency);
176
177 _last_n_subsets = n_subsets;
178 _last_n_new = 0;
179 _last_n = iv->length();
180
181 updateVectorLabels();
182
183 // should be updated by the update manager
184 //_sVector->update();
185 //_fVector->update();
186
187 unlockInputsAndOutputs();
188
189 return;
190 }
191
192
_adjustLengths()193 void PSD::_adjustLengths() {
194 int nPSDLen = PSDCalculator::calculateOutputVectorLength(_inputVectors[INVECTOR]->length(), _Average, _averageLength);
195
196 if (_PSDLength != nPSDLen) {
197 _sVector->resize(nPSDLen);
198 _fVector->resize(nPSDLen);
199
200 if ( (_sVector->length() == nPSDLen) && (_fVector->length() == nPSDLen) ) {
201 _PSDLength = nPSDLen;
202 } else {
203 Debug::self()->log(tr("Attempted to create a PSD that used all memory."), Debug::Error);
204 }
205
206 _last_n_subsets = 0;
207 _changed = true;
208 }
209 }
210
211
save(QXmlStreamWriter & s)212 void PSD::save(QXmlStreamWriter &s) {
213 s.writeStartElement(staticTypeTag);
214 s.writeAttribute("vector", _inputVectors[INVECTOR]->Name());
215 s.writeAttribute("samplerate", QString::number(_Frequency));
216 s.writeAttribute("gaussiansigma", QString::number(_gaussianSigma));
217 s.writeAttribute("average", QVariant(_Average).toString());
218 s.writeAttribute("fftlength", QString::number(int(ceil(log(double(_PSDLength*2)) / log(2.0)))));
219 s.writeAttribute("removemean", QVariant(_RemoveMean).toString());
220 s.writeAttribute("apodize", QVariant(_Apodize).toString());
221 s.writeAttribute("apodizefunction", QString::number(_apodizeFxn));
222 s.writeAttribute("vectorunits", _vectorUnits);
223 s.writeAttribute("rateunits", _rateUnits);
224 s.writeAttribute("outputtype", QString::number(_Output));
225 saveNameInfo(s, VECTORNUM|PSDNUM|SCALARNUM);
226
227 s.writeEndElement();
228 }
229
230
apodize() const231 bool PSD::apodize() const {
232 return _Apodize;
233 }
234
235
setApodize(bool in_apodize)236 void PSD::setApodize(bool in_apodize) {
237 _Apodize = in_apodize;
238 _changed = true;
239 }
240
241
removeMean() const242 bool PSD::removeMean() const {
243 return _RemoveMean;
244 }
245
246
setRemoveMean(bool in_removeMean)247 void PSD::setRemoveMean(bool in_removeMean) {
248 _RemoveMean = in_removeMean;
249 _changed = true;
250 }
251
252
average() const253 bool PSD::average() const {
254 return _Average;
255 }
256
257
setAverage(bool in_average)258 void PSD::setAverage(bool in_average) {
259 _Average = in_average;
260 _changed = true;
261 }
262
263
frequency() const264 double PSD::frequency() const {
265 return _Frequency;
266 }
267
268
setFrequency(double in_frequency)269 void PSD::setFrequency(double in_frequency) {
270 if (in_frequency > 0.0) {
271 _Frequency = in_frequency;
272 } else {
273 _Frequency = 1.0;
274 }
275 _changed = true;
276 }
277
278
length() const279 int PSD::length() const {
280 return _averageLength;
281 }
282
283
setLength(int in_length)284 void PSD::setLength(int in_length) {
285 if (in_length != _averageLength) {
286 _averageLength = in_length;
287 }
288 _changed = true;
289 }
290
291
output() const292 PSDType PSD::output() const {
293 return _Output;
294 }
295
296
setOutput(PSDType in_output)297 void PSD::setOutput(PSDType in_output) {
298 if (in_output != _Output) {
299 _Output = in_output;
300 }
301 _changed = true;
302 }
303
304
vector() const305 VectorPtr PSD::vector() const {
306 return _inputVectors[INVECTOR];
307 }
308
309
setVector(VectorPtr new_v)310 void PSD::setVector(VectorPtr new_v) {
311 Q_ASSERT(myLockStatus() == KstRWLock::WRITELOCKED);
312
313 VectorPtr v = _inputVectors[INVECTOR];
314 if (v) {
315 if (v == new_v) {
316 return;
317 }
318 v->unlock();
319 }
320
321 _inputVectors.remove(INVECTOR);
322 new_v->writeLock();
323 _inputVectors[INVECTOR] = new_v;
324 _changed = true;
325 }
326
327
slaveVectorsUsed() const328 bool PSD::slaveVectorsUsed() const {
329 return true;
330 }
331
332
propertyString() const333 QString PSD::propertyString() const {
334 return tr("PSD: %1").arg(_inputVectors[INVECTOR]->Name());
335 }
336
337
showNewDialog()338 void PSD::showNewDialog() {
339 DialogLauncher::self()->showPowerSpectrumDialog();
340 }
341
342
showEditDialog()343 void PSD::showEditDialog() {
344 DialogLauncher::self()->showPowerSpectrumDialog(this);
345 }
346
347
vectorUnits() const348 const QString& PSD::vectorUnits() const {
349 return _vectorUnits;
350 }
351
352
setVectorUnits(const QString & units)353 void PSD::setVectorUnits(const QString& units) {
354 _vectorUnits = units;
355 }
356
357
rateUnits() const358 const QString& PSD::rateUnits() const {
359 return _rateUnits;
360 }
361
362
setRateUnits(const QString & units)363 void PSD::setRateUnits(const QString& units) {
364 _rateUnits = units;
365 }
366
367
apodizeFxn() const368 ApodizeFunction PSD::apodizeFxn() const {
369 return _apodizeFxn;
370 }
371
372
setApodizeFxn(ApodizeFunction in_apodizeFxn)373 void PSD::setApodizeFxn(ApodizeFunction in_apodizeFxn) {
374 if (_apodizeFxn != in_apodizeFxn) {
375 _apodizeFxn = in_apodizeFxn;
376 }
377 _changed = true;
378 }
379
380
gaussianSigma() const381 double PSD::gaussianSigma() const {
382 return _gaussianSigma;
383 }
384
385
setGaussianSigma(double in_gaussianSigma)386 void PSD::setGaussianSigma(double in_gaussianSigma) {
387 if (_gaussianSigma != in_gaussianSigma) {
388 _gaussianSigma = in_gaussianSigma;
389 }
390 _changed = true;
391 }
392
393
makeDuplicate() const394 DataObjectPtr PSD::makeDuplicate() const {
395
396 PSDPtr powerspectrum = store()->createObject<PSD>();
397 Q_ASSERT(powerspectrum);
398
399 powerspectrum->writeLock();
400 powerspectrum->setVector(_inputVectors[INVECTOR]);
401 powerspectrum->setFrequency(_Frequency);
402 powerspectrum->setAverage(_Average);
403 powerspectrum->setLength(_averageLength);
404 powerspectrum->setApodize(_Apodize);
405 powerspectrum->setRemoveMean(_RemoveMean);
406 powerspectrum->setVectorUnits(_vectorUnits);
407 powerspectrum->setRateUnits(_rateUnits);
408 powerspectrum->setApodizeFxn(_apodizeFxn);
409 powerspectrum->setGaussianSigma(_gaussianSigma);
410 powerspectrum->setOutput(_Output);
411 if (descriptiveNameIsManual()) {
412 powerspectrum->setDescriptiveName(descriptiveName());
413 }
414 powerspectrum->registerChange();
415 powerspectrum->unlock();
416
417 return DataObjectPtr(powerspectrum);
418 }
419
420
updateVectorLabels()421 void PSD::updateVectorLabels() {
422 LabelInfo label_info;
423
424 switch (_Output) {
425 default:
426 case 0: // amplitude spectral density (default) [V/Hz^1/2]
427 label_info.quantity = tr("Spectral Density");
428 if (_vectorUnits.isEmpty() || _rateUnits.isEmpty()) {
429 label_info.units.clear();
430 } else {
431 label_info.units = QString("%1/%2^{1/2}").arg(_vectorUnits).arg(_rateUnits);
432 }
433 break;
434 case 1: // power spectral density [V^2/Hz]
435 label_info.quantity = tr("PSD");
436 if (_vectorUnits.isEmpty() || _rateUnits.isEmpty()) {
437 label_info.units.clear();
438 } else {
439 label_info.units = QString("%1^2/%2").arg(_vectorUnits).arg(_rateUnits);
440 }
441 break;
442 case 2: // amplitude spectrum [V]
443 label_info.quantity = tr("Amplitude Spectrum");
444 if (_vectorUnits.isEmpty() || _rateUnits.isEmpty()) {
445 label_info.units.clear();
446 } else {
447 label_info.units = QString("%1").arg(_vectorUnits);
448 }
449 break;
450 case 3: // power spectrum [V^2]
451 label_info.quantity = tr("Power Spectrum");
452 if (_vectorUnits.isEmpty() || _rateUnits.isEmpty()) {
453 label_info.units.clear();
454 } else {
455 label_info.units = QString("%1^2").arg(_vectorUnits);
456 }
457 break;
458 }
459 label_info.name.clear();
460 _sVector->setLabelInfo(label_info);
461
462 label_info.quantity = tr("Frequency");
463 label_info.units = _rateUnits;
464 _fVector->setLabelInfo(label_info);
465
466 label_info.quantity.clear();
467 label_info.units.clear();
468 label_info.name = _inputVectors[INVECTOR]->labelInfo().name;
469 label_info.file = _inputVectors[INVECTOR]->labelInfo().file;
470 _sVector->setTitleInfo(label_info);
471
472 }
473
_automaticDescriptiveName() const474 QString PSD::_automaticDescriptiveName() const {
475 return vector()->descriptiveName();
476 }
477
descriptionTip() const478 QString PSD::descriptionTip() const {
479 QString tip;
480
481 tip = tr("Spectrum: %1\n FFT Length: 2^%2").arg(Name()).arg(length());
482
483 if (average() || apodize() || removeMean()) {
484 tip += "\n ";
485 if (average()) tip += tr("Average; ");
486 if (apodize()) tip += tr("Apodize; ");
487 if (removeMean()) tip += tr("Remove Mean;");
488 }
489 tip += tr("\nInput: %1").arg(_inputVectors[INVECTOR]->descriptionTip());
490 return tip;
491 }
492 }
493 // vim: ts=2 sw=2 et
494