def phon2dB(loudness=None):
  """
  Loudness in phons to Sound Pressure Level (SPL) in dB using the
  ISO/FDIS 226:2003 model.

  This function needs Scipy, as ``scipy.interpolate.UnivariateSpline``
  objects are used as interpolators.

  Parameters
  ----------
  loudness :
    The loudness value in phons to be converted, or None (default) to get
    the threshold of hearing.

  Returns
  -------
  A callable that returns the SPL dB value for each given frequency in hertz.

  Note
  ----
  See ``phon2dB.iso226.schema`` and ``phon2dB.iso226.table`` to know the
  original frequency used for the result. The result for any other value is
  an interpolation (spline). Don't trust on values nor lower nor higher than
  the frequency limits there (20Hz and 12.5kHz) as they're not part of
  ISO226 and no value was collected to estimate them (they're just a spline
  interpolation to reach 1000dB at -30Hz and 32kHz). Likewise, the trustful
  loudness input range is from 20 to 90 phon, as written on ISO226, although
  other values aren't found by a spline interpolation but by using the
  formula on section 4.1 of ISO226.

  Hint
  ----
  The ``phon2dB.iso226.table`` also have other useful information, such as
  the threshold values in SPL dB.

  """
  from scipy.interpolate import UnivariateSpline

  table = phon2dB.iso226.table
  schema = phon2dB.iso226.schema
  freqs = [row[schema.index("freq")] for row in table]

  if loudness is None: # Threshold levels
    spl = [row[schema.index("threshold")] for row in table]

  else: # Curve for a specific phon value
    def get_pressure_level(freq, alpha, loudness_base, threshold):
      return 10 / alpha * math.log10(
        4.47e-3 * (10 ** (.025 * loudness) - 1.14) +
        (.4 * 10 ** ((threshold + loudness_base) / 10 - 9)) ** alpha
      ) - loudness_base + 94

    spl = [get_pressure_level(**dict(xzip(schema, args))) for args in table]

  interpolator = UnivariateSpline(freqs, spl, s=0)
  interpolator_low = UnivariateSpline([-30] + freqs, [1e3] + spl, s=0)
  interpolator_high = UnivariateSpline(freqs + [32000], spl + [1e3], s=0)

  @elementwise("freq", 0)
  def freq2dB_spl(freq):
    if freq < 20:
      return interpolator_low(freq).tolist()
    if freq > 12500:
      return interpolator_high(freq).tolist()
    return interpolator(freq).tolist()
  return freq2dB_spl