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 take from 'lodash-es/take.js'
// import reverse from 'lodash-es/reverse.js' //僅舊版hzs對稱建構使用, 已改為min(k,n-k)遮罩
import ispnum from 'wsemi/src/ispnum.mjs'
import isp0num from 'wsemi/src/isp0num.mjs'
import cdbl from 'wsemi/src/cdbl.mjs'
import { fft, ifft, complex } from 'mathjs'
//舊版: 取得>=n的最小2冪次, 供ml-fft補零至2冪次用; 改用mathjs後不需要, 保留供比對
/*
function get2n(n) {
let i = 1
let j = Math.pow(2, 52)
while (true) {
i *= 2
// console.log('n', n, 'i', i, 'n <= i', n <= i)
if (n <= i) {
break
}
// console.log('j', j, 'i >= j', i >= j)
if (i >= j) {
break
}
}
return i
}
*/
//舊版: 基於ml-fft, 先用get2n將n補零至2冪次(nCols)再做FFT; 非2冪次時實際為nCols點DFT(對DTFT的插值), 與真實n點DFT(MATLAB fft)不同, 保留供比對
/*
function _fft1d(arr, mode = 'norm') {
//check
if (mode !== 'norm' && mode !== 'inv') {
throw new Error(`invalid mode[${mode}]`)
}
//n
let n = size(arr)
// console.log('n', n)
//nCols
let nCols = get2n(n)
// console.log('nCols', nCols)
//init
FFT.init(nCols)
//fill
let re = new Array(nCols)
let im = new Array(nCols)
if (mode === 'norm') {
for (let i = 0; i < nCols; i++) {
let _i = get(arr, i, 0) //超過n至nCols之元素自動補0
let _j = 0
// console.log('_i', _i, '_j', _j)
re[i] = _i
im[i] = _j
}
}
else {
for (let i = 0; i < nCols; i++) {
let _i = get(arr, `${i}.0`, 0)
let _j = get(arr, `${i}.1`, 0)
// console.log('_i', _i, '_j', _j)
re[i] = _i
im[i] = _j
}
}
let res = []
if (mode === 'norm') {
FFT.fft(re, im)
for (let i = 0; i < nCols; i++) {
let _i = re[i]
let _j = im[i]
res.push([_i, _j])
}
}
else {
FFT.ifft(re, im)
for (let i = 0; i < nCols; i++) {
let _i = re[i]
// let _j = im[i]
// let _l = Math.sqrt(_i * _i, _j * _j)
res.push(_i)
}
}
return res
}
*/
//新版: 基於mathjs, 直接對任意n點做真實n點DFT(對齊MATLAB fft), 不再補零至2冪次
function _fft1d(arr, mode = 'norm') {
//check
if (mode !== 'norm' && mode !== 'inv') {
throw new Error(`invalid mode[${mode}]`)
}
//n
let n = size(arr)
// console.log('n', n)
//fft, mathjs對任意n: 2冪次走Cooley-Tukey, 其餘走Chirp-Z(CZT), 皆為真實n點DFT
let res = []
if (mode === 'norm') {
//fill, 實數輸入(虛部視為0)
let cs = new Array(n)
for (let i = 0; i < n; i++) {
cs[i] = get(arr, i, 0)
}
//fft, 回傳mathjs Complex物件陣列(含re, im)
let out = fft(cs)
//res
for (let i = 0; i < n; i++) {
res.push([out[i].re, out[i].im])
}
}
else {
//fill, 複數輸入[[re,im],...]
let cs = new Array(n)
for (let i = 0; i < n; i++) {
let _i = get(arr, `${i}.0`, 0)
let _j = get(arr, `${i}.1`, 0)
cs[i] = complex(_i, _j)
}
//ifft, mathjs之ifft已正規化(內部除以n), round-trip可還原原訊號
let out = ifft(cs)
//res, 只取實部(與舊版一致)
for (let i = 0; i < n; i++) {
res.push(out[i].re)
}
}
return res
}
/**
* FFT1D
*
* @param {Array} arr 輸入數據陣列
* @return {Array} 回傳轉換後數據陣列
* @example
*
* let arr
* let res
*
* arr = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15]
* res = wf.fft1d(arr)
* console.log(res)
* // => [
* // [ 120, 0 ],
* // [ -8, 40.21871593700678 ],
* // [ -8, 19.31370849898476 ],
* // [ -8.000000000000002, 11.972846101323912 ],
* // [ -8, 8 ],
* // [ -7.999999999999999, 5.345429103354391 ],
* // [ -8, 3.3137084989847594 ],
* // [ -8.000000000000002, 1.5912989390372623 ],
* // [ -8, 0 ],
* // [ -8, -1.5912989390372623 ],
* // [ -7.999999999999999, -3.3137084989847594 ],
* // [ -7.999999999999998, -5.345429103354393 ],
* // [ -7.999999999999999, -8 ],
* // [ -7.999999999999999, -11.97284610132391 ],
* // [ -7.999999999999997, -19.31370849898476 ],
* // [ -7.999999999999993, -40.21871593700678 ]
* // ]
*
*/
let fft1d = (arr) => {
return _fft1d(arr, 'norm')
}
/**
* iFFT1D
*
* @param {Array} arr 輸入數據陣列
* @return {Array} 回傳轉換後數據陣列
* @example
*
* let arr
* let res
*
* arr = [
* [120, 0],
* [-8, 40.21871593700678],
* [-8, 19.31370849898476],
* [-7.999999999999999, 11.972846101323913],
* [-8, 8],
* [-8, 5.345429103354391],
* [-8, 3.313708498984761],
* [-7.999999999999999, 1.5912989390372623],
* [-8, 0],
* [-7.999999999999999, -1.5912989390372623],
* [-8, -3.313708498984761],
* [-8, -5.345429103354391],
* [-8, -8],
* [-7.999999999999999, -11.972846101323913],
* [-8, -19.31370849898476],
* [-8, -40.21871593700678]
* ]
* res = wf.ifft1d(arr)
* console.log(res)
* // => [
* // 0, 1,
* // 2.0000000000000004, 3.0000000000000013,
* // 4, 5.000000000000002,
* // 6, 7.000000000000001,
* // 8, 9,
* // 10, 10.999999999999998,
* // 12, 12.999999999999998,
* // 14, 15
* // ]
*
*/
let ifft1d = (arr) => {
return _fft1d(arr, 'inv')
}
// function __fft2d(re, im) {
// let tre = []
// let tim = []
// let i = 0
// // x-axis
// for (let y = 0; y < _n; y++) {
// i = y * _n
// for (let x1 = 0; x1 < _n; x1++) {
// tre[x1] = re[x1 + i]
// tim[x1] = im[x1 + i]
// }
// fft1d(tre, tim)
// for (let x2 = 0; x2 < _n; x2++) {
// re[x2 + i] = tre[x2]
// im[x2 + i] = tim[x2]
// }
// }
// // y-axis
// for (let x = 0; x < _n; x++) {
// for (let y1 = 0; y1 < _n; y1++) {
// i = x + y1 * _n
// tre[y1] = re[i]
// tim[y1] = im[i]
// }
// fft1d(tre, tim)
// for (let y2 = 0; y2 < _n; y2++) {
// i = x + y2 * _n
// re[i] = tre[y2]
// im[i] = tim[y2]
// }
// }
// }
// function __ifft2d(re, im) {
// let tre = []
// let tim = []
// let i = 0
// // x-axis
// for (let y = 0; y < _n; y++) {
// i = y * _n
// for (let x1 = 0; x1 < _n; x1++) {
// tre[x1] = re[x1 + i]
// tim[x1] = im[x1 + i]
// }
// ifft1d(tre, tim)
// for (let x2 = 0; x2 < _n; x2++) {
// re[x2 + i] = tre[x2]
// im[x2 + i] = tim[x2]
// }
// }
// // y-axis
// for (let x = 0; x < _n; x++) {
// for (let y1 = 0; y1 < _n; y1++) {
// i = x + y1 * _n
// tre[y1] = re[i]
// tim[y1] = im[i]
// }
// ifft1d(tre, tim)
// for (let y2 = 0; y2 < _n; y2++) {
// i = x + y2 * _n
// re[i] = tre[y2]
// im[i] = tim[y2]
// }
// }
// }
// function fft2d(mat, mode = 'norm') {
// //check
// if (mode !== 'norm' && mode !== 'inv') {
// throw new Error(`invalid mode[${mode}]`)
// }
// //m, n
// let m = size(mat)
// let n = size(get(mat, 0, []))
// console.log('m', m)
// console.log('n', n)
// //nRows, nCols
// let nRows = get2n(m)
// let nCols = get2n(n)
// console.log('nRows', nRows)
// console.log('nCols', nCols)
// //data
// let data = new Array(nRows * nCols)
// for (let i = 0; i < nRows; i++) {
// for (let j = 0; j < nCols; j++) {
// data[i * nCols + j] = get(mat, `${i}.${j}`, 0) //超過元素自動補0
// }
// }
// console.log('data', data)
// let res = []
// if (mode === 'norm') {
// let ftData = FFTUtils.fft2DArray(data, nRows, nCols)
// console.log('ftData', ftData, size(ftData))
// }
// else {
// }
// // let ftData = FFTUtils.fft2DArray(data, nRows, nCols)
// // let ftRows = nRows * 2
// // let ftCols = nCols / 2 + 1
// // let iftData = FFTUtils.ifft2DArray(ftData, ftRows, ftCols)
// }
//基於mathjs, 對任意m×n矩陣做真實2D DFT(對齊MATLAB fft2); mathjs內部對每列每行各做1D(可分離), 各維皆支援任意長度(2冪次走Cooley-Tukey, 其餘走Chirp-Z)
function _fft2d(mat, mode = 'norm') {
//check
if (mode !== 'norm' && mode !== 'inv') {
throw new Error(`invalid mode[${mode}]`)
}
//m, n (列數, 行數); n以第0列長度為準, 缺項補0
let m = size(mat)
let n = size(get(mat, 0, []))
// console.log('m', m, 'n', n)
//fft 2D
let res = []
if (mode === 'norm') {
//fill, 實數輸入矩陣(虛部視為0)
let cs = []
for (let i = 0; i < m; i++) {
let row = new Array(n)
for (let j = 0; j < n; j++) {
row[j] = get(mat, `${i}.${j}`, 0)
}
cs.push(row)
}
//fft, 回傳mathjs Complex物件之巢狀陣列
let out = fft(cs)
//res, 每元素轉[re,im]
for (let i = 0; i < m; i++) {
let row = []
for (let j = 0; j < n; j++) {
row.push([out[i][j].re, out[i][j].im])
}
res.push(row)
}
}
else {
//fill, 複數輸入矩陣[[[re,im],...],...]
let cs = []
for (let i = 0; i < m; i++) {
let row = new Array(n)
for (let j = 0; j < n; j++) {
let _i = get(mat, `${i}.${j}.0`, 0)
let _j = get(mat, `${i}.${j}.1`, 0)
row[j] = complex(_i, _j)
}
cs.push(row)
}
//ifft, mathjs之ifft已正規化(內部除以m*n), round-trip可還原原矩陣
let out = ifft(cs)
//res, 只取實部(與1D一致)
for (let i = 0; i < m; i++) {
let row = []
for (let j = 0; j < n; j++) {
row.push(out[i][j].re)
}
res.push(row)
}
}
return res
}
/**
* FFT2D
*
* @param {Array} mat 輸入二維實數矩陣
* @return {Array} 回傳二維[re,im]複數矩陣
* @example
*
* let mat
* let res
*
* mat = [[1, 2, 3], [4, 5, 6]]
* res = wf.fft2d(mat)
* console.log(res)
* // => [
* // [ [ 21, 1.0445074572148558e-16 ], [ -2.9999999999999982, 1.7320508075688785 ], [ -3.0000000000000013, -1.7320508075688705 ] ],
* // [ [ -9, -2.326366143623307e-16 ], [ 0, -8.881784197001252e-16 ], [ 0, -2.7755575615628914e-15 ] ]
* // ]
*
*/
let fft2d = (mat) => {
return _fft2d(mat, 'norm')
}
/**
* iFFT2D
*
* @param {Array} mat 輸入二維[re,im]複數矩陣
* @return {Array} 回傳二維實數矩陣
* @example
*
* let mat
* let res
*
* mat = [
* [ [ 21, 1.0445074572148558e-16 ], [ -2.9999999999999982, 1.7320508075688785 ], [ -3.0000000000000013, -1.7320508075688705 ] ],
* [ [ -9, -2.326366143623307e-16 ], [ 0, -8.881784197001252e-16 ], [ 0, -2.7755575615628914e-15 ] ]
* ]
* res = wf.ifft2d(mat)
* console.log(res)
* // => [
* // [ 1.0000000000000002, 2.0000000000000004, 2.9999999999999987 ],
* // [ 4, 5.000000000000001, 5.999999999999999 ]
* // ]
*
*/
let ifft2d = (mat) => {
return _fft2d(mat, 'inv')
}
/**
* 振幅頻譜1D
*
* @param {Array} arr 輸入數據陣列
* @param {Number} dt 時間間隔，單位s
* @param {Object} [opt={}] 選項物件
* @param {Boolean} [opt.useOneTurn=true] true=輸入視為一個完整週期(頭尾重複)，週期(n-1)*dt；false=標準DFT，週期n*dt
* @return {Array} 回傳[{freq,real,imag,ampl},...]頻譜陣列
* @example
*
* let arr
* let res
*
* arr = [1, 2, 3, 4]
* res = wf.spectrum1d(arr, 1)
* console.log(res)
* // => [
* // { freq: 0, real: 10, imag: 0, ampl: 10 },
* // { freq: 0.3333333333333333, real: -2, imag: 2, ampl: 2.8284271247461903 }
* // ]
*
*/
function spectrum1d(arr, dt, opt = {}) {
//check dt
if (!ispnum(dt)) {
throw new Error(`dt[${dt}] is not a positive number`)
}
dt = cdbl(dt)
//useOneTurn: true(預設)=輸入視為一個完整週期(頭尾為同一相位、重複), 週期取(n-1)*dt; false=標準DFT, 週期取n*dt
let useOneTurn = get(opt, 'useOneTurn', true)
//fft1d
let rm = fft1d(arr)
// console.log('rm', rm)
// console.log('n1', size(arr))
//n
let n = size(rm)
// console.log('n', n)
//T, useOneTurn=true時週期(n-1)*dt, =false時標準DFT週期n*dt
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)
// console.log('F', F)
//hzs
// let hzs = range(df, F + df, df)
let hzs = range(0, F, df)
let hzsHalf = take(hzs, n / 2)
// hzs = [
// ...hzsHalf,
// ...reverse(hzsHalf),
// ]
hzs = hzsHalf
// console.log('hzs', JSON.stringify(hzs), size(hzs))
//rs
let rs = []
each(hzs, (v, k) => {
let freq = v
let real = rm[k][0]
let imag = rm[k][1]
let ampl = Math.sqrt(rm[k][0] ** 2 + rm[k][1] ** 2)
let r = {
freq,
real,
imag,
ampl,
}
rs.push(r)
})
return rs
}
/**
* 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
* @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)
//fft1d
let rm = fft1d(arr)
// 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
let res = ifft1d(rm)
// console.log('res', res)
return res
}
/**
* 振幅頻譜2D
*
* @param {Array} mat 輸入二維實數矩陣
* @param {Number} dt 時間間隔，單位s
* @param {Object} [opt={}] 選項物件
* @param {Boolean} [opt.useOneTurn=true] true=各軸視為一個完整週期(頭尾重複)，週期(len-1)*dt；false=標準DFT，週期len*dt
* @return {Array} 回傳二維[{freqRow,freqCol,real,imag,ampl},...]頻譜矩陣(僅正頻半象限)
* @example
*
* let mat
* let res
*
* mat = [[0,1,2,3],[1,2,3,4],[2,3,4,5],[3,4,5,6]]
* res = wf.spectrum2d(mat, 0.5)
* console.log(res)
* // => [
* // [
* // { freqRow: 0, freqCol: 0, real: 48, imag: 0, ampl: 48 },
* // { freqRow: 0, freqCol: 0.6666666666666666, real: -8, imag: 8, ampl: 11.313708498984761 }
* // ],
* // [
* // { freqRow: 0.6666666666666666, freqCol: 0, real: -8, imag: 8, ampl: 11.313708498984761 },
* // { freqRow: 0.6666666666666666, freqCol: 0.6666666666666666, real: 0, imag: 0, ampl: 0 }
* // ]
* // ]
*
*/
function spectrum2d(mat, dt, opt = {}) {
//check dt
if (!ispnum(dt)) {
throw new Error(`dt[${dt}] is not a positive number`)
}
dt = cdbl(dt)
//useOneTurn: true(預設)=各軸視為一個完整週期(頭尾重複), 週期(len-1)*dt; false=標準DFT, 週期len*dt
let useOneTurn = get(opt, 'useOneTurn', true)
//fft2d, 回傳m×n之[re,im]
let rm = fft2d(mat)
//m, n (列數, 行數)
let m = size(rm)
let n = size(get(rm, 0, []))
//df, 兩軸各自以該軸點數計算
let dfRow = 1 / ((useOneTurn ? m - 1 : m) * dt)
let dfCol = 1 / ((useOneTurn ? n - 1 : n) * dt)
//取正頻半象限 floor(m/2) × floor(n/2)(與1D取floor(n/2)一致)
let mHalf = Math.floor(m / 2)
let nHalf = Math.floor(n / 2)
//rs, 每格{freqRow, freqCol, real, imag, ampl}
let rs = []
for (let i = 0; i < mHalf; i++) {
let row = []
for (let j = 0; j < nHalf; j++) {
let real = rm[i][j][0]
let imag = rm[i][j][1]
let ampl = Math.sqrt(real ** 2 + imag ** 2)
row.push({
freqRow: i * dfRow,
freqCol: j * dfCol,
real,
imag,
ampl,
})
}
rs.push(row)
}
return rs
}
/**
* FFT2D 帶通濾波
*
* @param {Array} mat 輸入二維實數矩陣
* @param {Number} dt 時間間隔，單位s
* @param {Number} hzStart 帶通頻率下限，單位Hz（以徑向頻率 sqrt(fRow²+fCol²) 判斷）
* @param {Number} hzEnd 帶通頻率上限，單位Hz
* @param {Object} [opt={}] 選項物件
* @param {Boolean} [opt.useOneTurn=true] true=各軸視為一個完整週期(頭尾重複)，週期(len-1)*dt；false=標準DFT，週期len*dt
* @return {Array} 回傳帶通處理後二維實數矩陣
* @example
*
* let mat
* let res
*
* mat = [[0,1,2,3],[1,2,3,4],[2,3,4,5],[3,4,5,6]]
* res = wf.filter2d(mat, 0.5, 0.3, 0.8)
* console.log(res)
* // => [
* // [ -2, -2, 0, 0 ],
* // [ -2, -2, 0, 0 ],
* // [ 0, 0, 2, 2 ],
* // [ 0, 0, 2, 2 ]
* // ]
*
*/
function filter2d(mat, 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(預設)=各軸視為一個完整週期(頭尾重複), 週期(len-1)*dt; false=標準DFT, 週期len*dt
let useOneTurn = get(opt, 'useOneTurn', true)
//fft2d, 回傳m×n之[re,im]
let rm = fft2d(mat)
//m, n
let m = size(rm)
let n = size(get(rm, 0, []))
//df, 兩軸各自
let dfRow = 1 / ((useOneTurn ? m - 1 : m) * dt)
let dfCol = 1 / ((useOneTurn ? n - 1 : n) * dt)
//帶通, 以徑向頻率sqrt(fRow²+fCol²)判斷; bin(i,j)頻率大小取min(i,m-i)*dfRow與min(j,n-j)*dfCol(與共軛bin同值→一起遮罩, 輸出保持實數), 涵蓋全部bin
for (let i = 0; i < m; i++) {
for (let j = 0; j < n; j++) {
let fRow = Math.min(i, m - i) * dfRow
let fCol = Math.min(j, n - j) * dfCol
let fRad = Math.sqrt(fRow ** 2 + fCol ** 2)
let b = hzStart <= fRad && fRad <= hzEnd //允許通過
if (!b) {
rm[i][j][0] = 0
rm[i][j][1] = 0
}
}
}
//ifft2d
let res = ifft2d(rm)
return res
}
/**
* FFT與iFFT
*
* @returns {Array} 回傳數據陣列
* @example
*
* let arr
* let res
*
* arr = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15]
* res = wf.fft1d(arr)
* console.log(res)
* // => [
* // [ 120, 0 ],
* // [ -8, 40.21871593700678 ],
* // [ -8, 19.31370849898476 ],
* // [ -8.000000000000002, 11.972846101323912 ],
* // [ -8, 8 ],
* // [ -7.999999999999999, 5.345429103354391 ],
* // [ -8, 3.3137084989847594 ],
* // [ -8.000000000000002, 1.5912989390372623 ],
* // [ -8, 0 ],
* // [ -8, -1.5912989390372623 ],
* // [ -7.999999999999999, -3.3137084989847594 ],
* // [ -7.999999999999998, -5.345429103354393 ],
* // [ -7.999999999999999, -8 ],
* // [ -7.999999999999999, -11.97284610132391 ],
* // [ -7.999999999999997, -19.31370849898476 ],
* // [ -7.999999999999993, -40.21871593700678 ]
* // ]
*
* arr = [
* [120, 0],
* [-8, 40.21871593700678],
* [-8, 19.31370849898476],
* [-7.999999999999999, 11.972846101323913],
* [-8, 8],
* [-8, 5.345429103354391],
* [-8, 3.313708498984761],
* [-7.999999999999999, 1.5912989390372623],
* [-8, 0],
* [-7.999999999999999, -1.5912989390372623],
* [-8, -3.313708498984761],
* [-8, -5.345429103354391],
* [-8, -8],
* [-7.999999999999999, -11.972846101323913],
* [-8, -19.31370849898476],
* [-8, -40.21871593700678]
* ]
* res = wf.ifft1d(arr)
* console.log(res)
* // => [
* // 0, 1,
* // 2.0000000000000004, 3.0000000000000013,
* // 4, 5.000000000000002,
* // 6, 7.000000000000001,
* // 8, 9,
* // 10, 10.999999999999998,
* // 12, 12.999999999999998,
* // 14, 15
* // ]
*
*/
let WFft = {
fft1d,
ifft1d,
fft2d,
ifft2d,
spectrum1d,
spectrum2d,
filter1d,
filter2d,
}
export default WFft