508 lines
18 KiB
JavaScript
508 lines
18 KiB
JavaScript
import {
|
|
auditCovarianceEigenvalues,
|
|
computeComplexSampleCovariance,
|
|
selectWaxKailathSourceCount,
|
|
} from "/packages/acoustics-core/covariance.js";
|
|
import {
|
|
decomposeHermitian,
|
|
estimateHermitianEigenWorkload,
|
|
} from "/packages/acoustics-core/hermitian-eigen.js";
|
|
import {
|
|
computeMusicMap,
|
|
estimateMusicWorkload,
|
|
} from "/packages/acoustics-core/music.js";
|
|
|
|
function assertNear(actual, expected, tolerance, label) {
|
|
assert(
|
|
Math.abs(actual - expected) <= tolerance,
|
|
`${label}: expected ${expected} ± ${tolerance}, received ${actual}`,
|
|
);
|
|
}
|
|
|
|
function assertArrayNear(actual, expected, tolerance, label) {
|
|
assertEqual(actual.length, expected.length, `${label} length`);
|
|
for (let index = 0; index < actual.length; index += 1) {
|
|
assertNear(actual[index], expected[index], tolerance, `${label}[${index}]`);
|
|
}
|
|
}
|
|
|
|
function dftUnitary(order) {
|
|
const real = new Float64Array(order * order);
|
|
const imaginary = new Float64Array(order * order);
|
|
const scale = 1 / Math.sqrt(order);
|
|
for (let row = 0; row < order; row += 1) {
|
|
for (let column = 0; column < order; column += 1) {
|
|
const phase = (2 * Math.PI * row * column) / order;
|
|
const index = row * order + column;
|
|
real[index] = Math.cos(phase) * scale;
|
|
imaginary[index] = Math.sin(phase) * scale;
|
|
}
|
|
}
|
|
return { real, imaginary };
|
|
}
|
|
|
|
function hermitianFromEigenpairs(vectors, eigenvalues) {
|
|
const order = eigenvalues.length;
|
|
const real = new Float64Array(order * order);
|
|
const imaginary = new Float64Array(order * order);
|
|
for (let row = 0; row < order; row += 1) {
|
|
for (let column = row; column < order; column += 1) {
|
|
let valueReal = 0;
|
|
let valueImaginary = 0;
|
|
for (let eigenvector = 0; eigenvector < order; eigenvector += 1) {
|
|
const rowIndex = row * order + eigenvector;
|
|
const columnIndex = column * order + eigenvector;
|
|
const rowReal = vectors.real[rowIndex];
|
|
const rowImaginary = vectors.imaginary[rowIndex];
|
|
const columnReal = vectors.real[columnIndex];
|
|
const columnImaginary = vectors.imaginary[columnIndex];
|
|
valueReal += eigenvalues[eigenvector] *
|
|
(rowReal * columnReal + rowImaginary * columnImaginary);
|
|
valueImaginary += eigenvalues[eigenvector] *
|
|
(rowImaginary * columnReal - rowReal * columnImaginary);
|
|
}
|
|
const upper = row * order + column;
|
|
const lower = column * order + row;
|
|
real[upper] = valueReal;
|
|
imaginary[upper] = row === column ? 0 : valueImaginary;
|
|
real[lower] = valueReal;
|
|
imaginary[lower] = row === column ? 0 : -valueImaginary;
|
|
}
|
|
}
|
|
return { real, imaginary };
|
|
}
|
|
|
|
function projector(vectors, order, firstColumn, columnCount) {
|
|
const real = new Float64Array(order * order);
|
|
const imaginary = new Float64Array(order * order);
|
|
for (let row = 0; row < order; row += 1) {
|
|
for (let column = 0; column < order; column += 1) {
|
|
let valueReal = 0;
|
|
let valueImaginary = 0;
|
|
for (let vector = firstColumn; vector < firstColumn + columnCount; vector += 1) {
|
|
const rowIndex = row * order + vector;
|
|
const columnIndex = column * order + vector;
|
|
valueReal +=
|
|
vectors.real[rowIndex] * vectors.real[columnIndex] +
|
|
vectors.imaginary[rowIndex] * vectors.imaginary[columnIndex];
|
|
valueImaginary +=
|
|
vectors.imaginary[rowIndex] * vectors.real[columnIndex] -
|
|
vectors.real[rowIndex] * vectors.imaginary[columnIndex];
|
|
}
|
|
const index = row * order + column;
|
|
real[index] = valueReal;
|
|
imaginary[index] = valueImaginary;
|
|
}
|
|
}
|
|
return { real, imaginary };
|
|
}
|
|
|
|
function createAutomaticMusicFixture() {
|
|
const sampleRate = 3_200;
|
|
const fftSize = 32;
|
|
const snapshotCount = 64;
|
|
const frequencyHz = 500;
|
|
const speedOfSound = 343;
|
|
const source = [0.08, -0.04, 1];
|
|
const arrayDefinition = {
|
|
id: "automatic-music-audit-array",
|
|
microphones: [
|
|
{ id: "nw", position: [-0.1, 0.1, 0] },
|
|
{ id: "ne", position: [0.1, 0.1, 0] },
|
|
{ id: "sw", position: [-0.1, -0.1, 0] },
|
|
{ id: "se", position: [0.1, -0.1, 0] },
|
|
],
|
|
};
|
|
let state = 0x51f15e5d;
|
|
const random = () => {
|
|
state ^= state << 13;
|
|
state ^= state >>> 17;
|
|
state ^= state << 5;
|
|
return ((state >>> 0) / 0xffff_ffff) * 2 - 1;
|
|
};
|
|
const sourceRealBySnapshot = Float64Array.from(
|
|
{ length: snapshotCount },
|
|
random,
|
|
);
|
|
const sourceImaginaryBySnapshot = Float64Array.from(
|
|
{ length: snapshotCount },
|
|
random,
|
|
);
|
|
const channels = arrayDefinition.microphones.map((microphone) => {
|
|
const distance = Math.hypot(
|
|
source[0] - microphone.position[0],
|
|
source[1] - microphone.position[1],
|
|
source[2] - microphone.position[2],
|
|
);
|
|
const propagationPhase = -(2 * Math.PI * frequencyHz * distance) / speedOfSound;
|
|
const channel = new Float64Array(snapshotCount * fftSize);
|
|
for (let snapshot = 0; snapshot < snapshotCount; snapshot += 1) {
|
|
const sourceReal = sourceRealBySnapshot[snapshot];
|
|
const sourceImaginary = sourceImaginaryBySnapshot[snapshot];
|
|
for (let offset = 0; offset < fftSize; offset += 1) {
|
|
const phase = (2 * Math.PI * frequencyHz * offset) / sampleRate + propagationPhase;
|
|
channel[snapshot * fftSize + offset] =
|
|
4 * (sourceReal * Math.cos(phase) - sourceImaginary * Math.sin(phase)) /
|
|
distance +
|
|
0.08 * random();
|
|
}
|
|
}
|
|
return channel;
|
|
});
|
|
const experiment = {
|
|
stft: { fftSize, hopSize: fftSize, window: "hann" },
|
|
scan: {
|
|
z: source[2],
|
|
xMin: source[0],
|
|
xMax: source[0],
|
|
yMin: source[1],
|
|
yMax: source[1],
|
|
columns: 1,
|
|
rows: 1,
|
|
},
|
|
frequencyBands: [
|
|
{ id: "tone", minHz: 499, maxHz: 501, color: "#65f2bf" },
|
|
],
|
|
beamforming: {
|
|
algorithm: "music",
|
|
backend: "cpu",
|
|
normalization: "none",
|
|
compensateFrequencyResponse: true,
|
|
musicSourceCountMethod: "mdl",
|
|
},
|
|
};
|
|
return {
|
|
channels,
|
|
sampleRate,
|
|
arrayDefinition,
|
|
scene: { speedOfSound },
|
|
experiment,
|
|
};
|
|
}
|
|
|
|
test("complex Hermitian Jacobi resolves the [4, 1] two-by-two reference", () => {
|
|
const result = decomposeHermitian({
|
|
real: Float64Array.of(2, 1, 1, 3),
|
|
imaginary: Float64Array.of(0, 1, -1, 0),
|
|
});
|
|
assert(result.converged, "2x2 Jacobi convergence");
|
|
assertArrayNear(result.eigenvalues, [4, 1], 2e-14, "2x2 eigenvalues");
|
|
assertEqual(result.eigenvectors.layout, "row-major-columns", "vector layout");
|
|
assert(result.relativeResidual < 1e-14, "2x2 relative residual");
|
|
assert(result.orthogonalityError < 1e-14, "2x2 orthogonality");
|
|
assert(result.numericallyVerified, "2x2 numerical verification");
|
|
});
|
|
|
|
test("Hermitian Jacobi diagonalizes a constructed DFT4 matrix", () => {
|
|
const vectors = dftUnitary(4);
|
|
const matrix = hermitianFromEigenpairs(vectors, [9, 4, 2, 1]);
|
|
const result = decomposeHermitian(matrix);
|
|
assert(result.converged, "DFT4 convergence");
|
|
assertArrayNear(result.eigenvalues, [9, 4, 2, 1], 2e-12, "DFT4 eigenvalues");
|
|
assert(result.relativeResidual < 2e-13, "DFT4 relative residual");
|
|
assert(result.relativeOffDiagonal < result.relativeTolerance, "DFT4 off-diagonal bound");
|
|
});
|
|
|
|
test("repeated Hermitian eigenspaces preserve their invariant projector", () => {
|
|
const expectedVectors = dftUnitary(4);
|
|
const matrix = hermitianFromEigenpairs(expectedVectors, [5, 5, 1, 1]);
|
|
const result = decomposeHermitian(matrix);
|
|
assert(result.converged, "repeated-eigenspace convergence");
|
|
assertArrayNear(result.eigenvalues, [5, 5, 1, 1], 2e-12, "repeated eigenvalues");
|
|
const expectedProjector = projector(expectedVectors, 4, 0, 2);
|
|
const actualProjector = projector(result.eigenvectors, 4, 0, 2);
|
|
assertArrayNear(actualProjector.real, expectedProjector.real, 2e-12, "projector real");
|
|
assertArrayNear(
|
|
actualProjector.imaginary,
|
|
expectedProjector.imaginary,
|
|
2e-12,
|
|
"projector imaginary",
|
|
);
|
|
});
|
|
|
|
test("Hermitian eigenpairs are deterministic and phase-canonicalized", () => {
|
|
const matrix = hermitianFromEigenpairs(dftUnitary(4), [12, 7, 3, 0.5]);
|
|
const first = decomposeHermitian(matrix);
|
|
const second = decomposeHermitian(matrix);
|
|
assertArrayNear(second.eigenvalues, first.eigenvalues, 0, "deterministic eigenvalues");
|
|
assertArrayNear(second.eigenvectors.real, first.eigenvectors.real, 0, "deterministic real vectors");
|
|
assertArrayNear(
|
|
second.eigenvectors.imaginary,
|
|
first.eigenvectors.imaginary,
|
|
0,
|
|
"deterministic imaginary vectors",
|
|
);
|
|
for (let column = 0; column < first.order; column += 1) {
|
|
let maximumMagnitude = 0;
|
|
for (let row = 0; row < first.order; row += 1) {
|
|
const index = row * first.order + column;
|
|
const magnitude = Math.hypot(
|
|
first.eigenvectors.real[index],
|
|
first.eigenvectors.imaginary[index],
|
|
);
|
|
maximumMagnitude = Math.max(maximumMagnitude, magnitude);
|
|
}
|
|
let canonicalPivotFound = false;
|
|
for (let row = 0; row < first.order; row += 1) {
|
|
const index = row * first.order + column;
|
|
const magnitude = Math.hypot(
|
|
first.eigenvectors.real[index],
|
|
first.eigenvectors.imaginary[index],
|
|
);
|
|
if (
|
|
magnitude >= maximumMagnitude * (1 - 16 * Number.EPSILON) &&
|
|
first.eigenvectors.real[index] > 0 &&
|
|
first.eigenvectors.imaginary[index] === 0
|
|
) {
|
|
canonicalPivotFound = true;
|
|
}
|
|
}
|
|
assert(canonicalPivotFound, `column ${column} has a positive-real maximal pivot`);
|
|
}
|
|
});
|
|
|
|
test("Hermitian Jacobi reports non-convergence without returning it as success", () => {
|
|
const result = decomposeHermitian({
|
|
real: Float64Array.of(
|
|
4, 1, 0.5,
|
|
1, 3, 1.5,
|
|
0.5, 1.5, 2,
|
|
),
|
|
imaginary: Float64Array.of(
|
|
0, 0.5, -0.3,
|
|
-0.5, 0, 0.7,
|
|
0.3, -0.7, 0,
|
|
),
|
|
maximumSweeps: 1,
|
|
});
|
|
assert(!result.converged, "one-sweep result is explicitly unconverged");
|
|
assertEqual(result.sweeps, 1, "sweep limit is exact");
|
|
assert(!result.numericallyVerified, "partial eigenpairs are not verified");
|
|
assert(result.relativeResidual > result.relativeTolerance, "unconverged residual remains visible");
|
|
});
|
|
|
|
test("Hermitian workload rejects M=1024 before reading or allocating matrix data", () => {
|
|
const order = 1024;
|
|
const workload = estimateHermitianEigenWorkload({ order });
|
|
assert(workload.allocationBytes > 16 * 1024 * 1024, "M=1024 exceeds default bytes");
|
|
let numericRead = false;
|
|
const fakeMatrix = new Proxy(
|
|
{ length: order * order },
|
|
{
|
|
get(target, property) {
|
|
if (typeof property === "string" && /^\d+$/.test(property)) numericRead = true;
|
|
return target[property];
|
|
},
|
|
},
|
|
);
|
|
let error;
|
|
try {
|
|
decomposeHermitian({ real: fakeMatrix, imaginary: fakeMatrix, order });
|
|
} catch (caught) {
|
|
error = caught;
|
|
}
|
|
assert(error instanceof RangeError, "M=1024 is rejected by the workload cap");
|
|
assert(!numericRead, "workload cap fires before matrix values are read");
|
|
});
|
|
|
|
test("sample covariance uses compensated upper-triangle accumulation and exact conjugates", () => {
|
|
const cancellation = computeComplexSampleCovariance({
|
|
sensorCount: 2,
|
|
snapshotCount: 3,
|
|
layout: "snapshot-major",
|
|
real: Float64Array.of(1, 1e16, 1, 1, 1, -1e16),
|
|
imaginary: new Float64Array(6),
|
|
});
|
|
assertNear(cancellation.matrix.real[1], 1 / 3, 1e-16, "compensated cross covariance");
|
|
|
|
const complex = computeComplexSampleCovariance({
|
|
sensorCount: 2,
|
|
snapshotCount: 1,
|
|
real: Float64Array.of(1, 2),
|
|
imaginary: Float64Array.of(1, 3),
|
|
});
|
|
assertNear(complex.matrix.real[1], 5, 0, "upper real");
|
|
assertNear(complex.matrix.imaginary[1], -1, 0, "upper imaginary");
|
|
assertNear(complex.matrix.real[2], 5, 0, "lower real");
|
|
assertNear(complex.matrix.imaginary[2], 1, 0, "exact conjugate lower imaginary");
|
|
assertEqual(complex.matrix.imaginary[0], 0, "real diagonal 0");
|
|
assertEqual(complex.normalization, "1/T", "covariance normalization disclosed");
|
|
assertEqual(complex.centered, false, "covariance is not silently centered");
|
|
});
|
|
|
|
test("T<M rank-two covariance remains rank-deficient and source count is unavailable", () => {
|
|
const sensorCount = 4;
|
|
const snapshotCount = 2;
|
|
const covariance = computeComplexSampleCovariance({
|
|
sensorCount,
|
|
snapshotCount,
|
|
real: Float64Array.of(1, 2, 3, 4, 0, 1, 2, -1),
|
|
imaginary: Float64Array.of(0, 0, 0, 0, 1, -1, 0.5, 2),
|
|
});
|
|
const decomposition = decomposeHermitian({
|
|
...covariance.matrix,
|
|
order: sensorCount,
|
|
});
|
|
assert(decomposition.converged, "rank-two covariance decomposition");
|
|
const audit = auditCovarianceEigenvalues({
|
|
eigenvalues: decomposition.eigenvalues,
|
|
sensorCount,
|
|
matrixScale: covariance.matrixScale,
|
|
covarianceRoundoffBound: covariance.covarianceRoundoffBound,
|
|
eigensolverResidualFrobeniusNorm: decomposition.residualFrobeniusNorm,
|
|
});
|
|
assertEqual(
|
|
audit.roundoffDerivation,
|
|
"covariance-accumulation-plus-eigensolver-residual",
|
|
"PSD roundoff derivation",
|
|
);
|
|
assertNear(audit.roundoffClampedEigenvalues[2], 0, audit.roundoffTolerance, "rank eigenvalue 2");
|
|
assertNear(audit.roundoffClampedEigenvalues[3], 0, audit.roundoffTolerance, "rank eigenvalue 3");
|
|
const selected = selectWaxKailathSourceCount({
|
|
eigenvalues: audit.roundoffClampedEigenvalues,
|
|
snapshotCount,
|
|
independentSnapshots: true,
|
|
criterion: "mdl",
|
|
roundoffTolerance: audit.roundoffTolerance,
|
|
});
|
|
assert(!selected.available, "T<M MDL unavailable");
|
|
assertEqual(
|
|
selected.reason,
|
|
"snapshot-count-not-greater-than-sensors",
|
|
"T<M reason",
|
|
);
|
|
});
|
|
|
|
test("Wax-Kailath explicitly declines noiseless rank-deficient spectra", () => {
|
|
const result = selectWaxKailathSourceCount({
|
|
eigenvalues: Float64Array.of(5, 2, 0, 0),
|
|
snapshotCount: 100,
|
|
independentSnapshots: true,
|
|
criterion: "mdl",
|
|
roundoffTolerance: 1e-12,
|
|
});
|
|
assert(!result.available, "rank-deficient MDL unavailable");
|
|
assertEqual(
|
|
result.reason,
|
|
"noise-eigenvalue-not-strictly-above-roundoff",
|
|
"rank-deficient reason",
|
|
);
|
|
});
|
|
|
|
test("Wax-Kailath AIC and MDL select q=1 for [10,1,1,1]", () => {
|
|
for (const criterion of ["aic", "mdl"]) {
|
|
const result = selectWaxKailathSourceCount({
|
|
eigenvalues: Float64Array.of(10, 1, 1, 1),
|
|
snapshotCount: 100,
|
|
independentSnapshots: true,
|
|
criterion,
|
|
roundoffTolerance: 1e-14,
|
|
});
|
|
assert(result.available, `${criterion} available`);
|
|
assertEqual(result.sourceCount, 1, `${criterion} source count`);
|
|
}
|
|
});
|
|
|
|
test("Wax-Kailath AIC and MDL select q=0 for an identity covariance", () => {
|
|
for (const criterion of ["aic", "mdl"]) {
|
|
const result = selectWaxKailathSourceCount({
|
|
eigenvalues: Float64Array.of(1, 1, 1, 1),
|
|
snapshotCount: 100,
|
|
independentSnapshots: true,
|
|
criterion,
|
|
roundoffTolerance: 1e-14,
|
|
});
|
|
assert(result.available, `${criterion} available`);
|
|
assertEqual(result.sourceCount, 0, `${criterion} identity source count`);
|
|
}
|
|
});
|
|
|
|
test("PSD audit preserves tiny positives and rejects genuine negative eigenvalues", () => {
|
|
const tinyPositive = 1e-20;
|
|
const audit = auditCovarianceEigenvalues({
|
|
eigenvalues: Float64Array.of(10, 1, tinyPositive, -1e-14),
|
|
sensorCount: 4,
|
|
matrixScale: 10,
|
|
roundoffTolerance: 1e-12,
|
|
});
|
|
assertEqual(audit.rawEigenvalues[2], tinyPositive, "raw tiny positive");
|
|
assertEqual(
|
|
audit.roundoffClampedEigenvalues[2],
|
|
tinyPositive,
|
|
"tiny positive is not floored",
|
|
);
|
|
assertEqual(audit.roundoffClampedEigenvalues[3], 0, "roundoff negative clamps to zero");
|
|
let error;
|
|
try {
|
|
auditCovarianceEigenvalues({
|
|
eigenvalues: Float64Array.of(10, 1, 0.1, -1e-3),
|
|
sensorCount: 4,
|
|
roundoffTolerance: 1e-12,
|
|
});
|
|
} catch (caught) {
|
|
error = caught;
|
|
}
|
|
assert(error instanceof RangeError, "genuine negative eigenvalue is rejected");
|
|
});
|
|
|
|
test("MUSIC workload includes PSD-audit and automatic source-count typed arrays", () => {
|
|
const automatic = createAutomaticMusicFixture();
|
|
const automaticWorkload = estimateMusicWorkload(automatic);
|
|
const microphoneCount = automatic.arrayDefinition.microphones.length;
|
|
const fixed = {
|
|
...automatic,
|
|
experiment: structuredClone(automatic.experiment),
|
|
};
|
|
fixed.experiment.beamforming.musicSourceCountMethod = "fixed";
|
|
fixed.experiment.beamforming.musicSourceCount = 1;
|
|
const fixedWorkload = estimateMusicWorkload(fixed);
|
|
|
|
assertEqual(
|
|
automaticWorkload.covarianceAuditBytes,
|
|
2 * microphoneCount * Float64Array.BYTES_PER_ELEMENT,
|
|
"raw and roundoff-clamped PSD audit arrays",
|
|
);
|
|
assertEqual(fixedWorkload.sourceCountSelectionBytes, 0, "fixed K has no score array");
|
|
assertEqual(
|
|
automaticWorkload.sourceCountSelectionBytes,
|
|
microphoneCount * Float64Array.BYTES_PER_ELEMENT,
|
|
"automatic K criterion score array",
|
|
);
|
|
assertEqual(
|
|
automaticWorkload.perBinNumericalBytes - fixedWorkload.perBinNumericalBytes,
|
|
automaticWorkload.sourceCountSelectionBytes,
|
|
"automatic score bytes reach per-bin numerical peak",
|
|
);
|
|
assertEqual(
|
|
automaticWorkload.baseWorkingBytes - fixedWorkload.baseWorkingBytes,
|
|
automaticWorkload.sourceCountSelectionBytes,
|
|
"automatic score bytes reach total working peak",
|
|
);
|
|
});
|
|
|
|
test("MUSIC reports PSD derivation and labels MDL independence as an unverified model assumption", async () => {
|
|
const result = await computeMusicMap(createAutomaticMusicFixture());
|
|
const diagnostics = result.bands[0].diagnostics;
|
|
assertEqual(diagnostics.sourceCount.method, "mdl", "MDL selected model order");
|
|
assertEqual(diagnostics.sourceCount.value, 1, "one dominant source selected");
|
|
assertEqual(
|
|
diagnostics.covarianceAudit.roundoffDerivation,
|
|
"covariance-accumulation-plus-eigensolver-residual",
|
|
"PSD roundoff derivation is serialized from the audit API",
|
|
);
|
|
const independence =
|
|
diagnostics.sourceCount.assumptions.statisticalSnapshotIndependence;
|
|
assertEqual(independence.assumed, true, "statistical independence is an assumption");
|
|
assertEqual(
|
|
independence.evidence,
|
|
"non-overlapping-frames-only-not-statistical-proof",
|
|
"non-overlap is not mislabeled as statistical proof",
|
|
);
|
|
assertEqual(independence.modelSelectedByUser, true, "MDL model was selected by the user");
|
|
assertEqual(
|
|
diagnostics.sourceCount.assumptions.statisticalDataModel.verified,
|
|
false,
|
|
"Wax-Kailath data model is not claimed as verified",
|
|
);
|
|
});
|