Added power spectral density charts

index.html

@@ -458,8 +458,10 @@ <h4>Workspace</h4>
                                 <option value="0">Frequency</option>
                                 <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>
+                                <option value="5">PSD vs Throttle</option>
+                                <option value="6">PSD vs RPM</option>
+                                <option value="2">Error vs Setpoint</option>
                             </select>
                         </div>
 

src/graph_spectrum.js

@@ -112,6 +112,14 @@ export function FlightLogAnalyser(flightLog, canvas, analyserCanvas) {
           fftData = GraphSpectrumCalc.dataLoadFrequencyVsRpm();
           break;
 
+        case SPECTRUM_TYPE.PSD_VS_THROTTLE:
+          fftData = GraphSpectrumCalc.dataLoadFrequencyVsThrottle(true);
+          break;
+
+        case SPECTRUM_TYPE.PSD_VS_RPM:
+          fftData = GraphSpectrumCalc.dataLoadFrequencyVsRpm(true);
+          break;
+
         case SPECTRUM_TYPE.PIDERROR_VS_SETPOINT:
           fftData = GraphSpectrumCalc.dataLoadPidErrorVsSetpoint();
           break;

src/graph_spectrum_calc.js

@@ -79,7 +79,7 @@ GraphSpectrumCalc.setOutTime = function(time) {
   if ((time - this._analyserTimeRange.in) <= MAX_ANALYSER_LENGTH) {
     this._analyserTimeRange.out = time;
   } else {
-    this._analyserTimeRange.out = analyserTimeRange.in + MAX_ANALYSER_LENGTH;
+    this._analyserTimeRange.out = this._analyserTimeRange.in + MAX_ANALYSER_LENGTH;
   }
   return this._analyserTimeRange.out;
 };
@@ -110,11 +110,11 @@ GraphSpectrumCalc.dataLoadPSD = function(analyserZoomY) {
   const flightSamples = this._getFlightSamplesFreq(false);
 
   let pointsPerSegment = 512;
-  const multipiler = Math.floor(1 / analyserZoomY); // 0. ... 10
-  if (multipiler == 0) {
+  const multiplier = Math.floor(1 / analyserZoomY); // 0. ... 10
+  if (multiplier == 0) {
     pointsPerSegment = 256;
-  } else if (multipiler > 1) {
-    pointsPerSegment *= 2 ** Math.floor(multipiler / 2);
+  } else if (multiplier > 1) {
+    pointsPerSegment *= 2 ** Math.floor(multiplier / 2);
   }
   pointsPerSegment = Math.min(pointsPerSegment, flightSamples.samples.length);
   const overlapCount = Math.floor(pointsPerSegment / 2);
@@ -129,7 +129,7 @@ GraphSpectrumCalc.dataLoadPSD = function(analyserZoomY) {
     blackBoxRate : this._blackBoxRate,
     minimum: psd.min,
     maximum: psd.max,
-    maxNoiseIdx: psd.maxNoiseIdx,
+    maxNoiseFrequency: psd.maxNoiseFrequency,
   };
   return psdData;
 };
@@ -218,12 +218,85 @@ GraphSpectrumCalc._dataLoadFrequencyVsX = function(vsFieldNames, minValue = Infi
 
 };
 
-GraphSpectrumCalc.dataLoadFrequencyVsThrottle = function() {
-  return this._dataLoadFrequencyVsX(FIELD_THROTTLE_NAME, 0, 100);
+GraphSpectrumCalc._dataLoadPowerSpectralDensityVsX = function(vsFieldNames, minValue = Infinity, maxValue = -Infinity) {
+
+  const flightSamples = this._getFlightSamplesFreqVsX(vsFieldNames, minValue, maxValue, false);
+
+  // We divide it into FREQ_VS_THR_CHUNK_TIME_MS FFT chunks, we calculate the average throttle
+  // for each chunk. We use a moving window to get more chunks available.
+  const fftChunkLength = Math.round(this._blackBoxRate * FREQ_VS_THR_CHUNK_TIME_MS / 1000);
+  const fftChunkWindow = Math.round(fftChunkLength / FREQ_VS_THR_WINDOW_DIVISOR);
+
+  let maxNoise = 0; // Stores the maximum amplitude of the fft over all chunks
+  const psdLength = Math.floor(fftChunkLength / 2);
+   // Matrix where each row represents a bin of vs values, and the columns are amplitudes at frequencies
+  const matrixPsdOutput = new Array(NUM_VS_BINS)
+        .fill(null)
+        .map(() => (new Float64Array(psdLength)).fill(-100));
+
+  const numberSamples = new Uint32Array(NUM_VS_BINS); // Number of samples in each vs value, used to average them later.
+
+  for (let fftChunkIndex = 0; fftChunkIndex + fftChunkLength < flightSamples.samples.length; fftChunkIndex += fftChunkWindow) {
+
+    const fftInput = flightSamples.samples.slice(fftChunkIndex, fftChunkIndex + fftChunkLength);
+    const psd = this._psd(fftInput, fftChunkLength, 0, 'density');
+    maxNoise = Math.max(psd.max, maxNoise);
+    // calculate a bin index and put the fft value in that bin for each field (e.g. eRPM[0], eRPM[1]..) sepparately
+    for (const vsValueArray of flightSamples.vsValues) {
+      // Calculate average of the VS values in the chunk
+      let sumVsValues = 0;
+      for (let indexVs = fftChunkIndex; indexVs < fftChunkIndex + fftChunkLength; indexVs++) {
+        sumVsValues += vsValueArray[indexVs];
+      }
+      // Translate the average vs value to a bin index
+      const avgVsValue = sumVsValues / fftChunkLength;
+      let vsBinIndex = Math.floor(NUM_VS_BINS * (avgVsValue - flightSamples.minValue) / (flightSamples.maxValue - flightSamples.minValue));
+      // ensure that avgVsValue == flightSamples.maxValue does not result in an out of bounds access
+      if (vsBinIndex === NUM_VS_BINS) { vsBinIndex = NUM_VS_BINS - 1; }
+      numberSamples[vsBinIndex]++;
+
+      // add the output from the fft to the row given by the vs value bin index
+      for (let i = 0; i < psdLength; i++) {
+        matrixPsdOutput[vsBinIndex][i] += psd.psdOutput[i];
+      }
+    }
+  }
+
+  // Divide the values from the fft in each row (vs value bin) by the number of samples in the bin
+  for (let i = 0; i < NUM_VS_BINS; i++) {
+    if (numberSamples[i] > 1) {
+      for (let j = 0; j < psdLength; j++) {
+        matrixPsdOutput[i][j] /= numberSamples[i];
+      }
+    }
+  }
+
+  // The output data needs to be smoothed, the sampling is not perfect
+  // but after some tests we let the data as is, an we prefer to apply a
+  // blur algorithm to the heat map image
+
+  const psdData = {
+    fieldIndex   : this._dataBuffer.fieldIndex,
+    fieldName  : this._dataBuffer.fieldName,
+    fftLength  : psdLength,
+    fftOutput  : matrixPsdOutput,
+    maxNoise   : maxNoise,
+    blackBoxRate : this._blackBoxRate,
+    vsRange    : { min: flightSamples.minValue, max: flightSamples.maxValue},
+  };
+
+  return psdData;
+
+};
+
+GraphSpectrumCalc.dataLoadFrequencyVsThrottle = function(drawPSD = false) {
+  return drawPSD ? this._dataLoadPowerSpectralDensityVsX(FIELD_THROTTLE_NAME, 0, 100) :
+                   this._dataLoadFrequencyVsX(FIELD_THROTTLE_NAME, 0, 100);
 };
 
-GraphSpectrumCalc.dataLoadFrequencyVsRpm = function() {
-  const fftData = this._dataLoadFrequencyVsX(FIELD_RPM_NAMES, 0);
+GraphSpectrumCalc.dataLoadFrequencyVsRpm = function(drawPSD = false) {
+  const fftData = drawPSD ? this._dataLoadPowerSpectralDensityVsX(FIELD_RPM_NAMES, 0) :
+                            this._dataLoadFrequencyVsX(FIELD_RPM_NAMES, 0);
   fftData.vsRange.max *= 3.333 / this._motorPoles;
   fftData.vsRange.min *= 3.333 / this._motorPoles;
   return fftData;
@@ -345,7 +418,7 @@ GraphSpectrumCalc._getVsIndexes = function(vsFieldNames) {
   return fieldIndexes;
 };
 
-GraphSpectrumCalc._getFlightSamplesFreqVsX = function(vsFieldNames, minValue = Infinity, maxValue = -Infinity) {
+GraphSpectrumCalc._getFlightSamplesFreqVsX = function(vsFieldNames, minValue = Infinity, maxValue = -Infinity, scaled = true) {
 
   const allChunks = this._getFlightChunks();
   const vsIndexes = this._getVsIndexes(vsFieldNames);
@@ -356,7 +429,11 @@ GraphSpectrumCalc._getFlightSamplesFreqVsX = function(vsFieldNames, minValue = I
   let samplesCount = 0;
   for (const chunk of allChunks) {
     for (let frameIndex = 0; frameIndex < chunk.frames.length; frameIndex++) {
-      samples[samplesCount] = (this._dataBuffer.curve.lookupRaw(chunk.frames[frameIndex][this._dataBuffer.fieldIndex]));
+      if (scaled) {
+        samples[samplesCount] = (this._dataBuffer.curve.lookupRaw(chunk.frames[frameIndex][this._dataBuffer.fieldIndex]));
+      } else {
+        samples[samplesCount] = chunk.frames[frameIndex][this._dataBuffer.fieldIndex];
+      }
       for (let i = 0; i < vsIndexes.length; i++) {
         let vsFieldIx = vsIndexes[i];
         let value = chunk.frames[frameIndex][vsFieldIx];
@@ -503,14 +580,14 @@ GraphSpectrumCalc._normalizeFft = function(fftOutput, fftLength) {
     }
   }
 
-  maxNoiseIdx = maxNoiseIdx / fftLength * maxFrequency;
+  const maxNoiseFrequency = maxNoiseIdx / fftLength * maxFrequency;
 
   const fftData = {
     fieldIndex   : this._dataBuffer.fieldIndex,
     fieldName  : this._dataBuffer.fieldName,
     fftLength  : fftLength,
     fftOutput  : fftOutput,
-    maxNoiseIdx  : maxNoiseIdx,
+    maxNoiseFrequency  : maxNoiseFrequency,
     blackBoxRate : this._blackBoxRate,
   };
 
@@ -561,7 +638,7 @@ GraphSpectrumCalc._psd  = function(samples, pointsPerSegment, overlapCount, scal
       psdOutput: new Float64Array(0),
       min: 0,
       max: 0,
-      maxNoiseIdx: 0,
+      maxNoiseFrequency: 0,
     };
   }
   const maxFrequency = (this._blackBoxRate / 2.0);
@@ -599,16 +676,16 @@ GraphSpectrumCalc._psd  = function(samples, pointsPerSegment, overlapCount, scal
       psdOutput: psdOutput,
       min: min,
       max: max,
-      maxNoiseIdx: maxNoiseFrequency,
+      maxNoiseFrequency: 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;
+  const samplesCount = samples.length,
+        fftLength = Math.floor(pointsPerSegment / 2);
   let output = [];
   for (let i = 0; i <= samplesCount - pointsPerSegment; i += pointsPerSegment - overlapCount) {
     const fftInput = samples.slice(i, i + pointsPerSegment);
@@ -618,8 +695,8 @@ GraphSpectrumCalc._fft_segmented  = function(samples, pointsPerSegment, overlapC
     }
 
     const fftComplex = this._fft(fftInput);
-    const magnitudes = new Float64Array(pointsPerSegment / 2);
-    for (let i = 0; i < pointsPerSegment / 2; i++) {
+    const magnitudes = new Float64Array(fftLength);
+    for (let i = 0; i < fftLength; i++) {
       const re = fftComplex[2 * i];
       const im = fftComplex[2 * i + 1];
       magnitudes[i] = Math.hypot(re, im);