Added Power Spectral Density curves at the spectrum chart (Completion of #820 PR) (#826)

* Added computing of PSD - Power Spectral Density for spectrum chart

* Added PSD curves at spectrum chart

* Removed redundant check TODO code

* Code issues are resolved

* PSD caption is edited in index.html

Co-authored-by: MikeNomatter <[email protected]>

* Removed non execute condition. The FFT output is allways even value

* Powers scale computing code refactoring

* Draw the all points of PSD data

* Resolved issue of magnitude computing

* Added vertical grid lines

* Added grid with db value captions

* The PSD Y axises range is alligned by value 10db/Hz

* Added setup points per segment by using vertical slider (Y-scale slider)

* maximal num per segment count value is limited by data samples count

* The names variable are changed to JS code style

* Added the maximal noise frequency label

* The 'Max motor noise' label is renamed to 'Max noise'

* Added mouse marker at the PSD chart

* Zeroes frequency magnitudes removed from curves min-max range

* Decreased spectrums charts plot transparent

* Code style improvement

* Max noise definition low frequency limit is changed from 100Hz to 50Hz for PSD chart

* Added horizont mouse position marker line at PSD chart

* The max noise label is shifted down at PSD chart

* The PSD chart updates immediately after Y scale sliders change

* PSD unit label is changed at dBm/Hz

* The minimum PSD value is set as -70dBm

* Code style improvement

Co-authored-by: Mark Haslinghuis <[email protected]>

* Code style improvement

Co-authored-by: Mark Haslinghuis <[email protected]>

* Code style improvement

* Added checking for missing samples segment data

Co-authored-by: Mark Haslinghuis <[email protected]>

* Code style improvement

Co-authored-by: Mark Haslinghuis <[email protected]>

* Resolved wrong variable name

* Resolved missing coma

* Resolved "for" loop condition issue, when num per segment value is equal samples count

* Code style improvement

---------

Co-authored-by: MikeNomatter <[email protected]>
Co-authored-by: Mark Haslinghuis <[email protected]>

Author: 3 people

index.html

@@ -459,6 +459,7 @@ <h4>Workspace</h4>
                                 <option value="1">Freq. vs Throttle</option>
                                 <option value="3">Freq. vs RPM</option>
                                 <option value="2">Error vs Setpoint</option>
+                                <option value="4">Power spectral density</option>
                             </select>
                         </div>
 

src/graph_spectrum.js

@@ -116,6 +116,10 @@ export function FlightLogAnalyser(flightLog, canvas, analyserCanvas) {
           fftData = GraphSpectrumCalc.dataLoadPidErrorVsSetpoint();
           break;
 
+        case SPECTRUM_TYPE.POWER_SPECTRAL_DENSITY:
+          fftData = GraphSpectrumCalc.dataLoadPSD(analyserZoomY);
+          break;
+
         case SPECTRUM_TYPE.FREQUENCY:
         default:
           fftData = GraphSpectrumCalc.dataLoadFrequency();
@@ -187,6 +191,11 @@ export function FlightLogAnalyser(flightLog, canvas, analyserCanvas) {
         debounce(100, function () {
           analyserZoomY = 1 / (analyserZoomYElem.val() / 100);
           GraphSpectrumPlot.setZoom(analyserZoomX, analyserZoomY);
+          // Recalculate PSD with updated samples per segment count
+          if (userSettings.spectrumType == SPECTRUM_TYPE.POWER_SPECTRAL_DENSITY) {
+            dataLoad();
+            GraphSpectrumPlot.setData(fftData, userSettings.spectrumType);
+          }
           that.refresh();
         })
       )

src/graph_spectrum_calc.js

@@ -106,6 +106,33 @@ GraphSpectrumCalc.dataLoadFrequency = function() {
   return fftData;
 };
 
+GraphSpectrumCalc.dataLoadPSD = function(analyserZoomY) {
+  const flightSamples = this._getFlightSamplesFreq(false);
+
+  let pointsPerSegment = 512;
+  const multipiler = Math.floor(1 / analyserZoomY); // 0. ... 10
+  if (multipiler == 0) {
+    pointsPerSegment = 256;
+  } else if (multipiler > 1) {
+    pointsPerSegment *= 2 ** Math.floor(multipiler / 2);
+  }
+  pointsPerSegment = Math.min(pointsPerSegment, flightSamples.samples.length);
+  const overlapCount = Math.floor(pointsPerSegment / 2);
+
+  const psd =  this._psd(flightSamples.samples, pointsPerSegment, overlapCount);
+
+  const psdData = {
+    fieldIndex   : this._dataBuffer.fieldIndex,
+    fieldName  : this._dataBuffer.fieldName,
+    psdLength  : psd.psdOutput.length,
+    psdOutput  : psd.psdOutput,
+    blackBoxRate : this._blackBoxRate,
+    minimum: psd.min,
+    maximum: psd.max,
+    maxNoiseIdx: psd.maxNoiseIdx,
+  };
+  return psdData;
+};
 
 GraphSpectrumCalc._dataLoadFrequencyVsX = function(vsFieldNames, minValue = Infinity, maxValue = -Infinity) {
 
@@ -283,7 +310,7 @@ GraphSpectrumCalc._getFlightChunks = function() {
   return allChunks;
 };
 
-GraphSpectrumCalc._getFlightSamplesFreq = function() {
+GraphSpectrumCalc._getFlightSamplesFreq = function(scaled = true) {
 
   const allChunks = this._getFlightChunks();
 
@@ -293,7 +320,11 @@ GraphSpectrumCalc._getFlightSamplesFreq = function() {
   let samplesCount = 0;
   for (const chunk of allChunks) {
     for (const frame of chunk.frames) {
-      samples[samplesCount] = (this._dataBuffer.curve.lookupRaw(frame[this._dataBuffer.fieldIndex]));
+      if (scaled) {
+        samples[samplesCount] = this._dataBuffer.curve.lookupRaw(frame[this._dataBuffer.fieldIndex]);
+      } else {
+        samples[samplesCount] = frame[this._dataBuffer.fieldIndex];
+      }
       samplesCount++;
     }
   }
@@ -485,3 +516,116 @@ GraphSpectrumCalc._normalizeFft = function(fftOutput, fftLength) {
 
   return fftData;
 };
+
+/**
+ * Compute PSD for data samples by Welch method follow Python code
+ */
+GraphSpectrumCalc._psd  = function(samples, pointsPerSegment, overlapCount, scaling = 'density') {
+// Compute FFT for samples segments
+  const fftOutput = this._fft_segmented(samples, pointsPerSegment, overlapCount);
+
+  const dataCount = fftOutput[0].length;
+  const segmentsCount = fftOutput.length;
+  const psdOutput = new Float64Array(dataCount);
+
+// Compute power scale coef
+  let scale = 1;
+  if (userSettings.analyserHanning) {
+    const window = Array(pointsPerSegment).fill(1);
+    this._hanningWindow(window, pointsPerSegment);
+    if (scaling == 'density') {
+      let skSum = 0;
+      for (const value of window) {
+        skSum += value ** 2;
+      }
+      scale = 1 / (this._blackBoxRate * skSum);
+    } else if (scaling == 'spectrum') {
+      let sum = 0;
+      for (const value of window) {
+        sum += value;
+      }
+      scale = 1 / sum ** 2;
+    }
+  } else if (scaling == 'density') {
+    scale = 1 / pointsPerSegment;
+  } else if (scaling == 'spectrum') {
+    scale = 1 / pointsPerSegment ** 2;
+  }
+
+// Compute average for scaled power
+  let min = 1e6,
+      max = -1e6;
+      // Early exit if no segments were processed
+  if (segmentsCount === 0) {
+    return {
+      psdOutput: new Float64Array(0),
+      min: 0,
+      max: 0,
+      maxNoiseIdx: 0,
+    };
+  }
+  const maxFrequency = (this._blackBoxRate / 2.0);
+  const noise50HzIdx = 50 / maxFrequency * dataCount;
+  const noise3HzIdx = 3 / maxFrequency * dataCount;
+  let maxNoiseIdx = 0;
+  let maxNoise = -100;
+  for (let i = 0; i < dataCount; i++) {
+    psdOutput[i] = 0.0;
+    for (let j = 0; j < segmentsCount; j++) {
+      let p = scale * fftOutput[j][i] ** 2;
+      if (i != dataCount - 1) {
+        p *= 2;
+      }
+      psdOutput[i] += p;
+    }
+
+    const min_avg = 1e-7; // limit min value for -70db
+    let avg = psdOutput[i] / segmentsCount;
+    avg = Math.max(avg, min_avg);
+    psdOutput[i] = 10 * Math.log10(avg);
+    if (i > noise3HzIdx) {    // Miss big zero freq magnitude
+      min = Math.min(psdOutput[i], min);
+      max = Math.max(psdOutput[i], max);
+    }
+    if (i > noise50HzIdx && psdOutput[i] > maxNoise) {
+      maxNoise = psdOutput[i];
+      maxNoiseIdx = i;
+    }
+  }
+
+  const maxNoiseFrequency = maxNoiseIdx / dataCount * maxFrequency;
+
+  return {
+      psdOutput: psdOutput,
+      min: min,
+      max: max,
+      maxNoiseIdx: maxNoiseFrequency,
+    };
+};
+
+
+/**
+ * Compute FFT for samples segments by lenghts as pointsPerSegment with overlapCount overlap points count
+ */
+GraphSpectrumCalc._fft_segmented  = function(samples, pointsPerSegment, overlapCount) {
+  const samplesCount = samples.length;
+  let output = [];
+  for (let i = 0; i <= samplesCount - pointsPerSegment; i += pointsPerSegment - overlapCount) {
+    const fftInput = samples.slice(i, i + pointsPerSegment);
+
+    if (userSettings.analyserHanning) {
+      this._hanningWindow(fftInput, pointsPerSegment);
+    }
+
+    const fftComplex = this._fft(fftInput);
+    const magnitudes = new Float64Array(pointsPerSegment / 2);
+    for (let i = 0; i < pointsPerSegment / 2; i++) {
+      const re = fftComplex[2 * i];
+      const im = fftComplex[2 * i + 1];
+      magnitudes[i] = Math.hypot(re, im);
+    }
+    output.push(magnitudes);
+  }
+
+  return output;
+};