Skip to content

Apply the transfer function at the right FFT bins in amplify_motion() - #118

Open
jsh9 wants to merge 2 commits into
2026-10-10-parametrize-testsfrom
2026-10-11-fix-transfer-function-frequency-bins
Open

jsh9 wants to merge 2 commits into
2026-10-10-parametrize-testsfrom
2026-10-11-fix-transfer-function-frequency-bins

Conversation

@jsh9

@jsh9 jsh9 commented Oct 11, 2026

Copy link
Copy Markdown
Collaborator

Fixes #78. One of the four P0 fixes from #114, based on #74's branch (2026-10-10-parametrize-tests).

Problem

_get_freq_interval() built its frequency array as np.linspace(df, fmax, num=half_n), so element k was (k + 1)·df. amplify_motion() uses element k as the frequency of FFT bin k, where bin 0 is 0 Hz. So the transfer function was applied one bin too high.

  • A 1.0 Hz cosine through |TF| = f came out with amplitude 1.1 instead of 1.0.
  • This affects linear_site_resp(), Ground_Motion.amplify(), deconvolve(), amplify_by_tf() and Site_Effect_Adjustment.run().
  • Against an exact frequency-domain solution, linear_site_resp() had correlation 0.99791 (elastic) and 0.96969 (rigid). With this fix it is 1.00000 and 0.99992.

Changes

All in PySeismoSoil/helper_site_response.py:

  • _get_freq_interval(): f_array = np.arange(half_n) * df, so element k is the frequency of FFT bin k. This is the fix proposed in the issue.
  • fmax (beyond the issue): it is now (half_n - 1) * df, the last element of f_array (the Nyquist frequency when n is even). Before, it was one bin above that.
    • With extrap_tf=False, amplify_motion() halves the input until the transfer function reaches fmax. With the old value, a transfer function covering exactly 0–50 Hz for dt = 0.01 s still halved the output's sampling rate, although the docstring says that is enough.
    • With extrap_tf=True, the extrapolated point now falls exactly on the last bin.
    • The return-value docstrings are updated.
  • 0 Hz: the value used at 0 Hz is the same as before in every case, so tf_ss[0] = np.real(tf_ss[0]) is still right. The comments there and "Note (1)" are updated to describe how np.interp() fills 0 Hz.
  • show_fig plot: it now skips the 0 Hz point, which the logarithmic frequency axes cannot show.
  • CHANGELOG.md: a "Fixed" entry.
  • The example notebooks were re-run (second commit), as CONTRIBUTING.md requires. Only timings, folder names and object addresses changed in the printed outputs. The figures show the corrected results.

Tests

11 new or changed tests, all of which fail on the old code:

  • tests/test_helper_site_response.py:
    • test_amplify_motion__cosine_at_an_FFT_bin[amplify|deconvolve × even|odd length]: a cosine at FFT bin 10 through |TF| = 1 + f (zero phase) comes out multiplied (or divided) by |TF| at its frequency. Old code: 2.1 instead of 2.0.
    • test_amplify_motion__tf_up_to_the_Nyquist_frequency[even|odd]: with extrap_tf=False, a 0–50 Hz transfer function causes no downsampling. Old code: (500, 2) == (1000, 2) fails.
    • test_get_freq_interval[even|odd]: f_array equals np.fft.rfftfreq(n, dt), and fmax is its last value.
  • tests/test_helper_simulations.py::test_linear[elastic|rigid]:
    • The thresholds were 0.99 and 0.97; both are now 0.999. Actual correlations: 0.998177 / 0.972972 before, 0.999991 / 0.999860 now.
    • The comment that blamed sim.linear() for the difference is corrected.
  • tests/test_class_simulation.py::test_linear: the threshold goes from 0.99 to 0.999 (0.998177 before, 0.999991 now).

No benchmark values changed. The deconvolution round-trip tests passed before and still pass, because the shift cancels out there.

Testing

CI runs only on PRs into main, so it does not run here. These were run locally on Python 3.13:

  • python -m pytest tests: 249 passed.
  • pre-commit run -a and pydoclint PySeismoSoil: pass.
  • tox -e run-notebooks: all 14 notebooks ran without errors.

Notes for review

Two related off-by-one problems were found while fixing this, and are not changed here:

  • helper_signal_processing.fourier_transform() labels FFT bin 0 (0 Hz) as df (freq_array = np.arange(1, ...) / (N * dt)), so every frequency label is one bin too high. It feeds Ground_Motion.get_Fourier_spectrum(), calc_transfer_function(), compare_two_accel() and the GoF scores.
  • helper_simulations.linear() evaluates bin k at k·df + df/15, which is why it reaches only 0.99999 rather than exactly 1 against the exact solution.

The other P0 fixes (#77 in #117; #79 in #116; #60 in #115) are in separate PRs on the same base and change different library code. All of them re-run the notebooks and add a "Fixed" group to CHANGELOG.md, so whichever merges second will conflict there. To resolve: keep both entries and re-run the notebooks.

Because the base is not main, merging this PR will not close #78 automatically.

🤖 Generated with Claude Code

https://claude.ai/code/session_01W7fK8T2pPLf3tPCnMCKL1T


Generated by Claude Code

- `_get_freq_interval()` built its frequency array as
  `np.linspace(df, fmax, num=half_n)`, so element k was (k + 1) * df, but
  `amplify_motion()` uses element k as the frequency of FFT bin k (bin 0 is
  0 Hz). The transfer function was applied one bin too high. It now uses
  `np.arange(half_n) * df` (#78)
- `fmax` is now the last element of that array (the highest FFT bin, i.e.,
  the Nyquist frequency for an even length), instead of one bin above it.
  So with `extrap_tf=False`, a transfer function that reaches the Nyquist
  frequency no longer makes `amplify_motion()` downsample the motion
- The amplification and phase plots of `amplify_motion()` skip 0 Hz, which
  a logarithmic axis cannot show; update the comments on the 0 Hz value
- Add tests: a cosine at an FFT bin through a frequency-dependent transfer
  function (amplify and deconvolve), a transfer function up to the Nyquist
  frequency without `extrap_tf`, and the frequency array itself
- `sim.linear()` and `linear_site_resp()` now agree to a correlation of
  0.99999 (elastic) and 0.99986 (rigid), so raise the thresholds of the
  tests that compare them to 0.999 and fix the comment that blamed
  `sim.linear()` for the difference

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01W7fK8T2pPLf3tPCnMCKL1T
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01W7fK8T2pPLf3tPCnMCKL1T

This branch has not been deployed

No deployments
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants