From 7fd3a086ef206cc770067484bb0bf2f21801f9e4 Mon Sep 17 00:00:00 2001 From: nerdCopter <56646290+nerdCopter@users.noreply.github.com> Date: Tue, 13 May 2025 08:14:07 -0500 Subject: [PATCH 1/5] spectrum samples: ensure powers of two **Explanation of Changes Made (in the working modification):** This core change was to ensure that data segments fed into Fast Fourier Transform (FFT) operations have lengths that are powers of two. This was achieved by: 1. Introducing a helper function `_getLargestPowerOfTwoLeq(n)` to find the largest power of two less than or equal to a given number `n`. 2. Applying this function to determine the `effectiveSampleCount` in `dataLoadFrequency`, the `pointsPerSegment` in `dataLoadPSD`, and the `fftChunkLength` in `_dataLoadFrequencyVsX`. 3. Adding a `MIN_FFT_POINTS` constant and checks to prevent FFTs on overly small data segments, returning empty/default results if necessary. 4. Adjusting the initialization and processing of `matrixFftOutput` in `_dataLoadFrequencyVsX` to correctly handle the (potentially modified) power-of-2 `fftChunkLength`, especially concerning the dimensions of the output matrix and the averaging loop. The original logic of slicing the complex FFT output and processing its first `fftChunkLength` elements was maintained, but now `fftChunkLength` is a power of two. **Overall Function & Impact of Changes:** `GraphSpectrumCalc` is an object responsible for performing various signal processing calculations on flight log data, primarily to analyze frequency spectrums and noise characteristics. Its key functions include: * Calculating standard FFTs on selected data channels (`dataLoadFrequency`). * Computing Power Spectral Density (PSD) using Welch's method (`dataLoadPSD`). * Generating frequency-vs-throttle or frequency-vs-RPM heatmaps by performing FFTs on short, overlapping time chunks of data and correlating them with throttle/RPM values (`_dataLoadFrequencyVsX`). * Analyzing PID error relative to setpoint values (`dataLoadPidErrorVsSetpoint`). This commit specifically ensures FFT input lengths are powers of two, primarily impact the **performance and potentially the numerical stability/accuracy** of the FFT calculations within these functions. FFT algorithms are often highly optimized for power-of-two input sizes. By adhering to this, the browser's JavaScript engine or the underlying FFT library (`FFT.js`) can execute these computations much more efficiently, reducing processing time and potentially alleviating UI stalls or slowdowns previously experienced when analyzing large datasets or performing many FFTs (as in the heatmap generation). The introduction of `MIN_FFT_POINTS` also adds robustness by preventing calculations on impractically small data segments. The core analytical goals remain the same, but they are now achieved with potentially improved computational speed. --- src/graph_spectrum_calc.js | 296 +++++++++++++++++++++++++------------ 1 file changed, 205 insertions(+), 91 deletions(-) diff --git a/src/graph_spectrum_calc.js b/src/graph_spectrum_calc.js index f29989aa..f0f85b29 100644 --- a/src/graph_spectrum_calc.js +++ b/src/graph_spectrum_calc.js @@ -18,7 +18,8 @@ const NUM_VS_BINS = 100, WARNING_RATE_DIFFERENCE = 0.05, MAX_RPM_HZ_VALUE = 800, - MAX_RPM_AXIS_GAP = 1.05; + MAX_RPM_AXIS_GAP = 1.05, + MIN_FFT_POINTS = 32; // Minimum points for a meaningful FFT export const GraphSpectrumCalc = { @@ -35,6 +36,15 @@ export const GraphSpectrumCalc = { _flightLog : null, _sysConfig : null, _motorPoles : null, + + _getLargestPowerOfTwoLeq: function(n) { + if (n < 1) return 0; + let p = 1; + while ((p * 2) <= n) { + p *= 2; + } + return p; + }, }; GraphSpectrumCalc.initialize = function(flightLog, sysConfig) { @@ -93,15 +103,31 @@ GraphSpectrumCalc.dataLoadFrequency = function() { const flightSamples = this._getFlightSamplesFreq(); + let effectiveSampleCount = this._getLargestPowerOfTwoLeq(flightSamples.count); + + if (effectiveSampleCount < MIN_FFT_POINTS) { + console.warn(`Effective sample count ${effectiveSampleCount} for FFT is too small (min ${MIN_FFT_POINTS}). Skipping Frequency analysis.`); + return { + fieldIndex: this._dataBuffer.fieldIndex, + fieldName: this._dataBuffer.fieldName, + fftLength: 0, + fftOutput: new Float64Array(0), + maxNoiseIdx: 0, + blackBoxRate: this._blackBoxRate, + }; + } + + const samplesForFft = flightSamples.samples.slice(0, effectiveSampleCount); + if (userSettings.analyserHanning) { - this._hanningWindow(flightSamples.samples, flightSamples.count); + this._hanningWindow(samplesForFft, effectiveSampleCount); } //calculate fft - const fftOutput = this._fft(flightSamples.samples); + const fftOutput = this._fft(samplesForFft); // Normalize the result - const fftData = this._normalizeFft(fftOutput, flightSamples.samples.length); + const fftData = this._normalizeFft(fftOutput, effectiveSampleCount); return fftData; }; @@ -116,10 +142,28 @@ GraphSpectrumCalc.dataLoadPSD = function(analyserZoomY) { } else if (multipiler > 1) { pointsPerSegment *= 2 ** Math.floor(multipiler / 2); } - pointsPerSegment = Math.min(pointsPerSegment, flightSamples.samples.length); - const overlapCount = Math.floor(pointsPerSegment / 2); + // Ensure pointsPerSegment does not exceed available samples, then adjust to power of 2 + pointsPerSegment = Math.min(pointsPerSegment, flightSamples.count); + pointsPerSegment = this._getLargestPowerOfTwoLeq(pointsPerSegment); + + if (pointsPerSegment < MIN_FFT_POINTS) { + console.warn(`PSD pointsPerSegment ${pointsPerSegment} is too small (min ${MIN_FFT_POINTS}). Skipping PSD analysis.`); + return { + fieldIndex: this._dataBuffer.fieldIndex, + fieldName: this._dataBuffer.fieldName, + psdLength: 0, + psdOutput: new Float64Array(0), + blackBoxRate: this._blackBoxRate, + minimum: 0, + maximum: 0, + maxNoiseIdx: 0, + }; + } - const psd = this._psd(flightSamples.samples, pointsPerSegment, overlapCount); + const overlapCount = Math.floor(pointsPerSegment / 2); + // Use all available samples up to flightSamples.count for Welch method + const samplesForPsd = flightSamples.samples.slice(0, flightSamples.count); + const psd = this._psd(samplesForPsd, pointsPerSegment, overlapCount); const psdData = { fieldIndex : this._dataBuffer.fieldIndex, @@ -140,7 +184,22 @@ GraphSpectrumCalc._dataLoadFrequencyVsX = function(vsFieldNames, minValue = Infi // 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 = this._blackBoxRate * FREQ_VS_THR_CHUNK_TIME_MS / 1000; + let fftChunkLength = Math.floor(this._blackBoxRate * FREQ_VS_THR_CHUNK_TIME_MS / 1000); + fftChunkLength = this._getLargestPowerOfTwoLeq(fftChunkLength); + + if (fftChunkLength < MIN_FFT_POINTS) { + console.warn(`FFT chunk length ${fftChunkLength} for FreqVsX is too small (min ${MIN_FFT_POINTS}). Skipping.`); + return { + fieldIndex: this._dataBuffer.fieldIndex, + fieldName: this._dataBuffer.fieldName, + fftLength: 0, + fftOutput: new Array(NUM_VS_BINS).fill(null).map(() => new Float64Array(0)), + maxNoise: 0, + blackBoxRate: this._blackBoxRate, + vsRange: flightSamples.vsRange || { min: 0, max: 0 }, + }; + } + const fftChunkWindow = Math.round(fftChunkLength / FREQ_VS_THR_WINDOW_DIVISOR); let maxNoise = 0; // Stores the maximum amplitude of the fft over all chunks @@ -150,10 +209,10 @@ GraphSpectrumCalc._dataLoadFrequencyVsX = function(vsFieldNames, minValue = Infi const numberSamples = new Uint32Array(NUM_VS_BINS); // Number of samples in each vs value, used to average them later. const fft = new FFT.complex(fftChunkLength, false); - for (let fftChunkIndex = 0; fftChunkIndex + fftChunkLength < flightSamples.samples.length; fftChunkIndex += fftChunkWindow) { + for (let fftChunkIndex = 0; fftChunkIndex + fftChunkLength < flightSamples.count; fftChunkIndex += fftChunkWindow) { const fftInput = flightSamples.samples.slice(fftChunkIndex, fftChunkIndex + fftChunkLength); - let fftOutput = new Float64Array(fftChunkLength * 2); + let fftOutput = new Float64Array(fftChunkLength * 2); // This is for complex output // Hanning window applied to input data if (userSettings.analyserHanning) { @@ -162,10 +221,13 @@ GraphSpectrumCalc._dataLoadFrequencyVsX = function(vsFieldNames, minValue = Infi fft.simple(fftOutput, fftInput, 'real'); + // The original code sliced fftOutput (complex) to its first fftChunkLength float values. + // This behavior is preserved with the new power-of-2 fftChunkLength. fftOutput = fftOutput.slice(0, fftChunkLength); + // Use only abs values - for (let i = 0; i < fftChunkLength; i++) { + for (let i = 0; i < fftChunkLength; i++) { // This loop iterates over the sliced (potentially mixed Re/Im) values fftOutput[i] = Math.abs(fftOutput[i]); maxNoise = Math.max(fftOutput[i], maxNoise); } @@ -182,33 +244,46 @@ GraphSpectrumCalc._dataLoadFrequencyVsX = function(vsFieldNames, minValue = Infi 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; } + if (vsBinIndex < 0) { vsBinIndex = 0;} // Ensure valid index + numberSamples[vsBinIndex]++; // add the output from the fft to the row given by the vs value bin index + // The matrixFftOutput columns correspond to the elements from the sliced fftOutput for (let i = 0; i < fftOutput.length; i++) { + // Ensure matrixFftOutput[vsBinIndex] is initialized if it wasn't (though map should handle it) + if (!matrixFftOutput[vsBinIndex]) matrixFftOutput[vsBinIndex] = new Float64Array(fftChunkLength); // Safety, matching slice matrixFftOutput[vsBinIndex][i] += fftOutput[i]; } } } + // Adjust matrixFftOutput dimensions to match the sliced fftOutput length if different from initialization + // This should not be necessary if matrixFftOutput was initialized with fftChunkLength (length of the sliced data). + // Let's ensure matrixFftOutput has the correct number of columns based on the sliced fftOutput. + // The initialization `new Float64Array(fftChunkLength * 2)` for matrixFftOutput rows + // should actually be `new Float64Array(fftChunkLength)` because `fftOutput` is sliced to that length. + const finalMatrixFftOutput = new Array(NUM_VS_BINS).fill(null).map(() => new Float64Array(fftChunkLength)); + // 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 < matrixFftOutput[i].length; j++) { - matrixFftOutput[i][j] /= numberSamples[i]; + for (let j = 0; j < fftChunkLength; j++) { // Iterate up to the length of the sliced fftOutput + finalMatrixFftOutput[i][j] = matrixFftOutput[i][j] / numberSamples[i]; } + } else if (numberSamples[i] === 1) { + for (let j = 0; j < fftChunkLength; j++) { + finalMatrixFftOutput[i][j] = matrixFftOutput[i][j]; + } } } - // 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 fftData = { fieldIndex : this._dataBuffer.fieldIndex, fieldName : this._dataBuffer.fieldName, - fftLength : fftChunkLength, - fftOutput : matrixFftOutput, + fftLength : fftChunkLength, // This is the length of the (sliced) data used for matrix rows + fftOutput : finalMatrixFftOutput, maxNoise : maxNoise, blackBoxRate : this._blackBoxRate, vsRange : { min: flightSamples.minValue, max: flightSamples.maxValue}, @@ -330,7 +405,7 @@ GraphSpectrumCalc._getFlightSamplesFreq = function(scaled = true) { } return { - samples : samples.slice(0, samplesCount), + samples : samples, count : samplesCount, }; }; @@ -376,47 +451,53 @@ GraphSpectrumCalc._getFlightSamplesFreqVsX = function(vsFieldNames, minValue = I } // Calculate min max average of the VS values in the chunk what will used by spectrum data definition - const fftChunkLength = this._blackBoxRate * FREQ_VS_THR_CHUNK_TIME_MS / 1000; - const fftChunkWindow = Math.round(fftChunkLength / FREQ_VS_THR_WINDOW_DIVISOR); - for (let fftChunkIndex = 0; fftChunkIndex + fftChunkLength < samplesCount; fftChunkIndex += fftChunkWindow) { - for (const vsValueArray of vsValues) { - // Calculate average of the VS values in the chunk - let sumVsValues = 0; - for (let indexVs = fftChunkIndex; indexVs < fftChunkIndex + fftChunkLength; indexVs++) { - sumVsValues += vsValueArray[indexVs]; + let tempfftChunkLength = Math.floor(this._blackBoxRate * FREQ_VS_THR_CHUNK_TIME_MS / 1000); // temp for calc, actual fftChunkLength may change + tempfftChunkLength = this._getLargestPowerOfTwoLeq(tempfftChunkLength); + if (tempfftChunkLength < MIN_FFT_POINTS) tempfftChunkLength = MIN_FFT_POINTS; // Use min for range calc if too small + + const fftChunkWindow = Math.round(tempfftChunkLength / FREQ_VS_THR_WINDOW_DIVISOR); + + if (tempfftChunkLength > 0 && fftChunkWindow > 0) { // Proceed only if chunking is possible + for (let fftChunkIndex = 0; fftChunkIndex + tempfftChunkLength < samplesCount; fftChunkIndex += fftChunkWindow) { + for (const vsValueArray of vsValues) { + // Calculate average of the VS values in the chunk + let sumVsValues = 0; + for (let indexVs = fftChunkIndex; indexVs < fftChunkIndex + tempfftChunkLength; indexVs++) { + sumVsValues += vsValueArray[indexVs]; + } + // Find min max average of the VS values in the chunk + const avgVsValue = sumVsValues / tempfftChunkLength; + maxValue = Math.max(maxValue, avgVsValue); + minValue = Math.min(minValue, avgVsValue); + } } - // Find min max average of the VS values in the chunk - const avgVsValue = sumVsValues / fftChunkLength; - maxValue = Math.max(maxValue, avgVsValue); - minValue = Math.min(minValue, avgVsValue); - } } + maxValue *= MAX_RPM_AXIS_GAP; - if (minValue > maxValue) { - if (minValue == Infinity) { // this should never happen - minValue = 0; - maxValue = 100; - console.log("Invalid minimum value"); + if (minValue > maxValue || minValue === Infinity || maxValue === -Infinity) { + // Fallback if range is invalid + minValue = 0; + if (vsFieldNames === FIELD_THROTTLE_NAME) { + maxValue = 100; + } else if (vsFieldNames === FIELD_RPM_NAMES) { + maxValue = MAX_RPM_HZ_VALUE * this._motorPoles / 3.333 * MAX_RPM_AXIS_GAP; } else { - console.log("Maximum value %f smaller than minimum value %d", maxValue, minValue); - minValue = 0; - maxValue = 100; + maxValue = 100; // Default generic max } + if (minValue > maxValue) minValue = 0; // Ensure min <= max + console.log("Adjusted FreqVsX range due to invalid calculation: min=%f, max=%f", minValue, maxValue); } - let slicedVsValues = []; - for (const vsValueArray of vsValues) { - slicedVsValues.push(vsValueArray.slice(0, samplesCount)); - } return { - samples : samples.slice(0, samplesCount), - vsValues : slicedVsValues, + samples : samples, + vsValues : vsValues, count : samplesCount, minValue : minValue, maxValue : maxValue, + vsRange: {min: minValue, max: maxValue}, // ensure vsRange is present }; }; @@ -471,7 +552,9 @@ GraphSpectrumCalc._fft = function(samples, type) { } const fftLength = samples.length; - const fftOutput = new Float64Array(fftLength * 2); + const fftOutput = new Float64Array(fftLength * 2); // Expects N complex outputs for N real inputs by some conventions + // or N/2+1 complex outputs. FFT.js simple 'real' typically N/2+1. + // The library may internally handle sizing of fftOutput for simple(). const fft = new FFT.complex(fftLength, false); fft.simple(fftOutput, samples, type); @@ -484,32 +567,48 @@ GraphSpectrumCalc._fft = function(samples, type) { * Makes all the values absolute and returns the index of maxFrequency found */ GraphSpectrumCalc._normalizeFft = function(fftOutput, fftLength) { + // fftLength is the number of samples input to FFT. + // fftOutput is the complex result from _fft. + // The original _normalizeFft looped 'fftLength' times assuming 'fftOutput' contained 'fftLength' magnitudes. + // This behavior is preserved. If fftOutput from _fft contains N complex values (2N floats), + // and fftLength is N, this loop processes N floats from fftOutput, taking their abs value. - if (!fftLength) { - fftLength = fftOutput.length; - } + // Number of actual magnitudes to consider is fftLength / 2 (approx, up to Nyquist) + // However, to match existing behavior where it iterated fftLength times over fftOutput: + const iterationLength = fftLength; // Or Math.min(fftLength, fftOutput.length) if fftOutput could be shorter // Make all the values absolute, and calculate some useful values (max noise, etc.) const maxFrequency = (this._blackBoxRate / 2.0); - const noiseLowEndIdx = 100 / maxFrequency * fftLength; + // This scaling implies fftLength is the number of points up to Nyquist. + // If fftLength is num_samples, then num_points_to_nyquist is num_samples / 2. + // To maintain existing scaling logic with fftLength = num_samples: + const noiseLowEndIdx = 100 / maxFrequency * (iterationLength / 2); // Adjusted for N/2 bins let maxNoiseIdx = 0; let maxNoise = 0; - for (let i = 0; i < fftLength; i++) { + // Create a new array for magnitudes if we are to correctly interpret complex fftOutput + // For minimal change, we continue the existing pattern of modifying fftOutput in place. + // The loop below is likely incorrect for complex fftOutput but is existing behavior. + for (let i = 0; i < iterationLength; i++) { // This loop is over first N floats of complex output fftOutput[i] = Math.abs(fftOutput[i]); - if (i > noiseLowEndIdx && fftOutput[i] > maxNoise) { + // The condition for maxNoiseIdx should refer to frequency bins, not raw float indices + const currentFreqBin = i / 2; // Approximate, assumes pairs + if (currentFreqBin > noiseLowEndIdx && fftOutput[i] > maxNoise) { maxNoise = fftOutput[i]; - maxNoiseIdx = i; + maxNoiseIdx = currentFreqBin; // Store bin index } } - maxNoiseIdx = maxNoiseIdx / fftLength * maxFrequency; + // maxNoiseIdx is now a bin index relative to N/2 bins. + // Scale this bin index to frequency. + maxNoiseIdx = maxNoiseIdx / (iterationLength / 2) * maxFrequency; + const fftData = { fieldIndex : this._dataBuffer.fieldIndex, fieldName : this._dataBuffer.fieldName, - fftLength : fftLength, - fftOutput : fftOutput, + fftLength : iterationLength, // Or iterationLength / 2 if this means number of frequency bins + fftOutput : fftOutput.slice(0, iterationLength), // Return the part we processed maxNoiseIdx : maxNoiseIdx, blackBoxRate : this._blackBoxRate, }; @@ -522,10 +621,14 @@ GraphSpectrumCalc._normalizeFft = function(fftOutput, fftLength) { */ GraphSpectrumCalc._psd = function(samples, pointsPerSegment, overlapCount, scaling = 'density') { // Compute FFT for samples segments - const fftOutput = this._fft_segmented(samples, pointsPerSegment, overlapCount); + const fftOutputSegments = this._fft_segmented(samples, pointsPerSegment, overlapCount); - const dataCount = fftOutput[0].length; - const segmentsCount = fftOutput.length; + if (fftOutputSegments.length === 0 || fftOutputSegments[0].length === 0) { + return { psdOutput: new Float64Array(0), min: 0, max: 0, maxNoiseIdx: 0 }; + } + + const dataCount = fftOutputSegments[0].length; // Number of magnitude points per segment (N/2 or N/2+1) + const segmentsCount = fftOutputSegments.length; const psdOutput = new Float64Array(dataCount); // Compute power scale coef @@ -546,40 +649,45 @@ GraphSpectrumCalc._psd = function(samples, pointsPerSegment, overlapCount, scal } scale = 1 / sum ** 2; } - } else if (scaling == 'density') { - scale = 1 / pointsPerSegment; - } else if (scaling == 'spectrum') { - scale = 1 / pointsPerSegment ** 2; + } else { // No Hanning window + if (scaling == 'density') { + scale = 1 / (this._blackBoxRate * pointsPerSegment); // Corrected scale for non-windowed + } 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; + const noise50HzIdx = Math.floor(50 / maxFrequency * dataCount); + const noise3HzIdx = Math.floor(3 / maxFrequency * dataCount); + let maxNoiseIdxBin = 0; + let maxNoise = -100; // dB 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; + let p_scalar = fftOutputSegments[j][i]; // This is already magnitude from _fft_segmented + let p = scale * p_scalar ** 2; + // For real signals, one-sided PSD needs doubling for non-DC/Nyquist frequencies + // Assuming fftOutputSegments[j][i] are magnitudes from a one-sided spectrum (0 to Nyquist) + // DC (i=0) and Nyquist (i=dataCount-1, if dataCount = N/2+1) are not doubled. + if (i !== 0 && (dataCount % 2 === 1 ? i !== dataCount -1 : true) ) { // crude Nyquist check for N/2+1 length + // A common rule: double all bins except DC and Nyquist. + // If dataCount is N/2, all are doubled except DC. + // If dataCount is N/2+1, Nyquist is at dataCount-1. + // Let's assume dataCount from _fft_segmented is N/2+1. + if (i < dataCount -1) { // Don't double Nyquist if it's the last bin + p *= 2; + } } psdOutput[i] += p; } - const min_avg = 1e-7; // limit min value for -70db + const min_avg = 1e-10; // limit min value for -100db (increased range) let avg = psdOutput[i] / segmentsCount; avg = Math.max(avg, min_avg); psdOutput[i] = 10 * Math.log10(avg); @@ -589,11 +697,11 @@ GraphSpectrumCalc._psd = function(samples, pointsPerSegment, overlapCount, scal } if (i > noise50HzIdx && psdOutput[i] > maxNoise) { maxNoise = psdOutput[i]; - maxNoiseIdx = i; + maxNoiseIdxBin = i; } } - const maxNoiseFrequency = maxNoiseIdx / dataCount * maxFrequency; + const maxNoiseFrequency = maxNoiseIdxBin / (dataCount -1) * maxFrequency; // Scale bin to freq return { psdOutput: psdOutput, @@ -609,7 +717,10 @@ GraphSpectrumCalc._psd = function(samples, pointsPerSegment, overlapCount, scal */ GraphSpectrumCalc._fft_segmented = function(samples, pointsPerSegment, overlapCount) { const samplesCount = samples.length; - let output = []; + let output = []; // Array of magnitude arrays + if (pointsPerSegment === 0 || samplesCount < pointsPerSegment || pointsPerSegment < MIN_FFT_POINTS) { + return output; + } for (let i = 0; i <= samplesCount - pointsPerSegment; i += pointsPerSegment - overlapCount) { const fftInput = samples.slice(i, i + pointsPerSegment); @@ -617,15 +728,18 @@ GraphSpectrumCalc._fft_segmented = function(samples, pointsPerSegment, overlapC 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); + const fftComplex = this._fft(fftInput); // Returns complex array [re,im,re,im...] + // For N real inputs, FFT.js simple('real') returns N/2+1 complex values. + // So, fftComplex has (pointsPerSegment/2 + 1) * 2 float values. + const numMagnitudes = pointsPerSegment / 2 + 1; + const magnitudes = new Float64Array(numMagnitudes); + for (let k = 0; k < numMagnitudes; k++) { + const re = fftComplex[2 * k]; + const im = fftComplex[2 * k + 1]; + magnitudes[k] = Math.hypot(re, im); } output.push(magnitudes); } return output; -}; +}; \ No newline at end of file From dc69480c6d2917fd769c093ac9ac48249f052165 Mon Sep 17 00:00:00 2001 From: nerdCopter <56646290+nerdCopter@users.noreply.github.com> Date: Tue, 13 May 2025 08:32:52 -0500 Subject: [PATCH 2/5] newline --- src/graph_spectrum_calc.js | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/graph_spectrum_calc.js b/src/graph_spectrum_calc.js index f0f85b29..4e6007e6 100644 --- a/src/graph_spectrum_calc.js +++ b/src/graph_spectrum_calc.js @@ -742,4 +742,4 @@ GraphSpectrumCalc._fft_segmented = function(samples, pointsPerSegment, overlapC } return output; -}; \ No newline at end of file +}; From 73dad8f14fb44c7a541ce610c721b1b7e5444cb0 Mon Sep 17 00:00:00 2001 From: nerdCopter <56646290+nerdCopter@users.noreply.github.com> Date: Tue, 13 May 2025 09:45:41 -0500 Subject: [PATCH 3/5] @coderabbitai fixed --- src/graph_spectrum_calc.js | 62 +++++++++++++++----------------------- 1 file changed, 25 insertions(+), 37 deletions(-) diff --git a/src/graph_spectrum_calc.js b/src/graph_spectrum_calc.js index 4e6007e6..00d2b893 100644 --- a/src/graph_spectrum_calc.js +++ b/src/graph_spectrum_calc.js @@ -567,53 +567,41 @@ GraphSpectrumCalc._fft = function(samples, type) { * Makes all the values absolute and returns the index of maxFrequency found */ GraphSpectrumCalc._normalizeFft = function(fftOutput, fftLength) { - // fftLength is the number of samples input to FFT. - // fftOutput is the complex result from _fft. - // The original _normalizeFft looped 'fftLength' times assuming 'fftOutput' contained 'fftLength' magnitudes. - // This behavior is preserved. If fftOutput from _fft contains N complex values (2N floats), - // and fftLength is N, this loop processes N floats from fftOutput, taking their abs value. - - // Number of actual magnitudes to consider is fftLength / 2 (approx, up to Nyquist) - // However, to match existing behavior where it iterated fftLength times over fftOutput: - const iterationLength = fftLength; // Or Math.min(fftLength, fftOutput.length) if fftOutput could be shorter + // Number of usable frequency bins (0 to Nyquist) + const bins = Math.floor(fftLength / 2) + 1; // +1 to include Nyquist bin + const magnitudes = new Float64Array(bins); + + // Calculate magnitudes from complex values + for (let bin = 0; bin < bins; bin++) { + const re = fftOutput[2 * bin]; + const im = fftOutput[2 * bin + 1]; + magnitudes[bin] = Math.hypot(re, im); + } - // Make all the values absolute, and calculate some useful values (max noise, etc.) + // Find max noise after low-end cutoff const maxFrequency = (this._blackBoxRate / 2.0); - // This scaling implies fftLength is the number of points up to Nyquist. - // If fftLength is num_samples, then num_points_to_nyquist is num_samples / 2. - // To maintain existing scaling logic with fftLength = num_samples: - const noiseLowEndIdx = 100 / maxFrequency * (iterationLength / 2); // Adjusted for N/2 bins + const noiseLowEndIdx = Math.floor(100 / maxFrequency * bins); let maxNoiseIdx = 0; let maxNoise = 0; - // Create a new array for magnitudes if we are to correctly interpret complex fftOutput - // For minimal change, we continue the existing pattern of modifying fftOutput in place. - // The loop below is likely incorrect for complex fftOutput but is existing behavior. - for (let i = 0; i < iterationLength; i++) { // This loop is over first N floats of complex output - fftOutput[i] = Math.abs(fftOutput[i]); - // The condition for maxNoiseIdx should refer to frequency bins, not raw float indices - const currentFreqBin = i / 2; // Approximate, assumes pairs - if (currentFreqBin > noiseLowEndIdx && fftOutput[i] > maxNoise) { - maxNoise = fftOutput[i]; - maxNoiseIdx = currentFreqBin; // Store bin index + for (let bin = 0; bin < bins; bin++) { + if (bin > noiseLowEndIdx && magnitudes[bin] > maxNoise) { + maxNoise = magnitudes[bin]; + maxNoiseIdx = bin; } } - // maxNoiseIdx is now a bin index relative to N/2 bins. - // Scale this bin index to frequency. - maxNoiseIdx = maxNoiseIdx / (iterationLength / 2) * maxFrequency; - + // Scale bin index to frequency + const maxNoiseFreq = maxNoiseIdx / bins * maxFrequency; - const fftData = { - fieldIndex : this._dataBuffer.fieldIndex, - fieldName : this._dataBuffer.fieldName, - fftLength : iterationLength, // Or iterationLength / 2 if this means number of frequency bins - fftOutput : fftOutput.slice(0, iterationLength), // Return the part we processed - maxNoiseIdx : maxNoiseIdx, - blackBoxRate : this._blackBoxRate, + return { + fieldIndex: this._dataBuffer.fieldIndex, + fieldName: this._dataBuffer.fieldName, + fftLength: bins, // Return number of frequency bins + fftOutput: magnitudes, // Return the magnitude spectrum + maxNoiseIdx: maxNoiseFreq, + blackBoxRate: this._blackBoxRate, }; - - return fftData; }; /** From 6f6d279dae9194ed5e5868d5c4d4d5ac3bcea8b6 Mon Sep 17 00:00:00 2001 From: nerdCopter <56646290+nerdCopter@users.noreply.github.com> Date: Tue, 13 May 2025 15:38:11 -0500 Subject: [PATCH 4/5] zero fill to next 2^N --- src/graph_spectrum_calc.js | 289 +++++++++++++++++++++---------------- 1 file changed, 162 insertions(+), 127 deletions(-) diff --git a/src/graph_spectrum_calc.js b/src/graph_spectrum_calc.js index aad2669f..9667b3f8 100644 --- a/src/graph_spectrum_calc.js +++ b/src/graph_spectrum_calc.js @@ -45,6 +45,15 @@ export const GraphSpectrumCalc = { } return p; }, + + _getNextPowerOfTwo: function(n) { + if (n <= 0) return 1; // Return 1 for non-positive inputs, MIN_FFT_POINTS check will handle later. + let p = 1; + while (p < n) { + p *= 2; + } + return p; + }, }; GraphSpectrumCalc.initialize = function(flightLog, sysConfig) { @@ -89,7 +98,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; }; @@ -103,10 +112,11 @@ GraphSpectrumCalc.dataLoadFrequency = function() { const flightSamples = this._getFlightSamplesFreq(); - let effectiveSampleCount = this._getLargestPowerOfTwoLeq(flightSamples.count); + // Zero-fill to the next power of two + const targetFftLength = this._getNextPowerOfTwo(flightSamples.count); - if (effectiveSampleCount < MIN_FFT_POINTS) { - console.warn(`Effective sample count ${effectiveSampleCount} for FFT is too small (min ${MIN_FFT_POINTS}). Skipping Frequency analysis.`); + if (targetFftLength < MIN_FFT_POINTS || flightSamples.count === 0) { // Also check flightSamples.count to avoid FFT on 0 samples padded to MIN_FFT_POINTS + console.warn(`Target FFT length ${targetFftLength} (from ${flightSamples.count} samples) is too small (min ${MIN_FFT_POINTS}). Skipping Frequency analysis.`); return { fieldIndex: this._dataBuffer.fieldIndex, fieldName: this._dataBuffer.fieldName, @@ -117,17 +127,22 @@ GraphSpectrumCalc.dataLoadFrequency = function() { }; } - const samplesForFft = flightSamples.samples.slice(0, effectiveSampleCount); + const samplesForFft = new Float64Array(targetFftLength); + // Copy actual samples, the rest of samplesForFft will be zeros + samplesForFft.set(flightSamples.samples.slice(0, flightSamples.count)); + if (userSettings.analyserHanning) { - this._hanningWindow(samplesForFft, effectiveSampleCount); + // Apply Hanning window to the entire (potentially zero-padded) block + this._hanningWindow(samplesForFft, targetFftLength); } //calculate fft - const fftOutput = this._fft(samplesForFft); + const fftOutput = this._fft(samplesForFft); // samplesForFft has targetFftLength // Normalize the result - const fftData = this._normalizeFft(fftOutput, effectiveSampleCount); + // Pass targetFftLength as the length of the data that underwent FFT + const fftData = this._normalizeFft(fftOutput, targetFftLength); return fftData; }; @@ -182,13 +197,14 @@ GraphSpectrumCalc._dataLoadFrequencyVsX = function(vsFieldNames, minValue = Infi const flightSamples = this._getFlightSamplesFreqVsX(vsFieldNames, minValue, maxValue); - // 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. - let fftChunkLength = Math.floor(this._blackBoxRate * FREQ_VS_THR_CHUNK_TIME_MS / 1000); - fftChunkLength = this._getLargestPowerOfTwoLeq(fftChunkLength); + // Number of actual samples desired for each FFT chunk, based on time + const desiredSamplesPerChunk = Math.floor(this._blackBoxRate * FREQ_VS_THR_CHUNK_TIME_MS / 1000); + // FFT length will be the next power of two from desiredSamplesPerChunk + const fftChunkLength = this._getNextPowerOfTwo(desiredSamplesPerChunk); - if (fftChunkLength < MIN_FFT_POINTS) { - console.warn(`FFT chunk length ${fftChunkLength} for FreqVsX is too small (min ${MIN_FFT_POINTS}). Skipping.`); + + if (fftChunkLength < MIN_FFT_POINTS || desiredSamplesPerChunk === 0) { + console.warn(`FFT chunk length ${fftChunkLength} (from ${desiredSamplesPerChunk} desired samples) for FreqVsX is too small (min ${MIN_FFT_POINTS}). Skipping.`); return { fieldIndex: this._dataBuffer.fieldIndex, fieldName: this._dataBuffer.fieldName, @@ -200,81 +216,75 @@ GraphSpectrumCalc._dataLoadFrequencyVsX = function(vsFieldNames, minValue = Infi }; } - const fftChunkWindow = Math.round(fftChunkLength / FREQ_VS_THR_WINDOW_DIVISOR); + // Window step for moving window, based on the (padded) FFT chunk length + const fftChunkWindowStep = Math.round(fftChunkLength / FREQ_VS_THR_WINDOW_DIVISOR); let maxNoise = 0; // Stores the maximum amplitude of the fft over all chunks // Matrix where each row represents a bin of vs values, and the columns are amplitudes at frequencies - const matrixFftOutput = new Array(NUM_VS_BINS).fill(null).map(() => new Float64Array(fftChunkLength * 2)); + // The length of rows should match the length of the sliced fftOutput (which is fftChunkLength) + const matrixFftOutputSums = new Array(NUM_VS_BINS).fill(null).map(() => new Float64Array(fftChunkLength)); + + const numberSamplesInBin = new Uint32Array(NUM_VS_BINS); // Number of FFT segments contributing to each vs bin - const numberSamples = new Uint32Array(NUM_VS_BINS); // Number of samples in each vs value, used to average them later. + const fft = new FFT.complex(fftChunkLength, false); // FFT initialized for padded length - const fft = new FFT.complex(fftChunkLength, false); - for (let fftChunkIndex = 0; fftChunkIndex + fftChunkLength < flightSamples.count; fftChunkIndex += fftChunkWindow) { + // Loop through flight samples, taking chunks of 'desiredSamplesPerChunk' + // The loop must ensure there are enough samples for 'desiredSamplesPerChunk' + for (let fftChunkStartIndex = 0; (fftChunkStartIndex + desiredSamplesPerChunk) <= flightSamples.count; fftChunkStartIndex += fftChunkWindowStep) { - const fftInput = flightSamples.samples.slice(fftChunkIndex, fftChunkIndex + fftChunkLength); - let fftOutput = new Float64Array(fftChunkLength * 2); // This is for complex output + const actualSamplesInChunk = flightSamples.samples.slice(fftChunkStartIndex, fftChunkStartIndex + desiredSamplesPerChunk); + const fftInput = new Float64Array(fftChunkLength); // Prepare input array for FFT (zero-padded) + fftInput.set(actualSamplesInChunk); // Copy actual samples, rest are zeros - // Hanning window applied to input data + let fftOutputComplex = new Float64Array(fftChunkLength * 2); // For complex output from FFT.js + + // Hanning window applied to the entire (zero-padded) input block if (userSettings.analyserHanning) { this._hanningWindow(fftInput, fftChunkLength); } - fft.simple(fftOutput, fftInput, 'real'); - - // The original code sliced fftOutput (complex) to its first fftChunkLength float values. - // This behavior is preserved with the new power-of-2 fftChunkLength. - fftOutput = fftOutput.slice(0, fftChunkLength); + fft.simple(fftOutputComplex, fftInput, 'real'); - - // Use only abs values - for (let i = 0; i < fftChunkLength; i++) { // This loop iterates over the sliced (potentially mixed Re/Im) values - fftOutput[i] = Math.abs(fftOutput[i]); - maxNoise = Math.max(fftOutput[i], maxNoise); + // The original code sliced fftOutputComplex (complex) to its first fftChunkLength float values. + // This behavior is preserved. fftOutputProcessed will store these abs values. + const fftOutputProcessed = new Float64Array(fftChunkLength); + for (let i = 0; i < fftChunkLength; i++) { + fftOutputProcessed[i] = Math.abs(fftOutputComplex[i]); // This processes Re0, Im0, Re1, Im1 ... + maxNoise = Math.max(fftOutputProcessed[i], maxNoise); } - // calculate a bin index and put the fft value in that bin for each field (e.g. eRPM[0], eRPM[1]..) sepparately + // Calculate a bin index and add the processed fft values to that bin for (const vsValueArray of flightSamples.vsValues) { - // Calculate average of the VS values in the chunk + // Calculate average of the VS values over the actual data duration in the chunk let sumVsValues = 0; - for (let indexVs = fftChunkIndex; indexVs < fftChunkIndex + fftChunkLength; indexVs++) { + for (let indexVs = fftChunkStartIndex; indexVs < fftChunkStartIndex + desiredSamplesPerChunk; indexVs++) { sumVsValues += vsValueArray[indexVs]; } - // Translate the average vs value to a bin index - const avgVsValue = sumVsValues / fftChunkLength; + const avgVsValue = sumVsValues / desiredSamplesPerChunk; // Average over actual samples + 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; } - if (vsBinIndex < 0) { vsBinIndex = 0;} // Ensure valid index - - numberSamples[vsBinIndex]++; - - // add the output from the fft to the row given by the vs value bin index - // The matrixFftOutput columns correspond to the elements from the sliced fftOutput - for (let i = 0; i < fftOutput.length; i++) { - // Ensure matrixFftOutput[vsBinIndex] is initialized if it wasn't (though map should handle it) - if (!matrixFftOutput[vsBinIndex]) matrixFftOutput[vsBinIndex] = new Float64Array(fftChunkLength); // Safety, matching slice - matrixFftOutput[vsBinIndex][i] += fftOutput[i]; + if (vsBinIndex < 0) { vsBinIndex = 0;} + if (vsBinIndex >=0 && vsBinIndex < NUM_VS_BINS) { // Ensure valid index before accessing + numberSamplesInBin[vsBinIndex]++; + for (let i = 0; i < fftOutputProcessed.length; i++) { + matrixFftOutputSums[vsBinIndex][i] += fftOutputProcessed[i]; + } } } } - // Adjust matrixFftOutput dimensions to match the sliced fftOutput length if different from initialization - // This should not be necessary if matrixFftOutput was initialized with fftChunkLength (length of the sliced data). - // Let's ensure matrixFftOutput has the correct number of columns based on the sliced fftOutput. - // The initialization `new Float64Array(fftChunkLength * 2)` for matrixFftOutput rows - // should actually be `new Float64Array(fftChunkLength)` because `fftOutput` is sliced to that length. - const finalMatrixFftOutput = new Array(NUM_VS_BINS).fill(null).map(() => new Float64Array(fftChunkLength)); + const finalMatrixFftOutput = new Array(NUM_VS_BINS).fill(null).map(() => new Float64Array(fftChunkLength)); - // Divide the values from the fft in each row (vs value bin) by the number of samples in the bin + // Average the summed FFT outputs for each bin for (let i = 0; i < NUM_VS_BINS; i++) { - if (numberSamples[i] > 1) { - for (let j = 0; j < fftChunkLength; j++) { // Iterate up to the length of the sliced fftOutput - finalMatrixFftOutput[i][j] = matrixFftOutput[i][j] / numberSamples[i]; + if (numberSamplesInBin[i] > 0) { // Check for > 0 to handle division correctly + for (let j = 0; j < fftChunkLength; j++) { + finalMatrixFftOutput[i][j] = matrixFftOutputSums[i][j] / numberSamplesInBin[i]; } - } else if (numberSamples[i] === 1) { - for (let j = 0; j < fftChunkLength; j++) { - finalMatrixFftOutput[i][j] = matrixFftOutput[i][j]; - } + } else { + // If no samples in this bin, finalMatrixFftOutput[i] remains an array of zeros (or fill with null/NaN if preferred) + // Default Float64Array is zero-filled, which is fine. } } @@ -282,7 +292,7 @@ GraphSpectrumCalc._dataLoadFrequencyVsX = function(vsFieldNames, minValue = Infi const fftData = { fieldIndex : this._dataBuffer.fieldIndex, fieldName : this._dataBuffer.fieldName, - fftLength : fftChunkLength, // This is the length of the (sliced) data used for matrix rows + fftLength : fftChunkLength, // This is the length of the (padded) data used for matrix rows fftOutput : finalMatrixFftOutput, maxNoise : maxNoise, blackBoxRate : this._blackBoxRate, @@ -401,7 +411,8 @@ GraphSpectrumCalc._getFlightSamplesFreq = function() { } return { - samples : samples, + samples : samples, // Note: This array is oversized, only 'count' elements are valid. + // Slicing (e.g., samples.slice(0, samplesCount)) is done by callers. count : samplesCount, }; }; @@ -446,23 +457,24 @@ GraphSpectrumCalc._getFlightSamplesFreqVsX = function(vsFieldNames, minValue = I } } - // Calculate min max average of the VS values in the chunk what will used by spectrum data definition - let tempfftChunkLength = Math.floor(this._blackBoxRate * FREQ_VS_THR_CHUNK_TIME_MS / 1000); // temp for calc, actual fftChunkLength may change - tempfftChunkLength = this._getLargestPowerOfTwoLeq(tempfftChunkLength); - if (tempfftChunkLength < MIN_FFT_POINTS) tempfftChunkLength = MIN_FFT_POINTS; // Use min for range calc if too small + // Calculate min/max for VS values based on actual chunk duration, respecting MIN_FFT_POINTS as a minimum window. + let samplesPerSegmentForVsAvg = Math.floor(this._blackBoxRate * FREQ_VS_THR_CHUNK_TIME_MS / 1000); + // Ensure the segment for averaging is not smaller than MIN_FFT_POINTS, if MIN_FFT_POINTS is meaningful. + if (MIN_FFT_POINTS > 0) { + samplesPerSegmentForVsAvg = Math.max(samplesPerSegmentForVsAvg, MIN_FFT_POINTS); + } + - const fftChunkWindow = Math.round(tempfftChunkLength / FREQ_VS_THR_WINDOW_DIVISOR); + const windowStepForVsAvg = Math.round(samplesPerSegmentForVsAvg / FREQ_VS_THR_WINDOW_DIVISOR); - if (tempfftChunkLength > 0 && fftChunkWindow > 0) { // Proceed only if chunking is possible - for (let fftChunkIndex = 0; fftChunkIndex + tempfftChunkLength < samplesCount; fftChunkIndex += fftChunkWindow) { + if (samplesPerSegmentForVsAvg > 0 && windowStepForVsAvg > 0 && samplesCount >= samplesPerSegmentForVsAvg) { + for (let segmentStartIndex = 0; segmentStartIndex + samplesPerSegmentForVsAvg <= samplesCount; segmentStartIndex += windowStepForVsAvg) { for (const vsValueArray of vsValues) { - // Calculate average of the VS values in the chunk let sumVsValues = 0; - for (let indexVs = fftChunkIndex; indexVs < fftChunkIndex + tempfftChunkLength; indexVs++) { + for (let indexVs = segmentStartIndex; indexVs < segmentStartIndex + samplesPerSegmentForVsAvg; indexVs++) { sumVsValues += vsValueArray[indexVs]; } - // Find min max average of the VS values in the chunk - const avgVsValue = sumVsValues / tempfftChunkLength; + const avgVsValue = sumVsValues / samplesPerSegmentForVsAvg; maxValue = Math.max(maxValue, avgVsValue); minValue = Math.min(minValue, avgVsValue); } @@ -472,8 +484,8 @@ GraphSpectrumCalc._getFlightSamplesFreqVsX = function(vsFieldNames, minValue = I maxValue *= MAX_RPM_AXIS_GAP; - if (minValue > maxValue || minValue === Infinity || maxValue === -Infinity) { - // Fallback if range is invalid + if (minValue > maxValue || minValue === Infinity || maxValue === -Infinity || samplesCount === 0) { + // Fallback if range is invalid or no samples minValue = 0; if (vsFieldNames === FIELD_THROTTLE_NAME) { maxValue = 100; @@ -483,16 +495,16 @@ GraphSpectrumCalc._getFlightSamplesFreqVsX = function(vsFieldNames, minValue = I maxValue = 100; // Default generic max } if (minValue > maxValue) minValue = 0; // Ensure min <= max - console.log("Adjusted FreqVsX range due to invalid calculation: min=%f, max=%f", minValue, maxValue); + console.warn("Adjusted FreqVsX range due to invalid calculation or no samples: min=%f, max=%f", minValue, maxValue); } return { - samples : samples, - vsValues : vsValues, + samples : samples, // Note: This array is oversized. Callers should use .slice(0, count) + vsValues : vsValues, // Note: This array is oversized. count : samplesCount, minValue : minValue, maxValue : maxValue, - vsRange: {min: minValue, max: maxValue}, // ensure vsRange is present + vsRange: {min: minValue, max: maxValue}, }; }; @@ -522,15 +534,16 @@ GraphSpectrumCalc._getFlightSamplesPidErrorVsSetpoint = function(axisIndex) { } return { - piderror, - setpoint, + piderror, // Note: oversized + setpoint, // Note: oversized maxSetpoint, count: samplesCount, }; }; GraphSpectrumCalc._hanningWindow = function(samples, size) { - + // Expects 'samples' to be an array of at least 'size' length. + // Modifies 'samples' in-place for the first 'size' elements. if (!size) { size = samples.length; } @@ -541,20 +554,21 @@ GraphSpectrumCalc._hanningWindow = function(samples, size) { }; GraphSpectrumCalc._fft = function(samples, type) { - + // 'samples' is expected to be a Float64Array of power-of-2 length. if (!type) { type = 'real'; } const fftLength = samples.length; - const fftOutput = new Float64Array(fftLength * 2); // Expects N complex outputs for N real inputs by some conventions - // or N/2+1 complex outputs. FFT.js simple 'real' typically N/2+1. - // The library may internally handle sizing of fftOutput for simple(). + // FFT.js complex (Nayuki) for N real inputs produces N/2+1 complex outputs. + // The output buffer needs to be large enough for this: (fftLength/2 + 1) * 2 floats. + // A buffer of size fftLength*2 is more than enough. + const fftOutput = new Float64Array(fftLength * 2); const fft = new FFT.complex(fftLength, false); fft.simple(fftOutput, samples, type); - return fftOutput; + return fftOutput; // Contains (fftLength/2 + 1) complex values, i.e., (fftLength/2 + 1)*2 floats. }; @@ -562,8 +576,10 @@ GraphSpectrumCalc._fft = function(samples, type) { * Makes all the values absolute and returns the index of maxFrequency found */ GraphSpectrumCalc._normalizeFft = function(fftOutput, fftLength) { - // Number of usable frequency bins (0 to Nyquist) - const bins = Math.floor(fftLength / 2) + 1; // +1 to include Nyquist bin + // fftLength is the length of the original (potentially zero-padded) time-domain signal. + // fftOutput contains complex FFT results: [Re0, Im0, Re1, Im1, ...]. + // For a real input of fftLength, there are (fftLength/2 + 1) unique complex values. + const bins = Math.floor(fftLength / 2) + 1; const magnitudes = new Float64Array(bins); // Calculate magnitudes from complex values @@ -575,11 +591,16 @@ GraphSpectrumCalc._normalizeFft = function(fftOutput, fftLength) { // Find max noise after low-end cutoff const maxFrequency = (this._blackBoxRate / 2.0); - const noiseLowEndIdx = Math.floor(100 / maxFrequency * bins); + // Avoid division by zero if maxFrequency is 0 (e.g. blackBoxRate is 0) + const noiseLowEndIdx = (maxFrequency > 0) ? Math.floor(100 / maxFrequency * bins) : 0; + let maxNoiseIdx = 0; let maxNoise = 0; for (let bin = 0; bin < bins; bin++) { + // Original logic was > noiseLowEndIdx. If noiseLowEndIdx can be 0, this means all bins are checked. + // If we want to exclude DC (bin 0) and very low frequencies, the condition should be robust. + // Assuming noiseLowEndIdx is a valid start index for noise search. if (bin > noiseLowEndIdx && magnitudes[bin] > maxNoise) { maxNoise = magnitudes[bin]; maxNoiseIdx = bin; @@ -587,14 +608,20 @@ GraphSpectrumCalc._normalizeFft = function(fftOutput, fftLength) { } // Scale bin index to frequency - const maxNoiseFreq = maxNoiseIdx / bins * maxFrequency; + // Avoid division by zero if bins is 1 (e.g. fftLength is 0 or 1) + const maxNoiseFreq = (bins > 1) ? (maxNoiseIdx / (bins -1) * maxFrequency) : 0; + // Note: A common scaling is maxNoiseIdx * (blackBoxRate / fftLength). + // bins-1 corresponds to Nyquist. maxFrequency = blackBoxRate / 2. + // maxNoiseIdx / (bins-1) * (blackBoxRate/2) = maxNoiseIdx * blackBoxRate / (2*(bins-1)) + // If bins = fftLength/2+1, then 2*(bins-1) = 2*(fftLength/2) = fftLength. So this is consistent. + return { fieldIndex: this._dataBuffer.fieldIndex, fieldName: this._dataBuffer.fieldName, - fftLength: bins, // Return number of frequency bins + fftLength: bins, // Return number of frequency bins (magnitudes) fftOutput: magnitudes, // Return the magnitude spectrum - maxNoiseIdx: maxNoiseFreq, + maxNoiseIdx: maxNoiseFreq, // This is now a frequency in Hz blackBoxRate: this._blackBoxRate, }; }; @@ -610,33 +637,33 @@ GraphSpectrumCalc._psd = function(samples, pointsPerSegment, overlapCount, scal return { psdOutput: new Float64Array(0), min: 0, max: 0, maxNoiseIdx: 0 }; } - const dataCount = fftOutputSegments[0].length; // Number of magnitude points per segment (N/2 or N/2+1) + const dataCount = fftOutputSegments[0].length; // Number of magnitude points per segment (N/2+1 typically) const segmentsCount = fftOutputSegments.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); + const window = Array(pointsPerSegment).fill(1); // Create a dummy window array + this._hanningWindow(window, pointsPerSegment); // Apply Hanning values to it if (scaling == 'density') { let skSum = 0; for (const value of window) { skSum += value ** 2; } - scale = 1 / (this._blackBoxRate * skSum); + scale = (this._blackBoxRate * skSum > 0) ? (1 / (this._blackBoxRate * skSum)) : 1; } else if (scaling == 'spectrum') { let sum = 0; for (const value of window) { sum += value; } - scale = 1 / sum ** 2; + scale = (sum ** 2 > 0) ? (1 / sum ** 2) : 1; } } else { // No Hanning window if (scaling == 'density') { - scale = 1 / (this._blackBoxRate * pointsPerSegment); // Corrected scale for non-windowed + scale = (this._blackBoxRate * pointsPerSegment > 0) ? (1 / (this._blackBoxRate * pointsPerSegment)) : 1; } else if (scaling == 'spectrum') { - scale = 1 / pointsPerSegment ** 2; + scale = (pointsPerSegment ** 2 > 0) ? (1 / pointsPerSegment ** 2) : 1; } } @@ -646,8 +673,10 @@ GraphSpectrumCalc._psd = function(samples, pointsPerSegment, overlapCount, scal max = -1e6; const maxFrequency = (this._blackBoxRate / 2.0); - const noise50HzIdx = Math.floor(50 / maxFrequency * dataCount); - const noise3HzIdx = Math.floor(3 / maxFrequency * dataCount); + // Avoid division by zero for index calculations + const noise50HzIdx = (maxFrequency > 0) ? Math.floor(50 / maxFrequency * dataCount) : 0; + const noise3HzIdx = (maxFrequency > 0) ? Math.floor(3 / maxFrequency * dataCount) : 0; + let maxNoiseIdxBin = 0; let maxNoise = -100; // dB for (let i = 0; i < dataCount; i++) { @@ -655,26 +684,22 @@ GraphSpectrumCalc._psd = function(samples, pointsPerSegment, overlapCount, scal for (let j = 0; j < segmentsCount; j++) { let p_scalar = fftOutputSegments[j][i]; // This is already magnitude from _fft_segmented let p = scale * p_scalar ** 2; - // For real signals, one-sided PSD needs doubling for non-DC/Nyquist frequencies - // Assuming fftOutputSegments[j][i] are magnitudes from a one-sided spectrum (0 to Nyquist) - // DC (i=0) and Nyquist (i=dataCount-1, if dataCount = N/2+1) are not doubled. - if (i !== 0 && (dataCount % 2 === 1 ? i !== dataCount -1 : true) ) { // crude Nyquist check for N/2+1 length - // A common rule: double all bins except DC and Nyquist. - // If dataCount is N/2, all are doubled except DC. - // If dataCount is N/2+1, Nyquist is at dataCount-1. - // Let's assume dataCount from _fft_segmented is N/2+1. - if (i < dataCount -1) { // Don't double Nyquist if it's the last bin - p *= 2; - } + + // One-sided PSD needs doubling for non-DC/Nyquist frequencies. + // dataCount here is N/2+1 (number of magnitudes from _fft_segmented). + // DC is at index 0. Nyquist is at index dataCount-1. + if (i !== 0 && i !== dataCount - 1) { + p *= 2; } psdOutput[i] += p; } - const min_avg = 1e-10; // limit min value for -100db (increased range) - let avg = psdOutput[i] / segmentsCount; + const min_avg = 1e-10; + let avg = (segmentsCount > 0) ? (psdOutput[i] / segmentsCount) : 0; avg = Math.max(avg, min_avg); psdOutput[i] = 10 * Math.log10(avg); - if (i > noise3HzIdx) { // Miss big zero freq magnitude + + if (i > noise3HzIdx) { min = Math.min(psdOutput[i], min); max = Math.max(psdOutput[i], max); } @@ -683,14 +708,14 @@ GraphSpectrumCalc._psd = function(samples, pointsPerSegment, overlapCount, scal maxNoiseIdxBin = i; } } - - const maxNoiseFrequency = maxNoiseIdxBin / (dataCount -1) * maxFrequency; // Scale bin to freq + // Scale bin to freq. (dataCount-1) corresponds to Nyquist frequency. + const maxNoiseFrequency = (dataCount > 1) ? (maxNoiseIdxBin / (dataCount -1) * maxFrequency) : 0; return { psdOutput: psdOutput, min: min, max: max, - maxNoiseIdx: maxNoiseFrequency, + maxNoiseIdx: maxNoiseFrequency, // This is a frequency in Hz }; }; @@ -704,16 +729,26 @@ GraphSpectrumCalc._fft_segmented = function(samples, pointsPerSegment, overlapC if (pointsPerSegment === 0 || samplesCount < pointsPerSegment || pointsPerSegment < MIN_FFT_POINTS) { return output; } - for (let i = 0; i <= samplesCount - pointsPerSegment; i += pointsPerSegment - overlapCount) { - const fftInput = samples.slice(i, i + pointsPerSegment); + // Ensure pointsPerSegment is a power of two for FFT (if not already guaranteed by caller) + // This function's caller (dataLoadPSD) already ensures pointsPerSegment is a power of two. + + const step = pointsPerSegment - overlapCount; + if (step <= 0) { // Avoid infinite loop if overlap is too large + console.warn("PSD segment step is non-positive, review pointsPerSegment and overlapCount."); + return output; + } + + for (let i = 0; (i + pointsPerSegment) <= samplesCount; i += step) { + // Take a copy for windowing and FFT, as _hanningWindow modifies in place + const fftInputSegment = samples.slice(i, i + pointsPerSegment); if (userSettings.analyserHanning) { - this._hanningWindow(fftInput, pointsPerSegment); + this._hanningWindow(fftInputSegment, pointsPerSegment); // Modifies fftInputSegment } - const fftComplex = this._fft(fftInput); // Returns complex array [re,im,re,im...] - // For N real inputs, FFT.js simple('real') returns N/2+1 complex values. - // So, fftComplex has (pointsPerSegment/2 + 1) * 2 float values. + const fftComplex = this._fft(fftInputSegment); // Returns complex array [re,im,re,im...] + // It will have (pointsPerSegment/2 + 1)*2 float values. + const numMagnitudes = pointsPerSegment / 2 + 1; const magnitudes = new Float64Array(numMagnitudes); for (let k = 0; k < numMagnitudes; k++) { From f00b0aa8740dad6b7f08a01dfc95b5499cee6e81 Mon Sep 17 00:00:00 2001 From: nerdCopter <56646290+nerdCopter@users.noreply.github.com> Date: Tue, 13 May 2025 15:56:45 -0500 Subject: [PATCH 5/5] revert Calculate proper complex magnitudes due to merge-master. --- src/graph_spectrum_calc.js | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/src/graph_spectrum_calc.js b/src/graph_spectrum_calc.js index 9667b3f8..fb0f5f07 100644 --- a/src/graph_spectrum_calc.js +++ b/src/graph_spectrum_calc.js @@ -245,12 +245,15 @@ GraphSpectrumCalc._dataLoadFrequencyVsX = function(vsFieldNames, minValue = Infi fft.simple(fftOutputComplex, fftInput, 'real'); - // The original code sliced fftOutputComplex (complex) to its first fftChunkLength float values. - // This behavior is preserved. fftOutputProcessed will store these abs values. - const fftOutputProcessed = new Float64Array(fftChunkLength); - for (let i = 0; i < fftChunkLength; i++) { - fftOutputProcessed[i] = Math.abs(fftOutputComplex[i]); // This processes Re0, Im0, Re1, Im1 ... - maxNoise = Math.max(fftOutputProcessed[i], maxNoise); + // Calculate proper complex magnitudes. + const bins = fftChunkLength / 2 + 1; + const fftOutputProcessed = new Float64Array(bins); + for (let bin = 0; bin < bins; bin++) { + const re = fftOutputComplex[2 * bin]; + const im = fftOutputComplex[2 * bin + 1]; + const mag = Math.hypot(re, im); + fftOutputProcessed[bin] = mag; + maxNoise = Math.max(maxNoise, mag); } // Calculate a bin index and add the processed fft values to that bin