diff --git a/packages/dsp/src/fft.ts b/packages/dsp/src/fft.ts index acf2ff3cab..9abc7af30d 100644 --- a/packages/dsp/src/fft.ts +++ b/packages/dsp/src/fft.ts @@ -4,6 +4,78 @@ import { ComplexArray } from "./api"; const PI = Math.PI; +export const complexArray = (n: number): ComplexArray => [ + new Float64Array(n), + new Float64Array(n) +]; + +/** + * If given a {@link ComplexArray}, computes the complex conjugate, + * concatenates it in mirrored order to input (excluding bin 0) and + * returns it a new (complex) array. + * + * If given a {@link @thi.ng/api#NumericArray}, the `isImg` arg chooses + * from one of the following operations: If `true` (default), returns + * new array with negated values concatenated in reverse order. If + * `false`, returns new array with original values concatenated in + * reverse order. + * + * @example + * ```ts + * conjugate([0, 3, 2, 1], true) + * // Float64Array [ 0, 3, 2, 1, -1, -2, -3, -0 ] + * + * conjugate([0, 3, 2, 1], false) + * // Float64Array [ 0, 3, 2, 1, 1, 2, 3, 0 ] + * + * conjugate([[0, 1, 0, 1], [0, -0.5, 0, -0.25]]) + * [ + * Float64Array [ 0, 1, 0, 1, 1, 0, 1, 0 ], + * Float64Array [ 0, -0.5, 0, -0.25, 0.25, -0, 0.5, -0 ] + * ] + * ``` + * + * @example + * ```ts + * // generate 1Hz sine + * ifft(conjugate([0, 8, 0, 0, 0, 0, 0, 0]))[0] + * // [ + * // 0, 0.383, 0.707, 0.924, + * // 1, 0.924, 0.707, 0.383, + * // 0, -0.384, -0.707, -0.924, + * // -1, -0.924, -0.707, -0.383 + * // ] + * ``` + * + * @param src + * @param isImg + */ +export function conjugate(src: NumericArray, isImg?: boolean): NumericArray; +export function conjugate(complex: ComplexArray): ComplexArray; +export function conjugate(src: NumericArray | ComplexArray, isImg = true): any { + if (isNumber(src[0])) { + const n = src.length; + const dest = new Float64Array(n * 2); + dest.set(src); + for (let i = 1, j = n * 2 - 1; i < n; i++, j--) { + dest[j] = isImg ? -(src)[i] : (src)[i]; + } + return dest; + } else { + const n = src[0].length; + const res = complexArray(n * 2); + const [sreal, simg] = src; + const [dreal, dimg] = res; + (dreal).set(sreal); + (dimg).set(simg); + for (let i = 1, j = n * 2 - 1; i < n; i++, j--) { + dreal[j] = sreal[i]; + dimg[j] = -simg[i]; + } + return res; + } +} + const swapR = (real: NumericArray, n: number) => { const n2 = n >> 1; let ii: number; @@ -103,8 +175,7 @@ const transform = (real: NumericArray, img: NumericArray, n: number) => { * - https://betterexplained.com/articles/an-interactive-guide-to-the-fourier-transform/ * - http://toxi.co.uk/p5/fftDebug/fft4amit.pde * - * @param real - * @param img + * @param complex * @param window */ export const fft = ( @@ -145,10 +216,12 @@ export const fft = ( * * - https://www.dsprelated.com/showarticle/800.php (method #3) * - * @param real - * @param img + * @param complex */ -export const ifft = (complex: ComplexArray): ComplexArray => { +export const ifft = (src: NumericArray | ComplexArray): ComplexArray => { + let complex: ComplexArray = isNumber(src[0]) + ? [new Float64Array(src.length), src] + : src; fft([complex[1], complex[0]]); return scaleFFT(complex, 1 / complex[0].length); }; @@ -172,6 +245,14 @@ export const normalizeFFT = (complex: ComplexArray): ComplexArray => export const denormalizeFFT = (complex: ComplexArray): ComplexArray => scaleFFT(complex, Math.sqrt(complex[0].length)); +/** + * Computes magnitude spectrum for given FFT. By default only the first + * N/2 values are returned. + * + * @param complex - FFT result + * @param n - bin count + * @param out - result array + */ export const spectrumMag = ( complex: ComplexArray, n = complex[0].length / 2, @@ -185,38 +266,46 @@ export const spectrumMag = ( }; /** - * Computes power spectrum for the given FFT result arrays (length = N) - * and writes result to `out`. + * Computes power spectrum (optionally as dBFS) for the given raw, + * unnormalized FFT result arrays (length = N) and writes result to + * `out`. * * @remarks - * By default only the first N/2 values are returned. + * By default only the first N/2 values are returned. If `db` is true, + * the spectrum values are converted to dBFS. * - * Also by default, the FFT inputs are expected to be normalized (e.g. - * via {@link normalizeFFT}) and spectrum values converted to dB. + * - https://www.kvraudio.com/forum/viewtopic.php?t=276092 * - * @param real - * @param img - * @param normalized + * @param complex * @param db * @param n * @param out */ export const spectrumPow = ( complex: ComplexArray, - normalized = true, - db = true, + db = false, n = complex[0].length / 2, out: NumericArray = [] ) => { const [real, img] = complex; - const scale = normalized ? 1 : 1 / real.length; + const scale = (db ? 2 : 1) / real.length; for (let i = 0; i < n; i++) { - const p = (real[i] ** 2 + img[i] ** 2) * scale; - out[i] = db ? 20 * (Math.log(p) / Math.LN10) : p; + const p = real[i] ** 2 + img[i] ** 2; + out[i] = db + ? 20 * (Math.log(Math.sqrt(p) * scale) / Math.LN10) + : p * scale; } return out; }; +/** + * Computes phase spectrum for given FFT and writes results to `out`. By + * default only the first N/2 values are returned. + * + * @param complex - FFT result + * @param n - bin count + * @param out - result array + */ export const spectrumPhase = ( complex: ComplexArray, n = complex[0].length / 2, @@ -228,3 +317,39 @@ export const spectrumPhase = ( } return out; }; + +/** + * Returns FFT bin index for given frequency, sample rate and window + * size. See {@link binFreq} for reverse op. + * + * @param f - frequency + * @param fs - sample rate + * @param n - window size + */ +export const freqBin = (f: number, fs: number, n: number) => (f * n) / fs; + +/** + * Returns frequency for given FFT bin index, sample rate and window + * size. See {@link freqBin} for reverse op. + * + * @param bin - bin + * @param fs - sample rate + * @param n - window size + */ +export const binFreq = (bin: number, fs: number, n: number) => (bin * fs) / n; + +/** + * Returns array of bin center frequencies for given FFT window size and + * sample rate. By default only the first N/2 values are returned. + * + * @param n + * @param fs + * @param m + */ +export const fftFreq = (n: number, fs: number, m = n / 2) => { + const res = new Float64Array(m); + for (let i = 0; i < m; i++) { + res[i] = binFreq(i, fs, n); + } + return res; +}; diff --git a/packages/dsp/test/fft.ts b/packages/dsp/test/fft.ts index ec22e6fd20..980a8d9f1d 100644 --- a/packages/dsp/test/fft.ts +++ b/packages/dsp/test/fft.ts @@ -58,27 +58,17 @@ describe("fft", () => { ]) ); - assert.ok( - deltaEq(spectrumPow(res, false, false), [0, 3.414, 0, 0.586]) - ); + assert.ok(deltaEq(spectrumPow(res, false), [0, 3.414, 0, 0.586])); assert.ok( - deltaEq(spectrumPow(res, false, true), [ + deltaEq(spectrumPow(res, true), [ -Infinity, - 10.666, + 2.323, -Infinity, - -4.645 + -5.333 ]) ); - assert.ok( - deltaEq(spectrumPow(norm, true, false), [0, 3.414, 0, 0.586]) - ); - - assert.ok( - deltaEq(spectrumPow(norm), [-Infinity, 10.666, -Infinity, -4.645]) - ); - assert.ok(deltaEq(spectrumPhase(res), [0, 1.963, 0, 2.749])); }); });