Skip to content

Commit a49005c

Browse files
Add Discrete Fourier Transform (periodogram) view to XY output
Adds a frequency/period spectrum view to the XY output, toggled from a corner button. Each spectral component is drawn as a vertical bar (stem) on a real linear/logarithmic period axis, with a lin/log scale toggle. - frequency-utils.ts: radix-2 FFT + periodogram + period formatting - xy-output-component.tsx: frequency mode, bar rendering, axis toggles - output-components-style.css: toggle button styles Signed-off-by: Matthew Khouzam <matthew.khouzam@ericsson.com>
1 parent 2e36049 commit a49005c

4 files changed

Lines changed: 565 additions & 19 deletions

File tree

Lines changed: 107 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,107 @@
1+
import { nextPow2, fft, computePeriodogram, formatPeriod } from '../frequency-utils';
2+
3+
describe('nextPow2', () => {
4+
test('returns the next power of two', () => {
5+
expect(nextPow2(1)).toBe(1);
6+
expect(nextPow2(2)).toBe(2);
7+
expect(nextPow2(3)).toBe(4);
8+
expect(nextPow2(5)).toBe(8);
9+
expect(nextPow2(1000)).toBe(1024);
10+
});
11+
});
12+
13+
describe('fft', () => {
14+
test('DFT of a constant signal concentrates energy in the DC bin', () => {
15+
const re = [1, 1, 1, 1];
16+
const im = [0, 0, 0, 0];
17+
fft(re, im);
18+
// DC bin = sum of samples
19+
expect(re[0]).toBeCloseTo(4, 6);
20+
expect(im[0]).toBeCloseTo(0, 6);
21+
// All other bins ~ 0
22+
for (let k = 1; k < 4; k++) {
23+
expect(re[k]).toBeCloseTo(0, 6);
24+
expect(im[k]).toBeCloseTo(0, 6);
25+
}
26+
});
27+
28+
test('throws for non power-of-two lengths', () => {
29+
expect(() => fft([1, 2, 3], [0, 0, 0])).toThrow();
30+
});
31+
32+
test('matches naive DFT for a small signal', () => {
33+
const signal = [0, 1, 2, 3, 4, 5, 6, 7];
34+
const re = [...signal];
35+
const im = new Array(signal.length).fill(0);
36+
fft(re, im);
37+
38+
const n = signal.length;
39+
for (let k = 0; k < n; k++) {
40+
let expectedRe = 0;
41+
let expectedIm = 0;
42+
for (let t = 0; t < n; t++) {
43+
const angle = (-2 * Math.PI * k * t) / n;
44+
expectedRe += signal[t] * Math.cos(angle);
45+
expectedIm += signal[t] * Math.sin(angle);
46+
}
47+
expect(re[k]).toBeCloseTo(expectedRe, 6);
48+
expect(im[k]).toBeCloseTo(expectedIm, 6);
49+
}
50+
});
51+
});
52+
53+
describe('computePeriodogram', () => {
54+
test('returns empty for too-few samples or invalid span', () => {
55+
expect(computePeriodogram([1, 2, 3], 100)).toEqual([]);
56+
expect(computePeriodogram([1, 2, 3, 4], 0)).toEqual([]);
57+
expect(computePeriodogram([1, 2, 3, 4], -5)).toEqual([]);
58+
});
59+
60+
test('detects the dominant period of a pure sinusoid', () => {
61+
// Sample a sine wave with a known period over a known window.
62+
const n = 256;
63+
const timeSpan = 256; // ns; sample interval dt = 1 ns
64+
const cyclesInWindow = 8; // => period = timeSpan / cycles = 32 ns
65+
const values: number[] = [];
66+
for (let i = 0; i < n; i++) {
67+
values.push(Math.sin((2 * Math.PI * cyclesInWindow * i) / n));
68+
}
69+
70+
const spectrum = computePeriodogram(values, timeSpan);
71+
expect(spectrum.length).toBeGreaterThan(0);
72+
73+
// Results are sorted ascending by period.
74+
for (let i = 1; i < spectrum.length; i++) {
75+
expect(spectrum[i].period).toBeGreaterThanOrEqual(spectrum[i - 1].period);
76+
}
77+
78+
// The peak power should occur near the expected period (32 ns).
79+
const peak = spectrum.reduce((best, cur) => (cur.power > best.power ? cur : best), spectrum[0]);
80+
const expectedPeriod = timeSpan / cyclesInWindow;
81+
expect(peak.period).toBeCloseTo(expectedPeriod, 0);
82+
});
83+
84+
test('does not return periods longer than the observation window', () => {
85+
const n = 64;
86+
const timeSpan = 640;
87+
const values = Array.from({ length: n }, (_, i) => Math.sin((2 * Math.PI * 4 * i) / n));
88+
const spectrum = computePeriodogram(values, timeSpan);
89+
spectrum.forEach(pp => {
90+
expect(pp.period).toBeLessThanOrEqual(timeSpan);
91+
});
92+
});
93+
});
94+
95+
describe('formatPeriod', () => {
96+
test('formats using appropriate time units', () => {
97+
expect(formatPeriod(500)).toBe('500 ns');
98+
expect(formatPeriod(1500)).toBe('1.50 µs');
99+
expect(formatPeriod(2_000_000)).toBe('2.00 ms');
100+
expect(formatPeriod(3_000_000_000)).toBe('3.00 s');
101+
});
102+
103+
test('handles non-finite input', () => {
104+
expect(formatPeriod(NaN)).toBe('');
105+
expect(formatPeriod(Infinity)).toBe('');
106+
});
107+
});
Lines changed: 162 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,162 @@
1+
/**
2+
* Utilities to transform a uniformly-sampled time series into a
3+
* frequency/period spectrum (periodogram).
4+
*
5+
* The XY output samples each series uniformly across the displayed time
6+
* window, so we can estimate the power at each frequency using a Discrete
7+
* Fourier Transform and express the result as periods (period = 1 / frequency).
8+
*/
9+
10+
export interface PeriodPower {
11+
/** Period in the same time unit as the provided time span (e.g. nanoseconds). */
12+
period: number;
13+
/** Spectral power at that period. */
14+
power: number;
15+
}
16+
17+
/**
18+
* Returns the smallest power of two that is greater than or equal to `n`.
19+
*/
20+
export function nextPow2(n: number): number {
21+
let p = 1;
22+
while (p < n) {
23+
p <<= 1;
24+
}
25+
return p;
26+
}
27+
28+
/**
29+
* In-place iterative radix-2 Cooley-Tukey FFT.
30+
* `re` and `im` must have the same length, which must be a power of two.
31+
* On return, `re`/`im` hold the real/imaginary parts of the transform.
32+
*/
33+
export function fft(re: number[], im: number[]): void {
34+
const n = re.length;
35+
if (n <= 1) {
36+
return;
37+
}
38+
if ((n & (n - 1)) !== 0) {
39+
throw new Error('fft: input length must be a power of two');
40+
}
41+
42+
// Bit-reversal permutation.
43+
for (let i = 1, j = 0; i < n; i++) {
44+
let bit = n >> 1;
45+
for (; j & bit; bit >>= 1) {
46+
j ^= bit;
47+
}
48+
j ^= bit;
49+
if (i < j) {
50+
const tmpRe = re[i];
51+
re[i] = re[j];
52+
re[j] = tmpRe;
53+
const tmpIm = im[i];
54+
im[i] = im[j];
55+
im[j] = tmpIm;
56+
}
57+
}
58+
59+
// Butterfly stages.
60+
for (let len = 2; len <= n; len <<= 1) {
61+
const ang = (-2 * Math.PI) / len;
62+
const wRe = Math.cos(ang);
63+
const wIm = Math.sin(ang);
64+
for (let i = 0; i < n; i += len) {
65+
let curRe = 1;
66+
let curIm = 0;
67+
const half = len >> 1;
68+
for (let k = 0; k < half; k++) {
69+
const evenRe = re[i + k];
70+
const evenIm = im[i + k];
71+
const oddReRaw = re[i + k + half];
72+
const oddImRaw = im[i + k + half];
73+
const oddRe = oddReRaw * curRe - oddImRaw * curIm;
74+
const oddIm = oddReRaw * curIm + oddImRaw * curRe;
75+
76+
re[i + k] = evenRe + oddRe;
77+
im[i + k] = evenIm + oddIm;
78+
re[i + k + half] = evenRe - oddRe;
79+
im[i + k + half] = evenIm - oddIm;
80+
81+
const nextRe = curRe * wRe - curIm * wIm;
82+
curIm = curRe * wIm + curIm * wRe;
83+
curRe = nextRe;
84+
}
85+
}
86+
}
87+
}
88+
89+
/**
90+
* Computes a periodogram from a uniformly-sampled series.
91+
*
92+
* @param values The sampled y-values (uniformly spaced in time).
93+
* @param timeSpan The total duration covered by `values` (end - start), in any
94+
* time unit. The returned periods use the same unit.
95+
* @returns Array of {period, power}, sorted by ascending period, limited
96+
* to periods that fit inside the observation window and are at
97+
* least the Nyquist period (2 * sample interval).
98+
*/
99+
export function computePeriodogram(values: number[], timeSpan: number): PeriodPower[] {
100+
const n = values.length;
101+
if (n < 4 || !Number.isFinite(timeSpan) || timeSpan <= 0) {
102+
return [];
103+
}
104+
105+
// Detrend by removing the mean (DC component) so the constant offset does
106+
// not dominate the spectrum.
107+
let mean = 0;
108+
for (let i = 0; i < n; i++) {
109+
const v = Number(values[i]);
110+
mean += Number.isFinite(v) ? v : 0;
111+
}
112+
mean /= n;
113+
114+
const m = nextPow2(n);
115+
const re = new Array<number>(m).fill(0);
116+
const im = new Array<number>(m).fill(0);
117+
for (let i = 0; i < n; i++) {
118+
const v = Number(values[i]);
119+
re[i] = (Number.isFinite(v) ? v : 0) - mean;
120+
}
121+
122+
fft(re, im);
123+
124+
const dt = timeSpan / n; // sample interval
125+
const result: PeriodPower[] = [];
126+
for (let k = 1; k <= m >> 1; k++) {
127+
// With zero-padding to length m, bin k maps to frequency k / (m * dt),
128+
// hence period = (m * dt) / k.
129+
const period = (m * dt) / k;
130+
if (period > timeSpan) {
131+
// Longer than the observation window: not physically meaningful.
132+
continue;
133+
}
134+
const power = (re[k] * re[k] + im[k] * im[k]) / n;
135+
result.push({ period, power });
136+
}
137+
138+
// Ascending by period (shortest period / highest frequency first).
139+
result.sort((a, b) => a.period - b.period);
140+
return result;
141+
}
142+
143+
/**
144+
* Formats a period given in nanoseconds into a compact human-readable string
145+
* with an appropriate time unit.
146+
*/
147+
export function formatPeriod(ns: number): string {
148+
if (!Number.isFinite(ns)) {
149+
return '';
150+
}
151+
const abs = Math.abs(ns);
152+
if (abs >= 1e9) {
153+
return `${(ns / 1e9).toFixed(2)} s`;
154+
}
155+
if (abs >= 1e6) {
156+
return `${(ns / 1e6).toFixed(2)} ms`;
157+
}
158+
if (abs >= 1e3) {
159+
return `${(ns / 1e3).toFixed(2)} µs`;
160+
}
161+
return `${Math.round(ns)} ns`;
162+
}

0 commit comments

Comments
 (0)