Implementation of BEMA

Edit miro_to_struct to contain center measurement
Edit read_miro_struct to optionally load the center measurement
Edit ArraySignal named tuple to optionally contain center measurement

examples/miro_to_struct.m

@@ -25,6 +25,7 @@
 quadWeight = miro_data.quadWeight;
 scatterer  = miro_data.scatterer;
 avgAirTemp = miro_data.avgAirTemp; 
+irCenter   = miro_data.irCenter;
 
 clear miro_data path_name file_name;
 
@@ -53,7 +54,8 @@
 radius     = miro_data.radius;
 quadWeight = miro_data.quadWeight;
 scatterer  = miro_data.scatterer;
-avgAirTemp = miro_data.avgAirTemp; 
+avgAirTemp = miro_data.avgAirTemp;
+irCenter   = miro_data.irCenter;
 
 clear miro_data path_name file_name;
 
@@ -82,7 +84,8 @@
 radius     = miro_data.radius;
 quadWeight = miro_data.quadWeight;
 scatterer  = miro_data.scatterer;
-avgAirTemp = miro_data.avgAirTemp; 
+avgAirTemp = miro_data.avgAirTemp;
+irCenter   = miro_data.irCenter;
 
 clear miro_data path_name file_name;
 

sound_field_analysis/io.py

@@ -48,9 +48,9 @@ def __new__(cls, array_radius, array_type, transducer_type, scatter_radius=None,
 
     def __repr__(self):
         return 'ArrayConfiguration(\n' + ',\n'.join(
-                '    {0} = {1}'.format(name, repr(data).replace('\n', '\n      '))
-                for name, data in
-                zip(['array_radius', 'array_type', 'transducer_type', 'scatter_radius', 'dual_radius'], self)) + ')'
+            '    {0} = {1}'.format(name, repr(data).replace('\n', '\n      '))
+            for name, data in
+            zip(['array_radius', 'array_type', 'transducer_type', 'scatter_radius', 'dual_radius'], self)) + ')'
 
 
 class TimeSignal(namedtuple('TimeSignal', 'signal fs delay')):
@@ -84,8 +84,8 @@ def __new__(cls, signal, fs, delay=None):
 
     def __repr__(self):
         return 'TimeSignal(\n' + ',\n'.join(
-                '    {0} = {1}'.format(name, repr(data).replace('\n', '\n      '))
-                for name, data in zip(['signal', 'fs', 'delay'], self)) + ')'
+            '    {0} = {1}'.format(name, repr(data).replace('\n', '\n      '))
+            for name, data in zip(['signal', 'fs', 'delay'], self)) + ')'
 
 
 class SphericalGrid(namedtuple('SphericalGrid', 'azimuth colatitude radius weight')):
@@ -120,22 +120,24 @@ def __new__(cls, azimuth, colatitude, radius=None, weight=None):
 
     def __repr__(self):
         return 'SphericalGrid(\n' + ',\n'.join(
-                '    {0} = {1}'.format(name, repr(data).replace('\n', '\n      '))
-                for name, data in zip(['azimuth', 'colatitude', 'radius', 'weight'], self)) + ')'
+            '    {0} = {1}'.format(name, repr(data).replace('\n', '\n      '))
+            for name, data in zip(['azimuth', 'colatitude', 'radius', 'weight'], self)) + ')'
 
 
-class ArraySignal(namedtuple('ArraySignal', 'signal grid configuration temperature')):
+class ArraySignal(namedtuple('ArraySignal', 'signal grid center_signal configuration temperature')):
     """Named tuple ArraySignal"""
     __slots__ = ()
 
-    def __new__(cls, signal, grid, configuration=None, temperature=None):
+    def __new__(cls, signal, grid, center_signal=None, configuration=None, temperature=None):
         """
         Parameters
         ----------
         signal : TimeSignal
            Holds time domain signals and sampling frequency fs
         grid : SphericalGrid
            Location grid of all time domain signals
+        center_signal : TimeSignal
+           Center signal measurement
         configuration : ArrayConfiguration
            Information on array configuration
         temperature : array_like, optional
@@ -147,43 +149,61 @@ def __new__(cls, signal, grid, configuration=None, temperature=None):
             configuration = ArrayConfiguration(*configuration)
 
         # noinspection PyArgumentList
-        self = super(ArraySignal, cls).__new__(cls, signal, grid, configuration, temperature)
+        self = super(ArraySignal, cls).__new__(cls, signal, grid, center_signal, configuration, temperature)
         return self
 
     def __repr__(self):
         return 'ArraySignal(\n' + ',\n'.join(
-                '    {0} = {1}'.format(name, repr(data).replace('\n', '\n      '))
-                for name, data in zip(['signal', 'grid', 'configuration', 'temperature'], self)) + ')'
+            '    {0} = {1}'.format(name, repr(data).replace('\n', '\n      '))
+            for name, data in zip(['signal', 'grid', 'center_signal', 'configuration', 'temperature'], self)) + ')'
 
 
-def read_miro_struct(file_name, channel='irChOne', transducer_type='omni', scatter_radius=None):
+def read_miro_struct(file_name, channel='irChOne', transducer_type='omni', scatter_radius=None,
+                     get_center_signal=False):
     """ Reads miro matlab files.
 
-    Note: This function expects a slightly modified miro file in that it expects a field `colatitude` instead of
-    `elevation`. This is for avoiding confusion as may miro file contain colatitude data in the elevation field.
-
     Parameters
     ----------
     file_name : filepath
        Path to file that has been exported as a struct
     channel : string, optional
-       Channel that holds required signals. Default: 'irChOne'
+       Channel that holds required signals. [Default: 'irChOne']
     transducer_type : {omni, cardoid}, optional
-       Sets the type of transducer used in the recording. Default: omni
+       Sets the type of transducer used in the recording. [Default: 'omni']
     scatter_radius : float, option
-       Radius of the scatterer. Default: None
+       Radius of the scatterer. [Default: None]
+    get_center_signal : bool, optional
+        If center signal should be loaded. [Default: False]
 
     Returns
     -------
     array_signal : ArraySignal
-       Tuple containing a TimeSignal `signal`, SphericalGrid `grid`, ArrayConfiguration `configuration` and the air
-       temperature
+       Tuple containing a TimeSignal `signal`, SphericalGrid `grid`, TimeSignal 'center_signal',
+       ArrayConfiguration `configuration` and the air temperature
+
+    Notes
+    -----
+    This function expects a slightly modified miro file in that it expects a field `colatitude` instead of
+    `elevation`. This is for avoiding confusion as may miro file contain colatitude data in the elevation field.
+
+    To import center signal measurements the matlab method miro_to_struct has to be extended. Center measurements are
+    included in every measurement provided at http://audiogroup.web.th-koeln.de/.
+
     """
     current_data = sio.loadmat(file_name)
 
     time_signal = TimeSignal(signal=_np.squeeze(current_data[channel]).T,
                              fs=_np.squeeze(current_data['fs']))
 
+    center_signal = None
+    if get_center_signal:
+        try:
+            center_signal = TimeSignal(signal=_np.squeeze(current_data['irCenter']).T,
+                                       fs=_np.squeeze(current_data['fs']))
+        except KeyError:
+            print('WARNING: Center signal not included in miro struct, use extended miro_to_struct.m!')
+            center_signal = None
+
     mic_grid = SphericalGrid(azimuth=_np.squeeze(current_data['azimuth']),
                              colatitude=_np.squeeze(current_data['colatitude']),
                              radius=_np.squeeze(current_data['radius']),
@@ -199,7 +219,7 @@ def read_miro_struct(file_name, channel='irChOne', transducer_type='omni', scatt
         sphere_config = 'open'
     array_config = ArrayConfiguration(mic_grid.radius, sphere_config, transducer_type, scatter_radius)
 
-    return ArraySignal(time_signal, mic_grid, array_config, _np.squeeze(current_data['avgAirTemp']))
+    return ArraySignal(time_signal, mic_grid, center_signal, array_config, _np.squeeze(current_data['avgAirTemp']))
 
 
 def empty_time_signal(no_of_signals, signal_length):

sound_field_analysis/process.py

@@ -1,5 +1,7 @@
 """
 Functions that act on the Spatial Fourier Coefficients
+`BEMA`
+   Bandwidth Extension for Microphone Arrays
 
 `FFT`
    (Fast) Fourier Transform
@@ -31,39 +33,93 @@
 
 
 # noinspection PyUnusedLocal
-def BEMA(Pnm, ctSig, dn, transition, avgBandwidth, fade=True):
-    """BEMA Spatial Anti-Aliasing - NOT YET IMPLEMENTED
+def BEMA(Pnm, center_sig, dn, transition, avg_band_width=1, fade=True, max_order=None):
+    """BEMA Spatial Anti-Aliasing
 
     Parameters
     ----------
     Pnm : array_like
-       Spatial Fourier coefficients
-    ctSig : array_like
-       Signal of the center microphone
+       Sound field SH spatial Fourier coefficients
+    center_sig : array_like
+      center microphone in shape [0, NFFT]
     dn : array_like
        Radial filters for the current array configuration
     transition : int
-       Highest stable bin, approx: transition = (NFFT/FS+1) * (N*c)/(2*pi*r)
-    avgBandwidth : int
-       Averaging Bandwidth in oct
+       Highest stable bin, approx: transition = (NFFT/fs) * (N*c)/(2*pi*r)
+    avg_band_width : int, optional
+       Averaging Bandwidth in oct [Default: 1]
     fade : bool, optional
        Fade over if True, else hard cut [Default: True]
+    max_order : int, optional
+       Maximum transform order [Default: highest available order]
 
     Returns
     -------
     Pnm : array_like
-       Alias-free spatial Fourier coefficients
-
-    Note
-    ----
-    This was presented at the 2012 AES convention, see [1]_.
+       BEMA optimized sound field SH coefficients
 
     References
     ----------
     .. [1] B. Bernschütz, "Bandwidth Extension for Microphone Arrays",
        AES Convention 2012, Convention Paper 8751, 2012. http://www.aes.org/e-lib/browse.cfm?elib=16493
     """
-    print('!Warning, BEMA is not yet implemented. Continuing with initial coefficients!')
+
+    if not max_order:
+        max_order = int(_np.sqrt(Pnm.shape[0] - 1))  # SH order
+
+    transition = _np.floor(transition)
+
+    # computing spatiotemporal Image Imn
+    Imn = _np.zeros((Pnm.shape[0], 1), dtype=complex)  # preallocating the spectral image matrix
+    start_avg_bin = _np.floor(transition / (_np.power(2, avg_band_width)))  # first bin for averaging
+
+    modeCnt = 0
+    avgPower = 0
+    for n in range(0, max_order + 1):  # sh orders
+        for m in range(0, 2 * n + 1):  # modes
+            for bin in range(start_avg_bin - 1, transition):  # bins
+                synthBin = Pnm[modeCnt, bin] * dn[n, bin]
+                Imn[modeCnt, 0] += synthBin
+                avgPower += _np.abs(synthBin) ** 2
+            modeCnt += 1
+
+    # energy normalization (Eq. 9)
+    energySum = _np.sum(_np.power(_np.abs(Imn), 2))
+    normFactor = avgPower / (transition - start_avg_bin + 1)
+    Imn *= _np.sqrt(normFactor / energySum)
+
+    # normalize center signal to its own average level in the source band. (Eq. 13)
+    center_sig[_np.where(center_sig == 0)] = 1e-12
+    sq_avg = _np.sqrt(_np.mean(_np.power(_np.abs(center_sig[:, (start_avg_bin - 1):transition]), 2)))
+    center_sig = _np.multiply(center_sig, (1 / sq_avg))
+
+    Pnm_synth = _np.zeros(Pnm.shape, dtype=complex)
+    modeCnt = 0
+    for n in range(0, max_order + 1):
+        for m in range(0, 2 * n + 1):
+            for bin in range(start_avg_bin - 1, Pnm_synth.shape[1]):
+                Pnm_synth[modeCnt, bin] = Imn[modeCnt, 0] * (1 / dn[n, bin]) * center_sig[0, bin]
+            modeCnt += 1
+
+    # Phase correction (Eq. 16)
+    phaseOffset = _np.angle(Pnm[0, transition] / Pnm_synth[0, transition])
+    Pnm_synth *= _np.exp(1j * phaseOffset)
+
+    # Merge original and synthetic SH coefficients. (Eq. 14)
+    Pnm = _np.concatenate((Pnm[:, 0:(transition - 1)], Pnm_synth[:, (transition - 1):]), axis=1)
+
+    # Fade in of synthetic SH coefficients. (Eq. 17)
+    if fade:
+        Pnm_fade0 = Pnm[:, (start_avg_bin - 1):(transition - 1)]
+        Pnm_fadeS = Pnm_synth[:, (start_avg_bin - 1):(transition - 1)]
+
+        fadeUp = _np.linspace(0, 1, Pnm_fade0.shape[1])
+        fadeDown = _np.flip(fadeUp)
+
+        Pnm_fade0 *= fadeDown
+        Pnm_fadeS *= fadeUp
+        Pnm[:, start_avg_bin - 1:transition - 1] = Pnm_fade0 + Pnm_fadeS
+
     return Pnm