fair-acc / chart-fx

A scientific charting library focused on performance optimised real-time data visualisation at 25 Hz update rates for data sets with a few 10 thousand up to 5 million data points.
GNU Lesser General Public License v3.0
488 stars 90 forks source link

ShortTimeFourierTransform: Amplitude change with zero padding #607

Open Vado101 opened 10 months ago

Vado101 commented 10 months ago

Hi! I'm trying make spectrum analyzer using chart-fx library. When I use the real method of the chartfx-math/src/main/java/io/fair_acc/math/spectra/ShortTimeFourierTransform.java class with the FFT size greater than the data length and padding with zeros, then when applying different window functions, the signal amplitude noticeably different. I tried to implement my real method, where I placed the application of the window function above the code in which the padding with zeros occurs, and also calculated the amplitude compensation. In this case, when applying different window functions, the signal level is approximately the same. Please tell me if the real method works correctly and if so, why then the amplitude changes with different window functions.

Below is my modified real method

public static double[] real(double[] input, double[] output, int nFFT, int step, Apodization apodization,
                                ShortTimeFourierTransform.Padding padding, boolean dbScale, boolean truncateDCNy) {
      int nT = ceilDiv(input.length, step);
      double[] amplitudeData = output != null && output.length == nFFT / 2 * nT ? output : new double[nFFT / 2 * nT];
      double[] currentMagnitudeData = DoubleArrayCache.getInstance().getArray(nFFT / 2);
      DoubleFFT_1D fastFourierTransform = new DoubleFFT_1D((long)nFFT);
      double[] raw = DoubleArrayCache.getInstance().getArrayExact(nFFT);

      double[] temp = DoubleArrayCache.getInstance().getArrayExact(input.length);
      System.arraycopy(input, 0, temp, 0, temp.length);

      for(int i = 0; i < nT; ++i) {
          int offset = i * step;
          int validLength = temp.length - offset;
          if (validLength >= nFFT) {
              System.arraycopy(temp, offset, raw, 0, nFFT);
              apodization.apodize(raw);
          } else {
              apodization.apodize(temp);
              System.arraycopy(temp, offset, raw, 0, validLength);

              label32:
              switch(padding) {
                  case MIRROR:
                      int j = validLength;

                      while(true) {
                          if (j >= raw.length) {
                              break label32;
                          }

                          raw[j] = temp[temp.length - j + validLength - 1];
                          ++j;
                      }
                  case ZERO:
                      Arrays.fill(raw, validLength, raw.length, 0.0D);
                      break;
                  case ZOH:
                  default:
                      Arrays.fill(raw, validLength, raw.length, temp[temp.length - 1]);
              }
          }

          fastFourierTransform.realForward(raw);
          double k = (double) nFFT / (double) validLength;

          if (dbScale) {
              SpectrumTools.computeMagnitudeSpectrum_dB(raw, 0, nFFT, currentMagnitudeData, 0, truncateDCNy);
              k = Units.conversionLogUnits(k, true);

              // Amplitude compensation
              if (padding == ShortTimeFourierTransform.Padding.ZERO && nFFT > validLength) {
                  for (int j = 0; j < currentMagnitudeData.length; ++j) {
                      currentMagnitudeData[j] += k;
                  }
              }
          } else {
              SpectrumTools.computeMagnitudeSpectrum(raw, 0, nFFT, currentMagnitudeData, 0, truncateDCNy);

              // Amplitude compensation
              if (padding == ShortTimeFourierTransform.Padding.ZERO && nFFT > validLength) {
                  for (int j = 0; j < currentMagnitudeData.length; ++j) {
                      currentMagnitudeData[j] *= k;
                  }
              }
          }

          System.arraycopy(currentMagnitudeData, 0, amplitudeData, i * nFFT / 2, nFFT / 2);
      }

      DoubleArrayCache.getInstance().add(currentMagnitudeData);
      DoubleArrayCache.getInstance().add(raw);

      return amplitudeData;
 }

public static double conversionLogUnits(final double value, final boolean amplitude) {
    if (value == 0.0) {
        return value;
    }

    double temp = Math.log10(value);

    if (amplitude) {
        return 20.0 * temp;
    } else {
        return 10.0 * temp;
    }
}