93 lines
2.2 KiB
JavaScript
93 lines
2.2 KiB
JavaScript
// Helper: Next power of 2
|
|
function nextPow2(n) {
|
|
return 1 << (32 - Math.clz32(n - 1));
|
|
}
|
|
|
|
// Basic Cooley-Tukey FFT (real/imag arrays)
|
|
function fft1D(re, im) {
|
|
const N = re.length;
|
|
if (N <= 1) return;
|
|
|
|
// Bit-reversal permutation
|
|
const rev = new Uint32Array(N);
|
|
let logN = Math.log2(N);
|
|
for (let i = 0; i < N; i++) {
|
|
let x = i,
|
|
y = 0;
|
|
for (let j = 0; j < logN; j++) {
|
|
y <<= 1;
|
|
y |= x & 1;
|
|
x >>= 1;
|
|
}
|
|
rev[i] = y;
|
|
}
|
|
for (let i = 0; i < N; i++) {
|
|
if (i < rev[i]) {
|
|
[re[i], re[rev[i]]] = [re[rev[i]], re[i]];
|
|
[im[i], im[rev[i]]] = [im[rev[i]], im[i]];
|
|
}
|
|
}
|
|
|
|
// FFT
|
|
for (let s = 1; s <= logN; s++) {
|
|
const m = 1 << s;
|
|
const m2 = m >> 1;
|
|
const wAngle = (-2 * Math.PI) / m;
|
|
const cosW = Math.cos(wAngle);
|
|
const sinW = Math.sin(wAngle);
|
|
for (let k = 0; k < N; k += m) {
|
|
let wr = 1,
|
|
wi = 0;
|
|
for (let j = 0; j < m2; j++) {
|
|
const tRe = wr * re[k + j + m2] - wi * im[k + j + m2];
|
|
const tIm = wr * im[k + j + m2] + wi * re[k + j + m2];
|
|
const uRe = re[k + j];
|
|
const uIm = im[k + j];
|
|
re[k + j] = uRe + tRe;
|
|
im[k + j] = uIm + tIm;
|
|
re[k + j + m2] = uRe - tRe;
|
|
im[k + j + m2] = uIm - tIm;
|
|
const tempWr = wr;
|
|
wr = wr * cosW - wi * sinW;
|
|
wi = tempWr * sinW + wi * cosW;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// Inverse FFT
|
|
function ifft1D(re, im) {
|
|
// Conjugate
|
|
for (let i = 0; i < re.length; i++) im[i] = -im[i];
|
|
fft1D(re, im);
|
|
// Normalize and re-conjugate
|
|
const N = re.length;
|
|
for (let i = 0; i < N; i++) {
|
|
re[i] /= N;
|
|
im[i] = -im[i] / N;
|
|
}
|
|
}
|
|
|
|
// Apply ramp filter in frequency domain
|
|
export function applyRampFilter(projection) {
|
|
const N = nextPow2(projection.length);
|
|
const re = new Float32Array(N);
|
|
const im = new Float32Array(N);
|
|
for (let i = 0; i < projection.length; i++) {
|
|
re[i] = projection[i];
|
|
}
|
|
|
|
fft1D(re, im);
|
|
|
|
for (let i = 0; i < N / 2; i++) {
|
|
const freq = i / N;
|
|
re[i] *= freq;
|
|
im[i] *= freq;
|
|
re[N - i - 1] *= freq;
|
|
im[N - i - 1] *= freq;
|
|
}
|
|
|
|
ifft1D(re, im);
|
|
return Array.from(re.slice(0, projection.length));
|
|
}
|