Added power spectral density curves at the spectrum chart

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.PSD_VS_FREQUENCY:
+          fftData = GraphSpectrumCalc.dataLoadFrequencyPSD();
+          break;
+
         case SPECTRUM_TYPE.FREQUENCY:
         default:
           fftData = GraphSpectrumCalc.dataLoadFrequency();

src/graph_spectrum_calc.js

@@ -106,6 +106,29 @@ GraphSpectrumCalc.dataLoadFrequency = function() {
   return fftData;
 };
 
+GraphSpectrumCalc.dataLoadFrequencyPSD = function() {
+  const points_per_segment = 512,
+        overlap_count = 256;
+  const flightSamples = this._getFlightSamplesFreq(false);
+  const psd =  this._psd(flightSamples.samples, this._blackBoxRate, points_per_segment, overlap_count);
+  let min = 1e6,
+      max = -1e6;
+  for (const value of psd) {
+      min = Math.min(value, min);
+      max = Math.max(value, max);
+  }
+
+  const psdData = {
+    fieldIndex   : this._dataBuffer.fieldIndex,
+    fieldName  : this._dataBuffer.fieldName,
+    psdLength  : psd.length,
+    psdOutput  : psd,
+    blackBoxRate : this._blackBoxRate,
+    minimum: min,
+    maximum: max,
+  };
+  return psdData;
+};
 
 GraphSpectrumCalc._dataLoadFrequencyVsX = function(vsFieldNames, minValue = Infinity, maxValue = -Infinity) {
 
@@ -283,7 +306,7 @@ GraphSpectrumCalc._getFlightChunks = function() {
   return allChunks;
 };
 
-GraphSpectrumCalc._getFlightSamplesFreq = function() {
+GraphSpectrumCalc._getFlightSamplesFreq = function(scaled = true) {
 
   const allChunks = this._getFlightChunks();
 
@@ -293,7 +316,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 +512,83 @@ GraphSpectrumCalc._normalizeFft = function(fftOutput, fftLength) {
 
   return fftData;
 };
+
+/**
+ * Compute PSD for data samples by Welch method follow Python code
+ */
+GraphSpectrumCalc._psd  = function(samples, fs, n_per_seg, n_overlap, scaling = 'density') {
+// Compute FFT for samples segments
+  const fftOutput = this._fft_segmented(samples, n_per_seg, n_overlap);
+
+  const dataCount = fftOutput[0].length;
+  const segmentsCount = fftOutput.length;
+  const psdOutput = new Float64Array(dataCount);
+
+// Compute power scale coef
+  let scale = 1;
+  if (scaling == 'density') {
+    if (userSettings.analyserHanning) {
+      const window = Array(n_per_seg).fill(1);
+      this._hanningWindow(window, n_per_seg);
+      let skSum = 0;
+      for (const value of window) {
+        skSum += value ** 2;
+      }
+        scale = 1 / (fs * skSum);
+    } else {
+      scale = 1 / n_per_seg;
+    }
+  } else if (scaling == 'spectrum') {
+    if (userSettings.analyserHanning) {
+      const window = Array(n_per_seg).fill(1);
+      this._hanningWindow(window, n_per_seg);
+      let sum = 0;
+      for (const value of window) {
+        sum += value;
+      }
+      scale = 1 / sum ** 2;
+    } else {
+      scale = 1 / n_per_seg ** 2;
+    }
+  }
+
+// Compute average for scaled power
+  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-5; // limit min value for -50db
+    let avg = psdOutput[i] / segmentsCount;
+    avg = Math.max(avg, min_avg);
+    psdOutput[i] = 10 * Math.log10(avg);
+  }
+
+  return psdOutput;
+};
+
+
+/**
+ * Compute FFT for samples segments by lenghts as n_per_seg with n_overlap overlap points count
+ */
+GraphSpectrumCalc._fft_segmented  = function(samples, n_per_seg, n_overlap) {
+  const samplesCount = samples.length;
+  let output = [];
+  for (let i = 0; i < samplesCount - n_per_seg; i += n_per_seg - n_overlap) {
+    const fftInput = samples.slice(i, i + n_per_seg);
+
+    if (userSettings.analyserHanning) {
+      this._hanningWindow(fftInput, n_per_seg);
+    }
+
+    const fftOutput = this._fft(fftInput);
+    output.push(fftOutput.slice(0, n_per_seg));
+  }
+
+  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,
+  PSD_VS_FREQUENCY: 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.PSD_VS_FREQUENCY:
+      this._drawFrequencyPSDGraph(canvasCtx);
+      break;
   }
 };
 
@@ -294,7 +299,7 @@ GraphSpectrumPlot._drawFrequencyGraph = function (canvasCtx) {
     this._fftData.fieldName,
     WIDTH - 4,
     HEIGHT - 6,
-    "right"
+    "right",
   );
   this._drawHorizontalGridLines(
     canvasCtx,
@@ -304,7 +309,58 @@ GraphSpectrumPlot._drawFrequencyGraph = function (canvasCtx) {
     WIDTH,
     HEIGHT,
     MARGIN,
-    "Hz"
+    "Hz",
+  );
+};
+
+GraphSpectrumPlot._drawFrequencyPSDGraph = function (canvasCtx) {
+  const HEIGHT = canvasCtx.canvas.height - MARGIN;
+  const WIDTH = canvasCtx.canvas.width;
+  const LEFT = canvasCtx.canvas.left;
+  const TOP = canvasCtx.canvas.top;
+
+  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 = "orange";
+
+  canvasCtx.moveTo(0, 0);
+  const a1 = Math.abs(this._fftData.minimum),
+        a2 = Math.abs(this._fftData.maximum),
+        limit = Math.max(a1, a2);
+  const scaleY = HEIGHT / 2 / limit;
+  for (let pointNum = 0; pointNum < pointsCount; pointNum += 2) {
+    const freq = PLOTTED_BLACKBOX_RATE / 2 * pointNum / pointsCount;
+    const y = HEIGHT / 2 - this._fftData.psdOutput[pointNum] * scaleY;
+    canvasCtx.lineTo(freq * scaleX, y);
+  }
+  canvasCtx.stroke();
+
+  canvasCtx.restore();
+
+  this._drawAxisLabel(
+    canvasCtx,
+    this._fftData.fieldName,
+    WIDTH - 4,
+    HEIGHT - 6,
+    "right",
+  );
+  this._drawHorizontalGridLines(
+    canvasCtx,
+    PLOTTED_BLACKBOX_RATE / 2,
+    LEFT,
+    TOP,
+    WIDTH,
+    HEIGHT,
+    MARGIN,
+    "Hz",
   );
 };