Brownian noise of air

harris_sivian.py

"""
Created on Wed Nov 04 11:34:53 2015
"""
from numpy import logspace, log10, sqrt, pi
import matplotlib.pyplot as plt
from audio_electronics_shortcuts import dB

import pint
u = pint.UnitRegistry()


def air_sound_speed(T):
    """
    Speed of sound in air approximation from
    https://en.wikipedia.org/wiki/Speed_of_sound#Practical_formula_for_dry_air

    Parameters
    ----------
    T : pint.Quantity
        Absolute temperature of air

    Returns
    -------
    c : pint.Quantity
        Speed of sound in dry air

    Examples
    --------
    >>> air_sound_speed(295*u.K)
    344.370058948 meter / second
    """
    c = 20.05 * u.m/u.s * sqrt(T/u.K)
    return c.to('m/s')


def air_density(p, T):
    """
    Air density

    Parameters
    ----------
    p : pint.Quantity
        Air pressure
    T : pint.Quantity
        Absolute temperature of air

    Returns
    -------
    rho : pint.Quantity
        Density of dry air

    Examples
    --------
    >>> air_density(1*u.atm, 295*u.K)
    1.19653371887 kilogram / meter ** 3
    """
    R_specific = 287.058 * u.J / u.kg / u.K  # for dry air
    rho = p / (R_specific * T)
    return rho.to('kg/m^3')


def noise(rho, T, c, f1, f2):
    """
    Harris noise

    Parameters
    ----------
    rho : pint.Quantity
        Density of dry air
    T : pint.Quantity
        Absolute temperature of air
    c : pint.Quantity
        Speed of sound in dry air
    f1, f2 : pint.Quantity
        Lower and upper frequency bounds

    Returns
    -------
    P : pint.Quantity
        RMS pressure noise

    Examples
    --------
    >>> noise(rho, T, c, f1, f2)
    0.0395772367748 micropascal
    """
    k = u.boltzmann_constant
    P_rms = sqrt((8*pi*rho*T*k) / (3*c) * (f2**3 - f1**3))
    P_rms *= 1/sqrt(2)  # -3 dB (for diaphragm or something?)
    return P_rms.to('uPa')


if __name__ == '__main__':
    T = 300 * u.K  # Air temperature
    c = air_sound_speed(T)  # Speed of sound
    p = 100 * u.kPa  # Standard atmospheric pressure
    rho = air_density(p, T)

    for f1, f2 in ((1 * u.kHz, 6 * u.kHz),
                   (2.5 * u.kHz, 3.5 * u.kHz),
                   (20 * u.Hz, 20 * u.kHz),
                   (2999.5 * u.Hz, 3000.5 * u.Hz)):

        print(f'{f1:~P} to {f2:~P}:')
        P_rms = noise(rho, T, c, f1, f2)
        print(f'{P_rms:+.3~P}')
        print(f'{dB(P_rms/(20 * u.uPa)):+.1f} dB SPL')
        print('')

    f = logspace(log10(20), log10(20000), 100) * u.Hz
    k = u.boltzmann_constant
    # TODO: use noise() function for this
    P_rms = sqrt((8*pi*rho*T*k) / (3*c) * (((f + u.Hz/2))**3 - ((f - u.Hz/2))**3))
    P_rms *= 1/sqrt(2)  # -3 dB  (for diaphragm or something?)
    P_rms.ito('uPa')

    dBSPL = dB((P_rms/(20 * u.uPa)).astype(float))
    plt.plot(f, dBSPL)
    plt.xscale('log')
    plt.grid(True, color='0.7', linestyle='-', which='both', axis='both')
    plt.xlabel('Frequency [Hz]')
    plt.ylabel('Pressure in 1 Hz band [dB SPL]')
    plt.xlim(20, 20000)

From How loud is the thermal motion of air molecules?