diff --git a/VERSION b/VERSION index 524cb552..26aaba0e 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -1.1.1 +1.2.0 diff --git a/mesh2hrtf/Mesh2Input/mesh2input.py b/mesh2hrtf/Mesh2Input/mesh2input.py index feebd158..9dc105a6 100644 --- a/mesh2hrtf/Mesh2Input/mesh2input.py +++ b/mesh2hrtf/Mesh2Input/mesh2input.py @@ -15,7 +15,7 @@ bl_info = { "name": "Mesh2HRTF export add-on", "author": "The Mesh2HRTF developers", - "version": (1, 1, 1), + "version": (1, 2, 0), "blender": (2, 80, 0), "location": "File > Export", "description": "Export Blender scene as Mesh2HRTF project", @@ -169,6 +169,22 @@ class ExportMesh2HRTF(bpy.types.Operator, ExportHelper): min=0, max=24000, ) + bins_per_oct_decim: IntProperty( + name="Bins/octave", + description=("Value for distributing the frequencies in log scale," + "in bins/octave."), + default=18, + min=0, + max=36, + ) + num_octaves_decim: IntProperty( + name="Num. of octaves", + description=("Number of higher octaves in which the frequencies will" + "be distributed in log scale."), + default=2, + min=0, + max=6, + ) @classmethod def poll(cls, context): @@ -231,6 +247,10 @@ def draw(self, context): row.prop(self, "frequencyVectorType") row = layout.row() row.prop(self, "frequencyVectorValue") + row = layout.row() + row.prop(self, "bins_per_oct_decim") + row = layout.row() + row.prop(self, "num_octaves_decim") def save(operator, context, @@ -240,6 +260,8 @@ def save(operator, maxFrequency=20000, frequencyVectorType='Step size', frequencyVectorValue=100, + bins_per_oct_decim=18, + num_octaves_decim=2, pictures=True, evaluationGrids='ARI', materialSearchPaths='None', @@ -442,8 +464,9 @@ def save(operator, numFrequencySteps = int(frequencyVectorValue) frequencies, frequencyStepSize, numFrequencySteps = \ - _distribute_frequencies(minFrequency, maxFrequency, - frequencyStepSize, numFrequencySteps) + _distribute_frequencies_octaves(minFrequency, maxFrequency, + frequencyStepSize, numFrequencySteps, + bins_per_oct_decim, num_octaves_decim) # Write parameters.json (feedback for user, not used by NumCalc) -------------- _write_parameters_json( @@ -453,7 +476,8 @@ def save(operator, reference, computeHRIRs, sourceType, numSources, sourceXPosition, sourceYPosition, sourceZPosition, - frequencies, frequencyStepSize, numFrequencySteps) + frequencies, frequencyStepSize, numFrequencySteps, + bins_per_oct_decim, num_octaves_decim) # Render pictures of the model ------------------------------------------------ if pictures: @@ -870,7 +894,8 @@ def _write_parameters_json( reference, computeHRIRs, sourceType, numSources, sourceXPosition, sourceYPosition, sourceZPosition, - frequencies, frequencyStepSize, numFrequencySteps): + frequencies, frequencyStepSize, numFrequencySteps, + bins_per_oct_decim, num_octaves_decim): # calculate missing parameters if "ear" in sourceType: @@ -921,6 +946,8 @@ def _write_parameters_json( "frequencyStepSize": frequencyStepSize, "minFrequency": frequencies[0], "maxFrequency": frequencies[-1], + "bins_per_oct_decim": bins_per_oct_decim, + "num_octaves_decim": num_octaves_decim, "frequencies": frequencies } @@ -1031,8 +1058,9 @@ def _calculateReceiverProperties(obj, obj_data, unitFactor): return earCenter, earArea -def _distribute_frequencies(minFrequency, maxFrequency, - frequencyStepSize, numFrequencySteps): +def _distribute_frequencies_octaves(minFrequency, maxFrequency, + frequencyStepSize, numFrequencySteps, + bins_per_oct_decim, num_octaves_decim): """Calculate list of frequencies. Returns @@ -1041,11 +1069,11 @@ def _distribute_frequencies(minFrequency, maxFrequency, frequencies: list list that holds the frequencies in Hz that are calculated. frequencyStepSize: float - Step size between successive frequencies (written to parameters.json). - This is returned because it might be zero during the function call. + Step size between successive frequencies (written to Info.txt). This is + returned because it might be zero during the function call. numFrequencySteps: int - Number of frequencies to be simulated (written to parameters.json). - This is returned because it might be zero during the function call. + Number of frequencies to be simulated (written to Info.txt). This is + returned because it might be zero during the function call. """ @@ -1093,9 +1121,28 @@ def _distribute_frequencies(minFrequency, maxFrequency, else: frequencyStepSize = frequencies[1] - frequencies[0] + if bins_per_oct_decim > 0 and num_octaves_decim > 0: + + freqs_log = frequencies[-1]*np.power(2,-(1/bins_per_oct_decim)*(np.linspace( + bins_per_oct_decim*num_octaves_decim, 0, + int(1 + np.round(bins_per_oct_decim*num_octaves_decim))))) + + freqs_log = freqs_log[np.append(np.diff(freqs_log), frequencyStepSize + 1) > frequencyStepSize] + + frequencies = np.array(frequencies) + + freqs_lin = frequencies[frequencies < freqs_log[0]] + + frequencies = np.append(freqs_lin, freqs_log) + + frequencies = frequencies.tolist() + + numFrequencySteps = len(frequencies) + return frequencies, frequencyStepSize, numFrequencySteps + def _render_pictures(filepath1, unitFactor): """Render pictures of the 3D mesh and save to export folder.""" diff --git a/mesh2hrtf/Output2HRTF/output2hrtf.py b/mesh2hrtf/Output2HRTF/output2hrtf.py index 711e96ae..0b6fbb14 100644 --- a/mesh2hrtf/Output2HRTF/output2hrtf.py +++ b/mesh2hrtf/Output2HRTF/output2hrtf.py @@ -107,6 +107,25 @@ def output2hrtf(folder=None): sofa, params["sourceType"], params["sourceArea"], params["speedOfSound"], params["densityOfMedium"], mode="min") + # Resampling spectrum to regular sampling, if necessary + if num_octaves_decim > 0: + pressure = sofa.Data_Real + 1j*sofa.Data_Imag + frequencies = sofa.N + sampling_rate = round(2*sofa.N[-1]) + + num_octaves_decim = params["num_octaves_decim"] + + pressure, frequencies, inds_lin_freq = _convert_to_regular_samp_freq( + pressure, frequencies) + + inds_correct_phase = inds_lin_freq + + pressure = _correct_phase(pressure, inds_correct_phase) + + sofa.N = frequencies.flatten() + sofa.Data_Real = np.real(pressure) + sofa.Data_Imag = np.imag(pressure) + # write HRTF data to SOFA file sf.write_sofa(os.path.join( folder, 'Output2HRTF', f'HRTF_{grid}.sofa'), sofa) @@ -361,3 +380,73 @@ def _output_to_hrtf_load(foldername, filename, numFrequencies): data[ii, :] = tmpData if tmpData else np.nan return data, np.arange(start_index, numDatalines + start_index) + + + +def _convert_to_regular_samp_freq(pressure, frequencies): + """ + Converts irregularly sampled frequency axis into regular intervals via + interpolation. + """ + print('\n \n Irregular frequency sampling identified. \n \n' + 'Converting to regular frequency sampling... \n \n') + + freqs_decim = np.array(frequencies).flatten() + minFrequency = freqs_decim[0] + frequencyStepSize = freqs_decim[1] - freqs_decim[0] + frequencySteps = freqs_decim[-1]/frequencyStepSize + + # get all frequencies to be calculated + frequencies = np.array([ff*frequencyStepSize+minFrequency + for ff in range(int(frequencySteps))]) + pressure_decim = np.squeeze(pressure) + + # Interpolate spectra back to regular sampling in frequency + pressure_decim_mag = np.array([np.interp(frequencies, freqs_decim, np.abs( + pressure_decim)[i]) for i in range(pressure_decim.shape[0])]) + + pressure_decim_phase = np.array([np.interp(frequencies, freqs_decim, + np.angle(pressure_decim)[i]) + for i in range(pressure_decim.shape[0])]) + + # # Zeroing phase for high frequencies in the interpolated region + inds_lin_freq = ((np.diff(np.diff(np.concatenate((freqs_decim, + np.array([freqs_decim[-1]*2, freqs_decim[-1]*4])))))) == 0) + + pressure_decim_phase[:, np.sum(np.array(inds_lin_freq)) : ] = 0 + + pressure = np.zeros([pressure.shape[0], pressure.shape[1], + pressure_decim_mag.shape[1]], dtype = "complex_") + + pressure[:,0,:] = np.multiply(pressure_decim_mag, np.exp(1j*pressure_decim_phase)) + + inds_lin_freq = (frequencies < frequencies[np.sum(np.array(inds_lin_freq))]) + + return pressure, frequencies, inds_lin_freq + + +def _correct_phase(pressure, inds_correct_phase): + """ + Corrects the phase for high frequencies that were simulated using irregular + frequency sampling using the overall group delay of the linearly sampled region. + """ + new_phase_all = np.zeros(shape=pressure.shape) + + for Direction_ind in range(pressure.shape[0]): + for ear_ind in range(pressure.shape[1]): + + phase = np.unwrap(np.angle(pressure[Direction_ind, ear_ind, :])) + + phase_diff = np.diff(phase[inds_correct_phase]) + phase_diff_mean = np.average(phase_diff) + + phase[~inds_correct_phase] = np.linspace( + start = (phase[inds_correct_phase])[-1] + phase_diff_mean, + stop = (phase[inds_correct_phase])[-1] + phase_diff_mean*np.sum(~inds_correct_phase), + num = np.sum(~inds_correct_phase)) + + new_phase_all[Direction_ind, ear_ind, :] = phase + + pressure = np.multiply(np.abs(pressure), np.exp(1j*new_phase_all)) + + return pressure \ No newline at end of file