amplify_motion() applies the transfer function one frequency bin too high, because _get_freq_interval() starts its frequency array at df instead of 0 #78

import numpy as np
import PySeismoSoil.helper_simulations as sim
import PySeismoSoil.helper_site_response as sr

# 1. A 1.0 Hz cosine (exactly FFT bin 10) through |TF(f)| = f, zero phase
dt, n = 0.01, 1000  # df = 0.1 Hz
t = np.arange(n) * dt
x = np.cos(2 * np.pi * 1.0 * t)
f = np.linspace(0, 60, 6001)
y = sr.amplify_motion(np.column_stack((t, x)), (f, (f.copy(), np.zeros_like(f))))[0]
print(round(np.abs(y[:, 1]).max(), 4))  # 1.1    (expected 1.0 = |TF(1.0 Hz)|)

print(sr._get_freq_interval(np.column_stack((t, x)))[4][:3])
# [0.1 0.2 0.3]    (expected [0.  0.1 0.2], the frequencies of FFT bins 0, 1, 2)

# 2. Both linear solvers vs. the exact frequency-domain solution
vs = np.genfromtxt('tests/files/profile_FKSH14.txt')
acc = np.genfromtxt('tests/files/sample_accel.txt')
n = len(acc)
df = 1 / (n * (acc[1, 0] - acc[0, 0]))
h, Vs, xi, rho = vs[:, 0], vs[:, 1], vs[:, 2], vs[:, 3]
vs_star = Vs * np.sqrt(1 + 2j * xi)
alpha = rho[:-1] * vs_star[:-1] / (rho[1:] * vs_star[1:])

def exact_tf(freqs, boundary):  # surface / input, same formulas as linear_tf()
    out = []
    for fq in freqs:
        k = 2 * np.pi * fq / vs_star
        A = B = 1 + 0j
        for j in range(len(h) - 1):
            e = np.exp(1j * k[j] * h[j])
            A, B = (0.5 * A * (1 + alpha[j]) * e + 0.5 * B * (1 - alpha[j]) / e,
                    0.5 * A * (1 - alpha[j]) * e + 0.5 * B * (1 + alpha[j]) / e)
        out.append(1 / A if boundary == 'elastic' else 2 / (A + B))
    return np.array(out)

FFT = np.fft.fft(acc[:, 1])
f_bins = np.arange((n + 1) // 2) * df  # n is odd
for boundary in ['elastic', 'rigid']:
    H = exact_tf(f_bins, boundary)
    H[0] = 1.0
    exact = np.real(np.fft.ifft(FFT * np.concatenate([H, np.conj(H[1:][::-1])])))
    for name, r in [
        ('linear_site_resp', sr.linear_site_resp(vs, acc, boundary=boundary)[0][:, 1]),
        ('sim.linear', sim.linear(vs, acc, boundary=boundary)[3][:, 1]),
    ]:
        corr = np.corrcoef(r, exact)[0, 1]
        peak = np.abs(r).max() / np.abs(exact).max()
        print(f'{boundary:8s} {name:17s} corr={corr:.5f} peak ratio={peak:.4f}')
# elastic  linear_site_resp  corr=0.99791 peak ratio=1.0149
# elastic  sim.linear        corr=0.99999 peak ratio=1.0005
# rigid    linear_site_resp  corr=0.96969 peak ratio=1.0472
# rigid    sim.linear        corr=0.99988 peak ratio=1.0027