FFT & Signal Processing

1/**2 * Example 42: FFT & Signal Processing3 *4 * New in v1.0.0: Fast Fourier Transform (fft, ifft, rfft, irfft, fft2, ifft2, fftn)5 * for spectral analysis, filtering, and signal processing.6 */78import { fft, fft2, ifft, irfft, rfft, tensor } from "deepbox/ndarray";910console.log("=".repeat(60));11console.log("Example 42: FFT & Signal Processing");12console.log("=".repeat(60));1314// ============================================================================15// Part 1: Basic FFT — Time Domain to Frequency Domain16// ============================================================================17console.log("\n📊 Part 1: Basic FFT");18console.log("-".repeat(60));1920// A simple signal: sum of two sine waves at 5 Hz and 12 Hz21const sampleRate = 64;22const duration = 1; // 1 second23const n = sampleRate * duration;24const signalData: number[] = [];2526for (let i = 0; i < n; i++) {27  const t = i / sampleRate;28  // 5 Hz component (amplitude 1.0) + 12 Hz component (amplitude 0.5)29  signalData.push(Math.sin(2 * Math.PI * 5 * t) + 0.5 * Math.sin(2 * Math.PI * 12 * t));30}3132const signal = tensor(signalData);33console.log(`Signal: ${n} samples at ${sampleRate} Hz (1 second)`);34console.log(`  Components: 5 Hz (amplitude 1.0) + 12 Hz (amplitude 0.5)`);35console.log(`  Signal shape: ${signal.shape}`);3637// Compute FFT38const spectrum = fft(signal);39console.log(40  `  FFT output — real shape: ${spectrum.real.shape}, imag shape: ${spectrum.imag.shape}`41);4243// Compute magnitude spectrum44const realData = spectrum.real.data as Float64Array;45const imagData = spectrum.imag.data as Float64Array;46const magnitudes: number[] = [];47for (let i = 0; i < n; i++) {48  const re = Number(realData[spectrum.real.offset + i]);49  const im = Number(imagData[spectrum.imag.offset + i]);50  magnitudes.push(Math.sqrt(re * re + im * im) / n);51}5253// Show dominant frequencies (first half of spectrum)54console.log("\n  Frequency spectrum (top peaks):");55const halfN = Math.floor(n / 2);56const peaks: { freq: number; mag: number }[] = [];57for (let i = 1; i < halfN; i++) {58  const freq = (i * sampleRate) / n;59  const mag = magnitudes[i]! * 2; // multiply by 2 for single-sided spectrum60  if (mag > 0.1) {61    peaks.push({ freq, mag });62  }63}64peaks.sort((a, b) => b.mag - a.mag);65for (const peak of peaks.slice(0, 5)) {66  console.log(`    ${peak.freq.toFixed(1)} Hz — magnitude: ${peak.mag.toFixed(4)}`);67}6869// ============================================================================70// Part 2: Inverse FFT — Frequency Domain back to Time Domain71// ============================================================================72console.log("\n🔄 Part 2: Inverse FFT (Reconstruction)");73console.log("-".repeat(60));7475// Reconstruct the signal from its FFT76const reconstructed = ifft(spectrum.real, spectrum.imag);77console.log("Inverse FFT reconstruction:");78console.log(`  Reconstructed real shape: ${reconstructed.real.shape}`);7980// Check reconstruction error81const reconData = reconstructed.real.data as Float64Array;82let maxError = 0;83for (let i = 0; i < n; i++) {84  const original = Number((signal.data as Float64Array)[signal.offset + i]);85  const recon = Number(reconData[reconstructed.real.offset + i]);86  maxError = Math.max(maxError, Math.abs(original - recon));87}88console.log(`  Max reconstruction error: ${maxError.toExponential(4)}`);89console.log("  (Should be near machine epsilon ≈ 1e-15)");9091// ============================================================================92// Part 3: Real FFT (rfft) — Optimized for Real Signals93// ============================================================================94console.log("\n⚡ Part 3: Real FFT (rfft)");95console.log("-".repeat(60));9697// rfft is optimized for real-valued signals — returns only positive frequencies98const rspec = rfft(signal);99console.log("rfft — optimized for real signals:");100console.log(`  Input shape:  ${signal.shape} (${n} samples)`);101console.log(`  Output real:  ${rspec.real.shape} (${Math.floor(n / 2) + 1} frequencies)`);102console.log(`  Output imag:  ${rspec.imag.shape}`);103console.log("  Only positive frequencies (Nyquist symmetry exploited)");104105// Inverse real FFT106const irecon = irfft(rspec.real, rspec.imag);107console.log(`\n  irfft reconstruction shape: ${irecon.real.shape}`);108109// ============================================================================110// Part 4: 2D FFT — Image/Matrix Frequency Analysis111// ============================================================================112console.log("\n🖼️  Part 4: 2D FFT");113console.log("-".repeat(60));114115// Create a simple 2D pattern (checkerboard-like)116const size = 8;117const pattern2d: number[][] = [];118for (let i = 0; i < size; i++) {119  const row: number[] = [];120  for (let j = 0; j < size; j++) {121    row.push(Math.sin((2 * Math.PI * i) / size) + Math.cos((2 * Math.PI * 2 * j) / size));122  }123  pattern2d.push(row);124}125126const matrix = tensor(pattern2d);127console.log(`2D signal shape: ${matrix.shape}`);128129const spec2d = fft2(matrix);130console.log(`2D FFT output — real: ${spec2d.real.shape}, imag: ${spec2d.imag.shape}`);131132// Compute 2D magnitude133const real2d = spec2d.real.data as Float64Array;134const imag2d = spec2d.imag.data as Float64Array;135let maxMag2d = 0;136let maxI = 0;137let maxJ = 0;138for (let i = 0; i < size; i++) {139  for (let j = 0; j < size; j++) {140    const idx = spec2d.real.offset + i * size + j;141    const re = Number(real2d[idx]);142    const im = Number(imag2d[idx]);143    const mag = Math.sqrt(re * re + im * im);144    if (i + j > 0 && mag > maxMag2d) {145      maxMag2d = mag;146      maxI = i;147      maxJ = j;148    }149  }150}151console.log(`  Dominant 2D frequency: (${maxI}, ${maxJ}) with magnitude ${maxMag2d.toFixed(4)}`);152153// ============================================================================154// Part 5: Spectral Filtering — Removing High-Frequency Noise155// ============================================================================156console.log("\n🔇 Part 5: Spectral Filtering");157console.log("-".repeat(60));158159// Create a noisy signal160const cleanSignalData: number[] = [];161const noisySignalData: number[] = [];162for (let i = 0; i < n; i++) {163  const t = i / sampleRate;164  const clean = Math.sin(2 * Math.PI * 3 * t); // 3 Hz pure signal165  const noise = 0.3 * Math.sin(2 * Math.PI * 25 * t); // 25 Hz noise166  cleanSignalData.push(clean);167  noisySignalData.push(clean + noise);168}169170const noisySignal = tensor(noisySignalData);171console.log("Noisy signal: 3 Hz clean + 25 Hz noise");172173// FFT the noisy signal174const noisySpectrum = fft(noisySignal);175const noisyReal = new Float64Array(noisySpectrum.real.data as Float64Array);176const noisyImag = new Float64Array(noisySpectrum.imag.data as Float64Array);177178// Low-pass filter: zero out frequencies above 10 Hz179const cutoffBin = Math.floor((10 * n) / sampleRate);180for (let i = cutoffBin; i < n - cutoffBin; i++) {181  noisyReal[noisySpectrum.real.offset + i] = 0;182  noisyImag[noisySpectrum.imag.offset + i] = 0;183}184185console.log(`  Low-pass filter: cutoff at 10 Hz (bin ${cutoffBin})`);186console.log(`  Zeroed ${n - 2 * cutoffBin} frequency bins`);187188// Reconstruct filtered signal189const filteredReal = tensor(Array.from(noisyReal));190const filteredImag = tensor(Array.from(noisyImag));191const filtered = ifft(filteredReal, filteredImag);192193// Measure filtering quality194const filteredData = filtered.real.data as Float64Array;195let filterError = 0;196for (let i = 0; i < n; i++) {197  const diff = cleanSignalData[i]! - Number(filteredData[filtered.real.offset + i]);198  filterError += diff * diff;199}200const rmse = Math.sqrt(filterError / n);201console.log(`  RMSE after filtering: ${rmse.toFixed(6)} (lower is better)`);202203// ============================================================================204// Part 6: Parseval's Theorem — Energy Conservation205// ============================================================================206console.log("\n⚖️  Part 6: Parseval's Theorem");207console.log("-".repeat(60));208209// Parseval's theorem: energy in time domain equals energy in frequency domain210let timeEnergy = 0;211const sigData = signal.data as Float64Array;212for (let i = 0; i < n; i++) {213  const val = Number(sigData[signal.offset + i]);214  timeEnergy += val * val;215}216217let freqEnergy = 0;218for (let i = 0; i < n; i++) {219  const re = Number(realData[spectrum.real.offset + i]);220  const im = Number(imagData[spectrum.imag.offset + i]);221  freqEnergy += re * re + im * im;222}223freqEnergy /= n; // normalization224225console.log("Parseval's theorem: Σ|x[n]|² = (1/N) Σ|X[k]|²");226console.log(`  Time-domain energy:  ${timeEnergy.toFixed(6)}`);227console.log(`  Freq-domain energy:  ${freqEnergy.toFixed(6)}`);228console.log(`  Difference:          ${Math.abs(timeEnergy - freqEnergy).toExponential(4)}`);229230// ============================================================================231// Summary232// ============================================================================233console.log("\n💡 Key Takeaways");234console.log("-".repeat(60));235console.log("• fft/ifft: convert between time and frequency domains");236console.log("• rfft/irfft: optimized for real-valued signals (half the output)");237console.log("• fft2/ifft2: 2D FFT for images and matrices");238console.log("• Spectral filtering: modify frequency components then reconstruct");239console.log("• Parseval's theorem: energy is conserved between domains");240console.log("• FFT complexity: O(N log N) vs O(N²) for naive DFT");241242console.log("\n✅ FFT & Signal Processing Example Complete!");243console.log("=".repeat(60));244