import type { FnN3 , NumericArray } from "@thi.ng/api" ; import type { ComplexArray } from "./api.js" ; import { isComplex } from "./complex.js" ; import { magDb } from "./convert.js" ; import { invPowerScale , powerScale } from "./power.js" ; import { applyWindow } from "./window.js" ; param n export const complexArray = ( n ) : ComplexArray => [ new Float64Array ( n ) , new Float64Array ( n ) , ] ; /**
 * Creates a deep copy of given  ComplexArray .
 *
 *  param complex export const copyComplex = ( complex ) : ComplexArray => [ complex [ 0 ] . slice ( ) , complex [ 1 ] . slice ( ) , ] ; /**
 * If given a  ComplexArray , computes the complex conjugate,
 * concatenates it in mirrored order to input (excluding bin 0) and
 * returns it as new (complex) array.
 *
 *  remarks The length of the input buffer(s) is assumed to be a power of 2.
 *
 * 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 example param src param isImg export function conjugate ( src , isImg ? ) ; export function conjugate ( complex ) ; export function conjugate ( src , isImg = true ) { if ( isComplex ( src ) ) { const n = src [ 0 ] . length ; const res = complexArray ( n * 2 ) ; const [ sreal , simg ] = < ComplexArray > src ; const [ dreal , dimg ] = res ; ( < Float64Array > dreal ) . set ( sreal ) ; ( < Float64Array > 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 ; } else { const n = src . length ; const dest = new Float64Array ( n * 2 ) ; dest . set ( < NumericArray > src ) ; for ( let i = 1 , j = n * 2 - 1 ; i < n ; i ++ , j -- ) { dest [ j ] = isImg ? - ( < NumericArray > src ) [ i ] : ( < NumericArray > src ) [ i ] ; } return dest ; } } const swapR = ( real , n ) => { const n2 = n >> 1 ; let ii ; let jj ; let k ; let t ; for ( let i = 1 , j = 1 ; i < n ; i ++ ) { if ( i < j ) { ii = i - 1 ; jj = j - 1 ; t = real [ jj ] ; real [ jj ] = real [ ii ] ; real [ ii ] = t ; } k = n2 ; while ( k < j ) { j -= k ; k >>= 1 ; } j += k ; } } ; const swapRI = ( real , img , n ) => { const n2 = n >> 1 ; let ii ; let jj ; let k ; let t ; for ( let i = 1 , j = 1 ; i < n ; i ++ ) { if ( i < j ) { ii = i - 1 ; jj = j - 1 ; t = real [ jj ] ; real [ jj ] = real [ ii ] ; real [ ii ] = t ; t = img [ jj ] ; img [ jj ] = img [ ii ] ; img [ ii ] = t ; } k = n2 ; while ( k < j ) { j -= k ; k >>= 1 ; } j += k ; } } ; const transform = ( real , img , n ) => { let step = 1 ; let prevStep ; let i , j , ii , ip ; let tr , ti ; let ur , ui ; let wr , wi ; let t ; for ( let b = Math . ceil ( Math . log2 ( n ) ) ; b -- > 0 ; ) { prevStep = step ; step <<= 1 ; ur = 1 ; ui = 0 ; t = Math . PI / prevStep ; wr = Math . cos ( t ) ; wi = - Math . sin ( t ) ; for ( j = 1 ; j <= prevStep ; j ++ ) { for ( i = j ; i <= n ; i += step ) { ip = i + prevStep - 1 ; ii = i - 1 ; tr = real [ ip ] * ur - img [ ip ] * ui ; ti = real [ ip ] * ui + img [ ip ] * ur ; real [ ip ] = real [ ii ] - tr ; img [ ip ] = img [ ii ] - ti ; real [ ii ] += tr ; img [ ii ] += ti ; } t = ur ; ur = t * wr - ui * wi ; ui = t * wi + ui * wr ; } } } ; remarks If `window` is given, the `real` array will be pre-multiplied with
 * the `window`.
 *
 * No result scaling / normalization is performed. Use
 *  normalizeFFT for that purpose.
 *
 * - https://www.earlevel.com/main/2002/08/31/a-gentle-introduction-to-the-fft/
 * - https://www.earlevel.com/main/2019/04/30/waveutils-updated/
 * - https://betterexplained.com/articles/an-interactive-guide-to-the-fourier-transform/
 * - http://toxi.co.uk/p5/fftDebug/fft4amit.pde
 *
 *  param complex param window export const fft = ( complex , window ? ) : ComplexArray => { let real , img ; if ( isComplex ( complex ) ) { real = complex [ 0 ] ; img = < NumericArray > complex [ 1 ] ; } else { real = < NumericArray > complex ; } if ( window ) { applyWindow ( real , window ) ; } const n = real . length ; if ( img ) { swapRI ( real , img , n ) ; } else { swapR ( real , n ) ; img = new Float64Array ( n ) ; } transform ( real , img , n ) ; return [ real , img ] ; } ; /**
 * Inverse FFT via computing forward transform with swapped real/imaginary
 * components. Expects denormalized inputs (i.e. the same as the result of
 *  fft ).
 *
 *  remarks param complex export const ifft = ( src ) : ComplexArray => { let complex = isComplex ( src ) ? src : [ new Float64Array ( src . length ) , < NumericArray > src ] ; fft ( [ complex [ 1 ] , complex [ 0 ] ] ) ; return scaleFFT ( complex , 1 / complex [ 0 ] . length ) ; } ; export const scaleFFT = ( complex , scale ) : ComplexArray => { const [ real , img ] = complex ; const n = real . length ; for ( let i = 0 ; i < n ; i ++ ) { real [ i ] *= scale ; img [ i ] *= scale ; } return [ real , img ] ; } ; remarks param complex param window export const normalizeFFT = ( complex , window = 2 / complex [ 0 ] . length ) : ComplexArray => scaleFFT ( complex , powerScale ( window , 2 ) ) ; /**
 * Inverse operation of  normalizeFFT . De-normalizes the complex FFT
 * array by scaling each complex bin value with given scale factor (or, if given
 * as array, the scale factor derived from these window function samples).
 *
 *  remarks param complex param window export const denormalizeFFT = ( complex , window = complex [ 0 ] . length / 2 ) : ComplexArray => scaleFFT ( complex , invPowerScale ( window , 2 ) ) ; remarks It's recommended to apply this function prior computing
 *  spectrumPhase . The `eps` value might have to be adjusted and should
 * be approx. `max(spectrumMag(fft))/10000`.
 *
 * References:
 * - https://www.gaussianwaves.com/2015/11/interpreting-fft-results-obtaining-magnitude-and-phase-information/
 *
 *  param complex param eps export const thresholdFFT = ( complex , eps = 1e-12 ) => { const [ real , img ] = complex ; for ( let i = 0 , n = real . length ; i < n ; i ++ ) { if ( Math . hypot ( real [ i ] , img [ i ] ) < eps ) { real [ i ] = img [ i ] = 0 ; } } return complex ; } ; param complex param n param out export const spectrumMag = ( complex , n = complex [ 0 ] . length / 2 , out = [ ] ) => { const [ real , img ] = complex ; for ( let i = 0 ; i < n ; i ++ ) { out [ i ] = Math . hypot ( real [ i ] , img [ i ] ) ; } return out ; } ; remarks If `window` is given (scale factor or array), it will be used as (if number)
 * or to compute the scaling factor (if array) for each FFT bin's value. The
 * default (`window=1`) is the equivalent to a rectangular window (i.e. a
 * no-op). If windowing was used to compute the FFT, the same should be provided
 * to this function for correct results.
 *
 * **IMPORTANT:** If the FFT result has already been normalized using
 *  normalizeFFT , the scaling factor (`window` arg) MUST be set 1.0.
 *
 * By default only the first N/2 values are returned. If `db` is true, the
 * spectrum values are converted to dBFS.
 *
 * - https://holometer.fnal.gov/GH_FFT.pdf
 * - https://dsp.stackexchange.com/a/32080
 * - https://dsp.stackexchange.com/a/14935
 * - https://www.kvraudio.com/forum/viewtopic.php?t=276092
 *
 *  param complex param db param window param n param out export const spectrumPow = ( complex , db = false , window = 2 / complex [ 0 ] . length , n = complex [ 0 ] . length / 2 , out = [ ] ) => { const [ real , img ] = complex ; const scale = powerScale ( window , 2 ) ; for ( let i = 0 ; i < n ; i ++ ) { const p = real [ i ] ** 2 + img [ i ] ** 2 ; out [ i ] = db ? magDb ( Math . sqrt ( p ) * scale ) : p * scale ; } return out ; } ; remarks Consider applying  thresholdFFT prior to computing the phase spectrum
 * to avoid exploding floating point error magnitudes.
 *
 *  param complex param n param out export const spectrumPhase = ( complex , n = complex [ 0 ] . length / 2 , out = [ ] ) => { const [ real , img ] = complex ; for ( let i = 0 ; i < n ; i ++ ) { out [ i ] = Math . atan2 ( img [ i ] , real [ i ] ) ; } return out ; } ; /**
 * Returns FFT bin index for given frequency, sample rate and window
 * size. See  binFreq for reverse op.
 *
 *  param f param fs param n export const freqBin = ( f , fs , n ) => ( ( f * n ) / fs ) | 0 ; /**
 * Returns frequency for given FFT bin index, sample rate and window
 * size. See  freqBin for reverse op.
 *
 *  param bin param fs param n export const binFreq = ( bin , fs , n ) => ( bin * fs ) / n ; param n param fs param m export const fftFreq = ( n , fs , m = n / 2 ) => { const res = new Float64Array ( m ) ; for ( let i = 0 ; i <= m ; i ++ ) { res [ i ] = binFreq ( i , fs , n ) ; } return res ; } ;