diff --git a/examples/miro_to_struct.m b/examples/miro_to_struct.m index 41d701b..9979c19 100644 --- a/examples/miro_to_struct.m +++ b/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; diff --git a/sound_field_analysis/io.py b/sound_field_analysis/io.py index 6c69acd..873674a 100755 --- a/sound_field_analysis/io.py +++ b/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,15 +120,15 @@ 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 ---------- @@ -136,6 +136,8 @@ def __new__(cls, signal, grid, configuration=None, temperature=None): 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): diff --git a/sound_field_analysis/process.py b/sound_field_analysis/process.py index c9c4222..9cda682 100644 --- a/sound_field_analysis/process.py +++ b/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