diff --git a/index.html b/index.html index 5218a8d4..3f41083f 100644 --- a/index.html +++ b/index.html @@ -459,6 +459,7 @@

Workspace

+ diff --git a/src/graph_spectrum.js b/src/graph_spectrum.js index c89251a5..5d07c588 100644 --- a/src/graph_spectrum.js +++ b/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(); }) ) diff --git a/src/graph_spectrum_calc.js b/src/graph_spectrum_calc.js index b358dc22..f29989aa 100644 --- a/src/graph_spectrum_calc.js +++ b/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; +}; diff --git a/src/graph_spectrum_plot.js b/src/graph_spectrum_plot.js index b43d4aea..2ff35470 100644 --- a/src/graph_spectrum_plot.js +++ b/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,22 +1080,22 @@ 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; @@ -1016,7 +1103,7 @@ GraphSpectrumPlot._drawVerticalGridLines = function ( 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;