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

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;
+};

src/graph_spectrum_plot.js

@@ -17,6 +17,7 @@ export const SPECTRUM_TYPE = {
   FREQ_VS_THROTTLE: 1,
   PIDERROR_VS_SETPOINT: 2,
   FREQ_VS_RPM: 3,
+  POWER_SPECTRAL_DENSITY: 4,
 };
 
 export const SPECTRUM_OVERDRAW_TYPE = {
@@ -171,6 +172,10 @@ GraphSpectrumPlot._drawGraph = function (canvasCtx) {
     case SPECTRUM_TYPE.PIDERROR_VS_SETPOINT:
       this._drawPidErrorVsSetpointGraph(canvasCtx);
       break;
+
+    case SPECTRUM_TYPE.POWER_SPECTRAL_DENSITY:
+      this._drawPowerSpectralDensityGraph(canvasCtx);
+      break;
   }
 };
 
@@ -294,7 +299,7 @@ GraphSpectrumPlot._drawFrequencyGraph = function (canvasCtx) {
     this._fftData.fieldName,
     WIDTH - 4,
     HEIGHT - 6,
-    "right"
+    "right",
   );
   this._drawHorizontalGridLines(
     canvasCtx,
@@ -304,10 +309,92 @@ GraphSpectrumPlot._drawFrequencyGraph = function (canvasCtx) {
     WIDTH,
     HEIGHT,
     MARGIN,
-    "Hz"
+    "Hz",
   );
 };
 
+GraphSpectrumPlot._drawPowerSpectralDensityGraph = function (canvasCtx) {
+  const HEIGHT = canvasCtx.canvas.height - MARGIN;
+  const ACTUAL_MARGIN_LEFT = this._getActualMarginLeft();
+  const WIDTH = canvasCtx.canvas.width - ACTUAL_MARGIN_LEFT;
+  const LEFT = canvasCtx.canvas.offsetLeft + ACTUAL_MARGIN_LEFT;
+  const TOP = canvasCtx.canvas.offsetTop;
+
+  const PLOTTED_BLACKBOX_RATE = this._fftData.blackBoxRate / this._zoomX;
+
+  canvasCtx.save();
+  canvasCtx.translate(LEFT, TOP);
+  this._drawGradientBackground(canvasCtx, WIDTH, HEIGHT);
+
+  const pointsCount = this._fftData.psdLength;
+  const scaleX = 2 * WIDTH / PLOTTED_BLACKBOX_RATE * this._zoomX;
+  canvasCtx.beginPath();
+  canvasCtx.lineWidth = 1;
+  canvasCtx.strokeStyle = "white";
+
+  // Allign y axis range by 10db
+  const dbStep = 10;
+  const minY = Math.floor(this._fftData.minimum / dbStep) * dbStep;
+  const maxY = (Math.floor(this._fftData.maximum / dbStep) + 1) * dbStep;
+  const ticksCount = (maxY - minY) / dbStep;
+  const scaleY = HEIGHT / (maxY - minY);
+  //Store vsRange for _drawMousePosition
+  this._fftData.vsRange = {
+    min: minY,
+    max: maxY,
+  };
+  canvasCtx.moveTo(0, 0);
+  for (let pointNum = 0; pointNum < pointsCount; pointNum++) {
+    const freq = PLOTTED_BLACKBOX_RATE / 2 * pointNum / pointsCount;
+    const y = HEIGHT - (this._fftData.psdOutput[pointNum] - minY) * scaleY;
+    canvasCtx.lineTo(freq * scaleX, y);
+  }
+  canvasCtx.stroke();
+
+  this._drawAxisLabel(
+    canvasCtx,
+    this._fftData.fieldName,
+    WIDTH - 4,
+    HEIGHT - 6,
+    "right",
+  );
+  this._drawHorizontalGridLines(
+    canvasCtx,
+    PLOTTED_BLACKBOX_RATE / 2,
+    LEFT,
+    TOP,
+    WIDTH,
+    HEIGHT,
+    MARGIN,
+    "Hz",
+  );
+  this._drawVerticalGridLines(
+    canvasCtx,
+    LEFT,
+    TOP,
+    WIDTH,
+    HEIGHT,
+    minY,
+    maxY,
+    "dBm/Hz",
+    ticksCount,
+  );
+  const offset = 1;
+  this._drawInterestFrequency(
+    canvasCtx,
+    this._fftData.maxNoiseIdx,
+    PLOTTED_BLACKBOX_RATE,
+    "Max noise",
+    WIDTH,
+    HEIGHT,
+    15 * offset + MARGIN,
+    "rgba(255,0,0,0.50)",
+    3,
+  );
+
+  canvasCtx.restore();
+};
+
 GraphSpectrumPlot._drawFrequencyVsXGraph = function (canvasCtx) {
   const PLOTTED_BLACKBOX_RATE = this._fftData.blackBoxRate / this._zoomX;
 
@@ -862,7 +949,7 @@ GraphSpectrumPlot._drawFiltersAndMarkers = function (canvasCtx) {
       canvasCtx,
       this._fftData.maxNoiseIdx,
       PLOTTED_BLACKBOX_RATE,
-      "Max motor noise",
+      "Max noise",
       WIDTH,
       HEIGHT,
       15 * offset + MARGIN,
@@ -993,30 +1080,30 @@ GraphSpectrumPlot._drawVerticalGridLines = function (
   HEIGHT,
   minValue,
   maxValue,
-  label
+  label,
+  ticks = 5,
 ) {
-  const TICKS = 5;
 
-  for (let i = 0; i <= TICKS; i++) {
+  for (let i = 0; i <= ticks; i++) {
     canvasCtx.beginPath();
     canvasCtx.lineWidth = 1;
     canvasCtx.strokeStyle = "rgba(255,255,255,0.25)";
 
-    const verticalPosition = i * (HEIGHT / TICKS);
+    const verticalPosition = i * (HEIGHT / ticks);
     canvasCtx.moveTo(0, verticalPosition);
     canvasCtx.lineTo(WIDTH, verticalPosition);
 
     canvasCtx.stroke();
     const verticalAxisValue = (
-      (maxValue - minValue) * ((TICKS - i) / TICKS) +
+      (maxValue - minValue) * ((ticks - i) / ticks) +
       minValue
     ).toFixed(0);
     let textBaseline;
     switch (i) {
       case 0:
         textBaseline = "top";
         break;
-      case TICKS:
+      case ticks:
         textBaseline = "bottom";
         break;
       default:
@@ -1129,8 +1216,8 @@ GraphSpectrumPlot._drawGradientBackground = function (
   );
 
   if (this._isFullScreen) {
-    backgroundGradient.addColorStop(1, "rgba(0,0,0,0.9)");
-    backgroundGradient.addColorStop(0, "rgba(0,0,0,0.7)");
+    backgroundGradient.addColorStop(1, "rgba(0,0,0,1)");
+    backgroundGradient.addColorStop(0, "rgba(0,0,0,0.9)");
   } else {
     backgroundGradient.addColorStop(1, "rgba(255,255,255,0.25)");
     backgroundGradient.addColorStop(0, "rgba(255,255,255,0)");
@@ -1359,7 +1446,8 @@ GraphSpectrumPlot._drawMousePosition = function (
   if (
     this._spectrumType === SPECTRUM_TYPE.FREQUENCY ||
     this._spectrumType === SPECTRUM_TYPE.FREQ_VS_THROTTLE ||
-    this._spectrumType === SPECTRUM_TYPE.FREQ_VS_RPM
+    this._spectrumType === SPECTRUM_TYPE.FREQ_VS_RPM ||
+    this._spectrumType === SPECTRUM_TYPE.POWER_SPECTRAL_DENSITY
   ) {
     // Calculate frequency at mouse
     const sampleRate = this._fftData.blackBoxRate / this._zoomX;
@@ -1391,6 +1479,9 @@ GraphSpectrumPlot._drawMousePosition = function (
       case SPECTRUM_TYPE.FREQ_VS_RPM:
         unitLabel = "Hz";
         break;
+      case SPECTRUM_TYPE.POWER_SPECTRAL_DENSITY:
+        unitLabel = "dBm/Hz";
+        break;
       default:
         unitLabel = null;
         break;
@@ -1466,6 +1557,7 @@ GraphSpectrumPlot._getActualMarginLeft = function () {
   switch (this._spectrumType) {
     case SPECTRUM_TYPE.FREQ_VS_THROTTLE:
     case SPECTRUM_TYPE.FREQ_VS_RPM:
+    case SPECTRUM_TYPE.POWER_SPECTRAL_DENSITY:
       actualMarginLeft = this._isFullScreen
         ? MARGIN_LEFT_FULLSCREEN
         : MARGIN_LEFT;