From f07810e6da42e60c6f3d727fbac1527a8c19c7d8 Mon Sep 17 00:00:00 2001 From: thomaspmurphy Date: Mon, 15 Jun 2026 15:17:28 +0100 Subject: [PATCH 1/3] feat: add Livebook guides for convolution, windows, filters, peak finding, sound synthesis, and financial signals New guides covering: convolution theorem and cross-correlation, window functions and spectral leakage, FIR filtering, peak finding with argrelmax/ argrelmin, sound synthesis with waveforms and STFT, and signal processing for economic data. Also fixes the spectrogram guide for the new positional window-size API and updates firwin calls to require a list cutoff. --- guides/convolution.livemd | 320 +++++++++++++++++++++++ guides/fft.livemd | 4 +- guides/filtering.livemd | 4 +- guides/financial_signals.livemd | 267 +++++++++++++++++++ guides/peak_finding.livemd | 242 +++++++++++++++++ guides/sound_synthesis.livemd | 447 ++++++++++++++++++++++++++++++++ guides/spectrogram.livemd | 20 +- guides/windows.livemd | 255 ++++++++++++++++++ 8 files changed, 1541 insertions(+), 18 deletions(-) create mode 100644 guides/convolution.livemd create mode 100644 guides/financial_signals.livemd create mode 100644 guides/peak_finding.livemd create mode 100644 guides/sound_synthesis.livemd create mode 100644 guides/windows.livemd diff --git a/guides/convolution.livemd b/guides/convolution.livemd new file mode 100644 index 0000000..8fd803b --- /dev/null +++ b/guides/convolution.livemd @@ -0,0 +1,320 @@ +# The Convolution Theorem + +```elixir +Mix.install([ + {:nx_signal, path: __DIR__ |> Path.join("..") |> Path.expand()}, + {:kino_vega_lite, "~> 0.1"}, + {:tucan, "~> 0.5"} +]) +``` + +## The most important theorem in DSP + +Convolution in the time domain is equivalent to pointwise multiplication in +the frequency domain, and vice versa: + +$$ +x[n] * h[n] \;\longleftrightarrow\; X(f) \cdot H(f) +$$ + +$$ +x[n] \cdot h[n] \;\longleftrightarrow\; X(f) * H(f) +$$ + +This single fact underpins almost everything in signal processing: + +* **Filtering**: applying a filter $h$ to signal $x$ is just multiplication + in the frequency domain, then transform back. +* **Fast convolution**: direct convolution costs $O(N^2)$; via FFT it costs + $O(N \log N)$. +* **Spectral leakage**: multiplying by a window in time convolves your + spectrum with the window's own Fourier transform (the dual). + +## What is convolution? + +The cleanest way to understand convolution is through the **impulse response**: +a linear time-invariant (LTI) system is completely described by what it does to +a unit impulse $\delta[n]$. Call that response $h[n]$. For any input $x[n]$, +the output is: + +$$ +y[n] = \sum_{k} x[k] \, h[n - k] = (x * h)[n] +$$ + +Every sample of $x$ contributes a scaled, delayed copy of $h$; the output is +their sum. + +```elixir +# A simple 5-tap box (moving average) filter +h = Nx.broadcast(1.0 / 5, {5}) |> Nx.as_type(:f32) + +# An input with impulses at positions 3, 10, and 16 (amplitudes 1, 0.75, 0.5) +x = + Nx.broadcast(0.0, {30}) + |> Nx.indexed_put(Nx.tensor([[3], [10], [16]]), Nx.tensor([1.0, 0.75, 0.5])) + +y = NxSignal.Convolution.convolve(x, h, mode: :full) + +samples_x = Enum.map(0..29, & &1) +samples_y = Enum.map(0..(Nx.size(y) - 1), & &1) + +impulse_data = + Enum.zip(samples_x, Nx.to_flat_list(x)) + |> Enum.map(fn {i, v} -> %{sample: i, value: v, signal: "Input x[n]"} end) + +filter_data = + Enum.zip(Enum.map(0..4, & &1), Nx.to_flat_list(h)) + |> Enum.map(fn {i, v} -> %{sample: i, value: v, signal: "Filter h[n]"} end) + +output_data = + Enum.zip(samples_y, Nx.to_flat_list(y)) + |> Enum.map(fn {i, v} -> %{sample: i, value: v, signal: "Output y[n] = x * h"} end) + +all_data = impulse_data ++ filter_data ++ output_data + +Tucan.lineplot(all_data, "sample", "value") +|> Tucan.color_by("signal") +|> Tucan.facet_by(:row, "signal") +|> Tucan.Axes.set_x_title("Sample") +|> Tucan.Axes.set_y_title("Amplitude") +|> Tucan.set_title("Convolution as scaled, delayed copies of the impulse response") +|> Tucan.set_width(640) +|> Tucan.set_height(80) +``` + +Each impulse in $x$ produces a copy of $h$ at that position, scaled by the +impulse's amplitude. The output $y$ is their superposition. + +## The theorem in action + +The theorem states that two routes exist to compute $y = x * h$: + +**Route A - time domain:** direct convolution (flip-and-slide). + +**Route B - frequency domain:** $Y = \text{IFFT}(X \cdot H)$. + +Both must give the same answer. Let's verify with a real filter. + +```elixir +fs = 8_000 +dur = 0.5 +n = trunc(fs * dur) +t = Nx.linspace(0, dur, n: n, endpoint: false, type: :f32) + +# Signal: 440 Hz tone + 1 800 Hz tone + Gaussian noise +key = Nx.Random.key(42) +{noise, _key} = Nx.Random.normal(key, shape: {n}, type: :f32) + +signal = + Nx.sin(Nx.multiply(t, 2 * :math.pi() * 440)) + |> Nx.add(Nx.multiply(0.6, Nx.sin(Nx.multiply(t, 2 * :math.pi() * 1800)))) + |> Nx.add(Nx.multiply(0.2, noise)) + +# 51-tap low-pass FIR filter, cutoff at 1 000 Hz +cutoff_norm = 1000.0 / (fs / 2.0) +h_fir = NxSignal.Filters.firwin(51, [cutoff_norm]) + +# Route A: direct time-domain convolution +route_a = NxSignal.Convolution.convolve(signal, h_fir, mode: :same, method: :direct) + +# Route B: manual frequency-domain multiplication +fft_len = n + Nx.size(h_fir) - 1 +# Next power of two for efficiency +fft_len_padded = Integer.pow(2, ceil(:math.log2(fft_len))) + +x_fft = Nx.fft(Nx.as_type(signal, {:c, 64}), length: fft_len_padded) +h_fft = Nx.fft(Nx.as_type(h_fir, {:c, 64}), length: fft_len_padded) +y_full = Nx.ifft(Nx.multiply(x_fft, h_fft)) + +# Trim to :same length (centre of the full convolution) +trim = div(Nx.size(h_fir), 2) +route_b = Nx.real(y_full[trim..(trim + n - 1)]) + +# Confirm both routes agree +max_diff = + Nx.subtract(route_a, route_b) + |> Nx.abs() + |> Nx.reduce_max() + +IO.inspect(Nx.to_number(max_diff), label: "max |Route A − Route B|") +``` + +```elixir +# Visualise: input spectrum, filter response, and output spectrum side by side +half = div(n, 2) +freqs_hz = Nx.linspace(0, fs / 2.0, n: half, type: :f32) |> Nx.to_flat_list() + +to_amp_list = fn sig -> + sig + |> Nx.as_type({:c, 64}) + |> Nx.fft() + |> Nx.abs() + |> Nx.slice([0], [half]) + |> Nx.to_flat_list() +end + +h_padded = Nx.pad(Nx.as_type(h_fir, {:c, 64}), Nx.as_type(0, {:c, 64}), [{0, n - Nx.size(h_fir), 0}]) +h_spectrum = h_padded |> Nx.fft() |> Nx.abs() |> Nx.slice([0], [half]) |> Nx.to_flat_list() + +spectrum_data = + Enum.zip(freqs_hz, to_amp_list.(signal)) + |> Enum.map(fn {f, a} -> %{frequency: f, amplitude: a, panel: "Input |X(f)|"} end) + |> Kernel.++( + Enum.zip(freqs_hz, h_spectrum) + |> Enum.map(fn {f, a} -> %{frequency: f, amplitude: a, panel: "Filter |H(f)|"} end) + ) + |> Kernel.++( + Enum.zip(freqs_hz, to_amp_list.(route_a)) + |> Enum.map(fn {f, a} -> %{frequency: f, amplitude: a, panel: "Output |X(f)·H(f)|"} end) + ) + +Tucan.lineplot(spectrum_data, "frequency", "amplitude") +|> Tucan.facet_by(:row, "panel") +|> Tucan.Axes.set_x_title("Frequency (Hz)") +|> Tucan.Axes.set_y_title("Amplitude") +|> Tucan.set_title("The convolution theorem: three views of filtering") +|> Tucan.set_width(640) +|> Tucan.set_height(100) +``` + +The 1 800 Hz component visible in the input spectrum is absent from the +output. This is because the filter zeroed it in the frequency domain. + +## Cross-correlation and delay estimation + +Cross-correlation is like convolution but without the time-reversal of $h$: + +$$ +(x \star h)[n] = \sum_k x[k] \, h[k - n] +$$ + +It measures how similar $x$ and $h$ are as a function of **lag** $n$. The +peak of the cross-correlation identifies the delay between two signals. This +technique is used in sonar, seismology, and audio alignment. + +```elixir +fs_corr = 8_000 +burst_len = trunc(0.04 * fs_corr) # 40 ms reference burst +total_len = trunc(0.4 * fs_corr) # 400 ms total window +true_delay = 80 # samples ≈ 10 ms + +t_burst = Nx.linspace(0, burst_len / fs_corr, n: burst_len, endpoint: false, type: :f32) +reference = Nx.sin(Nx.multiply(t_burst, 2 * :math.pi() * 600)) + +# Received signal: noise + attenuated, delayed copy of the reference +key2 = Nx.Random.key(7) +{noise2, _} = Nx.Random.normal(key2, shape: {total_len}, type: :f32) +noise2 = Nx.multiply(noise2, 0.4) + +rx = + Nx.pad(reference, 0.0, [{true_delay, total_len - burst_len - true_delay, 0}]) + |> Nx.multiply(0.7) + |> Nx.add(noise2) + +# Cross-correlate received signal with reference +corr = NxSignal.Convolution.correlate(rx, reference, mode: :full) +corr_amps = Nx.abs(corr) + +# Lag axis: negative lags first, zero at index burst_len - 1 +n_corr = Nx.size(corr_amps) +lags = Enum.map(0..(n_corr - 1), fn i -> i - (burst_len - 1) end) + +{_val, peak_idx} = Nx.top_k(corr_amps, k: 1) +detected_delay = Nx.to_number(peak_idx[0]) - (burst_len - 1) + +IO.puts("True delay: #{true_delay} samples") +IO.puts("Detected delay: #{detected_delay} samples") +``` + +```elixir +corr_data = + Enum.zip(lags, Nx.to_flat_list(corr_amps)) + |> Enum.filter(fn {lag, _} -> lag >= -20 and lag <= 200 end) + |> Enum.map(fn {lag, a} -> %{lag: lag, correlation: a} end) + +peak_data = [%{lag: detected_delay, correlation: Nx.to_number(corr_amps[detected_delay + burst_len - 1])}] + +VegaLite.new(width: 680, height: 220, title: "Cross-correlation: peak at delay = #{detected_delay} samples") +|> VegaLite.layers([ + VegaLite.new() + |> VegaLite.data_from_values(corr_data) + |> VegaLite.mark(:line, color: "steelblue") + |> VegaLite.encode_field(:x, "lag", type: :quantitative, title: "Lag (samples)") + |> VegaLite.encode_field(:y, "correlation", type: :quantitative, title: "|Correlation|"), + VegaLite.new() + |> VegaLite.data_from_values(peak_data) + |> VegaLite.mark(:point, color: "tomato", size: 120, filled: true) + |> VegaLite.encode_field(:x, "lag", type: :quantitative) + |> VegaLite.encode_field(:y, "correlation", type: :quantitative) +]) +``` + +The correlation peak (red dot) is at exactly the true delay, even though the +echo is buried in noise and invisible to the eye in the raw received signal. + +## Windowing and spectral leakage + +The second form of the theorem states that **multiplication in time equals +convolution in frequency**. This is why window functions matter. + +When we multiply a signal by a rectangular window (i.e. we simply observe +$N$ samples and discard the rest), we convolve its spectrum with a $\text{sinc}$ +function, spreading energy from each spectral line into neighbouring bins. +A smooth window has a more compact Fourier transform, so the spectral smearing +is less severe. + +```elixir +fs_win = 200 +n_win = 200 +f_tone2 = 10.5 + +t_win = Nx.linspace(0, n_win / fs_win, n: n_win, endpoint: false, type: :f32) +x_win = Nx.sin(Nx.multiply(t_win, 2 * :math.pi() * f_tone2)) + +rect_w = NxSignal.Windows.rectangular(n_win, type: :f32) +hann_w = NxSignal.Windows.hann(n_win, is_periodic: false) + +half_win = div(n_win, 2) +freqs_win = NxSignal.fft_frequencies(fs_win, fft_length: n_win)[0..half_win] |> Nx.to_flat_list() + +dual_data = + [{rect_w, "Rectangular (multiply by 1)"}, {hann_w, "Hann (smooth taper)"}] + |> Enum.flat_map(fn {w, name} -> + amps = + Nx.multiply(x_win, w) + |> Nx.as_type({:c, 64}) + |> Nx.fft() + |> Nx.abs() + |> Nx.slice([0], [half_win + 1]) + peak = Nx.reduce_max(amps) + db = + Nx.divide(amps, peak) + |> Nx.log10() + |> Nx.multiply(20) + |> Nx.max(-100.0) + |> Nx.to_flat_list() + Enum.zip(freqs_win, db) + |> Enum.map(fn {f, d} -> %{frequency: f, amplitude_db: d, window: name} end) + end) + +Tucan.lineplot(dual_data, "frequency", "amplitude_db") +|> Tucan.color_by("window") +|> VegaLite.encode_field(:y, "amplitude_db", type: :quantitative, title: "Amplitude (dB)", scale: [domain: [-100, 5]]) +|> Tucan.Axes.set_x_title("Frequency (Hz)") +|> Tucan.set_title("Dual theorem: multiplication in time = convolution in frequency (spectral leakage)") +|> Tucan.set_width(680) +|> Tucan.set_height(260) +``` + +## Summary + +| Function | Method | When to use | +| --------------------------------------------- | -------------- | ----------------------------------- | +| `Convolution.convolve(x, h, method: :direct)` | $O(N \cdot M)$ | Short kernels ($M \lesssim 50$) | +| `Convolution.convolve(x, h, method: :fft)` | $O(N \log N)$ | Long kernels or large signals | +| `Convolution.fftconvolve(x, h)` | $O(N \log N)$ | Convenience alias for FFT method | +| `Convolution.correlate(x, h)` | Direct or FFT | Delay estimation, matched filtering | + +The convolution theorem is the reason `method: :fft` exists: instead of +$O(N \cdot M)$ multiply-accumulates, you pay three FFTs at $O(N \log N)$ each. +For any kernel longer than roughly 50 taps the FFT method wins. diff --git a/guides/fft.livemd b/guides/fft.livemd index 561dbab..9f50d77 100644 --- a/guides/fft.livemd +++ b/guides/fft.livemd @@ -2,7 +2,7 @@ ```elixir Mix.install([ - {:nx_signal, "~> 0.2"}, + {:nx_signal, "~> 0.3"}, {:vega_lite, "~> 0.1"}, {:kino_vega_lite, "~> 0.1"} ]) @@ -120,7 +120,7 @@ We can confirm this visual inspection with a peek into our data. We use `Nx.top_ {values, indices} = Nx.top_k(amplitudes, k: 5) { - values, + values, frequencies[indices] } ``` diff --git a/guides/filtering.livemd b/guides/filtering.livemd index 0635ca0..b894f43 100644 --- a/guides/filtering.livemd +++ b/guides/filtering.livemd @@ -2,7 +2,7 @@ ```elixir Mix.install([ - {:nx_signal, "~> 0.2"}, + {:nx_signal, "~> 0.3"}, {:vega_lite, "~> 0.1"}, {:kino_vega_lite, "~> 0.1"} ]) @@ -55,7 +55,7 @@ VegaLite.new(width: 600, height: 400, title: "Signal sample") fc = 600 # Design the FIR filter coefficients using the window method -h = NxSignal.Filters.firwin(window_length, fc, sampling_rate: fs, window: :hann) +h = NxSignal.Filters.firwin(window_length, [fc], sampling_rate: fs, window: :hann) # Separate periodic Hann window for STFT analysis stft_window = NxSignal.Windows.hann(window_length) diff --git a/guides/financial_signals.livemd b/guides/financial_signals.livemd new file mode 100644 index 0000000..569087f --- /dev/null +++ b/guides/financial_signals.livemd @@ -0,0 +1,267 @@ +# Signal Processing for Economic Data + +```elixir +Mix.install([ + {:nx_signal, path: __DIR__ |> Path.join("..") |> Path.expand()}, + {:kino, "~> 0.13"}, + {:kino_vega_lite, "~> 0.1"}, + {:tucan, "~> 0.5"} +]) +``` + +## Introduction + +Financial time series are signals: sequences of numbers indexed by time. Every tool from digital signal processing applies to them directly. This notebook uses a constructed S&P 500-like daily price series to demonstrate: + +1. The **FFT** for revealing hidden cycles in market data. +2. A **FIR low-pass filter** for smoothing out noise, reproducing the familiar "moving average" used by traders. +3. The **Zoom FFT** for zooming in on the annual-cycle frequency band with much finer resolution than the full FFT allows. + +## Generating a price series + +We construct a synthetic S&P 500-like price series from first principles: +a linear log-price trend, a 252-day (annual) seasonal cycle, and a Gaussian +random walk for daily noise. Building the series this way lets us control +exactly which frequency components are present, so the FFT results below +are easy to interpret. + +```elixir +n_days = 3780 # ~15 years of trading days + +t = Nx.iota({n_days}, type: :f32) +{noise, _} = Nx.Random.normal(Nx.Random.key(42), shape: {n_days}, type: :f32) + +# Log price = linear growth trend + annual cycle + random walk +trend_per_day = :math.log(4.8) / n_days # roughly 4.8× over 15 years +annual_cycle = Nx.multiply(0.04, Nx.sin(Nx.multiply(t, 2 * :math.pi() / 252.0))) + +random_walk = + noise + |> Nx.multiply(0.01) + |> Nx.to_flat_list() + |> Enum.scan(0.0, &(&1 + &2)) + |> Nx.tensor(type: :f32) + +log_prices = + Nx.add(Nx.multiply(t, trend_per_day), annual_cycle) + |> Nx.add(random_walk) + |> Nx.add(:math.log(1000.0)) + +closes = Nx.exp(log_prices) + +# Date axis: Mon-Fri trading days starting 2010-01-04 +dates = + Stream.iterate(~D[2010-01-04], &Date.add(&1, 1)) + |> Stream.filter(fn d -> Date.day_of_week(d) not in [6, 7] end) + |> Enum.take(n_days) + |> Enum.map(&Date.to_string/1) + +IO.puts("#{n_days} synthetic trading days") +IO.puts("From #{List.first(dates)} to #{List.last(dates)}") +``` + +## The price series and log returns + +Raw closing prices form a non-stationary signal: the mean drifts upward over +time. For frequency analysis we instead work with **log returns** +$r_t = \log(P_t / P_{t-1})$, which are approximately stationary. + +```elixir +n = Nx.size(closes) + +# Log returns: r[t] = log(P[t] / P[t-1]) +log_returns = + Nx.log(closes[1..(n - 1)]) |> Nx.subtract(Nx.log(closes[0..(n - 2)])) + +price_data = + Enum.zip(dates, Nx.to_flat_list(closes)) + |> Enum.map(fn {d, c} -> %{date: d, close: c} end) + +Tucan.lineplot(price_data, "date", "close") +|> VegaLite.encode_field(:x, "date", type: :temporal, title: "Date") +|> Tucan.set_title("S&P 500 Daily Close") +|> Tucan.set_width(680) +|> Tucan.set_height(240) +``` + +```elixir +return_data = + Enum.zip(Enum.drop(dates, 1), Nx.to_flat_list(log_returns)) + |> Enum.map(fn {d, r} -> %{date: d, return: r} end) + +Tucan.lineplot(return_data, "date", "return") +|> VegaLite.encode_field(:x, "date", type: :temporal, title: "Date") +|> Tucan.set_title("S&P 500 Daily Log Returns") +|> Tucan.set_width(680) +|> Tucan.set_height(200) +``` + +## Frequency analysis with the FFT + +`Nx.fft` decomposes the log-return signal into its constituent frequencies. +The x-axis below is expressed in **trading days per cycle**: a period of 252 +corresponds to the annual cycle, 63 to the quarterly cycle, 21 to the monthly +cycle. + +```elixir +m = Nx.size(log_returns) + +# Subtract mean to suppress the DC (zero-frequency) component +returns_ac = Nx.subtract(log_returns, Nx.mean(log_returns)) + +# Apply a Hann window to reduce spectral leakage +window = NxSignal.Windows.hann(m, is_periodic: true) +windowed = Nx.multiply(returns_ac, window) + +spectrum = Nx.fft(Nx.as_type(windowed, {:c, 64})) +power = spectrum |> Nx.abs() |> Nx.pow(2) + +# Frequency axis: cycles per trading day +# Only the positive frequencies (first half) are meaningful for real signals +half = div(m, 2) +freq = Nx.linspace(0, 0.5, n: half, endpoint: false, type: :f32) + +# Convert to period (trading days per cycle); skip DC bin (index 0) +period_days = Nx.divide(1.0, freq[1..(half - 1)]) +power_half = power[1..(half - 1)] + +# Clip period axis to a readable range: 5 to 500 trading days +{periods_roi, powers_roi} = + Enum.zip(Nx.to_flat_list(period_days), Nx.to_flat_list(power_half)) + |> Enum.filter(fn {p, _} -> p >= 5 and p <= 500 end) + |> Enum.unzip() + +spectrum_data = + Enum.zip(periods_roi, powers_roi) + |> Enum.map(fn {p, pw} -> %{period: p, power: pw} end) + +VegaLite.new(width: 680, height: 260, title: "Power Spectrum of S&P 500 Log Returns") +|> VegaLite.data_from_values(spectrum_data) +|> VegaLite.mark(:bar, tooltip: true) +|> VegaLite.encode_field(:x, "period", + type: :quantitative, + title: "Period (trading days)", + scale: [type: "log"] + ) +|> VegaLite.encode_field(:y, "power", + type: :quantitative, + title: "Power", + scale: [type: "log"] + ) +``` + +Vertical lines around period ≈ 252 (annual) and 63 (quarterly) are expected +artefacts of well-known seasonal patterns in equity markets. + +## Smoothing with a low-pass FIR filter + +A **moving average** is, mathematically speaking, a low-pass FIR filter: it suppresses +high-frequency fluctuations (day-to-day noise) and preserves the underlying +trend. + +`NxSignal.Filters.firwin/3` designs an FIR filter using the window method. +Below we design a linear-phase filter whose cutoff is set to 1/63 cycles per +trading day. This is equivalent to a 63-day (roughly quarterly) moving average, but +with a sharper roll-off than the simple boxcar. + +```elixir +alias NxSignal.Convolution + +# Normalised cutoff: (1/63 cycles/day) / (Nyquist = 0.5 cycles/day) = 2/63 +cutoff_norm = 2 / 63 + +# 127-tap filter (odd tap count → linear phase Type I) +coeffs = NxSignal.Filters.firwin(127, [cutoff_norm], window: :hamming) + +smoothed = + Convolution.fftconvolve(closes, coeffs, mode: :same) + +# Trim the edge transient introduced by the filter delay (63 samples each side) +trim = 63 + +trim_dates = Enum.slice(dates, trim, n - 2 * trim) +trim_close = Nx.to_flat_list(closes[trim..(n - trim - 1)]) +trim_smooth = Nx.to_flat_list(smoothed[trim..(n - trim - 1)]) + +# Combine raw and smoothed into a single long-form dataset +chart_data_long = + Enum.zip([trim_dates, trim_close, trim_smooth]) + |> Enum.flat_map(fn {d, raw, filt} -> + [%{date: d, close: raw, series: "Raw"}, %{date: d, close: filt, series: "FIR smoothed (63-day)"}] + end) + +Tucan.lineplot(chart_data_long, "date", "close") +|> Tucan.color_by("series") +|> VegaLite.encode_field(:x, "date", type: :temporal, title: "Date") +|> Tucan.set_title("S&P 500: Raw vs FIR Smoothed (63-day)") +|> Tucan.set_width(680) +|> Tucan.set_height(260) +``` + +The smoothed line tracks the broad trend while the raw line shows the daily +volatility that the filter removes. This is exactly the "50-day moving average" +chart familiar to technical analysts, produced here with a principled FIR +design rather than a simple average. + +## Zoom FFT: the annual cycle in detail + +The standard FFT spaces its bins uniformly across $[0, 0.5]$ cycles/day. +Near the annual frequency ($1/252 \approx 0.00397$ cycles/day) the bin +spacing is only $1/N$ cycles/day. With $N \approx 3\,700$ daily samples that +is about $2.7 \times 10^{-4}$ cycles/day, one bin per 8 calendar months. +That makes it hard to see whether the annual peak splits into sub-peaks or +drifts between years. + +`NxSignal.zoom_fft/4` concentrates all $M$ output bins on a narrow band, giving +resolution proportional to $M$ rather than $N$. + +```elixir +m_zoom = 1024 + +# Normalised frequency band around the annual cycle +# Annual: 1/252 ≈ 0.00397 cycles/day; examine ± 50% either side +f1_norm = (1 / 400) / 0.5 # lower bound (400-day period) +f2_norm = (1 / 150) / 0.5 # upper bound (150-day period) + +zoom = NxSignal.zoom_fft(returns_ac, f1_norm, f2_norm, m: m_zoom) +zoom_power = zoom |> Nx.abs() |> Nx.pow(2) |> Nx.to_flat_list() + +# Frequency axis for the zoom output (in cycles/day) +f_lo = 1 / 400 / 1.0 +f_hi = 1 / 150 / 1.0 + +zoom_periods = + Enum.map(0..(m_zoom - 1), fn k -> + freq = f_lo + k * (f_hi - f_lo) / m_zoom + 1 / freq + end) + +zoom_data = + Enum.zip(zoom_periods, zoom_power) + |> Enum.map(fn {p, pw} -> %{period: p, power: pw} end) + +Tucan.lineplot(zoom_data, "period", "power", tooltip: true) +|> VegaLite.encode_field(:x, "period", + type: :quantitative, + title: "Period (trading days)", + scale: [domain: [150, 400]] +) +|> Tucan.Axes.set_y_title("Power") +|> Tucan.set_title("Zoom FFT: Annual Cycle Band (150-400 trading-day periods)") +|> Tucan.set_width(680) +|> Tucan.set_height(260) +``` + +The Zoom FFT provides $1024/N$ times finer resolution than the standard FFT in +this band, making it straightforward to see how spectral energy is distributed +around the annual cycle frequency and whether it is concentrated (a clean +annual periodicity) or spread across a wider range of periods. + +## Summary + +| Technique | NxSignal function | Finance application | +| --------------- | -------------------------------------------- | ---------------------------------------------------------------------------- | +| FFT | `Nx.fft` | Identify dominant cycles and seasonal patterns | +| FIR low-pass | `Filters.firwin` + `Convolution.fftconvolve` | Noise removal; reproduces the moving average used in technical analysis | +| Zoom FFT | `NxSignal.zoom_fft` | High-resolution examination of a specific frequency band (e.g. annual cycle) | +| Window function | `Windows.hann` | Reduce spectral leakage before FFT | diff --git a/guides/peak_finding.livemd b/guides/peak_finding.livemd new file mode 100644 index 0000000..f240e6c --- /dev/null +++ b/guides/peak_finding.livemd @@ -0,0 +1,242 @@ +# Peak Finding + +```elixir +Mix.install([ + {:nx_signal, path: __DIR__ |> Path.join("..") |> Path.expand()}, + {:kino_vega_lite, "~> 0.1"}, + {:tucan, "~> 0.5"} +]) +``` + +## Introduction + +A local maximum is a sample that is strictly greater than all of its +neighbours within a given window. `NxSignal.PeakFinding` provides three +functions for finding such extrema: + +* `argrelmax(signal, order)` - local maxima +* `argrelmin(signal, order)` - local minima +* `argrelextrema(signal, comparator, order)` - generalised version + +The `order` parameter controls the half-width of the comparison window. +A point must beat all neighbours within `order` samples on each side to +qualify. Small `order` → sensitive (finds every little wiggle); large +`order` → selective (finds only broad peaks). + +All three functions return `%{indices: tensor, valid_indices: scalar}`. The +first `valid_indices` entries of `indices` are the detected positions; the +remainder are padded with −1. + +```elixir +# Helper: extract valid peak indices as a flat Elixir list +defmodule Peaks do + def indices(result) do + n = Nx.to_number(result.valid_indices) + if n == 0, do: [], else: result.indices |> Nx.slice([0], [n]) |> Nx.to_flat_list() + end +end +``` + +## Basic usage: effect of the order parameter + +A noisy signal consisting of two sinusoids illustrates how `order` controls +which peaks are returned. + +```elixir +fs = 200 +dur = 2.0 +t = Nx.linspace(0, dur, n: trunc(fs * dur), endpoint: false, type: :f32) + +# Signal: 5 Hz + 13 Hz + noise +key = Nx.Random.key(0) +{noise, _} = Nx.Random.normal(key, shape: {trunc(fs * dur)}, type: :f32) + +signal = + Nx.sin(Nx.multiply(t, 2 * :math.pi() * 5)) + |> Nx.add(Nx.multiply(0.4, Nx.sin(Nx.multiply(t, 2 * :math.pi() * 13)))) + |> Nx.add(Nx.multiply(0.15, noise)) + +t_list = Nx.to_flat_list(t) +s_list = Nx.to_flat_list(signal) +``` + +```elixir +# Show how the detected peaks change with order +for order <- [3, 10, 20] do + peak_idx = NxSignal.PeakFinding.argrelmax(signal, order: order) |> Peaks.indices() + + signal_data = Enum.zip(t_list, s_list) |> Enum.map(fn {ti, v} -> %{time: ti, amplitude: v} end) + peak_data = Enum.map(peak_idx, fn i -> %{time: Enum.at(t_list, i), amplitude: Enum.at(s_list, i)} end) + + chart = + VegaLite.new(width: 640, height: 160, + title: "argrelmax (order=#{order}): #{length(peak_idx)} peaks found" + ) + |> VegaLite.layers([ + VegaLite.new() + |> VegaLite.data_from_values(signal_data) + |> VegaLite.mark(:line, color: "steelblue") + |> VegaLite.encode_field(:x, "time", type: :quantitative, title: "Time (s)") + |> VegaLite.encode_field(:y, "amplitude", type: :quantitative, title: "Amplitude"), + VegaLite.new() + |> VegaLite.data_from_values(peak_data) + |> VegaLite.mark(:point, color: "tomato", size: 80, filled: true) + |> VegaLite.encode_field(:x, "time", type: :quantitative) + |> VegaLite.encode_field(:y, "amplitude", type: :quantitative) + ]) + + Kino.render(chart) +end + +:ok +``` + +With `order = 3` every small ripple is detected. With `order = 20`, roughly +half the period of the 5 Hz component, only the broad structural peaks +survive. + +**Rule of thumb**: set `order` to roughly half the period of the feature you +want to detect, measured in samples. + +## Spectral peak detection + +One of the most important applications is finding the dominant frequency +components of a signal automatically from its FFT magnitude. + +```elixir +fs_spec = 1_000 +dur_spec = 2.0 +n_spec = trunc(fs_spec * dur_spec) +t_spec = Nx.linspace(0, dur_spec, n: n_spec, endpoint: false, type: :f32) + +# Three tones with different amplitudes +x_spec = + Nx.sin(Nx.multiply(t_spec, 2 * :math.pi() * 50)) + |> Nx.add(Nx.multiply(0.6, Nx.sin(Nx.multiply(t_spec, 2 * :math.pi() * 130)))) + |> Nx.add(Nx.multiply(0.3, Nx.sin(Nx.multiply(t_spec, 2 * :math.pi() * 210)))) + +# FFT magnitude (positive frequencies only) +half_spec = div(n_spec, 2) +freqs_spec = NxSignal.fft_frequencies(fs_spec, fft_length: n_spec)[0..half_spec] |> Nx.to_flat_list() + +mag = + x_spec + |> Nx.as_type({:c, 64}) + |> Nx.fft() + |> Nx.abs() + |> Nx.slice([0], [half_spec + 1]) + +# Find spectral peaks: order=3 works well for FFT magnitudes +peak_result = NxSignal.PeakFinding.argrelmax(mag, order: 3) +peak_indices = Peaks.indices(peak_result) +mag_list = Nx.to_flat_list(mag) + +# Sort by amplitude, keep top 5 +peak_freqs = + peak_indices + |> Enum.map(fn i -> {Enum.at(freqs_spec, i), Enum.at(mag_list, i)} end) + |> Enum.sort_by(fn {_, a} -> -a end) + |> Enum.take(5) + +IO.inspect(Enum.map(peak_freqs, fn {f, _} -> Float.round(f, 1) end), label: "Detected frequencies (Hz)") +``` + +```elixir +spec_data = Enum.zip(freqs_spec, mag_list) |> Enum.map(fn {f, a} -> %{frequency: f, amplitude: a} end) +peak_pts = Enum.map(peak_indices, fn i -> %{frequency: Enum.at(freqs_spec, i), amplitude: Enum.at(mag_list, i)} end) + +VegaLite.new(width: 680, height: 240, title: "Spectral peak detection: 50, 130, 210 Hz") +|> VegaLite.layers([ + VegaLite.new() + |> VegaLite.data_from_values(spec_data) + |> VegaLite.mark(:line, color: "steelblue") + |> VegaLite.encode_field(:x, "frequency", type: :quantitative, title: "Frequency (Hz)") + |> VegaLite.encode_field(:y, "amplitude", type: :quantitative, title: "Amplitude"), + VegaLite.new() + |> VegaLite.data_from_values(peak_pts) + |> VegaLite.mark(:point, color: "tomato", size: 100, filled: true) + |> VegaLite.encode_field(:x, "frequency", type: :quantitative) + |> VegaLite.encode_field(:y, "amplitude", type: :quantitative) +]) +``` + +## Envelope detection + +An amplitude-modulated (AM) signal's positive envelope is traced by its +local maxima. Finding those peaks recovers the modulating waveform. + +```elixir +fs_am = 2_000 +dur_am = 1.0 +n_am = trunc(fs_am * dur_am) +t_am = Nx.linspace(0, dur_am, n: n_am, endpoint: false, type: :f32) + +# Carrier: 200 Hz; modulation: 5 Hz (one envelope cycle every 200 ms) +carrier = Nx.sin(Nx.multiply(t_am, 2 * :math.pi() * 200)) +modulation = Nx.add(0.5, Nx.multiply(0.5, Nx.sin(Nx.multiply(t_am, 2 * :math.pi() * 5)))) +am_signal = Nx.multiply(carrier, modulation) + +# Order ≈ half the carrier period in samples: 2000 / 200 / 2 = 5 +env_result = NxSignal.PeakFinding.argrelmax(am_signal, order: 5) +env_indices = Peaks.indices(env_result) + +t_am_list = Nx.to_flat_list(t_am) +am_list = Nx.to_flat_list(am_signal) + +am_data = Enum.zip(t_am_list, am_list) |> Enum.map(fn {ti, v} -> %{time: ti, amplitude: v} end) +env_data = Enum.map(env_indices, fn i -> %{time: Enum.at(t_am_list, i), amplitude: Enum.at(am_list, i)} end) + +VegaLite.new(width: 680, height: 220, + title: "Envelope detection: AM signal (carrier 200 Hz, modulation 5 Hz)" +) +|> VegaLite.layers([ + VegaLite.new() + |> VegaLite.data_from_values(am_data) + |> VegaLite.mark(:line, color: "lightsteelblue", stroke_width: 0.8) + |> VegaLite.encode_field(:x, "time", type: :quantitative, title: "Time (s)") + |> VegaLite.encode_field(:y, "amplitude", type: :quantitative, title: "Amplitude"), + VegaLite.new() + |> VegaLite.data_from_values(env_data) + |> VegaLite.mark(:point, color: "tomato", size: 30, filled: true) + |> VegaLite.encode_field(:x, "time", type: :quantitative) + |> VegaLite.encode_field(:y, "amplitude", type: :quantitative) +]) +``` + +The detected peaks (red) trace the sinusoidal envelope at 5 Hz. Connecting +them would reproduce the modulating signal. + +## Local minima + +`argrelmin` finds troughs, the negative envelope of the same signal: + +```elixir +min_result = NxSignal.PeakFinding.argrelmin(am_signal, order: 5) +min_indices = Peaks.indices(min_result) +min_data = Enum.map(min_indices, fn i -> %{time: Enum.at(t_am_list, i), amplitude: Enum.at(am_list, i)} end) + +VegaLite.new(width: 680, height: 220, title: "Local minima: negative envelope") +|> VegaLite.layers([ + VegaLite.new() + |> VegaLite.data_from_values(am_data) + |> VegaLite.mark(:line, color: "lightsteelblue", stroke_width: 0.8) + |> VegaLite.encode_field(:x, "time", type: :quantitative, title: "Time (s)") + |> VegaLite.encode_field(:y, "amplitude", type: :quantitative, title: "Amplitude"), + VegaLite.new() + |> VegaLite.data_from_values(min_data) + |> VegaLite.mark(:point, color: "darkorange", size: 30, filled: true) + |> VegaLite.encode_field(:x, "time", type: :quantitative) + |> VegaLite.encode_field(:y, "amplitude", type: :quantitative) +]) +``` + +## Summary + +| Function | Returns | Typical use | +| ------------------------------------------ | -------------- | ----------------------------------------------------------------- | +| `argrelmax(signal, order)` | Local maxima | Peak frequency identification, positive envelope, ridge detection | +| `argrelmin(signal, order)` | Local minima | Trough detection, negative envelope | +| `argrelextrema(signal, comparator, order)` | Custom extrema | Any comparison function (e.g. `>=` for plateaus) | + +All return `%{indices: tensor, valid_indices: scalar}`. Always use +`valid_indices` to slice off the padding before processing results. diff --git a/guides/sound_synthesis.livemd b/guides/sound_synthesis.livemd new file mode 100644 index 0000000..9062297 --- /dev/null +++ b/guides/sound_synthesis.livemd @@ -0,0 +1,447 @@ +# Sound Synthesis with NxSignal + +```elixir +Mix.install([ + {:nx_signal, path: __DIR__ |> Path.join("..") |> Path.expand()}, + {:kino, "~> 0.13"}, + {:kino_vega_lite, "~> 0.1"}, + {:tucan, "~> 0.5"} +]) +``` + +## Audio helper + +This module converts a 1-D f32 Nx tensor of normalised PCM samples +(values in $[-1, 1]$) into a WAV binary that `Kino.Audio` can play directly +in the browser. + +```elixir +defmodule Audio do + @moduledoc """ + Utilities for rendering Nx tensors as playable audio in Livebook. + """ + + @doc """ + Encode a 1-D f32 tensor as a 16-bit mono WAV binary. + + ## Examples + + Audio.play(signal, 44_100) + """ + def to_wav(samples, sample_rate \\ 44_100) do + pcm = + samples + |> Nx.as_type(:f32) + |> Nx.clip(-1.0, 1.0) + |> Nx.multiply(32_767) + |> Nx.as_type({:s, 16}) + |> Nx.to_binary() + + data_size = byte_size(pcm) + byte_rate = sample_rate * 2 + chunk_size = 36 + data_size + + << + "RIFF", chunk_size::little-32, + "WAVE", + "fmt ", 16::little-32, + 1::little-16, # PCM + 1::little-16, # mono + sample_rate::little-32, + byte_rate::little-32, + 2::little-16, # block align (1 channel × 2 bytes) + 16::little-16, # bits per sample + "data", data_size::little-32, + pcm::binary + >> + end + + @doc "Render a tensor as an inline audio player." + def play(samples, sample_rate \\ 44_100, opts \\ []) do + samples |> to_wav(sample_rate) |> Kino.Audio.new(:wav, opts) + end + + @doc "Normalise a tensor so its peak is at ±1." + def normalise(t) do + peak = t |> Nx.abs() |> Nx.reduce_max() + Nx.divide(t, Nx.max(peak, 1.0e-9)) + end +end +``` + +## Your first sound + +A 440 Hz sine wave: this is concert pitch A4. + +```elixir +fs = 44_100 +t = Nx.linspace(0, 1.0, n: fs, type: :f32) + +signal = Nx.sin(Nx.multiply(t, 2 * :math.pi() * 440)) + +Audio.play(signal, fs) +``` + +## Waveform shapes + +The *shape* of a wave determines its timbre. Let's listen to three classic +waveforms all tuned to the same pitch (220 Hz, A3) and compare their +character. + +Sine waves are pure: they are just the fundamental frequency, no harmonics. This is why they sound a bit cold to the human ear. Square and sawtooth waves are rich in harmonics, which is why they sound buzzy or bright. + +```elixir +fs = 44_100 +dur = 0.8 +t = Nx.linspace(0, dur, n: trunc(fs * dur), type: :f32) +phi = Nx.multiply(t, 2 * :math.pi() * 220) + +sine = Nx.sin(phi) +square = NxSignal.Waveforms.square(phi) |> Nx.as_type(:f32) +sawtooth = NxSignal.Waveforms.sawtooth(phi) + +Kino.Layout.grid([ + Kino.Layout.grid([Kino.Text.new("Sine"), Audio.play(sine, fs)], columns: 2), + Kino.Layout.grid([Kino.Text.new("Square"), Audio.play(square, fs)], columns: 2), + Kino.Layout.grid([Kino.Text.new("Sawtooth"), Audio.play(sawtooth, fs)], columns: 2) +]) +``` + +Showing the first two cycles (≈ 9 ms) makes the structural differences between +the waveforms immediately visible. + +```elixir +# 2 cycles at 220 Hz = 2/220 s ≈ 401 samples at 44 100 Hz +n_show = 401 +t_show = t |> Nx.slice([0], [n_show]) |> Nx.to_flat_list() + +wave_data = + [{sine, "Sine"}, {square, "Square"}, {sawtooth, "Sawtooth"}] + |> Enum.flat_map(fn {wave, name} -> + wave + |> Nx.slice([0], [n_show]) + |> Nx.to_flat_list() + |> Enum.zip(t_show) + |> Enum.map(fn {v, ti} -> %{time: ti, amplitude: v, waveform: name} end) + end) + +Tucan.lineplot(wave_data, "time", "amplitude") +|> Tucan.color_by("waveform") +|> Tucan.facet_by(:row, "waveform") +|> Tucan.set_title("Waveform shapes: first 2 cycles at 220 Hz") +|> Tucan.set_width(640) +|> Tucan.set_height(80) +``` + +## Musical intervals and chords + +A musical interval is a ratio of frequencies. The octave is 2:1, a perfect +fifth is 3:2, and a major third is 5:4. Stack them and you get a major chord. + +```elixir +# Helper: generate a sine tone at frequency f for dur seconds +tone = fn f, dur, fs -> + t = Nx.linspace(0, dur, n: trunc(fs * dur), type: :f32) + Nx.sin(Nx.multiply(t, 2 * :math.pi() * f)) +end + +fs = 44_100 +dur = 1.5 + +# C major chord: C4 (261.6 Hz), E4 (329.6 Hz), G4 (392.0 Hz) +c = tone.(261.6, dur, fs) +e = tone.(329.6, dur, fs) +g = tone.(392.0, dur, fs) + +chord = + Nx.add(c, Nx.add(e, g)) + |> Audio.normalise() + +Kino.Layout.grid([ + Kino.Layout.grid([Kino.Text.new("C4"), Audio.play(c, fs)], columns: 2), + Kino.Layout.grid([Kino.Text.new("E4"), Audio.play(e, fs)], columns: 2), + Kino.Layout.grid([Kino.Text.new("G4"), Audio.play(g, fs)], columns: 2), + Kino.Layout.grid([Kino.Text.new("C major"), Audio.play(chord, fs)], columns: 2) +]) +``` + +## Chirp sweep + +`NxSignal.Waveforms.chirp/5` generates a sinusoid whose instantaneous +frequency sweeps continuously from `f0` to `f1` over the interval `[0, t1]`. +This is the same mathematical structure used in radar and sonar pulses, and +also in the bat echolocation clicks that inspired the name "chirp z-transform". + +```elixir +fs = 44_100 +dur = 3.0 +t = Nx.linspace(0, dur, n: trunc(fs * dur), type: :f32) + +# Linear sweep: 80 Hz → 3 400 Hz over 3 seconds +chirp_linear = NxSignal.Waveforms.chirp(t, 80, dur, 3400, method: :linear) + +# Logarithmic sweep: same range but exponential spacing — +# sounds more even to the ear because pitch perception is logarithmic +chirp_log = NxSignal.Waveforms.chirp(t, 80, dur, 3400, method: :logarithmic) + +Kino.Layout.grid([ + Kino.Layout.grid([Kino.Text.new("Linear sweep"), Audio.play(chirp_linear, fs)], columns: 2), + Kino.Layout.grid([Kino.Text.new("Logarithmic sweep"), Audio.play(chirp_log, fs)], columns: 2) +]) +``` + +A spectrogram computed via `NxSignal.stft/3` makes the sweep visible: each +column is the short-time power spectrum at that moment, so a rising tone +appears as a diagonal line sweeping upwards. + +```elixir +win_len = 2048 + +{spectrum, times, freqs} = + NxSignal.stft( + chirp_linear, + NxSignal.Windows.hann(win_len, is_periodic: true), + sampling_rate: fs, + overlap_length: win_len - 512 + ) + +# Restrict to positive frequencies below 4 000 Hz +max_bin = trunc(4000 / (fs / win_len)) + 1 + +magnitude_db = + spectrum + |> Nx.slice_along_axis(0, max_bin, axis: 1) + |> Nx.abs() + |> Nx.log10() + |> Nx.multiply(20) + +{n_frames, n_bins} = Nx.shape(magnitude_db) + +# Scaling factors for axis label expressions +hz_per_bin = Float.round(fs / win_len, 2) +times_list = Nx.to_flat_list(times) +s_per_frame = Float.round( + (List.last(times_list) - List.first(times_list)) / (n_frames - 1), + 5 +) + +mag_list = Nx.to_flat_list(magnitude_db) + +# Integer frame/bin indices: VegaLite quantitative y-axis puts bin 0 (= 0 Hz) +# at the bottom, giving the correct upward-sweeping diagonal. +spec_data = + Enum.with_index(mag_list, fn val, idx -> + %{frame: div(idx, n_bins), bin: rem(idx, n_bins), power_db: val} + end) + +VegaLite.new(width: 680, height: 300, + title: "Spectrogram: linear chirp 80 → 3 400 Hz" +) +|> VegaLite.data_from_values(spec_data) +|> VegaLite.mark(:rect) +|> VegaLite.encode_field(:x, "frame", + type: :quantitative, + bin: [step: 1], + title: "Time (s)", + axis: [tick_count: 6, label_expr: "format(datum.value * #{s_per_frame}, '.1f')"] + ) +|> VegaLite.encode_field(:y, "bin", + type: :quantitative, + bin: [step: 1], + title: "Frequency (Hz)", + axis: [tick_count: 5, label_expr: "format(datum.value * #{hz_per_bin}, '.0f')"] + ) +|> VegaLite.encode_field(:color, "power_db", + type: :quantitative, + aggregate: :mean, + scale: [scheme: "viridis"], + title: "Power (dB)" + ) +``` + +## Envelope shaping + +A raw oscillator plays at full volume for its entire duration, which sounds +unnatural. Real instruments have an **ADSR envelope**: + +* **A**ttack: ramp up from silence +* **D**ecay: settle to the sustain level +* **S**ustain: hold while the note is held +* **R**elease: fade to silence + +```elixir +defmodule ADSR do + @doc """ + Build a linear ADSR envelope tensor. + + Times are in seconds; `sustain` is a level in [0, 1]. + """ + def envelope(attack, decay, sustain, release, fs) do + make = fn n_samples, from, to -> + Nx.linspace(from, to, n: n_samples, type: :f32) + end + + n_a = trunc(attack * fs) + n_d = trunc(decay * fs) + n_r = trunc(release * fs) + # sustain phase: fixed at sustain level for 1 second + n_s = fs + + Nx.concatenate([ + make.(n_a, 0.0, 1.0), + make.(n_d, 1.0, sustain), + Nx.broadcast(Nx.tensor(sustain, type: :f32), {n_s}), + make.(n_r, sustain, 0.0) + ]) + end +end + +fs = 44_100 + +env = ADSR.envelope(1, 0.1, 0.6, 1, fs) + +t = Nx.linspace(0, Nx.size(env) / fs, n: Nx.size(env), type: :f32) +note = Nx.sin(Nx.multiply(t, 2 * :math.pi() * 440)) |> Nx.multiply(env) + +Kino.Layout.grid([ + Kino.Layout.grid([Kino.Text.new("Raw sine"), Audio.play(Nx.sin(Nx.multiply(t, 2 * :math.pi() * 440)), fs)], columns: 2), + Kino.Layout.grid([Kino.Text.new("With ADSR"), Audio.play(note, fs)], columns: 2) +]) +``` + +```elixir +env_data = + Enum.zip(Nx.to_flat_list(t), Nx.to_flat_list(env)) + |> Enum.map(fn {ti, v} -> %{time: ti, level: v} end) + +Tucan.lineplot(env_data, "time", "level") +|> Tucan.Axes.set_x_title("Time (s)") +|> Tucan.Axes.set_y_title("Amplitude") +|> Tucan.set_title("ADSR envelope: attack 10 ms, decay 100 ms, sustain 0.6, release 400 ms") +|> Tucan.set_width(640) +|> Tucan.set_height(180) +``` + +## Additive synthesis + +Any periodic timbre can be constructed by summing sine waves at the +fundamental frequency and its harmonics (integer multiples). Choosing which +harmonics to include, and how loud each one is, directly sculpts the sound. + +Below we build three timbres all on the same fundamental (G2, 98 Hz): + +* **Pure**: just the fundamental +* **Warm**: fundamental + a few quiet even harmonics +* **Bright**: fundamental + many harmonics with equal weighting (approximates sawtooth) + +```elixir +fs = 44_100 +dur = 1.5 +f0 = 98.0 # G2 + +t = Nx.linspace(0, dur, n: trunc(fs * dur), type: :f32) + +# Build a tone as a sum of harmonics with given amplitudes +additive = fn partials -> + partials + |> Enum.reduce(Nx.broadcast(0.0, {trunc(fs * dur)}), fn {n, amp}, acc -> + harmonic = Nx.sin(Nx.multiply(t, 2 * :math.pi() * f0 * n)) + Nx.add(acc, Nx.multiply(amp, harmonic)) + end) + |> Audio.normalise() +end + +pure = additive.([{1, 1.0}]) +warm = additive.([{1, 1.0}, {2, 0.5}, {3, 0.25}, {4, 0.12}]) +bright = additive.(for n <- 1..16, do: {n, 1 / n}) + +Kino.Layout.grid([ + Kino.Layout.grid([Kino.Text.new("Pure (1 partial)"), Audio.play(pure, fs)], columns: 2), + Kino.Layout.grid([Kino.Text.new("Warm (4 partials)"), Audio.play(warm, fs)], columns: 2), + Kino.Layout.grid([Kino.Text.new("Bright (16 partials)"), Audio.play(bright, fs)], columns: 2) +]) +``` + +The spectrum shows clearly why these tones sound different: each harmonic +partial appears as a discrete peak. + +```elixir +n = Nx.size(pure) +half = div(n, 2) +freqs = Nx.linspace(0, fs / 2, n: half, type: :f32) + +spectrum_data = + [{pure, "Pure"}, {warm, "Warm"}, {bright, "Bright"}] + |> Enum.flat_map(fn {sig, name} -> + amps = + sig + |> Nx.as_type({:c, 64}) + |> Nx.fft() + |> Nx.abs() + |> Nx.slice([0], [half]) + + Enum.zip(Nx.to_flat_list(freqs), Nx.to_flat_list(amps)) + |> Enum.filter(fn {f, _} -> f <= 2000 end) + |> Enum.map(fn {f, a} -> %{frequency: f, amplitude: a, timbre: name} end) + end) + +Tucan.lineplot(spectrum_data, "frequency", "amplitude") +|> Tucan.color_by("timbre") +|> Tucan.Axes.set_x_title("Frequency (Hz)") +|> Tucan.Axes.set_y_title("Amplitude") +|> Tucan.set_title("Spectra of additive synthesis tones (G2, 98 Hz, up to 2 kHz)") +|> Tucan.set_width(640) +|> Tucan.set_height(220) +``` + +## FIR low-pass filter + +`NxSignal.Filters.firwin/3` designs an FIR filter using the window method. +Here we apply a low-pass filter to the bright (16-partial) tone from above and +listen to the harmonics being progressively removed. + +```elixir +alias NxSignal.Convolution + +cutoff_hz = 500.0 +nyquist = fs / 2.0 + +# Design a 101-tap linear-phase low-pass filter +coeffs = NxSignal.Filters.firwin(101, [cutoff_hz / nyquist]) + +filtered = + Convolution.convolve(bright, coeffs, mode: :same, method: :fft) + |> Audio.normalise() + +Kino.Layout.grid([ + Kino.Layout.grid([Kino.Text.new("Unfiltered"), Audio.play(bright, fs)], columns: 2), + Kino.Layout.grid([Kino.Text.new("Low-pass 500 Hz"), Audio.play(filtered, fs)], columns: 2) +]) +``` + +The filter removes all harmonics above 500 Hz. Only five of the sixteen partials pass through. The attenuation at the cutoff is clearly visible in the spectrum. + +```elixir +filter_spec_data = + [{bright, "Unfiltered"}, {filtered, "Low-pass 500 Hz"}] + |> Enum.flat_map(fn {sig, label} -> + amps = + sig + |> Nx.as_type({:c, 64}) + |> Nx.fft() + |> Nx.abs() + |> Nx.slice([0], [half]) + + Enum.zip(Nx.to_flat_list(freqs), Nx.to_flat_list(amps)) + |> Enum.filter(fn {f, _} -> f <= 2000 end) + |> Enum.map(fn {f, a} -> %{frequency: f, amplitude: a, signal: label} end) + end) + +Tucan.lineplot(filter_spec_data, "frequency", "amplitude") +|> Tucan.color_by("signal") +|> Tucan.Axes.set_x_title("Frequency (Hz)") +|> Tucan.Axes.set_y_title("Amplitude") +|> Tucan.set_title("Effect of FIR low-pass filter on spectrum (cutoff 500 Hz)") +|> Tucan.set_width(640) +|> Tucan.set_height(220) +``` diff --git a/guides/spectrogram.livemd b/guides/spectrogram.livemd index fba659b..d6b5953 100644 --- a/guides/spectrogram.livemd +++ b/guides/spectrogram.livemd @@ -38,7 +38,7 @@ n = Nx.iota({full_n}) d = %{data: Nx.to_flat_list(data[[1000..1250]]), n: Nx.to_flat_list(n[[1000..1250]])} ``` - + ```elixir VegaLite.new(width: 600, height: 600, title: "Audio Sample") @@ -51,7 +51,6 @@ VegaLite.new(width: 600, height: 600, title: "Audio Sample") ```elixir defmodule Spectrogram do alias VegaLite, as: Vl - import Nx.Defn def calculate_stft_and_plot_spectrogram( input, @@ -60,7 +59,8 @@ defmodule Spectrogram do plot_cutoff_frequency \\ 4000 ) do n_window = ceil(fs * window_duration_ms) - {spectrogram, f, t, max_f} = stft(input, fs: fs, n_window: n_window) + window = NxSignal.Windows.hann(n_window, is_periodic: true) + {spectrogram, f, t, max_f} = stft(input, window, fs) max_f = Nx.to_number(max_f) spectrogram = Nx.slice(spectrogram, [0, 0], [Nx.size(t), max_f]) @@ -71,23 +71,15 @@ defmodule Spectrogram do |> plot() end - defn stft(input, opts) do - fs = opts[:fs] - n_window = opts[:n_window] - - # ms to samples - window = NxSignal.Windows.hann(n: n_window, is_periodic: true) - - # use the default overlap of 50% + defp stft(input, window, fs) do {s, t, f} = NxSignal.stft(input, window, sampling_rate: fs, fft_length: 1024) max_f = - Nx.select(f >= fs / 2, Nx.iota(f.shape), Nx.size(f) + 1) + Nx.select(Nx.greater_equal(f, fs / 2), Nx.iota(Nx.shape(f)), Nx.size(f) + 1) |> Nx.argmin() spectrogram = Nx.abs(s) - # to dBFS - spectrogram = 20 * Nx.log(spectrogram / Nx.reduce_max(spectrogram)) / Nx.log(10) + spectrogram = Nx.multiply(20, Nx.divide(Nx.log(Nx.divide(spectrogram, Nx.reduce_max(spectrogram))), Nx.log(10))) {spectrogram, f, t, max_f} end diff --git a/guides/windows.livemd b/guides/windows.livemd new file mode 100644 index 0000000..517172c --- /dev/null +++ b/guides/windows.livemd @@ -0,0 +1,255 @@ +# Window Functions + +```elixir +Mix.install([ + {:nx_signal, path: __DIR__ |> Path.join("..") |> Path.expand()}, + {:kino_vega_lite, "~> 0.1"}, + {:tucan, "~> 0.5"} +]) +``` + +## Introduction + +When we compute a DFT we analyse only a finite slice of a signal. If the +slice does not contain an exact integer number of cycles of a frequency +component, that component appears to "jump" discontinuously at the edges of +the window. The Fourier transform sees this as energy spread across many +bins, not just the true one. This is **spectral leakage**. + +A **window function** tapers the signal smoothly to zero at both ends before +the FFT, eliminating the discontinuity. The cost is a slightly wider main lobe +around each true frequency; the benefit is dramatically lower sidelobe levels +elsewhere. Every window represents a different point on this trade-off. + +Formally, multiplying the signal by a window $w[n]$ in the time domain is +equivalent to convolving the spectrum with $W(f)$—the Fourier transform of +the window—in the frequency domain. A narrow $W(f)$ preserves frequency +resolution; low sidelobes on $W(f)$ suppress leakage. + +## The leakage problem in practice + +A tone at exactly a DFT bin frequency leaks nothing; one between bins leaks +into every other bin. 10.5 Hz sits halfway between the 10 Hz and 11 Hz bins. + +```elixir +fs = 200 +n = 200 +f_tone = 10.5 # halfway between bins; worst case for leakage + +t = Nx.linspace(0, n / fs, n: n, endpoint: false, type: :f32) +x = Nx.sin(Nx.multiply(t, 2 * :math.pi() * f_tone)) + +# Rectangular window is just the raw signal (multiply by all-ones) +rect_window = NxSignal.Windows.rectangular(n, type: :f32) +hann_window = NxSignal.Windows.hann(n, is_periodic: false) + +half = div(n, 2) +freqs = NxSignal.fft_frequencies(fs, fft_length: n)[0..half] |> Nx.to_flat_list() + +to_db = fn window, name -> + amps = + Nx.multiply(x, window) + |> Nx.as_type({:c, 64}) + |> Nx.fft() + |> Nx.abs() + |> Nx.slice([0], [half + 1]) + + peak = Nx.reduce_max(amps) + amps_norm = Nx.divide(amps, peak) + + # Floor at -120 dB to avoid -inf for zero bins + db = + amps_norm + |> Nx.log10() + |> Nx.multiply(20) + |> Nx.max(-120.0) + |> Nx.to_flat_list() + + Enum.zip(freqs, db) + |> Enum.map(fn {f, d} -> %{frequency: f, amplitude_db: d, window: name} end) +end + +leakage_data = + to_db.(rect_window, "Rectangular") ++ + to_db.(hann_window, "Hann") + +Tucan.lineplot(leakage_data, "frequency", "amplitude_db") +|> Tucan.color_by("window") +|> Tucan.Axes.set_x_title("Frequency (Hz)") +|> Tucan.Axes.set_y_title("Amplitude (dB)") +|> Tucan.set_title("Spectral leakage: 10.5 Hz tone, rectangular vs Hann window") +|> Tucan.set_width(680) +|> Tucan.set_height(280) +``` + +The rectangular window leaks energy across the entire spectrum. The Hann +window concentrates it near 10.5 Hz, with the remainder buried in the noise +floor. + +## Window shapes + +NxSignal provides seven windows. Here is how they look in the time domain: + +```elixir +n_shape = 64 + +windows_list = [ + {"Rectangular", NxSignal.Windows.rectangular(n_shape, type: :f32)}, + {"Bartlett", NxSignal.Windows.bartlett(n_shape)}, + {"Triangular", NxSignal.Windows.triangular(n_shape)}, + {"Hann", NxSignal.Windows.hann(n_shape, is_periodic: false)}, + {"Hamming", NxSignal.Windows.hamming(n_shape, is_periodic: false)}, + {"Blackman", NxSignal.Windows.blackman(n_shape, is_periodic: false)}, + {"Kaiser β=14", NxSignal.Windows.kaiser(n_shape, beta: 14)} +] + +samples = Nx.iota({n_shape}, type: :f32) |> Nx.to_flat_list() + +shape_data = + Enum.flat_map(windows_list, fn {name, w} -> + Enum.zip(samples, Nx.to_flat_list(w)) + |> Enum.map(fn {i, v} -> %{sample: i, amplitude: v, window: name} end) + end) + +Tucan.lineplot(shape_data, "sample", "amplitude") +|> Tucan.color_by("window") +|> Tucan.Axes.set_x_title("Sample") +|> Tucan.Axes.set_y_title("Amplitude") +|> Tucan.set_title("Window functions in the time domain (N = 64)") +|> Tucan.set_width(680) +|> Tucan.set_height(260) +``` + +## Frequency responses: the sidelobe comparison + +This is the chart that matters most for choosing a window. Zero-padding to +1024 points interpolates the frequency axis so the sidelobe structure is +clearly visible. + +```elixir +pad_to = 1024 + +freq_resp_data = + Enum.flat_map(windows_list, fn {name, w} -> + padded = Nx.pad(Nx.as_type(w, {:c, 64}), Nx.as_type(0, {:c, 64}), [{0, pad_to - n_shape, 0}]) + + amps = Nx.fft(padded) |> Nx.abs() |> Nx.slice([0], [div(pad_to, 2)]) + peak = Nx.reduce_max(amps) + + db = + Nx.divide(amps, peak) + |> Nx.log10() + |> Nx.multiply(20) + |> Nx.max(-120.0) + |> Nx.to_flat_list() + + # Normalised frequency 0..0.5 + bins = Enum.map(0..(div(pad_to, 2) - 1), fn k -> k / pad_to end) + + Enum.zip(bins, db) + |> Enum.map(fn {f, d} -> %{norm_freq: f, amplitude_db: d, window: name} end) + end) + +Tucan.lineplot(freq_resp_data, "norm_freq", "amplitude_db") +|> Tucan.color_by("window") +|> VegaLite.encode_field(:y, "amplitude_db", + type: :quantitative, + title: "Amplitude (dB)", + scale: [domain: [-120, 5]] + ) +|> Tucan.Axes.set_x_title("Normalised frequency (cycles/sample)") +|> Tucan.set_title("Window frequency responses: sidelobe levels (N = 64, zero-padded to 1024)") +|> Tucan.set_width(680) +|> Tucan.set_height(300) +``` + +The trade-off is plain: the rectangular window has the narrowest main lobe +but the highest sidelobes (−13 dB). Blackman and Kaiser push sidelobes below +−60 dB at the cost of a wider main lobe. + +## Key parameters + +| Window | Peak sidelobe | Typical use | +| -------------- | ------------- | --------------------------------------------- | +| Rectangular | −13 dB | Transients; when leakage is not a concern | +| Bartlett | −25 dB | General purpose; linear taper | +| Triangular | −27 dB | Similar to Bartlett | +| Hann | −31 dB | **Default choice** for spectrum analysis | +| Hamming | −41 dB | FIR filter design (non-zero endpoints) | +| Blackman | −57 dB | High dynamic range measurements | +| Kaiser (β=14) | −60 dB+ | Configurable; use when requirements are known | + +## The Kaiser β parameter + +Kaiser is unique in that a single parameter $\beta$ continuously controls the +sidelobe level vs main-lobe width trade-off. As $\beta \to 0$ the window +approaches rectangular; $\beta = 1$ is already nearly flat across most of the +window. + +```elixir +kaiser_data = + [{"β = 1", NxSignal.Windows.kaiser(n_shape, beta: 1)}, + {"β = 5", NxSignal.Windows.kaiser(n_shape, beta: 5)}, + {"β = 10", NxSignal.Windows.kaiser(n_shape, beta: 10)}, + {"β = 14", NxSignal.Windows.kaiser(n_shape, beta: 14)}] + |> Enum.flat_map(fn {name, w} -> + padded = Nx.pad(Nx.as_type(w, {:c, 64}), Nx.as_type(0, {:c, 64}), [{0, pad_to - n_shape, 0}]) + amps = Nx.fft(padded) |> Nx.abs() |> Nx.slice([0], [div(pad_to, 2)]) + peak = Nx.reduce_max(amps) + db = + Nx.divide(amps, peak) + |> Nx.log10() + |> Nx.multiply(20) + |> Nx.max(-120.0) + |> Nx.to_flat_list() + bins = Enum.map(0..(div(pad_to, 2) - 1), fn k -> k / pad_to end) + Enum.zip(bins, db) + |> Enum.map(fn {f, d} -> %{norm_freq: f, amplitude_db: d, window: name} end) + end) + +Tucan.lineplot(kaiser_data, "norm_freq", "amplitude_db") +|> Tucan.color_by("window") +|> VegaLite.encode_field(:y, "amplitude_db", + type: :quantitative, + title: "Amplitude (dB)", + scale: [domain: [-120, 5]] + ) +|> Tucan.Axes.set_x_title("Normalised frequency (cycles/sample)") +|> Tucan.set_title("Kaiser window: effect of β on sidelobe level (N = 64)") +|> Tucan.set_width(680) +|> Tucan.set_height(280) +``` + +Larger $\beta$ widens the main lobe but drives sidelobes lower. For +`NxSignal.Filters.firwin`, the `:kaiser` window is available via the +`window:` option and its $\beta$ is derived automatically from the stopband +attenuation specification. + +## Choosing a window + +**If in doubt, use Hann.** It is the best all-round choice for general +spectrum analysis: moderate main-lobe width and sidelobes low enough for most +practical SNR ranges. + +Use **Hamming** when designing FIR filters because its non-zero endpoints give +slightly better passband ripple. + +Use **Blackman** or **Kaiser** when you need to measure a weak component +sitting close to a strong one. The −57 dB / −60 dB sidelobe floor allows you +to see signals 600× weaker than a nearby dominant tone. + +Use **Rectangular** for transient or impulse signals where the signal is +already zero at the edges, or for coherent averaging where you know the signal +is periodic with exactly $N$ samples. + +## Connection to filtering + +Windowing is multiplication in the time domain. By the convolution theorem, +multiplication in time equals convolution in frequency: applying window +$w[n]$ to a signal smears each spectral line by convolving it with $W(f)$. + +A rectangular window has a $\text{sinc}$-shaped $W(f)$ with slow-decaying +sidelobes. Smooth windows have more concentrated $W(f)$, so nearby spectral +lines stay separated. This is also why `NxSignal.Filters.firwin` is named for +the **win**dow method: it multiplies the ideal (sinc-shaped) impulse response +by a window to control the filter's transition band and stopband behaviour. From d521c27ed9648755c8e7243c6674f18266eabb77 Mon Sep 17 00:00:00 2001 From: thomaspmurphy Date: Tue, 16 Jun 2026 14:06:21 +0100 Subject: [PATCH 2/3] refactor(guides): zip_with, rfft, cumulative_sum per review feedback --- guides/convolution.livemd | 86 ++++++++++++++++----------------- guides/financial_signals.livemd | 22 ++++----- 2 files changed, 50 insertions(+), 58 deletions(-) diff --git a/guides/convolution.livemd b/guides/convolution.livemd index 8fd803b..50987b4 100644 --- a/guides/convolution.livemd +++ b/guides/convolution.livemd @@ -55,20 +55,20 @@ x = y = NxSignal.Convolution.convolve(x, h, mode: :full) -samples_x = Enum.map(0..29, & &1) -samples_y = Enum.map(0..(Nx.size(y) - 1), & &1) - impulse_data = - Enum.zip(samples_x, Nx.to_flat_list(x)) - |> Enum.map(fn {i, v} -> %{sample: i, value: v, signal: "Input x[n]"} end) + Enum.zip_with(Enum.to_list(0..29), Nx.to_flat_list(x), fn i, v -> + %{sample: i, value: v, signal: "Input x[n]"} + end) filter_data = - Enum.zip(Enum.map(0..4, & &1), Nx.to_flat_list(h)) - |> Enum.map(fn {i, v} -> %{sample: i, value: v, signal: "Filter h[n]"} end) + Enum.zip_with(Enum.to_list(0..4), Nx.to_flat_list(h), fn i, v -> + %{sample: i, value: v, signal: "Filter h[n]"} + end) output_data = - Enum.zip(samples_y, Nx.to_flat_list(y)) - |> Enum.map(fn {i, v} -> %{sample: i, value: v, signal: "Output y[n] = x * h"} end) + Enum.zip_with(Enum.to_list(0..(Nx.size(y) - 1)), Nx.to_flat_list(y), fn i, v -> + %{sample: i, value: v, signal: "Output y[n] = x * h"} + end) all_data = impulse_data ++ filter_data ++ output_data @@ -117,18 +117,13 @@ h_fir = NxSignal.Filters.firwin(51, [cutoff_norm]) # Route A: direct time-domain convolution route_a = NxSignal.Convolution.convolve(signal, h_fir, mode: :same, method: :direct) -# Route B: manual frequency-domain multiplication -fft_len = n + Nx.size(h_fir) - 1 -# Next power of two for efficiency -fft_len_padded = Integer.pow(2, ceil(:math.log2(fft_len))) - -x_fft = Nx.fft(Nx.as_type(signal, {:c, 64}), length: fft_len_padded) -h_fft = Nx.fft(Nx.as_type(h_fir, {:c, 64}), length: fft_len_padded) -y_full = Nx.ifft(Nx.multiply(x_fft, h_fft)) - -# Trim to :same length (centre of the full convolution) -trim = div(Nx.size(h_fir), 2) -route_b = Nx.real(y_full[trim..(trim + n - 1)]) +# Route B: FFT-based convolution +# Under the hood NxSignal zero-pads both signals to the next power of two, +# multiplies their FFTs, IFFTs the result, and trims to the requested length: +# +# Y = IFFT(FFT(x, L) · FFT(h, L)), L = next_pow2(N + M - 1) +# +route_b = NxSignal.Convolution.convolve(signal, h_fir, mode: :same, method: :fft) # Confirm both routes agree max_diff = @@ -141,32 +136,36 @@ IO.inspect(Nx.to_number(max_diff), label: "max |Route A − Route B|") ```elixir # Visualise: input spectrum, filter response, and output spectrum side by side -half = div(n, 2) -freqs_hz = Nx.linspace(0, fs / 2.0, n: half, type: :f32) |> Nx.to_flat_list() +# rfft returns div(n, 2) + 1 bins (DC through Nyquist inclusive) +half = div(n, 2) + 1 +freqs_hz = Nx.linspace(0, fs / 2.0, n: half, endpoint: true, type: :f32) |> Nx.to_flat_list() to_amp_list = fn sig -> sig - |> Nx.as_type({:c, 64}) - |> Nx.fft() + |> Nx.rfft() |> Nx.abs() - |> Nx.slice([0], [half]) |> Nx.to_flat_list() end -h_padded = Nx.pad(Nx.as_type(h_fir, {:c, 64}), Nx.as_type(0, {:c, 64}), [{0, n - Nx.size(h_fir), 0}]) -h_spectrum = h_padded |> Nx.fft() |> Nx.abs() |> Nx.slice([0], [half]) |> Nx.to_flat_list() - -spectrum_data = - Enum.zip(freqs_hz, to_amp_list.(signal)) - |> Enum.map(fn {f, a} -> %{frequency: f, amplitude: a, panel: "Input |X(f)|"} end) - |> Kernel.++( - Enum.zip(freqs_hz, h_spectrum) - |> Enum.map(fn {f, a} -> %{frequency: f, amplitude: a, panel: "Filter |H(f)|"} end) - ) - |> Kernel.++( - Enum.zip(freqs_hz, to_amp_list.(route_a)) - |> Enum.map(fn {f, a} -> %{frequency: f, amplitude: a, panel: "Output |X(f)·H(f)|"} end) - ) +h_padded = Nx.pad(h_fir, 0.0, [{0, n - Nx.size(h_fir), 0}]) +h_spectrum = h_padded |> Nx.rfft() |> Nx.abs() |> Nx.to_flat_list() + +input_spec = + Enum.zip_with(freqs_hz, to_amp_list.(signal), fn f, a -> + %{frequency: f, amplitude: a, panel: "Input |X(f)|"} + end) + +filter_spec = + Enum.zip_with(freqs_hz, h_spectrum, fn f, a -> + %{frequency: f, amplitude: a, panel: "Filter |H(f)|"} + end) + +output_spec = + Enum.zip_with(freqs_hz, to_amp_list.(route_a), fn f, a -> + %{frequency: f, amplitude: a, panel: "Output |X(f)·H(f)|"} + end) + +spectrum_data = input_spec ++ filter_spec ++ output_spec Tucan.lineplot(spectrum_data, "frequency", "amplitude") |> Tucan.facet_by(:row, "panel") @@ -282,10 +281,8 @@ dual_data = |> Enum.flat_map(fn {w, name} -> amps = Nx.multiply(x_win, w) - |> Nx.as_type({:c, 64}) - |> Nx.fft() + |> Nx.rfft() |> Nx.abs() - |> Nx.slice([0], [half_win + 1]) peak = Nx.reduce_max(amps) db = Nx.divide(amps, peak) @@ -293,8 +290,7 @@ dual_data = |> Nx.multiply(20) |> Nx.max(-100.0) |> Nx.to_flat_list() - Enum.zip(freqs_win, db) - |> Enum.map(fn {f, d} -> %{frequency: f, amplitude_db: d, window: name} end) + Enum.zip_with(freqs_win, db, fn f, d -> %{frequency: f, amplitude_db: d, window: name} end) end) Tucan.lineplot(dual_data, "frequency", "amplitude_db") diff --git a/guides/financial_signals.livemd b/guides/financial_signals.livemd index 569087f..9b85358 100644 --- a/guides/financial_signals.livemd +++ b/guides/financial_signals.livemd @@ -38,9 +38,7 @@ annual_cycle = Nx.multiply(0.04, Nx.sin(Nx.multiply(t, 2 * :math.pi() / 252.0)) random_walk = noise |> Nx.multiply(0.01) - |> Nx.to_flat_list() - |> Enum.scan(0.0, &(&1 + &2)) - |> Nx.tensor(type: :f32) + |> Nx.cumulative_sum(axis: 0) log_prices = Nx.add(Nx.multiply(t, trend_per_day), annual_cycle) @@ -74,8 +72,7 @@ log_returns = Nx.log(closes[1..(n - 1)]) |> Nx.subtract(Nx.log(closes[0..(n - 2)])) price_data = - Enum.zip(dates, Nx.to_flat_list(closes)) - |> Enum.map(fn {d, c} -> %{date: d, close: c} end) + Enum.zip_with(dates, Nx.to_flat_list(closes), fn d, c -> %{date: d, close: c} end) Tucan.lineplot(price_data, "date", "close") |> VegaLite.encode_field(:x, "date", type: :temporal, title: "Date") @@ -86,8 +83,9 @@ Tucan.lineplot(price_data, "date", "close") ```elixir return_data = - Enum.zip(Enum.drop(dates, 1), Nx.to_flat_list(log_returns)) - |> Enum.map(fn {d, r} -> %{date: d, return: r} end) + Enum.zip_with(Enum.drop(dates, 1), Nx.to_flat_list(log_returns), fn d, r -> + %{date: d, return: r} + end) Tucan.lineplot(return_data, "date", "return") |> VegaLite.encode_field(:x, "date", type: :temporal, title: "Date") @@ -132,8 +130,7 @@ power_half = power[1..(half - 1)] |> Enum.unzip() spectrum_data = - Enum.zip(periods_roi, powers_roi) - |> Enum.map(fn {p, pw} -> %{period: p, power: pw} end) + Enum.zip_with(periods_roi, powers_roi, fn p, pw -> %{period: p, power: pw} end) VegaLite.new(width: 680, height: 260, title: "Power Spectrum of S&P 500 Log Returns") |> VegaLite.data_from_values(spectrum_data) @@ -185,10 +182,10 @@ trim_smooth = Nx.to_flat_list(smoothed[trim..(n - trim - 1)]) # Combine raw and smoothed into a single long-form dataset chart_data_long = - Enum.zip([trim_dates, trim_close, trim_smooth]) - |> Enum.flat_map(fn {d, raw, filt} -> + Enum.zip_with([trim_dates, trim_close, trim_smooth], fn [d, raw, filt] -> [%{date: d, close: raw, series: "Raw"}, %{date: d, close: filt, series: "FIR smoothed (63-day)"}] end) + |> List.flatten() Tucan.lineplot(chart_data_long, "date", "close") |> Tucan.color_by("series") @@ -237,8 +234,7 @@ zoom_periods = end) zoom_data = - Enum.zip(zoom_periods, zoom_power) - |> Enum.map(fn {p, pw} -> %{period: p, power: pw} end) + Enum.zip_with(zoom_periods, zoom_power, fn p, pw -> %{period: p, power: pw} end) Tucan.lineplot(zoom_data, "period", "power", tooltip: true) |> VegaLite.encode_field(:x, "period", From 1eafd5bf7570247943df80b11c3b0d98a85c4d40 Mon Sep 17 00:00:00 2001 From: thomaspmurphy Date: Tue, 16 Jun 2026 14:21:21 +0100 Subject: [PATCH 3/3] fix(guides): lower audio volume, fix peak indices slice for rank-2 tensor --- guides/peak_finding.livemd | 2 +- guides/sound_synthesis.livemd | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/guides/peak_finding.livemd b/guides/peak_finding.livemd index f240e6c..9fa9c28 100644 --- a/guides/peak_finding.livemd +++ b/guides/peak_finding.livemd @@ -32,7 +32,7 @@ remainder are padded with −1. defmodule Peaks do def indices(result) do n = Nx.to_number(result.valid_indices) - if n == 0, do: [], else: result.indices |> Nx.slice([0], [n]) |> Nx.to_flat_list() + if n == 0, do: [], else: result.indices |> Nx.slice([0, 0], [n, 1]) |> Nx.to_flat_list() end end ``` diff --git a/guides/sound_synthesis.livemd b/guides/sound_synthesis.livemd index 9062297..5b4ab13 100644 --- a/guides/sound_synthesis.livemd +++ b/guides/sound_synthesis.livemd @@ -33,7 +33,7 @@ defmodule Audio do samples |> Nx.as_type(:f32) |> Nx.clip(-1.0, 1.0) - |> Nx.multiply(32_767) + |> Nx.multiply(22_000) |> Nx.as_type({:s, 16}) |> Nx.to_binary()