lasp_imptube.py

Go to the documentation of this file.

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""!
Author: J.A. de Jong - ASCEE
 
Description: Two-microphone impedance tube methods
"""
__all__ = ['TwoMicImpedanceTube']
# from lrftubes import Air
from .lasp_measurement import Measurement
from numpy import pi, sqrt, exp
import numpy as np
from scipy.interpolate import UnivariateSpline
# from lrftubes import PrsDuct
from functools import lru_cache
 
class TwoMicImpedanceTube:
    def __init__(self, mnormal: Measurement,
                       mswitched: Measurement,
                 s: float,
                 d1: float,
                 d2: float,
                 fl: float = None,
                 fu: float = None,
                 periodic_method=False,
                 # mat= Air(),
                 D_imptube = 50e-3,
                 thermoviscous = True,
                 **kwargs):
 
        """
        Initialize two-microphone impedance tube methods
 
        Args:
            mnormal: Measurement in normal configuration
            mswitched: Measurement in normal configuration
            s: Microphone distance
            fl: Lower evaluation frequency
            fu: Lower evaluation frequency
 
            kwargs: tuple with extra arguments, of which:
                N: Period length of periodic excitation *obligatory*
                chan0: Measurement channel index of mic 0
                chan1: Measurement channel index of mic 1
 
        """
        self.mnormal = mnormal
        self.mswitched = mswitched
 
        self.mat = mat
        self.s = s
        self.d1 = d1
        self.d2 = d2
 
        self.fl = fl
        if fl is None:
            ksmin = 0.1*pi
            kmin = ksmin/s
            self.fl = kmin*mat.c0/2/pi
 
        self.fu = fu
        if fu is None:
            ksmax = 0.8*pi
            kmax = ksmax/s
            self.fu = kmax*mat.c0/2/pi
 
        self.thermoviscous = thermoviscous
        self.D_imptube = D_imptube
 
        self.periodic_method = periodic_method
        self.channels = [kwargs.pop('chan0', 0), kwargs.pop('chan1', 1)]
        # Compute calibration correction
        if periodic_method:
            self.N = kwargs.pop('N')
            freq, C1 = mnormal.periodicCPS(self.N, channels=self.channels)
            freq, C2 = mswitched.periodicCPS(self.N, channels=self.channels)
        else:
            self.nfft = kwargs.pop('nfft', 16000)
            freq, C1 = mnormal.CPS(nfft=self.nfft, channels=self.channels)
            freq, C2 = mswitched.CPS(nfft=self.nfft, channels=self.channels)
 
        # Store this, as it is often used for just a single sample.
        self.C1 = C1 
        self.freq = freq
 
        self.il = np.where(self.freq<= self.fl)[0][-1]
        self.ul = np.where(self.freq > self.fu)[0][0]
 
        # Calibration correction factor
        # self.K = 0*self.freq + 1.0
        K = sqrt(C2[:,0,1]*C1[:,0,0]/(C2[:,1,1]*C1[:,1,0]))
        # self.K = UnivariateSpline(self.freq, K.real)(self.freq) +\
        #         1j*UnivariateSpline(self.freq, K.imag)(self.freq)
        self.K = K
 
    def cut_to_limits(self, ar):
        return ar[self.il:self.ul]
 
    def getFreq(self):
        """
        Array of frequencies, cut to limits of validity
        """
        return self.cut_to_limits(self.freq)
    
    @lru_cache
    def G_AB(self, meas):
        if meas is self.mnormal:
            C = self.C1
            freq = self.freq
        else:
            if self.periodic_method:
                freq, C = meas.periodicCPS(self.N, self.channels)
            else:
                freq, C = meas.CPS(nfft=self.nfft, channels=self.channels)
 
        # Microphone transfer function
        G_AB = self.K*C[:,1,0]/C[:,0,0]
        return self.getFreq(), self.cut_to_limits(G_AB)
    
    def k(self, freq):
        """
        Wave number, or thermoviscous wave number
        """
        if self.thermoviscous:
            D = self.D_imptube
            S = pi/4*D**2
            d = PrsDuct(0, S=S, rh=D/4, cs='circ')
            d.mat = self.mat
            omg = 2*pi*freq
            G, Z = d.GammaZc(omg)
            return G
        else:
            return 2*pi*freq/self.mat.c0
 
 
    def R(self, meas):
        freq, G_AB = self.G_AB(meas)
        s = self.s
        k = self.k(freq)
        d1 = self.d1
        RpA = (G_AB - exp(-1j*k*s))/(exp(1j*k*s)-G_AB)
 
        R = RpA*exp(2*1j*k*(s+d1))
 
        return freq, R
 
    def alpha(self, meas):
        """
        Acoustic absorption coefficient
        """
        freq, R = self.R(meas)
        return freq, 1 - np.abs(R)**2
    
    def z(self, meas):
        """
        Acoustic impedance at the position of the sample, in front of the sample
        """
        freq, R = self.R(meas)
        return freq, self.mat.z0*(1+R)/(1-R)
 
    def zs(self, meas):
        """
        Sample impedance jump, assuming a cavity behind the sample with
        thickness d2
        """
        freq, R = self.R(meas)
        z0 = self.mat.z0
        k = 2*pi*freq/self.mat.c0
        d2 = self.d2
 
        zs = 2*z0*(1-R*exp(2*1j*k*d2))/((R-1)*(exp(2*d2*k*1j)-1))
        return freq, zs