import os
import numpy
import sofar as sf #
https://pypi.org/project/sofar/ pip install sofar
import subprocess
import shutil
# ---------------------------------------- User Settings: --------------------------------------------
Input1 = r"C:\mesh2hrtf-tools\SOFA_Merge_Folder\left"
Input2 = r"C:\mesh2hrtf-tools\SOFA_Merge_Folder\right"
# Input1 = r"C:\hrtf\SOFA_Merge_Folder" # example (use r"..." raw strings without "\" in the end!)
SavePath = os.path.dirname(os.path.realpath(__file__)) # Automatically detect the working directory
# Result files will be saved in a sub-folder created in "SavePath"
# SavePath = r"C:\hrtf\SOFA_Merge_Folder" # alternatively specify a fixed folder (use r"..." raw strings!)
Convert_to_spherical = True # keep this True, unless you have very specific demand for cartesian coordinates.
# Some software only works with "spherical" coordinates of "SourcePosition" "in SOFA files
Diagnostic_plot = True # shows and saves to disk the plot of the merged HRIR.sofa file: great for visual quality check.
# requires: --- pyfar --- and --- matplotlib ---
# ---------------------------------------- End of User Settings -----------------------------------------
#
# "join_sofa_files.py" merges SOFA data from two separate Mesh2HRTF simulations for LEFT & RIGHT ear.
# it uses a lot of example code from Mesh2HRTF "Output2HRTF_Main.py" & "NumCalcManager.py"
#
# HOW TO USE:
# Main tutorial is here:
https://sourceforge.net/p/mesh2hrtf/wiki/Basic_HRTF_post_processing
#
# mode A0 - no inputs = (recommended usage mode) same as A1, but actually executes "Output2HRTF.py" as well to run full
# pre-processing in one go and applies fixes to the possible "Non-convergence issue".
#
# mode A1 - no inputs = scan start folder for 2 projects to merge.
# 1- run this .py file from a dedicated folder (for example use folder "SOFA_Merge_Folder")
# 2- Move into the "SOFA_Merge_Folder" exactly the 2 project folders that need to be merged.
# 3- "open/run" this "join_sofa_files.py" file with Python. (Python3 with "sofar" & "numpy" must be installed)
# (an easy Drag&Drop way to run this script is from Spyder which can be installed by Anaconda)
# Without input arguments this script automatically merges the SOFA files of the 2 projects it finds.
# 4- Merged SOFA files will be saved in a folder next to the project that was found. DONE.
#
# mode A2 - Input1 only = scan Input1 folder for 2 projects to merge.
# 1- Move into any folder exactly the 2 project folders that need to be merged.
# 2- "run" this "join_sofa_files.py" file and specify "Input1" path to the folder that contains projects to merge.
# This script searches and merges the SOFA files from the 2 projects it finds inside "input1" path.
# 3- Merged SOFA files will be saved in a folder next to the project that was found. DONE.
#
# NOT_IMPLEMENTED: mode B - 1 input = Use the given path to identify a pair of projects to merge
# where folder-name for LEFT ends on '_L' and RIGHT ends on '_R'.
#
# mode C - 2 inputs = Merge the 2 projects that were given as input.
#
# made for Python 3.7+ see license details at the end of the script.
# v0.50 2022-01-23 First working version. Tested on Windows 10. (by Sergejs Dombrovskis)
# v0.90 2022-01-24 Mode-A1 is ready + most descriptions in place. (by S.D.)
# v0.93 2022-02-06 added Mode-A2 and Mode-C. (by S.D.)
# v1.00 2022-03-15 added Mode-A0 with automatic execution of "Output2HRTF.py" + small improvements +
# added multi-SamplingRate HRIR outputs + extra sanity checks (by S.D.)
# v1.10 2022-03-17 added diagnostic plotting of the merged HRIR.sofa file (by S.D.)
# v1.20 2022-03-20 added effective workaround to produce usable data even if simulation encountered the
# "Non-Convergence issue" - only in A0 mode (by S.D.)
#
def merge_write_sofa(sofa_left, sofa_right, basepath, filename, sofa_type='HRTF'):
"""Write complex pressure or impulse responses as SOFA file."""
# adapted from "write_to_sofa()" in "Output2HRTF_Main.py"
# create empty SOFA object
if sofa_type == 'HRTF':
convention = "SimpleFreeFieldHRTF" # if numSources == 2 else "GeneralTF"
elif sofa_type == 'HRIR':
convention = "SimpleFreeFieldHRIR" # if numSources == 2 else "GeneralFIR"
else:
convention = 'Error - Unrecognized sofa_type provided.'
sofa = sf.Sofa(convention)
# write meta data
sofa.GLOBAL_ApplicationName = sofa_left.GLOBAL_ApplicationName # 'Mesh2HRTF'
sofa.GLOBAL_ApplicationVersion = sofa_left.GLOBAL_ApplicationVersion
sofa.GLOBAL_History = sofa_left.GLOBAL_History # "numerically simulated data"
if len(sofa_left.GLOBAL_Title) == 0: # title added to improve compatibility with Matlab SOFA API 1.1.3
sofa.GLOBAL_Title = 'untitled ' + sofa_type + ' SOFA data'
else:
sofa.GLOBAL_Title = sofa_left.GLOBAL_Title
# Source and receiver data
if Convert_to_spherical and \
sofa_left.SourcePosition_Type == 'cartesian' and sofa_left.SourcePosition_Units == 'meter':
sofa.SourcePosition = cart2sph_in_deg(sofa_left.SourcePosition) # see sub-function
sofa.SourcePosition_Units = 'degree, degree, metre' # alternatively 'degrees'
sofa.SourcePosition_Type = 'spherical'
else:
sofa.SourcePosition = sofa_left.SourcePosition
sofa.SourcePosition_Units = sofa_left.SourcePosition_Units # "meter"
sofa.SourcePosition_Type = sofa_left.SourcePosition_Type # "cartesian"
sofa.ReceiverPosition = [sofa_left.ReceiverPosition.reshape(3, 1),
sofa_right.ReceiverPosition.reshape(3, 1)]
sofa.ReceiverPosition_Units = sofa_left.ReceiverPosition_Units # "meter"
sofa.ReceiverPosition_Type = sofa_left.ReceiverPosition_Type # "cartesian"
if sofa_type == 'HRTF':
tmp_data_real = numpy.zeros((sofa_left.Data_Real.shape[0], 2,
sofa_left.Data_Real.shape[2])) # init
tmp_data_imag = tmp_data_real # init
for jj in range(sofa_left.Data_Real.shape[2]): # merge arrays
tmp_data_real[:, 0, jj] = sofa_left.Data_Real[:, 0, jj]
tmp_data_real[:, 1, jj] = sofa_right.Data_Real[:, 0, jj]
tmp_data_imag[:, 0, jj] = sofa_left.Data_Imag[:, 0, jj]
tmp_data_imag[:, 1, jj] = sofa_right.Data_Imag[:, 0, jj]
sofa.Data_Real = tmp_data_real
sofa.Data_Imag = tmp_data_imag
sofa.N = sofa_left.N # frequencies
# init defaults for HRTF:
output_file_sampling_rates = [int(2 * sofa.N[-1])]
elif sofa_type == 'HRIR':
tmp_data_ir = numpy.zeros((sofa_left.Data_IR.shape[0], 2, sofa_left.Data_IR.shape[2])) # init
for jj in range(sofa_left.Data_IR.shape[2]): # merge arrays
tmp_data_ir[:, 0, jj] = sofa_left.Data_IR[:, 0, jj]
tmp_data_ir[:, 1, jj] = sofa_right.Data_IR[:, 0, jj]
sofa.Data_IR = tmp_data_ir
sofa.Data_SamplingRate = sofa_left.Data_SamplingRate # fs
sofa.Data_Delay = numpy.ones((1, 2)) * sofa_left.Data_Delay
# init defaults for HRTF:
output_file_sampling_rates = [int(sofa.Data_SamplingRate)]
output_data_length = [sofa.Data_IR.shape[2]]
# extra "multi-sample rate HRIR" output processing
if numpy.round(sofa.Data_SamplingRate / sofa.Data_IR.shape[2]) == 150:
frequency_step = 150 # detected a valid multi-sample rate frequency step
elif numpy.round(sofa.Data_SamplingRate / sofa.Data_IR.shape[2]) == 75:
frequency_step = 75 # detected a valid multi-sample rate frequency step
else:
print('\n --- Note: the Frequency Step of this dataset is not equal to 150 or 75 Hz, therefore ')
print(' "multi-SamplingRate" HRIR saving is not possible. The only available SamplingRate is ' +
str(sofa.Data_SamplingRate) + ' ' + sofa.Data_SamplingRate_Units + '\n')
frequency_step = 0
if frequency_step > 0.5: # "multi-sample rate HRIR" output is possible
common_sampling_rates = [192000, 96000, 88200, 48000, 44100] # edit this DESCENDING list if you have to
for rate in common_sampling_rates:
if sofa.Data_IR.shape[2] > rate / frequency_step:
# add this sampling rate to output:
output_file_sampling_rates.append(rate)
output_data_length.append(int(rate / frequency_step))
else:
raise Exception("impossible
")
# save the merged files
for f_nr in range(len(output_file_sampling_rates)):
if f_nr > 0: # extra "multi-SamplingRate HRIR" output processing
# trim data down to correct sampling rate:
# noinspection PyUnboundLocalVariable
sofa.Data_IR = sofa.Data_IR[:, :,
utput_data_length[f_nr]]
sofa.Data_SamplingRate = float(output_file_sampling_rates[f_nr])
# https://sofa_right.readthedocs.io/en/latest/sofa_right.html#sofa_right-documentation
sofa.verify()
# Save as SOFA file
path = os.path.join(basepath, filename[:-5] + '_' + str(output_file_sampling_rates[f_nr]) + '.sofa')
# Need to delete old FILE first if file already exists
if os.path.exists(path):
os.remove(path)
# Save
sf.write_sofa(path, sofa, version='latest', compression=9) # max compression
if f_nr == 0:
full_reference_path = path
# noinspection PyUnboundLocalVariable
return full_reference_path
# cart2sph conversion with DEGREE output
def cart2sph_in_deg(xyz):
theta_phi_r = numpy.empty(xyz.shape)
# based on pyfar example
https://github.com/pyfar/pyfar/blob/main/pyfar/classes/coordinates.py
theta_phi_r[:, 2] = numpy.sqrt(xyz[:, 0] ** 2 + xyz[:, 1] ** 2 + xyz[:, 2] ** 2) # radius
z_div_r = numpy.where(theta_phi_r[:, 2] != 0, xyz[:, 2] / theta_phi_r[:, 2], 0)
theta_phi_r[:, 1] = numpy.rad2deg(numpy.arccos(z_div_r)) # colatitude
theta_phi_r[:, 0] = numpy.rad2deg(numpy.mod(numpy.arctan2(xyz[:, 1], xyz[:, 0]), 2 * numpy.pi)) # azimuth
# added fixes, including rounding to 2-3 decimal places:
theta_phi_r[:, 2] = numpy.around(theta_phi_r[:, 2], 3) # this one round to 3 decimal places (1mm precision)
theta_phi_r[:, 1] = numpy.around(- (theta_phi_r[:, 1] - 90), 2) # to match what "SOFAplotHRTF.m" does
theta_phi_r[:, 0] = numpy.around(theta_phi_r[:, 0], 3)
# theta_phi_r[:, 0] = numpy.around(numpy.where(theta_phi_r[:, 0] > 180,
# theta_phi_r[:, 0] - 360, theta_phi_r[:, 0]), 2)
return theta_phi_r
def scan_for_errors(project_path): # find any broken data in this project
# NOTE: made only to work with 1 source!!!
num_sources = 1
lowest_corrupt_frq = 0 # init
for source in range(num_sources):
full_p = os.path.join(project_path, 'NumCalc', 'source_' + str(source + 1))
# list all NumCalc logs
all_nc_list = [f for f in os.listdir(full_p) if
f.endswith('.out') and f.startswith('NC') and os.path.isfile(os.path.join(full_p, f))]
for file in all_nc_list:
# read the file
with open(os.path.join(full_p, file)) as ff:
lines = ff.readlines()
ff.close()
frequency_key = ', Frequency = '
non_convergence_key = 'Warning: Maximum number of iterations is reached!'
frequencies = [] # init
non_convergence_freq = [] # init
# find line containing the number of frequency steps
for line_nr in range(len(lines)):
if lines[line_nr].__contains__(frequency_key):
where_is_it = lines[line_nr].find(frequency_key)
# note down what is the frequency of this log and on which line
# (because more than 1 frequency may be in the same worker/log-file)
frequencies.append(int(lines[line_nr].strip()[where_is_it + len(frequency_key):-3]))
elif lines[line_nr].startswith(non_convergence_key):
if frequencies[-1] > 24000: # affects very high frequencies only
print('Warning - Non-Convergence detected at over 24kHz (at ' + str(frequencies[-1]) + 'Hz)')
else:
print('PROBLEM!! - Non-Convergence issue in important range: @ ' + str(frequencies[-1]) + 'Hz')
non_convergence_freq.append(frequencies[-1]) # note down which frequency failed
if len(non_convergence_freq) > 0: # we have non-convergence
if lowest_corrupt_frq == 0:
lowest_corrupt_frq = min(non_convergence_freq)
elif lowest_corrupt_frq > min(non_convergence_freq):
lowest_corrupt_frq = min(non_convergence_freq) # keep only the lowest problematic frequency
return lowest_corrupt_frq
def write_new_numFrequencies(path1, path2, corrupt_frq1, corrupt_frq2):
if min(corrupt_frq1, corrupt_frq2) == 0: # ignore the Zero init value
lowest_corrupt_frq = max(corrupt_frq1, corrupt_frq2)
else:
lowest_corrupt_frq = min(corrupt_frq1, corrupt_frq2)
for pro in [1, 2]: # execute for BOTH PROJECTS
if pro == 1:
project_path = path1
else:
project_path = path2
info_path = os.path.join(project_path, 'Info.txt')
# 1 Back up the Info.txt
info_bckp_path = os.path.join(project_path, 'Info_BACKUP_original.txt')
if not os.path.isfile(info_bckp_path): # backup does not yet exist
shutil.copy2(info_path, info_bckp_path)
# 2 read the original file
with open(info_path) as ff:
lines = ff.readlines()
ff.close()
# 3 Detect important information:
for line_nr in range(len(lines)):
if lines[line_nr].startswith('Minimum evaluated Frequency: '):
min_frequency = float(lines[line_nr].strip()[len('Minimum evaluated Frequency: '):])
elif lines[line_nr].startswith('Highest evaluated Frequency: '):
line_of_max_frequency = line_nr
orig_max_frequency = float(lines[line_nr].strip()[len('Highest evaluated Frequency: '):])
elif lines[line_nr].startswith('Frequency Stepsize: '):
frequency_step = float(lines[line_nr].strip()[len('Frequency Stepsize: '):])
elif lines[line_nr].startswith('Frequency Steps: '):
line_of_frequency_steps_nr = line_nr
orig_frequency_steps_nr = int(lines[line_nr].strip()[len('Frequency Steps: '):])
# 4 MODIFY the number of simulated frequency steps in the Info.txt
new_max_frequency = lowest_corrupt_frq - frequency_step
new_frequency_steps_nr = 1 + round((new_max_frequency - min_frequency) / frequency_step)
lines[line_of_frequency_steps_nr] = ('Frequency Steps: ' + str(int(new_frequency_steps_nr)) + '\n')
lines[line_of_max_frequency] = ('Highest evaluated Frequency: ' + str(new_max_frequency) + '\n')
# 5 add extra modification info for the Info.txt
lines.append('\n\n\n--------------------------------------------------------------------------------\n\n')
lines.append('NOTE: this file was MODIFIED by "join_sofa_files.py" to reduce maximum Frequency\n')
lines.append(' because simulation issue(s) were detected at higher frequencies.\n')
lines.append('\n Initial "Highest evaluated Frequency" was ' + str(orig_max_frequency) + '\n')
# noinspection PyUnboundLocalVariable
lines.append(' Initial "Frequency Steps" was ' + str(orig_frequency_steps_nr) + '\n')
# 6 OVER-WRITE the new Info.txt file
f3 = open(info_path, "w")
for li in lines:
f3.write(li)
f3.close()
def plot_hrir_pair(sofa_data_path, noise_floor=-70):
"""
Plots both left and right ear HRIRs in an easy to diagnose format
# based on tutorial:
https://sofar.readthedocs.io/en/latest/working_with_sofa_files.html#working-with-sofa
Params:
sofa_data_path : full path to the sofa HRIR file to load and plot
noise_floor : were to limit colormap to provide easy to read image
= -50 # noise floor value used in "SOFAplotHRTF.m" from SOFA Matlab API
= -70 # value that visually matches the look of "SOFAplotHRTF.m" from SOFA MatlabAPI
"""
# read the SOFA file using pyfar
data_ir, source_coordinates, receiver_coordinates = pf.io.read_sofa(sofa_data_path)
# find all source positions on or in the vicinity of the HORIZONTAL plane using the find_slice method
_, mask = source_coordinates.find_slice('elevation', unit='deg', value=0, show=False) # do not Show the plot
angles = source_coordinates.get_sph('top_elev', 'deg')[mask, 0]
# noinspection PyTypeChecker
axes_subplot = plt.subplots(2, 1, figsize=(8, 6), sharex=True, label=sofa_data_path)
axes_subplot[0].set_size_inches(8, 10, forward=True)
for Ear in [0, 1]: # Ear 0 == LEFT and 1 == RIGHT ear
# main plot function
ax_h, qm_h, cb_h = pf.plot.time_2d(data_ir[mask, Ear], indices=angles, dB=True, log_prefix=20,
log_reference=1, unit=None, ax=axes_subplot[1][Ear],
cmap=mpl.cm.get_cmap(name='hot_r'),
orientation='horizontal') # YlOrBr hot hot_r RdBu
natural_limits = qm_h.get_clim() # get the full range of data
qm_h.set_clim([natural_limits[1] + noise_floor, natural_limits[1]]) # set limits according to noise floor
ax_h[0].set_ylabel("Azimuth angle in degrees")
# ax_h.set_ylim(0, 3)
if Ear == 0:
ax_h[0].set_title("Left ear HRIR (Horizontal plane)")
elif Ear == 1:
ax_h[0].set_title("Right ear HRIR (Horizontal plane)")
axes_subplot[0].savefig((sofa_data_path[:-5] + '.png'), dpi=100)
# ToDo - accept different input scenarios:
# A (DONE) - no inputs = scan start folder for 2 projects to merge.
# B - 1 input = Use the given path to identify a pair of projects to merge
# where folder-name for LEFT ends on '_L' and RIGHT ends on '_R'.
# C - 2 inputs = Merge the 2 projects that were given as input.
# # Identify starting conditions
print('join_sofa_files started with SavePath = "' + SavePath + '"')
if len(Input1) > 0 and (os.path.basename(Input1)[-2:] == "_L" or os.path.basename(Input1)[-2:] == "_R"):
ProjectName_withL = True # we can search for a mathing pairs of "_L" and "_R" projects
else:
ProjectName_withL = False # unless "Input2" is specified, search in this folder for 2 merge-able projects
if len(Input1) == 0 or (not ProjectName_withL and len(Input2) == 0):
# mode A = scan folder for 2 projects to merge.
counter = 0 # init
if len(Input1) == 0:
Search_Folder = SavePath # no inputs - search folder where the script is (or any custom SavePath)
else:
Search_Folder = Input1 # search in the first Input1 path
Path1 = '' # init
Path2 = '' # init
for subdir in os.listdir(Search_Folder):
# we only search folders that do not start with "." (extra scripts or files are allowed)
if os.path.isdir(os.path.join(Search_Folder, subdir)) and not subdir.startswith('.'):
counter += 1
if counter == 1:
Path1 = os.path.join(Search_Folder, subdir) # project 1
elif counter == 2:
Path2 = os.path.join(Search_Folder, subdir) # project 2
else:
raise Exception('more than 2 folders to merge in the Search_Folder = "' + Search_Folder + '"')
if len(Path2) == 0:
raise Exception('less than 2 folders to merge in the Search_Folder = "' + Search_Folder + '"')
elif len(Input2) == 0: # mode B - 1 input = Use the given path to identify a pair of projects to merge
# where folder-name for LEFT ends on '_L' and RIGHT ends on '_R'.
raise Exception("not implemented")
else: # mode C - 2 inputs = Merge the 2 projects that were given as input.
Path1 = Input1
Path2 = Input2
# DONE - autodetect details of the input scenario:
# a - folder with Mesh2HRTF project
# b - folder with .sofa files
# c - complete file path to "HRTF_<something>.sofa" of "HRIR_<something>.sofa"
if os.path.isdir(Path1):
# check if on Path1 is a Mesh2HRTF project folder that is not yet post-processed:
Lowest_corrupt_frq_1 = 0 # init
Lowest_corrupt_frq_2 = 0 # init
if not os.path.isdir(os.path.join(Path1, 'Output2HRTF')) and os.path.isfile(os.path.join(Path1, 'Output2HRTF.py')):
# Mode-A0 detected:
print('\nChecking for errors in "' + Path1 + '"...\n')
Lowest_corrupt_frq_1 = scan_for_errors(Path1) # find any broken data on Path1
print('\nChecking for errors in "' + Path2 + '"...\n')
Lowest_corrupt_frq_2 = scan_for_errors(Path2) # find any broken data on Path1
if Lowest_corrupt_frq_1 > 1 or Lowest_corrupt_frq_2 > 1: # there was coppupt data in the simulation
# exclude all frequencies that are corrupt by specifying the last Corrupt one.
write_new_numFrequencies(Path1, Path2, Lowest_corrupt_frq_1, Lowest_corrupt_frq_2)
print('\nRunning Mesh2HRTF post-processing in "' + Path1 + '"...\n')
os.chdir(Path1) # necessary for the 'Output2HRTF.py'
subprocess.call("Output2HRTF.py", shell=True) # run the normal Mesh2HRTF post-processing
# noinspection PyUnboundLocalVariable
os.chdir(Search_Folder) # change back to Search_Folder (just in case)
# check if on Path2 is a Mesh2HRTF project folder that is not yet post-processed:
if not os.path.isdir(os.path.join(Path2, 'Output2HRTF')) and os.path.isfile(os.path.join(Path2, 'Output2HRTF.py')):
print('\nRunning Mesh2HRTF post-processing in "' + Path2 + '"...\n')
os.chdir(Path2) # necessary for the 'Output2HRTF.py'
subprocess.call("Output2HRTF.py", shell=True) # run the normal Mesh2HRTF post-processing
os.chdir(Search_Folder) # change back to Search_Folder (just in case)
if os.path.isdir(os.path.join(Path1, 'Output2HRTF')): # a - folder with Mesh2HRTF project
SofaFLDR_1 = os.path.join(Path1, 'Output2HRTF')
SofaFLDR_2 = os.path.join(Path2, 'Output2HRTF')
else: # b - folder with .sofa files
SofaFLDR_1 = Path1
SofaFLDR_2 = Path2
AllSofaFiles1 = os.listdir(SofaFLDR_1) # searching for '*.sofa' files
AllSofaFiles2 = os.listdir(SofaFLDR_2)
FolderName = os.path.basename(Path1)
else: # os.path.isfile(Path1): # c - complete file path to "HRTF_<something>.sofa" of "HRIR_<something>.sofa"
AllSofaFiles1 = [os.path.basename(Path1)]
AllSofaFiles2 = [os.path.basename(Path2)]
SofaFLDR_1 = os.path.dirname(Path1)
SofaFLDR_2 = os.path.dirname(Path2)
FolderName = os.path.basename(SofaFLDR_1)
if FolderName.endswith('_L') or FolderName.endswith('_R'):
FolderName = FolderName[:-2] # drop the "_L" or "_R" from the folder name (tidy up)
BasePath = os.path.join(SavePath, 'SOFA_' + FolderName + '_merged')
# autodetect which is Left and which is Right based on source coord.
HRTFsofa_read_1 = sf.read_sofa(os.path.join(SofaFLDR_1, AllSofaFiles1[0]))
if HRTFsofa_read_1.ReceiverPosition.shape[0] > 1: # more than one receiver in this SOFA file
raise Exception('The SOFA file is not for a single ear! File in: "' + SofaFLDR_1 + '"')
if HRTFsofa_read_1.ReceiverPosition[0, 1] > 0: # sofa file is for Left ear (positive Y coordinate)
AllSofaFiles_L = AllSofaFiles1
AllSofaFiles_R = AllSofaFiles2
SofaFLDR_L = SofaFLDR_1
SofaFLDR_R = SofaFLDR_2
else: # flip inputs because the 2nd folder contains Left ear
AllSofaFiles_R = AllSofaFiles1
AllSofaFiles_L = AllSofaFiles2
SofaFLDR_R = SofaFLDR_1
SofaFLDR_L = SofaFLDR_2
# Save complex pressure as SOFA file
print('\n Merging HRTF SOFA files from: "' + os.path.basename(Path1) + '" and "' +
os.path.basename(Path2) + '"')
del SofaFLDR_2, SofaFLDR_1, AllSofaFiles2, AllSofaFiles1, Path1, Path2
for ii in range(2): # range(len(AllSofaFiles_L))
# # Load the data
sofa_read_L = sf.read_sofa(os.path.join(SofaFLDR_L, AllSofaFiles_L[ii]))
sofa_read_R = sf.read_sofa(os.path.join(SofaFLDR_R, AllSofaFiles_R[ii]))
# detect data type:
if sofa_read_L.GLOBAL_DataType == 'FIR':
SOFA_Type = 'HRIR'
# sanity checks that the files can be merged
if sofa_read_L.Data_SamplingRate != sofa_read_R.Data_SamplingRate:
raise Exception("The SamplingRates of the two input SOFA files do not match!")
elif sofa_read_L.Data_IR.shape[1] != 1 or sofa_read_R.Data_IR.shape[1] != 1:
raise Exception("The input SOFA files contain more than one ear - impossible to merge")
elif sofa_read_L.Data_IR.shape[2] != sofa_read_R.Data_IR.shape[2]:
raise Exception("The number of steps in the two input SOFA files do not match!")
elif sofa_read_L.ReceiverPosition_Units != sofa_read_R.ReceiverPosition_Units:
raise Exception("The ReceiverPosition_Units in the two input SOFA files do not match!")
else:
SOFA_Type = 'HRTF'
# sanity checks that the files can be merged
if sofa_read_L.N[0] != sofa_read_R.N[0] or sofa_read_L.N[-1] != sofa_read_R.N[-1]:
raise Exception("The sampling Frequencies of the two input SOFA files do not match!")
if ii == 0 and not os.path.isdir(BasePath):
os.mkdir(BasePath) # create new output folder
# noinspection PyTypeChecker
Full_Ref_Path = merge_write_sofa(sofa_read_L, sofa_read_R, BasePath, AllSofaFiles_L[ii], SOFA_Type)
if SOFA_Type == 'HRIR': # note down what should be plotted in the end
SOFA_to_plot = Full_Ref_Path
# make a simple copy of the Info.txt
if os.path.isfile(os.path.join(os.path.dirname(SofaFLDR_L), 'Info.txt')):
shutil.copy2(os.path.join(os.path.dirname(SofaFLDR_L), 'Info.txt'), os.path.join(BasePath, 'Info.txt'))
print('\n Done: Merged SOFA files are saved in:\n "' + BasePath + '"\n')
if Diagnostic_plot:
# noinspection PyBroadException
try:
#############################################################################
# plotting part: it is OK if plotting crashes, because the merging is already complete.
# these imports are here so that key functionality works even without "pyfar" and "matplotlib"
import pyfar as pf #
https://pypi.org/project/sofar/ pip install pyfar
import matplotlib as mpl
import matplotlib.pyplot as plt
# noinspection PyUnboundLocalVariable
plot_hrir_pair(SOFA_to_plot, -50)
#############################################################################
except Exception:
print('\n WARNING: plotting command crashed. Likely reasons are:')
print(' 1 - missing "pyfar" or missing "matplotlib"')
print(' (to install "pyfar" in command line type "pip install pyfar" ')
print(' 2 - obsolete python packages (especially numpy) - try running: "pip list --outdated"')
print(' (to update a package use for example: "pip install --upgrade numpy"')
print(' 3 - actual bug in this script (debug, or open an issue if you need this plot)\n')
input('Hit Enter to close - Plotting is complete
')
else:
input('Hit Enter to exit
')
#
# made for the Mesh2HRTF codebase
#
# mesh2hrtf.sourceforge.net
#
# Mesh2HRTF is licensed under the GNU Lesser General Public License as
# published by the Free Software Foundation, either version 3 of the License,
# or (at your option) any later version. Mesh2HRTF is distributed in the hope
# that it will be useful, but WITHOUT ANY WARRANTY; without even the implied
# warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# Lesser General Public License for more details. You should have received a
# copy of the GNU LesserGeneral Public License along with Mesh2HRTF. If not,
# see <
http://www.gnu.org/licenses/lgpl.html>.
#
# If you use Mesh2HRTF:
# - Provide credits:
# "Mesh2HRTF, H. Ziegelwanger, ARI, OEAW (mesh2hrtf.sourceforge.net)"
# - In your publication, cite both articles:
# [1] Ziegelwanger, H., Kreuzer, W., and Majdak, P. (2015). "Mesh2HRTF:
# Open-source software package for the numerical calculation of head-
# related transfer functions," in Proceedings of the 22nd ICSV,
# Florence, IT.
# [2] Ziegelwanger, H., Majdak, P., and Kreuzer, W. (2015). "Numerical
# calculation of listener-specific head-related transfer functions and
# sound localization: Microphone model and mesh discretization," The
# Journal of the Acoustical Society of America, 138, 208-222.
). While Blender is the Frankenstein monster itself, where motivated people kept adding tools over tools over tools just because they could. Decades of that did not result in something very intuitive IMO.
