[Kst] kdeextragear-2/kst/kst

[Barth Netterfield]

CVS commit by netterfield: 

BUG:
Fix psd normalization and special case bugs.
The normalization is now such that an input time stream with a 1 hz sample rate of white
Gaussian noise with unit rms produces a psd with an rms amplitude of sqrt(2)/rtHz, from
0 to 0.5Hz, independent of apodization, mean removal, or averaging.

Approved by George


  M +36 -23    kstpsd.cpp   1.12
  M +3 -3      kstpsd.h   1.5


--- kdeextragear-2/kst/kst/kstpsd.cpp  #1.11:1.12
@@ -139,5 +139,5 @@ void KstPSD::commonConstructor(const QSt
   _a = new double[_ALen];
   _w = new double[_ALen];
-  GenW();
+  _WLen = 0;
 
   _last_f0 = 0;
@@ -174,21 +174,23 @@ KstPSD::~KstPSD() {
 }
 
-void KstPSD::GenW() {
+void KstPSD::GenW(int len) {
   int i;
   double sW = 0;
 
-  for (i = 0; i < _ALen; i++) {
-    //w[i] = 4;
-    _w[i] = (1 - cos(2*M_PI*double(i)/double(_ALen)));
+  if (len != _WLen) {
+    for (i = 0; i < len; i++) {
+      _w[i] = (1 - cos(2*M_PI*double(i)/double(len)));
     sW += _w[i] * _w[i];
   }
-  sW = sqrt(sW/double(_ALen));
-  for (i = 0; i < _ALen; i++) {
+    sW = sqrt(sW/double(len));
+    for (i = 0; i < len; i++) {
     _w[i] /= sW;
   }
+    _WLen = len;
+  }
 }
 
-inline double cabs(double r, double i) {
-  return sqrt(r*r + i*i);
+inline double cabs2(double r, double i) {
+  return (r*r + i*i);
 }
 
@@ -204,6 +206,8 @@ KstObject::UpdateType KstPSD::update(int
   double *psd, *f;
   double mean;
+  double nf;
   bool force = false;
   KstVectorPtr iv = _inputVectors[INVECTOR];
+  bool done;
 
   if (update_counter <= 0) {
@@ -217,4 +221,6 @@ KstObject::UpdateType KstPSD::update(int
   v_len = iv->sampleCount();
   n_subsets = v_len/_PSDLen+1;
+  if (v_len%_PSDLen==0) n_subsets--;
+
   _last_n_new += iv->numNew();
 
@@ -234,11 +240,16 @@ KstObject::UpdateType KstPSD::update(int
   }
 
-  for (i_subset = 0; i_subset < n_subsets; i_subset++) {
+  nf = 0;
+  done = false;
+
+  for (i_subset = 0; !done; i_subset++) {
     /* copy each chunk into a[] and find mean */
-    if (i_subset*_PSDLen + _ALen <= v_len) {
+    if (i_subset*_PSDLen + _ALen < v_len) {
       copyLen = _ALen;
     } else {
       copyLen = v_len - i_subset*_PSDLen;
+      done = true;
     }
+
     mean = 0;
     for (i_samp = 0; i_samp < copyLen; i_samp++) {
@@ -252,4 +263,5 @@ KstObject::UpdateType KstPSD::update(int
     }
 
+    if (apodize()) GenW(copyLen);
     /* Remove Mean and apodize */
     if (removeMean() && apodize()) {
@@ -266,4 +278,8 @@ KstObject::UpdateType KstPSD::update(int
       }
     }
+
+    nf += copyLen;
+
+
     for (;i_samp < _ALen; i_samp++) {
       _a[i_samp] = 0.0;
@@ -272,9 +288,8 @@ KstObject::UpdateType KstPSD::update(int
     /* fft a */
     rdft(_ALen, 1, _a);
-    /* sum each bin into psd[] */
-    psd[0] += _a[0];
-    psd[_PSDLen-1] += fabs(_a[1]);
+    psd[0] += _a[0]*_a[0];
+    psd[_PSDLen-1] += _a[1]*_a[1];
     for (i_samp=1; i_samp<_PSDLen-1; i_samp++) {
-      psd[i_samp]+= cabs(_a[i_samp*2], _a[i_samp*2+1]);
+      psd[i_samp]+= cabs2(_a[i_samp*2], _a[i_samp*2+1]);
     }
   }
@@ -284,8 +299,7 @@ KstObject::UpdateType KstPSD::update(int
   _last_n_new = 0;
 
-  _norm_factor = 1.0/(sqrt(double(_Freq)*double(_PSDLen))*double(n_subsets));
-
+  nf = 1.0/(double(_Freq)*double(nf/2.0));
   for ( i_samp = 0; i_samp<_PSDLen; i_samp++ ) {
-    psd[i_samp] *= _norm_factor;
+    psd[i_samp] = sqrt(psd[i_samp]*nf);
   }
 
@@ -373,5 +387,4 @@ void KstPSD::_adjustLengths() {
     delete[] _w;
     _w = new double[_ALen];
-    GenW();
 
     _last_f0 = 0;

--- kdeextragear-2/kst/kst/kstpsd.h  #1.4:1.5
@@ -71,5 +71,5 @@ class KstPSD : public KstDataObject {
 
   private:
-    void GenW();
+    void GenW(int len);
     void commonConstructor(const QString& in_tag, KstVectorPtr in_V,
         double freq, bool average, int len,
@@ -86,9 +86,9 @@ class KstPSD : public KstDataObject {
     double _Freq;
     int _Len;
-    double _norm_factor;
     double *_w;
     int _PSDLen;
     double *_a;
     int _ALen;
+    int _WLen;
 
     KstVectorMap::Iterator _sVector, _fVector;