diff --git a/benchmark/reimMatrixFFTByColumnsComparison.ts b/benchmark/reimMatrixFFTByColumnsComparison.ts
new file mode 100644
index 00000000..dc4dcddb
--- /dev/null
+++ b/benchmark/reimMatrixFFTByColumnsComparison.ts
@@ -0,0 +1,108 @@
+/* eslint-disable no-console */
+import { reimMatrixFFT } from '../src/reimMatrix/reimMatrixFFT.ts';
+import { reimMatrixFFTByColumns } from '../src/reimMatrix/reimMatrixFFTByColumns.ts';
+import type { DataReImMatrix } from '../src/types/index.ts';
+
+/**
+ * Transpose a complex matrix (with separate re and im arrays)
+ * @param data
+ */
+function transposeComplexMatrix(data: DataReImMatrix): DataReImMatrix {
+  const { re, im } = data;
+  const numRows = re.length;
+  if (numRows === 0) return { re: [], im: [] };
+
+  const numCols = re[0].length;
+  const transRe: Float64Array[] = [];
+  const transIm: Float64Array[] = [];
+
+  for (let col = 0; col < numCols; col++) {
+    const newRow = new Float64Array(numRows);
+    const newIm = new Float64Array(numRows);
+    for (let row = 0; row < numRows; row++) {
+      newRow[row] = re[row][col];
+      newIm[row] = im[row][col];
+    }
+    transRe.push(newRow);
+    transIm.push(newIm);
+  }
+
+  return { re: transRe, im: transIm };
+}
+
+// Parameters
+const numRows = 512; // number of columns to transform
+const numCols = 4096; // length of each column (size of FFT)
+const iterations = 20;
+
+console.log(
+  `Matrix size: ${numRows} rows × ${numCols} columns (${numRows * numCols} complex numbers)`,
+);
+console.log(`Iterations: ${iterations}`);
+console.log('');
+
+// Build input data
+const inputData: DataReImMatrix = {
+  re: Array.from({ length: numRows }, () => {
+    const arr = new Float64Array(numCols);
+    for (let i = 0; i < numCols; i++) {
+      arr[i] = Math.random();
+    }
+    return arr;
+  }),
+  im: Array.from({ length: numRows }, () => {
+    const arr = new Float64Array(numCols);
+    for (let i = 0; i < numCols; i++) {
+      arr[i] = Math.random();
+    }
+    return arr;
+  }),
+};
+
+// Warmup runs
+reimMatrixFFTByColumns(inputData);
+const transposed = transposeComplexMatrix(inputData);
+reimMatrixFFT(transposed);
+
+console.log(
+  '--- Approach 1: reimMatrixFFTByColumns (Transform columns directly) ---',
+);
+{
+  let totalTime = 0;
+  for (let i = 0; i < iterations; i++) {
+    const start = performance.now();
+    reimMatrixFFTByColumns(inputData, { inPlace: true });
+    const elapsed = performance.now() - start;
+    totalTime += elapsed;
+  }
+  const avgTime = totalTime / iterations;
+  console.log(`Total time: ${totalTime.toFixed(2)} ms`);
+  console.log(`Average time per iteration: ${avgTime.toFixed(4)} ms`);
+  console.log(`Time per column: ${(avgTime / numRows).toFixed(6)} ms`);
+}
+
+console.log('');
+console.log(
+  '--- Approach 2: Transpose → reimMatrixFFT → Transpose (Transform rows via transposition) ---',
+);
+{
+  let totalTime = 0;
+  for (let i = 0; i < iterations; i++) {
+    const start = performance.now();
+    const transposed = transposeComplexMatrix(inputData);
+    const fftResult = reimMatrixFFT(transposed, { inPlace: true });
+    transposeComplexMatrix(fftResult);
+    const elapsed = performance.now() - start;
+    totalTime += elapsed;
+  }
+  const avgTime = totalTime / iterations;
+  console.log(`Total time: ${totalTime.toFixed(2)} ms`);
+  console.log(`Average time per iteration: ${avgTime.toFixed(4)} ms`);
+  console.log(`Time per column: ${(avgTime / numRows).toFixed(6)} ms`);
+}
+
+console.log('');
+console.log('Notes:');
+console.log('- Approach 1 is optimized for column-wise FFT transformation');
+console.log('- Approach 2 transposes twice, which may have overhead');
+console.log('- Both should produce mathematically equivalent results');
diff --git a/src/__tests__/__snapshots__/index.test.ts.snap b/src/__tests__/__snapshots__/index.test.ts.snap
index d3b02335..2d70898d 100644
--- a/src/__tests__/__snapshots__/index.test.ts.snap
+++ b/src/__tests__/__snapshots__/index.test.ts.snap
@@ -9,6 +9,8 @@ exports[`existence of exported functions 1`] = `
   "reimZeroFilling",
   "reimArrayFFT",
   "reimMatrixFFT",
+  "reimMatrixFFTByColumns",
+  "reimMatrixPhaseCorrection",
   "getOutputArray",
   "xAbsolute",
   "xAbsoluteMedian",
diff --git a/src/reim/__tests__/reimAbsolute.test.ts b/src/reim/__tests__/reimAbsolute.test.ts
index d7fd18fe..146aeca5 100644
--- a/src/reim/__tests__/reimAbsolute.test.ts
+++ b/src/reim/__tests__/reimAbsolute.test.ts
@@ -9,3 +9,12 @@ test('reimAbsolute', () => {
 
   expect(result).toStrictEqual(Float64Array.from([0, 5, 10]));
 });
+
+test('reimAbsolute inplace on re', () => {
+  const re = [0, 3, 6];
+  const im = [0, 4, 8];
+  const result = reimAbsolute({ re, im }, { output: re });
+
+  expect(result).toStrictEqual([0, 5, 10]);
+  expect(re).toStrictEqual([0, 5, 10]);
+});
diff --git a/src/reim/__tests__/reimPhaseCorrection.test.ts b/src/reim/__tests__/reimPhaseCorrection.test.ts
index 1ca29fef..aaf55e56 100644
--- a/src/reim/__tests__/reimPhaseCorrection.test.ts
+++ b/src/reim/__tests__/reimPhaseCorrection.test.ts
@@ -125,3 +125,105 @@ test('180 zero order phasing', () => {
 
   expect(diff[index]).toBeCloseTo(0, 4);
 });
+
+test('inPlace option: false (default) should create new arrays', () => {
+  const re = new Float64Array([0, 1, 2, 3]);
+  const im = new Float64Array([0, 1, 2, 3]);
+  const originalRe = Array.from(re);
+  const originalIm = Array.from(im);
+
+  const result = reimPhaseCorrection({ re, im }, Math.PI / 4, 0, {
+    inPlace: false,
+  });
+
+  // Input arrays should not be modified
+  expect(Array.from(re)).toStrictEqual(originalRe);
+  expect(Array.from(im)).toStrictEqual(originalIm);
+
+  // Result should be different objects
+  expect(result.re).not.toBe(re);
+  expect(result.im).not.toBe(im);
+
+  // Result should have correct values
+  expect(result.re).toHaveLength(4);
+  expect(result.im).toHaveLength(4);
+});
+
+test('inPlace option: true should modify input arrays', () => {
+  const re = new Float64Array([0, 1, 2, 3]);
+  const im = new Float64Array([0, 1, 2, 3]);
+
+  const result = reimPhaseCorrection({ re, im }, Math.PI / 4, 0, {
+    inPlace: true,
+  });
+
+  // Result should reference the same arrays
+  expect(result.re).toBe(re);
+  expect(result.im).toBe(im);
+
+  // Input arrays should be modified
+  expect(Array.from(re)).not.toStrictEqual([0, 1, 2, 3]);
+  expect(Array.from(im)).not.toStrictEqual([0, 1, 2, 3]);
+});
+
+test('inPlace option: true and false should produce same values', () => {
+  const phi0 = Math.PI / 3;
+  const phi1 = Math.PI / 6;
+
+  const re1 = new Float64Array([1, 2, 3, 4, 5]);
+  const im1 = new Float64Array([1, 2, 3, 4, 5]);
+  const result1 = reimPhaseCorrection({ re: re1, im: im1 }, phi0, phi1, {
+    inPlace: false,
+  });
+
+  const re2 = new Float64Array([1, 2, 3, 4, 5]);
+  const im2 = new Float64Array([1, 2, 3, 4, 5]);
+  const result2 = reimPhaseCorrection({ re: re2, im: im2 }, phi0, phi1, {
+    inPlace: true,
+  });
+
+  // Values should be the same regardless of inPlace option
+  for (let i = 0; i < result1.re.length; i++) {
+    expect(result1.re[i]).toBeCloseTo(result2.re[i], 10);
+    expect(result1.im[i]).toBeCloseTo(result2.im[i], 10);
+  }
+});
+
+test('inPlace option with reverse flag', () => {
+  const re1 = new Float64Array([1, 2, 3, 4]);
+  const im1 = new Float64Array([1, 2, 3, 4]);
+  const result1 = reimPhaseCorrection({ re: re1, im: im1 }, Math.PI / 4, 0.5, {
+    inPlace: false,
+    reverse: true,
+  });
+
+  const re2 = new Float64Array([1, 2, 3, 4]);
+  const im2 = new Float64Array([1, 2, 3, 4]);
+  const result2 = reimPhaseCorrection({ re: re2, im: im2 }, Math.PI / 4, 0.5, {
+    inPlace: true,
+    reverse: true,
+  });
+
+  // Result should reference the same arrays when inPlace: true
+  expect(result2.re).toBe(re2);
+  expect(result2.im).toBe(im2);
+
+  // Values should be identical
+  for (let i = 0; i < result1.re.length; i++) {
+    expect(result1.re[i]).toBeCloseTo(result2.re[i], 10);
+    expect(result1.im[i]).toBeCloseTo(result2.im[i], 10);
+  }
+});
+
+test('inPlace option with zero phase correction', () => {
+  const re = new Float64Array([1, 2, 3]);
+  const im = new Float64Array([1, 2, 3]);
+
+  const result = reimPhaseCorrection({ re, im }, 0, 0, { inPlace: true });
+
+  // With zero phase correction, values should remain unchanged
+  expect(result.re).toBe(re);
+  expect(result.im).toBe(im);
+  expect(Array.from(result.re)).toStrictEqual([1, 2, 3]);
+  expect(Array.from(result.im)).toStrictEqual([1, 2, 3]);
+});
diff --git a/src/reim/reimAbsolute.ts b/src/reim/reimAbsolute.ts
index 4cfea882..b9420d60 100644
--- a/src/reim/reimAbsolute.ts
+++ b/src/reim/reimAbsolute.ts
@@ -1,18 +1,25 @@
+import type { DoubleArray } from 'cheminfo-types';
+
 import type { DataReIm } from '../types/index.ts';
 
 /**
  * Calculates reimAbsolute value of a complex spectrum.
  * @param data - complex spectrum
+ * @param options
+ * @param options.output
  * @returns - reimAbsolute value
  */
-export function reimAbsolute(data: DataReIm): Float64Array<ArrayBuffer> {
+export function reimAbsolute(
+  data: DataReIm,
+  options: { output?: DoubleArray } = {},
+): DoubleArray {
   const length = data.re.length;
+  const { output = new Float64Array(length) } = options;
   const re = data.re;
   const im = data.im;
-  const newArray = new Float64Array(length);
   for (let i = 0; i < length; i++) {
-    newArray[i] = Math.hypot(re[i], im[i]);
+    output[i] = Math.hypot(re[i], im[i]);
   }
 
-  return newArray;
+  return output;
 }
diff --git a/src/reim/reimAutoPhaseCorrection.ts b/src/reim/reimAutoPhaseCorrection.ts
index 8c5dda5d..8af93ab9 100644
--- a/src/reim/reimAutoPhaseCorrection.ts
+++ b/src/reim/reimAutoPhaseCorrection.ts
@@ -36,6 +36,11 @@ export interface AutoPhaseCorrectionOptions {
    * @default false
    */
   reverse?: boolean;
+  /**
+   * Apply the phase correction directly in the input data
+   * @default false
+   */
+  inPlace?: boolean;
 }
 
 /**
@@ -48,13 +53,14 @@ export interface AutoPhaseCorrectionOptions {
 export function reimAutoPhaseCorrection(
   data: DataReIm,
   options: AutoPhaseCorrectionOptions = {},
-): { data: DataReIm<Float64Array>; ph0: number; ph1: number } {
+): { data: DataReIm; ph0: number; ph1: number } {
   const {
     magnitudeMode = true,
     minRegSize = 30,
     factorNoise = 3,
     maxDistanceToJoin = 256,
     reverse = false,
+    inPlace = false,
   } = options;
 
   const finalPeaks = detectBaselineRegions(data, {
@@ -96,7 +102,7 @@ export function reimAutoPhaseCorrection(
     { re, im },
     toRadians(ph0),
     toRadians(ph1),
-    { reverse },
+    { reverse, inPlace },
   );
 
   return { data: phased, ph0, ph1 };
diff --git a/src/reim/reimPhaseCorrection.ts b/src/reim/reimPhaseCorrection.ts
index 1ab7a722..a7fd49b5 100644
--- a/src/reim/reimPhaseCorrection.ts
+++ b/src/reim/reimPhaseCorrection.ts
@@ -1,7 +1,16 @@
 import type { DataReIm } from '../types/index.ts';
 
 export interface ReimPhaseCorrectionOptions {
+  /**
+   * Apply the phase correction from the last element of the array
+   * @default false
+   */
   reverse?: boolean;
+  /**
+   * Apply the phase correction directly in the input data
+   * @default false
+   */
+  inPlace?: boolean;
 }
 
 /**
@@ -10,25 +19,31 @@ export interface ReimPhaseCorrectionOptions {
  * @param phi0 - Angle in radians for zero order phase correction
  * @param phi1 - Angle in radians for first order phase correction
  * @param options
- * @returns - returns a new object {re:[], im:[]}
+ * @returns - returns a new object {re:[], im:[]} unless inPlace=true
  */
 export function reimPhaseCorrection(
   data: DataReIm,
   phi0 = 0,
   phi1 = 0,
   options: ReimPhaseCorrectionOptions = {},
-): DataReIm<Float64Array> {
-  const { reverse = false } = options;
+): DataReIm {
+  const { reverse = false, inPlace = false } = options;
 
   phi0 = Number.isFinite(phi0) ? phi0 : 0;
   phi1 = Number.isFinite(phi1) ? phi1 : 0;
 
+  const length = data.re.length;
+
+  // Decide target arrays
   const re = data.re;
   const im = data.im;
-  const length = data.re.length;
+
+  const outRe = inPlace ? re : new Float64Array(length);
+  const outIm = inPlace ? im : new Float64Array(length);
 
   let firstAngle = phi0;
   let delta = phi1 / length;
+
   if (reverse) {
     delta *= -1;
     firstAngle += phi1;
@@ -36,20 +51,24 @@ export function reimPhaseCorrection(
 
   const alpha = 2 * Math.sin(delta / 2) ** 2;
   const beta = Math.sin(delta);
+
   let cosTheta = Math.cos(firstAngle);
   let sinTheta = Math.sin(firstAngle);
 
-  const newRe = new Float64Array(length);
-  const newIm = new Float64Array(length);
   for (let i = 0; i < length; i++) {
-    newRe[i] = re[i] * cosTheta - im[i] * sinTheta;
-    newIm[i] = im[i] * cosTheta + re[i] * sinTheta;
-    // calculate angles i+1 from i
+    const r = re[i];
+    const ii = im[i];
+
+    outRe[i] = r * cosTheta - ii * sinTheta;
+    outIm[i] = ii * cosTheta + r * sinTheta;
+
+    // Recursive angle update (stable incremental rotation)
     const newCosTheta = cosTheta - (alpha * cosTheta + beta * sinTheta);
     const newSinTheta = sinTheta - (alpha * sinTheta - beta * cosTheta);
+
     cosTheta = newCosTheta;
     sinTheta = newSinTheta;
   }
 
-  return { re: newRe, im: newIm };
+  return { re: outRe, im: outIm };
 }
diff --git a/src/reimMatrix/__tests__/reimMatrixFFTByColumns.test.ts b/src/reimMatrix/__tests__/reimMatrixFFTByColumns.test.ts
new file mode 100644
index 00000000..017774e8
--- /dev/null
+++ b/src/reimMatrix/__tests__/reimMatrixFFTByColumns.test.ts
@@ -0,0 +1,321 @@
+import { expect, test } from 'vitest';
+
+import { reimMatrixFFTByColumns } from '../reimMatrixFFTByColumns.ts';
+
+test('reimMatrixFFTByColumns: round-trip on a single column', () => {
+  const data = {
+    re: [
+      Float64Array.from([0, 3]), // row 0
+      Float64Array.from([0, 6]), // row 1
+      Float64Array.from([0, 5]), // row 2
+      Float64Array.from([0, 4]), // row 3
+    ],
+    im: [
+      Float64Array.from([0, 4]),
+      Float64Array.from([0, 8]),
+      Float64Array.from([0, 3]),
+      Float64Array.from([0, 2]),
+    ],
+  };
+
+  const transformed = reimMatrixFFTByColumns(data, { applyZeroShift: true });
+  const restored = reimMatrixFFTByColumns(transformed, {
+    inverse: true,
+    applyZeroShift: true,
+  });
+
+  for (let i = 0; i < data.re.length; i++) {
+    expect(restored.re[i]).toStrictEqual(data.re[i]);
+    expect(restored.im[i]).toStrictEqual(data.im[i]);
+  }
+});
+
+test('reimMatrixFFTByColumns: round-trip on multiple columns', () => {
+  const data = {
+    re: [
+      Float64Array.from([1, 0, 2, 3]),
+      Float64Array.from([0, 1, 4, 5]),
+      Float64Array.from([0, 3, 6, 7]),
+      Float64Array.from([1, 2, 3, 4]),
+    ],
+    im: [
+      Float64Array.from([0, 0, 1, 2]),
+      Float64Array.from([0, 0, 2, 3]),
+      Float64Array.from([0, 4, 8, 1]),
+      Float64Array.from([0, 1, 2, 3]),
+    ],
+  };
+
+  const transformed = reimMatrixFFTByColumns(data);
+  const restored = reimMatrixFFTByColumns(transformed, { inverse: true });
+
+  for (let i = 0; i < data.re.length; i++) {
+    expect(restored.re[i]).toStrictEqual(data.re[i]);
+    expect(restored.im[i]).toStrictEqual(data.im[i]);
+  }
+});
+
+test('reimMatrixFFTByColumns: applyZeroShift round-trip on multiple columns', () => {
+  const data = {
+    re: [
+      Float64Array.from([0, 3, 6, 5]),
+      Float64Array.from([1, 2, 3, 4]),
+      Float64Array.from([2, 1, 5, 3]),
+      Float64Array.from([0, 4, 2, 1]),
+    ],
+    im: [
+      Float64Array.from([0, 4, 8, 3]),
+      Float64Array.from([0, 1, 0, 1]),
+      Float64Array.from([1, 2, 3, 0]),
+      Float64Array.from([2, 0, 1, 2]),
+    ],
+  };
+
+  const transformed = reimMatrixFFTByColumns(data, { applyZeroShift: true });
+  const restored = reimMatrixFFTByColumns(transformed, {
+    inverse: true,
+    applyZeroShift: true,
+  });
+
+  for (let i = 0; i < data.re.length; i++) {
+    expect(restored.re[i]).toStrictEqual(data.re[i]);
+    expect(restored.im[i]).toStrictEqual(data.im[i]);
+  }
+});
+
+test('reimMatrixFFTByColumns: output arrays are independent (no shared buffers)', () => {
+  const data = {
+    re: [
+      Float64Array.from([1, 0]),
+      Float64Array.from([0, 1]),
+      Float64Array.from([1, 1]),
+      Float64Array.from([0, 0]),
+    ],
+    im: [
+      Float64Array.from([0, 0]),
+      Float64Array.from([0, 0]),
+      Float64Array.from([0, 0]),
+      Float64Array.from([0, 0]),
+    ],
+  };
+
+  const result = reimMatrixFFTByColumns(data);
+
+  expect(result.re[0]).not.toBe(result.re[1]);
+  expect(result.im[0]).not.toBe(result.im[1]);
+  expect(result.re[0]).not.toBe(data.re[0]);
+  expect(result.im[0]).not.toBe(data.im[0]);
+});
+
+test('reimMatrixFFTByColumns: returns empty matrix for empty input', () => {
+  expect(reimMatrixFFTByColumns({ re: [], im: [] })).toStrictEqual({
+    re: [],
+    im: [],
+  });
+});
+
+test('reimMatrixFFTByColumns: throws RangeError when rows have different lengths', () => {
+  const data = {
+    re: [
+      Float64Array.from([1, 0, 0, 0]),
+      Float64Array.from([0, 1, 0, 0, 0, 0, 0, 0]),
+      Float64Array.from([1, 0, 0, 0]),
+      Float64Array.from([0, 1, 0, 0]),
+    ],
+    im: [
+      Float64Array.from([0, 0, 0, 0]),
+      Float64Array.from([0, 0, 0, 0, 0, 0, 0, 0]),
+      Float64Array.from([0, 0, 0, 0]),
+      Float64Array.from([0, 0, 0, 0]),
+    ],
+  };
+
+  expect(() => reimMatrixFFTByColumns(data)).toThrowError(RangeError);
+});
+
+test('reimMatrixFFTByColumns: throws RangeError indicating which row has the wrong length', () => {
+  const data = {
+    re: [
+      Float64Array.from([1, 0, 0, 0]),
+      Float64Array.from([0, 1, 0, 0]),
+      Float64Array.from([0, 0, 1, 0, 0, 0, 0, 0]),
+      Float64Array.from([1, 1, 0, 0]),
+    ],
+    im: [
+      Float64Array.from([0, 0, 0, 0]),
+      Float64Array.from([0, 0, 0, 0]),
+      Float64Array.from([0, 0, 0, 0, 0, 0, 0, 0]),
+      Float64Array.from([0, 0, 0, 0]),
+    ],
+  };
+
+  expect(() => reimMatrixFFTByColumns(data)).toThrowError(/row 2/);
+});
+
+test('reimMatrixFFTByColumns inPlace: result shares re/im references with input', () => {
+  const data = {
+    re: [
+      Float64Array.from([1, 0]),
+      Float64Array.from([0, 3]),
+      Float64Array.from([2, 1]),
+      Float64Array.from([1, 2]),
+    ],
+    im: [
+      Float64Array.from([0, 0]),
+      Float64Array.from([0, 4]),
+      Float64Array.from([1, 0]),
+      Float64Array.from([0, 1]),
+    ],
+  };
+
+  const result = reimMatrixFFTByColumns(data, { inPlace: true });
+
+  expect(result.re[0]).toBe(data.re[0]);
+  expect(result.im[0]).toBe(data.im[0]);
+  expect(result.re[1]).toBe(data.re[1]);
+  expect(result.im[1]).toBe(data.im[1]);
+});
+
+test('reimMatrixFFTByColumns inPlace: round-trip restores original values', () => {
+  const data = {
+    re: [
+      Float64Array.from([0, 3, 6, 5]),
+      Float64Array.from([1, 2, 3, 4]),
+      Float64Array.from([2, 1, 5, 3]),
+      Float64Array.from([0, 4, 2, 1]),
+    ],
+    im: [
+      Float64Array.from([0, 4, 8, 3]),
+      Float64Array.from([0, 1, 0, 1]),
+      Float64Array.from([1, 2, 3, 0]),
+      Float64Array.from([2, 0, 1, 2]),
+    ],
+  };
+  const originals = {
+    re: data.re.map((row) => Float64Array.from(row)),
+    im: data.im.map((row) => Float64Array.from(row)),
+  };
+
+  reimMatrixFFTByColumns(data, { inPlace: true, applyZeroShift: true });
+  reimMatrixFFTByColumns(data, {
+    inPlace: true,
+    inverse: true,
+    applyZeroShift: true,
+  });
+
+  for (let i = 0; i < data.re.length; i++) {
+    expect(data.re[i]).toStrictEqual(originals.re[i]);
+    expect(data.im[i]).toStrictEqual(originals.im[i]);
+  }
+});
+
+test('reimMatrixFFTByColumns: processes each column independently', () => {
+  // Create a matrix where each column has distinct, recognizable patterns
+  const data = {
+    re: [
+      Float64Array.from([1, 2, 3, 4]), // column 0: [1, 0, 0, 0], column 1: [2, 0, 0, 0], etc.
+      Float64Array.from([0, 0, 0, 0]),
+      Float64Array.from([0, 0, 0, 0]),
+      Float64Array.from([0, 0, 0, 0]),
+    ],
+    im: [
+      Float64Array.from([0, 0, 0, 0]),
+      Float64Array.from([0, 0, 0, 0]),
+      Float64Array.from([0, 0, 0, 0]),
+      Float64Array.from([0, 0, 0, 0]),
+    ],
+  };
+
+  const transformed = reimMatrixFFTByColumns(data);
+
+  // Column 0 contains [1, 0, 0, 0]
+  const col0Re = transformed.re.map((row) => row[0]);
+  const col0Im = transformed.im.map((row) => row[0]);
+
+  // Column 1 contains [2, 0, 0, 0]
+  const col1Re = transformed.re.map((row) => row[1]);
+  const col1Im = transformed.im.map((row) => row[1]);
+
+  // Column 2 contains [3, 0, 0, 0]
+  const col2Re = transformed.re.map((row) => row[2]);
+  const col2Im = transformed.im.map((row) => row[2]);
+
+  // The FFT of [1, 0, 0, 0] scaled by amplitude
+  // FFT of [2, 0, 0, 0] should be the FFT of [1, 0, 0, 0] scaled by 2
+  // FFT of [3, 0, 0, 0] should be the FFT of [1, 0, 0, 0] scaled by 3
+
+  // Verify scaling relationship: col1 = col0 * 2, col2 = col0 * 3
+  for (let i = 0; i < 4; i++) {
+    expect(col1Re[i]).toBeCloseTo(col0Re[i] * 2, 10);
+    expect(col1Im[i]).toBeCloseTo(col0Im[i] * 2, 10);
+    expect(col2Re[i]).toBeCloseTo(col0Re[i] * 3, 10);
+    expect(col2Im[i]).toBeCloseTo(col0Im[i] * 3, 10);
+  }
+});
+
+test('reimMatrixFFTByColumns: handles matrix with single column', () => {
+  const data = {
+    re: [
+      Float64Array.from([1]),
+      Float64Array.from([2]),
+      Float64Array.from([3]),
+      Float64Array.from([4]),
+    ],
+    im: [
+      Float64Array.from([0]),
+      Float64Array.from([0]),
+      Float64Array.from([0]),
+      Float64Array.from([0]),
+    ],
+  };
+
+  const transformed = reimMatrixFFTByColumns(data);
+  const restored = reimMatrixFFTByColumns(transformed, { inverse: true });
+
+  for (let i = 0; i < data.re.length; i++) {
+    expect(restored.re[i]).toStrictEqual(data.re[i]);
+    expect(restored.im[i]).toStrictEqual(data.im[i]);
+  }
+});
+
+test('reimMatrixFFTByColumns: handles multiple row matrix with few columns', () => {
+  const data = {
+    re: [
+      Float64Array.from([1, 2]),
+      Float64Array.from([3, 4]),
+      Float64Array.from([5, 6]),
+      Float64Array.from([7, 8]),
+    ],
+    im: [
+      Float64Array.from([0, 1]),
+      Float64Array.from([2, 3]),
+      Float64Array.from([1, 0]),
+      Float64Array.from([0, 2]),
+    ],
+  };
+
+  const transformed = reimMatrixFFTByColumns(data);
+  const restored = reimMatrixFFTByColumns(transformed, { inverse: true });
+
+  expect(restored.re[0]).toStrictEqual(data.re[0]);
+  expect(restored.im[0]).toStrictEqual(data.im[0]);
+});
+
+test('reimMatrixFFTByColumns: throws error for matrices with mismatched column counts', () => {
+  const data = {
+    re: [
+      Float64Array.from([1, 2]),
+      Float64Array.from([3]), // different column count
+      Float64Array.from([5, 6]),
+      Float64Array.from([7, 8]),
+    ],
+    im: [
+      Float64Array.from([0, 0]),
+      Float64Array.from([0]),
+      Float64Array.from([0, 0]),
+      Float64Array.from([0, 0]),
+    ],
+  };
+
+  expect(() => reimMatrixFFTByColumns(data)).toThrowError(RangeError);
+});
diff --git a/src/reimMatrix/__tests__/reimMatrixPhaseCorrection.test.ts b/src/reimMatrix/__tests__/reimMatrixPhaseCorrection.test.ts
new file mode 100644
index 00000000..60c4edbe
--- /dev/null
+++ b/src/reimMatrix/__tests__/reimMatrixPhaseCorrection.test.ts
@@ -0,0 +1,390 @@
+import { expect, test } from 'vitest';
+
+import { reimMatrixPhaseCorrection } from '../reimMatrixPhaseCorrection.ts';
+
+test('reimMatrixPhaseCorrection: no correction when phi0 and phi1 are 0', () => {
+  const data = {
+    re: [Float64Array.from([0, 1, 2, 3]), Float64Array.from([4, 5, 6, 7])],
+    im: [Float64Array.from([0, 1, 2, 3]), Float64Array.from([4, 5, 6, 7])],
+  };
+
+  const result = reimMatrixPhaseCorrection(data, 0, 0);
+
+  expect(Array.from(result.re[0])).toStrictEqual(Array.from(data.re[0]));
+  expect(Array.from(result.re[1])).toStrictEqual(Array.from(data.re[1]));
+  expect(Array.from(result.im[0])).toStrictEqual(Array.from(data.im[0]));
+  expect(Array.from(result.im[1])).toStrictEqual(Array.from(data.im[1]));
+});
+
+test('reimMatrixPhaseCorrection: direction rows by default', () => {
+  const data = {
+    re: [Float64Array.from([1, 0, 0, 0])],
+    im: [Float64Array.from([0, 0, 0, 0])],
+  };
+
+  const resultDefault = reimMatrixPhaseCorrection(data, Math.PI / 4, 0);
+  const resultRows = reimMatrixPhaseCorrection(data, Math.PI / 4, 0, {
+    direction: 'rows',
+  });
+
+  expect(Array.from(resultDefault.re[0])).toStrictEqual(
+    Array.from(resultRows.re[0]),
+  );
+  expect(Array.from(resultDefault.im[0])).toStrictEqual(
+    Array.from(resultRows.im[0]),
+  );
+});
+
+test('reimMatrixPhaseCorrection: zero order phase correction on rows', () => {
+  const data = {
+    re: [Float64Array.from([1, 0, 0, 0])],
+    im: [Float64Array.from([0, 0, 0, 0])],
+  };
+
+  const result = reimMatrixPhaseCorrection(data, Math.PI, 0, {
+    direction: 'rows',
+  });
+
+  // After applying pi phase, the real and imaginary parts should be negated
+  expect(result.re[0][0]).toBeCloseTo(-1, 5);
+  expect(result.im[0][0]).toBeCloseTo(0, 5);
+});
+
+test('reimMatrixPhaseCorrection: first order phase correction on rows', () => {
+  const data = {
+    re: [Float64Array.from([1, 1, 1, 1])],
+    im: [Float64Array.from([0, 0, 0, 0])],
+  };
+
+  const result = reimMatrixPhaseCorrection(data, 0, Math.PI / 2, {
+    direction: 'rows',
+  });
+
+  // With first order phase correction, different elements get different phases
+  expect(result.re[0][0]).toBeCloseTo(1, 5); // First element gets phase 0
+  expect(Math.hypot(result.re[0][3], result.im[0][3])).toBeCloseTo(1, 5); // Magnitude preserved
+});
+
+test('reimMatrixPhaseCorrection: combined zero and first order correction on rows', () => {
+  const data = {
+    re: [Float64Array.from([1, 1, 1, 1])],
+    im: [Float64Array.from([0, 0, 0, 0])],
+  };
+
+  const result = reimMatrixPhaseCorrection(data, Math.PI / 4, Math.PI / 2, {
+    direction: 'rows',
+  });
+
+  // All magnitudes should be preserved
+  for (let i = 0; i < 4; i++) {
+    const magnitude = Math.hypot(result.re[0][i], result.im[0][i]);
+
+    expect(magnitude).toBeCloseTo(1, 5);
+  }
+});
+
+test('reimMatrixPhaseCorrection: phase correction on columns', () => {
+  const data = {
+    re: [
+      Float64Array.from([1, 0]),
+      Float64Array.from([1, 0]),
+      Float64Array.from([1, 0]),
+      Float64Array.from([1, 0]),
+    ],
+    im: [
+      Float64Array.from([0, 0]),
+      Float64Array.from([0, 0]),
+      Float64Array.from([0, 0]),
+      Float64Array.from([0, 0]),
+    ],
+  };
+
+  const result = reimMatrixPhaseCorrection(data, Math.PI / 4, 0, {
+    direction: 'columns',
+  });
+
+  // All rows in the same column should be corrected the same way
+  for (let i = 0; i < 4; i++) {
+    expect(result.re[i][0]).toBeCloseTo(result.re[0][0], 5);
+    expect(result.im[i][0]).toBeCloseTo(result.im[0][0], 5);
+  }
+});
+
+test('reimMatrixPhaseCorrection: reverse option on rows', () => {
+  const data = {
+    re: [Float64Array.from([1, 1, 1, 1])],
+    im: [Float64Array.from([0, 0, 0, 0])],
+  };
+
+  const resultForward = reimMatrixPhaseCorrection(
+    structuredClone(data),
+    0,
+    Math.PI / 2,
+    {
+      direction: 'rows',
+      reverse: false,
+    },
+  );
+
+  const resultReverse = reimMatrixPhaseCorrection(data, 0, Math.PI / 2, {
+    direction: 'rows',
+    reverse: true,
+  });
+
+  // Reverse should apply phase in opposite direction
+  // Both should have same magnitude but different phases for later elements
+  for (let i = 0; i < 4; i++) {
+    const magForward = Math.hypot(
+      resultForward.re[0][i],
+      resultForward.im[0][i],
+    );
+    const magReverse = Math.hypot(
+      resultReverse.re[0][i],
+      resultReverse.im[0][i],
+    );
+
+    expect(Math.abs(magForward - magReverse)).toBeLessThan(1e-14);
+  }
+
+  // But later elements should have different real parts
+  expect(resultForward.re[0][2]).not.toBeCloseTo(resultReverse.re[0][2], 2);
+});
+
+test('reimMatrixPhaseCorrection: inPlace option creates new arrays when false', () => {
+  const data = {
+    re: [Float64Array.from([1, 2, 3])],
+    im: [Float64Array.from([4, 5, 6])],
+  };
+
+  const result = reimMatrixPhaseCorrection(data, Math.PI / 4, 0, {
+    inPlace: false,
+  });
+
+  // Results should not share array references with input
+  expect(result.re[0]).not.toStrictEqual(data.re[0]);
+  expect(result.im[0]).not.toStrictEqual(data.im[0]);
+});
+
+test('reimMatrixPhaseCorrection: inPlace option modifies original arrays when true', () => {
+  const data = {
+    re: [Float64Array.from([1, 0, 0, 0])],
+    im: [Float64Array.from([0, 0, 0, 0])],
+  };
+
+  const originalReBuffer = data.re[0].buffer;
+  const originalImBuffer = data.im[0].buffer;
+
+  const result = reimMatrixPhaseCorrection(data, Math.PI / 2, 0, {
+    inPlace: true,
+  });
+
+  // Results should share the same array buffers
+  expect(result.re[0].buffer).toStrictEqual(originalReBuffer);
+  expect(result.im[0].buffer).toStrictEqual(originalImBuffer);
+});
+
+test('reimMatrixPhaseCorrection: multiple rows with row-wise correction', () => {
+  const data = {
+    re: [
+      Float64Array.from([1, 1, 1, 1]),
+      Float64Array.from([2, 2, 2, 2]),
+      Float64Array.from([3, 3, 3, 3]),
+    ],
+    im: [
+      Float64Array.from([0, 0, 0, 0]),
+      Float64Array.from([0, 0, 0, 0]),
+      Float64Array.from([0, 0, 0, 0]),
+    ],
+  };
+
+  const result = reimMatrixPhaseCorrection(data, Math.PI / 4, Math.PI / 6, {
+    direction: 'rows',
+  });
+
+  // Each row should be corrected independently
+  expect(result.re).toHaveLength(3);
+  expect(result.im).toHaveLength(3);
+
+  // Magnitude of each row should be preserved
+  for (let row = 0; row < 3; row++) {
+    for (let col = 0; col < 4; col++) {
+      const expectedMagnitude = row + 1; // Original magnitude
+      const actualMagnitude = Math.hypot(
+        result.re[row][col],
+        result.im[row][col],
+      );
+
+      expect(actualMagnitude).toBeCloseTo(expectedMagnitude, 5);
+    }
+  }
+});
+
+test('reimMatrixPhaseCorrection: multiple columns with column-wise correction', () => {
+  const data = {
+    re: [
+      Float64Array.from([1, 2, 3]),
+      Float64Array.from([1, 2, 3]),
+      Float64Array.from([1, 2, 3]),
+      Float64Array.from([1, 2, 3]),
+    ],
+    im: [
+      Float64Array.from([0, 0, 0]),
+      Float64Array.from([0, 0, 0]),
+      Float64Array.from([0, 0, 0]),
+      Float64Array.from([0, 0, 0]),
+    ],
+  };
+
+  const result = reimMatrixPhaseCorrection(data, Math.PI / 4, Math.PI / 6, {
+    direction: 'columns',
+  });
+
+  // Each column should be corrected independently
+  for (let col = 0; col < 3; col++) {
+    for (let row = 0; row < 4; row++) {
+      const expectedMagnitude = col + 1;
+      const actualMagnitude = Math.hypot(
+        result.re[row][col],
+        result.im[row][col],
+      );
+
+      expect(actualMagnitude).toBeCloseTo(expectedMagnitude, 5);
+    }
+  }
+});
+
+test('reimMatrixPhaseCorrection: empty matrix', () => {
+  const data = { re: [], im: [] };
+
+  const result = reimMatrixPhaseCorrection(data, Math.PI / 4, Math.PI / 6);
+
+  expect(result.re).toStrictEqual([]);
+  expect(result.im).toStrictEqual([]);
+});
+
+test('reimMatrixPhaseCorrection: single element', () => {
+  const data = {
+    re: [Float64Array.from([5])],
+    im: [Float64Array.from([0])],
+  };
+
+  const result = reimMatrixPhaseCorrection(data, Math.PI / 4, 0);
+
+  expect(result.re[0][0]).toBeCloseTo(5 * Math.cos(Math.PI / 4), 5);
+  expect(result.im[0][0]).toBeCloseTo(5 * Math.sin(Math.PI / 4), 5);
+});
+
+test('reimMatrixPhaseCorrection: throws RangeError when rows have different lengths', () => {
+  const data = {
+    re: [
+      Float64Array.from([1, 0, 0, 0]),
+      Float64Array.from([0, 1, 0, 0, 0, 0, 0, 0]),
+    ],
+    im: [
+      Float64Array.from([0, 0, 0, 0]),
+      Float64Array.from([0, 0, 0, 0, 0, 0, 0, 0]),
+    ],
+  };
+
+  expect(() => reimMatrixPhaseCorrection(data, Math.PI / 4, 0)).toThrowError(
+    RangeError,
+  );
+});
+
+test('reimMatrixPhaseCorrection: throws RangeError indicating which row has the wrong length', () => {
+  const data = {
+    re: [
+      Float64Array.from([1, 0, 0, 0]),
+      Float64Array.from([0, 1, 0, 0]),
+      Float64Array.from([0, 0, 1, 0, 0, 0, 0, 0]),
+    ],
+    im: [
+      Float64Array.from([0, 0, 0, 0]),
+      Float64Array.from([0, 0, 0, 0]),
+      Float64Array.from([0, 0, 0, 0, 0, 0, 0, 0]),
+    ],
+  };
+
+  expect(() => reimMatrixPhaseCorrection(data, Math.PI / 4, 0)).toThrowError(
+    /row 2/,
+  );
+});
+
+test('reimMatrixPhaseCorrection: handles undefined options gracefully', () => {
+  const data = {
+    re: [Float64Array.from([1, 0, 0, 0])],
+    im: [Float64Array.from([0, 0, 0, 0])],
+  };
+
+  const result = reimMatrixPhaseCorrection(data);
+
+  expect(result.re).toHaveLength(1);
+  expect(result.im).toHaveLength(1);
+});
+
+test('reimMatrixPhaseCorrection: handles NaN and Infinity gracefully', () => {
+  const data = {
+    re: [Float64Array.from([1, 0, 0, 0])],
+    im: [Float64Array.from([0, 0, 0, 0])],
+  };
+
+  const result = reimMatrixPhaseCorrection(data, Number.NaN, Infinity);
+
+  // NaN and Infinity should be converted to 0
+  expect(Array.from(result.re[0])).toStrictEqual(Array.from(data.re[0]));
+  expect(Array.from(result.im[0])).toStrictEqual(Array.from(data.im[0]));
+});
+
+test('reimMatrixPhaseCorrection: all options combined', () => {
+  const data = {
+    re: [Float64Array.from([1, 1, 1, 1]), Float64Array.from([2, 2, 2, 2])],
+    im: [Float64Array.from([0, 0, 0, 0]), Float64Array.from([0, 0, 0, 0])],
+  };
+
+  const result = reimMatrixPhaseCorrection(data, Math.PI / 6, Math.PI / 4, {
+    reverse: true,
+    inPlace: false,
+    direction: 'rows',
+  });
+
+  // Verify it produces a valid result
+  expect(result.re).toHaveLength(2);
+  expect(result.im).toHaveLength(2);
+
+  // Magnitudes should be preserved
+  for (let row = 0; row < 2; row++) {
+    for (let col = 0; col < 4; col++) {
+      const expectedMagnitude = row + 1;
+      const actualMagnitude = Math.hypot(
+        result.re[row][col],
+        result.im[row][col],
+      );
+
+      expect(actualMagnitude).toBeCloseTo(expectedMagnitude, 5);
+    }
+  }
+});
+
+test('reimMatrixPhaseCorrection: preserves magnitude along all directions', () => {
+  const data = {
+    re: [Float64Array.from([3, 4, 5]), Float64Array.from([6, 8, 10])],
+    im: [Float64Array.from([4, 3, 12]), Float64Array.from([8, 6, 24])],
+  };
+
+  const originalMagnitudes = data.re.map((row, i) =>
+    row.map((r, j) => Math.hypot(r, data.im[i][j])),
+  );
+
+  const result = reimMatrixPhaseCorrection(data, Math.PI / 3, Math.PI / 6, {
+    direction: 'rows',
+  });
+
+  // Verify magnitudes are preserved
+  for (let row = 0; row < data.re.length; row++) {
+    for (let col = 0; col < data.re[row].length; col++) {
+      const magnitude = Math.hypot(result.re[row][col], result.im[row][col]);
+
+      expect(magnitude).toBeCloseTo(originalMagnitudes[row][col], 5);
+    }
+  }
+});
diff --git a/src/reimMatrix/index.ts b/src/reimMatrix/index.ts
index 5ab93eda..16444afc 100644
--- a/src/reimMatrix/index.ts
+++ b/src/reimMatrix/index.ts
@@ -1 +1,3 @@
 export * from './reimMatrixFFT.ts';
+export * from './reimMatrixFFTByColumns.ts';
+export * from './reimMatrixPhaseCorrection.ts';
diff --git a/src/reimMatrix/reimMatrixFFTByColumns.ts b/src/reimMatrix/reimMatrixFFTByColumns.ts
new file mode 100644
index 00000000..4c29dda1
--- /dev/null
+++ b/src/reimMatrix/reimMatrixFFTByColumns.ts
@@ -0,0 +1,102 @@
+import FFT from 'fft.js';
+
+import { zeroShift } from '../reim/zeroShift.ts';
+import type { DataReImMatrix } from '../types/index.ts';
+
+export interface ReimMatrixFFTByColumnsOptions {
+  inverse?: boolean;
+  applyZeroShift?: boolean;
+  /**
+   * Write the result back into the input arrays instead of allocating new ones.
+   * @default false
+   */
+  inPlace?: boolean;
+}
+
+/**
+ * Apply FFT to each column of a complex matrix, reusing a single FFT instance for
+ * all columns to avoid repeated instantiation overhead.
+ * All columns must have the same length (all rows must have the same number of columns).
+ * @param data - complex matrix with re and im arrays of Float64Array rows
+ * @param options - options
+ * @returns complex matrix with transformed columns
+ */
+export function reimMatrixFFTByColumns(
+  data: DataReImMatrix,
+  options: ReimMatrixFFTByColumnsOptions = {},
+): DataReImMatrix {
+  const { re, im } = data;
+  const numRows = re.length;
+
+  if (numRows === 0) return { re: [], im: [] };
+
+  const { inverse = false, applyZeroShift = false, inPlace = false } = options;
+
+  const numColumns = re[0].length;
+  const csize = numRows << 1;
+
+  // Validate all rows have the same number of columns
+  for (let j = 0; j < numRows; j++) {
+    if (re[j].length !== numColumns || im[j].length !== numColumns) {
+      throw new RangeError(
+        `All rows must have the same length. Expected ${numColumns} but row ${j} has length ${re[j].length}.`,
+      );
+    }
+  }
+
+  if (numColumns === 0) {
+    return { re: new Array(numRows), im: new Array(numRows) };
+  }
+
+  // Single FFT instance and working buffers reused across all columns
+  const fft = new FFT(numRows);
+  const complexArray = new Float64Array(csize);
+  const output = new Float64Array(csize);
+
+  const resultRe = new Array<Float64Array>(numRows);
+  const resultIm = new Array<Float64Array>(numRows);
+
+  // Initialize result arrays
+  for (let i = 0; i < numRows; i++) {
+    resultRe[i] = inPlace ? re[i] : new Float64Array(numColumns);
+    resultIm[i] = inPlace ? im[i] : new Float64Array(numColumns);
+  }
+
+  // Process each column
+  for (let col = 0; col < numColumns; col++) {
+    // Extract column data into complex array
+    for (let row = 0; row < numRows; row++) {
+      complexArray[row << 1] = re[row][col];
+      complexArray[(row << 1) + 1] = im[row][col];
+    }
+
+    const colRe = new Float64Array(numRows);
+    const colIm = new Float64Array(numRows);
+
+    if (inverse) {
+      const input = applyZeroShift
+        ? zeroShift(complexArray, true)
+        : complexArray;
+      fft.inverseTransform(output, input);
+      for (let i = 0; i < csize; i += 2) {
+        colRe[i >>> 1] = output[i];
+        colIm[i >>> 1] = output[i + 1];
+      }
+    } else {
+      fft.transform(output, complexArray);
+      const source = applyZeroShift ? zeroShift(output) : output;
+      for (let i = 0; i < csize; i += 2) {
+        colRe[i >>> 1] = source[i];
+        colIm[i >>> 1] = source[i + 1];
+      }
+    }
+
+    // Store transformed column back into result matrix
+    for (let row = 0; row < numRows; row++) {
+      resultRe[row][col] = colRe[row];
+      resultIm[row][col] = colIm[row];
+    }
+  }
+
+  return { re: resultRe, im: resultIm };
+}
diff --git a/src/reimMatrix/reimMatrixPhaseCorrection.ts b/src/reimMatrix/reimMatrixPhaseCorrection.ts
new file mode 100644
index 00000000..db30a9ae
--- /dev/null
+++ b/src/reimMatrix/reimMatrixPhaseCorrection.ts
@@ -0,0 +1,119 @@
+import type { DataReImMatrix } from '../types/index.ts';
+
+export interface ReimMatrixPhaseCorrectionOptions {
+  reverse?: boolean;
+  inPlace?: boolean;
+  /**
+   * Direction for phase correction: 'rows' or 'columns'
+   * @default 'rows'
+   */
+  direction?: 'rows' | 'columns';
+}
+
+/**
+ * Apply phase correction to a complex matrix along rows or columns.
+ * All rows must have the same length.
+ * @param data - complex matrix with re and im arrays of Float64Array rows
+ * @param phi0 - Angle in radians for zero order phase correction
+ * @param phi1 - Angle in radians for first order phase correction
+ * @param options - options including direction ('rows' or 'columns', default 'rows')
+ * @returns - complex matrix with corrected rows or columns
+ */
+export function reimMatrixPhaseCorrection(
+  data: DataReImMatrix,
+  phi0 = 0,
+  phi1 = 0,
+  options: ReimMatrixPhaseCorrectionOptions = {},
+): DataReImMatrix {
+  const { reverse = false, inPlace = false, direction = 'rows' } = options;
+
+  phi0 = Number.isFinite(phi0) ? phi0 : 0;
+  phi1 = Number.isFinite(phi1) ? phi1 : 0;
+
+  const { re, im } = data;
+  const numRows = re.length;
+
+  if (numRows === 0) return { re: [], im: [] };
+
+  const numColumns = re[0].length;
+
+  for (let j = 0; j < numRows; j++) {
+    if (re[j].length !== numColumns || im[j].length !== numColumns) {
+      throw new RangeError(
+        `All rows must have the same length. Expected ${numColumns} but row ${j} has length ${re[j].length}.`,
+      );
+    }
+  }
+
+  const resultRe = new Array<Float64Array>(numRows);
+  const resultIm = new Array<Float64Array>(numRows);
+
+  // Initialize result arrays
+  for (let i = 0; i < numRows; i++) {
+    resultRe[i] = inPlace ? re[i] : new Float64Array(numColumns);
+    resultIm[i] = inPlace ? im[i] : new Float64Array(numColumns);
+  }
+
+  const size = direction === 'rows' ? numColumns : numRows;
+
+  let delta = phi1 / (size - 1);
+
+  if (reverse) {
+    delta = -delta;
+  }
+
+  const alpha = 2 * Math.sin(delta / 2) ** 2;
+  const beta = Math.sin(delta);
+
+  if (direction === 'rows') {
+    for (let j = 0; j < numRows; j++) {
+      let firstAngle = phi0;
+      if (reverse) {
+        firstAngle += phi1;
+      }
+
+      let cosTheta = Math.cos(firstAngle);
+      let sinTheta = Math.sin(firstAngle);
+
+      for (let i = 0; i < numColumns; i++) {
+        const r = re[j][i];
+        const ii = im[j][i];
+
+        resultRe[j][i] = r * cosTheta - ii * sinTheta;
+        resultIm[j][i] = ii * cosTheta + r * sinTheta;
+
+        const newCosTheta = cosTheta - (alpha * cosTheta + beta * sinTheta);
+        const newSinTheta = sinTheta - (alpha * sinTheta - beta * cosTheta);
+
+        cosTheta = newCosTheta;
+        sinTheta = newSinTheta;
+      }
+    }
+  } else {
+    for (let col = 0; col < numColumns; col++) {
+      let firstAngle = phi0;
+      if (reverse) {
+        firstAngle += phi1;
+      }
+
+      let cosTheta = Math.cos(firstAngle);
+      let sinTheta = Math.sin(firstAngle);
+
+      for (let row = 0; row < numRows; row++) {
+        const r = re[row][col];
+        const ii = im[row][col];
+
+        resultRe[row][col] = r * cosTheta - ii * sinTheta;
+        resultIm[row][col] = ii * cosTheta + r * sinTheta;
+
+        const newCosTheta = cosTheta - (alpha * cosTheta + beta * sinTheta);
+        const newSinTheta = sinTheta - (alpha * sinTheta - beta * cosTheta);
+
+        cosTheta = newCosTheta;
+        sinTheta = newSinTheta;
+      }
+    }
+  }
+
+  return { re: resultRe, im: resultIm };
+}