import get from 'lodash-es/get.js'
import size from 'lodash-es/size.js'
import each from 'lodash-es/each.js'
import range from 'lodash-es/range.js'
import ispnum from 'wsemi/src/ispnum.mjs'
import isp0num from 'wsemi/src/isp0num.mjs'
import cdbl from 'wsemi/src/cdbl.mjs'
import fft1d from './fft1d.mjs'
import ifft1d from './ifft1d.mjs'
/**
* FFT1D Filter
*
* @param {Array} arr 輸入數據陣列
* @param {Number} dt 輸入數據點時間間隔數字，單位s
* @param {Number} hzStart 輸入過濾用帶通頻率下限數字，單位Hz
* @param {Number} hzEnd 輸入過濾用帶通頻率上限數字，單位Hz
* @param {Object} [opt={}] 選項物件
* @param {Boolean} [opt.useOneTurn=true] true=輸入視為一個完整週期(頭尾重複)，週期(n-1)*dt；false=標準DFT，週期n*dt
* @param {String} [opt.type='dft'] 輸入計算方式字串，'dft'為使用mathjs對任意n點做真實n點DFT(2冪次走Cooley-Tukey、其餘走Chirp-Z)，數據品質最佳但非2冪次時較慢；'pow2'為先補零至2冪次(最少4點)再使用ml-fft之radix-2 FFT，速度極快適合前端即時繪圖，但輸出點數與頻率解析度df係以補零後之2冪次點數計算，預設'dft'
* @return {Array} 回傳帶通處理後數據陣列
* @example
*
* let dt
* let arr
* let res
*
* //dt=0.0078125s, 3+6hz
* dt = 0.0078125
*
* arr = [
* 0,
* 0.437015151709824,
* 0.845854910274064,
* 1.20056554679302,
* 1.47944976553089,
* 1.66674368151922,
* 1.75379573376597,
* 1.73964987434863,
* 1.63098631369783,
* 1.44142799002054,
* 1.19027504868833,
* 0.900778315875612,
* 0.598101848038141,
* 0.307150781019375,
* 0.0504516520458098,
* -0.153732804251563,
* -0.292893218813452,
* -0.361241031239776,
* -0.360072875476548,
* -0.297503430771426,
* -0.187593110348962,
* -0.0489494660021425,
* 0.0970731816865672,
* 0.228416556922734,
* 0.324423348821458,
* 0.367818520155133,
* 0.346391996239585,
* 0.254233601317238,
* 0.0924099202087415,
* -0.130978839760706,
* -0.401370102712605,
* -0.698891832710318,
* -1,
* -1.27946118721924,
* -1.51251056875181,
* -1.67699974648618,
* -1.75534914481383,
* -1.73613585202716,
* -1.61517856456688,
* -1.39602400854158,
* -1.08979021355164,
* -0.714376916729265,
* -0.293107462345689,
* 0.147084814656977,
* 0.577773754381215,
* 0.971283137555866,
* 1.30286634912854,
* 1.55263964022464,
* 1.70710678118655,
* 1.76014786721285,
* 1.7133908766509,
* 1.57593734934667,
* 1.36346871276832,
* 1.09681259653473,
* 0.80009440465607,
* 0.498634516368545,
* 0.216772751324739,
* -0.0241926543480825,
* -0.207774827040493,
* -0.323625771825179,
* -0.368309299491684,
* -0.345455359932455,
* -0.265285555765141,
* -0.1435542027991,
* -3.67544536472586E-16,
* 0.1435542027991,
* 0.26528555576514,
* 0.345455359932455,
* 0.368309299491684,
* 0.323625771825179,
* 0.207774827040493,
* 0.0241926543480835,
* -0.216772751324738,
* -0.498634516368544,
* -0.800094404656069,
* -1.09681259653473,
* -1.36346871276832,
* -1.57593734934667,
* -1.7133908766509,
* -1.76014786721285,
* -1.70710678118655,
* -1.55263964022464,
* -1.30286634912855,
* -0.971283137555864,
* -0.577773754381218,
* -0.14708481465698,
* 0.293107462345686,
* 0.714376916729258,
* 1.08979021355163,
* 1.39602400854158,
* 1.61517856456688,
* 1.73613585202716,
* 1.75534914481383,
* 1.67699974648618,
* 1.51251056875181,
* 1.27946118721924,
* 1,
* 0.698891832710321,
* 0.40137010271261,
* 0.130978839760704,
* -0.0924099202087405,
* -0.254233601317238,
* -0.346391996239584,
* -0.367818520155133,
* -0.324423348821457,
* -0.228416556922734,
* -0.0970731816865679,
* 0.0489494660021418,
* 0.18759311034896,
* 0.297503430771423,
* 0.360072875476548,
* 0.361241031239776,
* 0.292893218813452,
* 0.153732804251566,
* -0.0504516520458086,
* -0.307150781019376,
* -0.598101848038137,
* -0.900778315875611,
* -1.19027504868833,
* -1.44142799002054,
* -1.63098631369783,
* -1.73964987434863,
* -1.75379573376597,
* -1.66674368151922,
* -1.47944976553089,
* -1.20056554679302,
* -0.845854910274064,
* -0.43701515170983,
* ]
*
* res = wf.filter1d(arr, dt, 0, 2)
* console.log(res)
* // => [
* // -1.967983691641977e-17,
* // -2.307740440357446e-17,
* // -2.6914726362424797e-17,
* // -3.118255834826216e-17,
* // -3.587061878001477e-17,
* // -4.0967613709476845e-17,
* // -4.646126402941945e-17,
* // -5.2338335055034716e-17,
* // -5.858466840745145e-17,
* // -6.518521612250948e-17,
* // -7.212407690262343e-17,
* // -7.938453442440174e-17,
* // -8.694909760973329e-17,
* // -9.479954276332726e-17,
* // -1.029169574751909e-16,
* // -1.112817861822817e-16,
* // ... (112 more items)
* // ]
*
* res = wf.filter1d(arr, dt, 0, 4)
* console.log(res)
* // => [
* // -4.560690164444523e-16,
* // 0.14673047445536136,
* // 0.290284677254462,
* // 0.4275550934302818,
* // 0.5555702330196021,
* // 0.6715589548470182,
* // 0.7730104533627369,
* // 0.8577286100002721,
* // 0.9238795325112868,
* // 0.970031253194544,
* // 0.9951847266721972,
* // 0.9987954562051725,
* // 0.9807852804032308,
* // 0.9415440651830213,
* // 0.8819212643483555,
* // 0.8032075314806454,
* // ... (112 more items)
* // ]
*
* res = wf.filter1d(arr, dt, 4, 8)
* console.log(res)
* // => [
* // -4.482834651769285e-16,
* // 0.2902846772544623,
* // 0.5555702330196025,
* // 0.7730104533627374,
* // 0.9238795325112875,
* // 0.995184726672198,
* // 0.9807852804032318,
* // 0.8819212643483564,
* // 0.7071067811865488,
* // 0.47139673682599875,
* // 0.19509032201612916,
* // -0.0980171403295599,
* // -0.3826834323650895,
* // -0.6343932841636456,
* // -0.8314696123025458,
* // -0.9569403357322097,
* // ... (112 more items)
* // ]
*
*/
function filter1d(arr, dt, hzStart, hzEnd, opt = {}) {
//check dt
if (!ispnum(dt)) {
throw new Error(`dt[${dt}] is not a positive number`)
}
dt = cdbl(dt)
//check hzStart
if (!isp0num(hzStart)) {
throw new Error(`hzStart[${hzStart}] is not a positive number`)
}
hzStart = cdbl(hzStart)
//check hzEnd
if (!isp0num(hzEnd)) {
throw new Error(`hzEnd[${hzEnd}] is not a positive number`)
}
hzEnd = cdbl(hzEnd)
//useOneTurn: true(預設)=輸入視為一個完整週期(頭尾重複), 週期(n-1)*dt; false=標準DFT, 週期n*dt
let useOneTurn = get(opt, 'useOneTurn', true)
//type: 'dft'(預設)=mathjs真實n點DFT; 'pow2'=補零至2冪次後用ml-fft, 速度極快
let type = get(opt, 'type', 'dft')
//fft1d
let rm = fft1d(arr, { type })
// console.log('rm', rm)
// console.log('n1', size(arr))
//n
let n = size(rm)
// console.log('n', n)
//T
let T = useOneTurn ? dt * (n - 1) : dt * n
// console.log('T', T)
//df
let df = 1 / T
// console.log('df', df)
//F
// let F = df * (n - 1) //舊版hzs建構用, 新遮罩改用min(k,n-k)不需要
// console.log('F', F)
//hzs
//舊版: 對稱建構hzs; 對奇數n時 take(hzs, n/2) 取整少一個 -> mirror後長度n-1, 末尾bin永遠漏遮罩, 且負頻半段與正頻不共軛對齊 -> 輸出非實數(被ifft取實部截斷), 保留供比對
// let hzs = range(df, F + df, df)
// let hzs = range(0, F, df)
// let hzsHalf = take(hzs, n / 2)
// hzs = [
// ...hzsHalf,
// ...reverse(hzsHalf),
// ]
//新版: 每個bin k的頻率大小=min(k,n-k)*df(bin k 與其共軛 bin n-k 同頻率→一起遮罩, 輸出保持實數), 涵蓋全部n個bin, 奇偶n皆正確
let hzs = range(0, n).map((k) => Math.min(k, n - k) * df)
// console.log('hzs', JSON.stringify(hzs), size(hzs))
//帶通, 依照指定起訖頻率清除rm
each(hzs, (v, k) => {
let b = hzStart <= v && v <= hzEnd //允許通過
if (!b) {
rm[k][0] = 0
rm[k][1] = 0
}
})
// console.log('rm(clean)', JSON.stringify(rm))
//ifft1d, rm之點數已為2冪次, 'pow2'不會再補零
let res = ifft1d(rm, { type })
// console.log('res', res)
return res
}
export default filter1d