1 /***************************************************************************
2 psdcalculator.cpp: Power Spectra Calculator for KST
3 -------------------
4 begin : 2006
5 copyright : (C) 2006 by Kst
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 utility class for calculating power spectra
19 */
20
21 #include "psdcalculator.h"
22
23 #include <assert.h>
24
25
26
27 #include "debug.h"
28 #include "vector.h"
29
30 #include <qnamespace.h>
31 #include <math_kst.h>
32 #include "measuretime.h"
33
34 extern "C" void rdft(int n, int isgn, double *a);
35
36 #define PSDMINLEN 2
37 #define PSDMAXLEN 27
38
cabs2(double r,double i)39 inline double PSDCalculator::cabs2(double r, double i) {
40 return r*r + i*i;
41 }
42
PSDCalculator()43 PSDCalculator::PSDCalculator()
44 {
45 _a = 0L;
46 _b = 0L;
47 _w = 0L;
48
49 _fft_len = 0;
50
51 _prev_apodize_function = WindowUndefined;
52 _prev_gaussian_sigma = 1.0;
53 _prev_output_len = 0;
54 _prev_cross_spec = false;
55 }
56
57
~PSDCalculator()58 PSDCalculator::~PSDCalculator() {
59 delete[] _w;
60 _w = 0L;
61 delete[] _a;
62 _a = 0L;
63 delete[] _b;
64 _b = 0L;
65 }
66
updateWindowFxn(ApodizeFunction apodizeFxn,double gaussianSigma)67 void PSDCalculator::updateWindowFxn(ApodizeFunction apodizeFxn, double gaussianSigma) {
68 const double a = double(_fft_len) / 2.0;
69 double x;
70 double sW = 0.0;
71
72 switch (apodizeFxn) {
73 case WindowOriginal:
74 for (int i = 0; i < _fft_len; ++i) {
75 _w[i] = 1.0 - cos(2.0 * M_PI * double(i) / double(_fft_len));
76 sW += _w[i] * _w[i];
77 }
78 break;
79
80 case WindowBartlett:
81 for (int i = 0; i < _fft_len; ++i) {
82 x = i - a;
83 _w[i] = 1.0 - fabs(x) / a;
84 sW += _w[i] * _w[i];
85 }
86 break;
87
88 case WindowBlackman:
89 for (int i = 0; i < _fft_len; ++i) {
90 x = i - a;
91 _w[i] = 0.42 + 0.5 * cos(M_PI * x / a) + 0.08 * cos(2 * M_PI * x/a);
92 sW += _w[i] * _w[i];
93 }
94 break;
95
96 case WindowConnes:
97 for (int i = 0; i < _fft_len; ++i) {
98 x = i - a;
99 _w[i] = pow(static_cast<double>(1.0 - (x * x) / (a * a)), 2);
100 sW += _w[i] * _w[i];
101 }
102 break;
103
104 case WindowCosine:
105 for (int i = 0; i < _fft_len; ++i) {
106 x = i - a;
107 _w[i] = cos(M_PI * x / (2.0 * a));
108 sW += _w[i] * _w[i];
109 }
110 break;
111
112 case WindowGaussian:
113 for (int i = 0; i < _fft_len; ++i) {
114 x = i - a;
115 _w[i] = exp(-1.0 * x * x/(2.0 * gaussianSigma * gaussianSigma));
116 }
117 break;
118
119 case WindowHamming:
120 for (int i = 0; i < _fft_len; ++i) {
121 x = i - a;
122 _w[i] = 0.53836 + 0.46164 * cos(M_PI * x / a);
123 sW += _w[i] * _w[i];
124 }
125 break;
126
127 case WindowHann:
128 for (int i = 0; i < _fft_len; ++i) {
129 x = i - a;
130 _w[i] = pow(static_cast<double>(cos(M_PI * x/(2.0 * a))), 2);
131 sW += _w[i] * _w[i];
132 }
133 break;
134
135 case WindowWelch:
136 for (int i = 0; i < _fft_len; ++i) {
137 x = i - a;
138 _w[i] = 1.0 - x * x / (a * a);
139 sW += _w[i] * _w[i];
140 }
141 break;
142
143 case WindowUniform:
144 default:
145 for (int i = 0; i < _fft_len; ++i) {
146 _w[i] = 1.0;
147 }
148 sW = _fft_len;
149 break;
150 }
151
152 double norm = sqrt((double)_fft_len/sW); // normalization constant s.t. sum over (w^2) is _awLen
153
154 for (int i = 0; i < _fft_len; ++i) {
155 _w[i] *= norm;
156 }
157
158 _prev_apodize_function = apodizeFxn;
159 _prev_gaussian_sigma = gaussianSigma;
160 }
161
162
calculatePowerSpectrum(double const * input,int input_len,double * output,int output_len,bool removeMean,bool average,int average_len,bool apodize,ApodizeFunction apodize_function,double gaussian_sigma,PSDType output_type,double sampling_freq,double const * input2,int input2_len,double * output2)163 int PSDCalculator::calculatePowerSpectrum(
164 double const *input, int input_len,
165 double *output, int output_len,
166 bool removeMean,
167 bool average, int average_len,
168 bool apodize, ApodizeFunction apodize_function, double gaussian_sigma,
169 PSDType output_type, double sampling_freq,
170 double const *input2, int input2_len, double *output2) {
171
172 bool cross_spectra = false;
173
174 if ((input2_len == input_len) && input2) {
175 cross_spectra = true;
176 }
177
178 if ((input2) && (input2_len == input_len)) {
179 cross_spectra = true;
180 }
181
182 if (output_len != calculateOutputVectorLength(input_len, average, average_len)) {
183 Kst::Debug::self()->log(Kst::Debug::tr("in PSDCalculator::calculatePowerSpectrum: received output array with wrong length."), Kst::Debug::Error);
184 return -1;
185 }
186
187 if (output_len != _prev_output_len) {
188 delete[] _a;
189 delete[] _w;
190 delete[] _b;
191
192 _fft_len = output_len*2;
193 _prev_output_len = output_len;
194
195 _a = new double[_fft_len];
196 _b = new double[_fft_len];
197 _w = new double[_fft_len];
198
199 updateWindowFxn(apodize_function, gaussian_sigma);
200 }
201
202 if ( (_prev_apodize_function != apodize_function) || (_prev_gaussian_sigma != gaussian_sigma) ) {
203 updateWindowFxn(apodize_function, gaussian_sigma);
204 }
205
206 int currentCopyLen;
207 int nsamples = 0;
208 int i_samp;
209 int ioffset;
210
211 memset(output, 0, sizeof(double)*output_len); // initialize output.
212 if (cross_spectra) {
213 memset(output2, 0, sizeof(double)*output_len); // initialize complex output for xspectra.
214 }
215
216 // Mingw build could be 10 times slower (Gaussian apod, mostly 0 then?)
217 //MeasureTime time_in_rfdt("rdft()");
218
219 bool done = false;
220 for (int i_subset = 0; !done; i_subset++) {
221 ioffset = i_subset*output_len; //overlapping average => i_subset*outputLen
222
223 // only zero pad if we really have to. It is better to adjust the last chunk's
224 // overlap.
225 if (ioffset + _fft_len*5/4 < input_len) {
226 currentCopyLen = _fft_len; //will copy a complete window.
227 } else if (_fft_len<input_len) { // count the last one from the end.
228 ioffset = input_len-_fft_len - 1;
229 currentCopyLen = _fft_len; //will copy a complete window.
230 done = true;
231 } else {
232 currentCopyLen = input_len - ioffset; //will copy a partial window.
233 memset(&_a[currentCopyLen], 0, sizeof(double)*(_fft_len - currentCopyLen)); //zero the leftovers.
234 if (cross_spectra) {
235 memset(&_b[currentCopyLen], 0, sizeof(double)*(_fft_len - currentCopyLen)); //zero the leftovers.
236 }
237 done = true;
238 }
239
240 double mean = 0.0;
241 double mean2 = 0.0;
242
243 if (removeMean) {
244 for (i_samp = 0; i_samp < currentCopyLen; ++i_samp) {
245 mean += input[i_samp + ioffset];
246 }
247 mean /= (double)currentCopyLen;
248 if (cross_spectra) {
249 for (i_samp = 0; i_samp < currentCopyLen; ++i_samp) {
250 mean2 += input2[i_samp + ioffset];
251 }
252 mean2 /= (double)currentCopyLen;
253 }
254 }
255
256 // apply the PSD options (removeMean, apodize, etc.)
257 // separate cases for speed- although this shouldn't really matter- the rdft should be the most time consuming step by far for any large data set.
258 if (removeMean && apodize) {
259 for (i_samp = 0; i_samp < currentCopyLen; ++i_samp) {
260 _a[i_samp] = (input[i_samp + ioffset] - mean)*_w[i_samp];
261 }
262 } else if (removeMean) {
263 for (i_samp = 0; i_samp < currentCopyLen; ++i_samp) {
264 _a[i_samp] = input[i_samp + ioffset] - mean;
265 }
266 } else if (apodize) {
267 for (i_samp = 0; i_samp < currentCopyLen; ++i_samp) {
268 _a[i_samp] = input[i_samp + ioffset]*_w[i_samp];
269 }
270 } else {
271 for (i_samp = 0; i_samp < currentCopyLen; ++i_samp) {
272 _a[i_samp] = input[i_samp + ioffset];
273 }
274 }
275
276 if (cross_spectra) {
277 if (removeMean && apodize) {
278 for (i_samp = 0; i_samp < currentCopyLen; ++i_samp) {
279 _b[i_samp] = (input2[i_samp + ioffset] - mean2)*_w[i_samp];
280 }
281 } else if (removeMean) {
282 for (i_samp = 0; i_samp < currentCopyLen; ++i_samp) {
283 _b[i_samp] = input2[i_samp + ioffset] - mean2;
284 }
285 } else if (apodize) {
286 for (i_samp = 0; i_samp < currentCopyLen; ++i_samp) {
287 _b[i_samp] = input2[i_samp + ioffset]*_w[i_samp];
288 }
289 } else {
290 for (i_samp = 0; i_samp < currentCopyLen; ++i_samp) {
291 _b[i_samp] = input2[i_samp + ioffset];
292 }
293 }
294 }
295
296 nsamples += currentCopyLen;
297
298 #if !defined(__QNX__)
299 rdft(_fft_len, 1, _a); //real discrete fourier transorm on _a.
300 if (cross_spectra) {
301 rdft(_fft_len, 1, _b); //real discrete fourier transorm on _b.
302 }
303 #else
304 Q_ASSERT(0); // there is a linking problem when not compling with pch. . .
305 #endif
306
307 if (cross_spectra) {
308 output[0] += _a[0] * _b[0];
309 output2[0] = 0;
310 output[output_len-1] += _a[1] * _b[1];
311 output2[output_len-1] = 0;
312 for (i_samp = 1; i_samp < output_len - 1; ++i_samp) {
313 output[i_samp] += _a[i_samp*2] * _b[i_samp*2] +
314 _a[i_samp*2+1] * _b[i_samp*2+1];
315 output2[i_samp] = -_a[i_samp*2] * _b[i_samp*2+1] +
316 _a[i_samp*2+1] * _b[i_samp*2];
317 }
318 } else {
319 output[0] += _a[0] * _a[0];
320 output[output_len-1] += _a[1] * _a[1];
321 for (i_samp = 1; i_samp < output_len - 1; ++i_samp) {
322 output[i_samp] += cabs2(_a[i_samp * 2], _a[i_samp * 2 + 1]);
323 }
324 }
325 }
326
327 // FIXME: NORMALIZATION.
328 /* This normalization doesn't give the same results as the original KstPSD.
329
330 double frequencyStep = .5*(double)inputSamplingFreq/(double)(outputLen-1);
331
332 //normalization factors which were left out earlier for speed.
333 // - 2.0 for the negative frequencies which were neglected by the rdft //FIXME: double check.
334 // - /(_awLen*_awLen) for the constant Wss from numerical recipes in C. (ensure that the window function agrees with this.)
335 // - /i_subset to average the powers in all the subsets.
336 double norm = 2.0/(double)_awLen/(double)_awLen/(double)i_subset;
337 */
338
339 // original normalization
340 double frequencyStep = 2.0*(double)sampling_freq/(double)nsamples; //OLD value for frequencyStep.
341 double norm = 2.0/(double)nsamples*2.0/(double)nsamples; //OLD value for norm.
342
343 switch (output_type) {
344 default:
345 case PSDAmplitudeSpectralDensity: // amplitude spectral density (default) [V/Hz^1/2]
346 norm /= frequencyStep;
347 for (i_samp = 0; i_samp < output_len; ++i_samp) {
348 output[i_samp] = sqrt(output[i_samp]*norm);
349 }
350 break;
351
352 case PSDPowerSpectralDensity: // power spectral density [V^2/Hz]
353 norm /= frequencyStep;
354 for (i_samp = 0; i_samp < output_len; ++i_samp) {
355 output[i_samp] *= norm;
356 }
357 if (cross_spectra) {
358 for (i_samp = 0; i_samp < output_len; ++i_samp) {
359 output2[i_samp] *= norm;
360 }
361 }
362 break;
363
364 case PSDAmplitudeSpectrum: // amplitude spectrum [V]
365 for (i_samp = 0; i_samp < output_len; ++i_samp) {
366 output[i_samp] = sqrt(output[i_samp]*norm);
367 }
368 break;
369
370 case PSDPowerSpectrum: // power spectrum [V^2]
371 for (i_samp = 0; i_samp < output_len; ++i_samp) {
372 output[i_samp] *= norm;
373 }
374 if (cross_spectra) {
375 for (i_samp = 0; i_samp < output_len; ++i_samp) {
376 output2[i_samp] *= norm;
377 }
378 }
379 break;
380 }
381
382 return 0;
383 }
384
385
calculateOutputVectorLength(int inputLen,bool average,int averageLen)386 int PSDCalculator::calculateOutputVectorLength(int inputLen, bool average, int averageLen) {
387 int psdloglen;
388
389 if (average && pow(2.0, averageLen) < inputLen) {
390 psdloglen = averageLen;
391 } else {
392 psdloglen = int(ceil(log(double(inputLen)) / log(2.0)));
393 }
394
395 if (psdloglen < PSDMINLEN) {
396 psdloglen = PSDMINLEN;
397 }
398
399 if (psdloglen > PSDMAXLEN) {
400 psdloglen = PSDMAXLEN;
401 }
402
403 return int(pow(2.0, psdloglen - 1));
404 }
405
406 // vim: ts=2 sw=2 et
407