diff --git a/lingvodoc/schema/gql_cognate.py b/lingvodoc/schema/gql_cognate.py index 0680ac16..2b32b412 100644 --- a/lingvodoc/schema/gql_cognate.py +++ b/lingvodoc/schema/gql_cognate.py @@ -1964,7 +1964,7 @@ def f(): for tier_result in tier_result_list: - for (interval_str, interval_r_length, j_local, p_mean, i_list, f_list, + for (interval_str, interval_r_length, j_local, s_local, p_mean, i_list, f_list, sign_longest, sign_highest, source_index) in tier_result.interval_data_list: interval_str_list = interval_str.split() diff --git a/lingvodoc/schema/query.py b/lingvodoc/schema/query.py index f5eaddb3..6c5a419d 100644 --- a/lingvodoc/schema/query.py +++ b/lingvodoc/schema/query.py @@ -5181,7 +5181,7 @@ def perform_phonological_statistical_distance( # ...for all intervals. - for (interval_str, interval_r_length, j_local, p_mean, i_list, f_list, + for (interval_str, interval_r_length, j_local, s_local, p_mean, i_list, f_list, sign_longest, sign_highest, source_index) in tier_result.interval_data_list: formant_list.append(tuple(map(float, f_list[:2]))) diff --git a/lingvodoc/views/v2/jitter.py b/lingvodoc/views/v2/jitter.py index 012af57e..768f5087 100644 --- a/lingvodoc/views/v2/jitter.py +++ b/lingvodoc/views/v2/jitter.py @@ -75,11 +75,15 @@ def unidirectional_autowindow(pulse, tmin, tmax): tmax = pulse['xmax'] return tmin, tmax +def get_window_points(pulse, tmin, tmax): + left, right = bisect.bisect(pulse['t'], tmin), bisect.bisect(pulse['t'], tmax) + return left, right - bool(right) # decrease right if not zero + def get_mean_period(pulse, tmin, tmax): tmin, tmax = unidirectional_autowindow(pulse, tmin, tmax) - left, right = bisect.bisect(pulse['t'], tmin), bisect.bisect(pulse['t'], tmax) - first, last = left, right - bool(right) # decrease right if it is not zero + first, last = get_window_points(pulse, tmin, tmax) + numberOfPeriods = 0 dsum = 0.0 for ipoint in range(first, last): @@ -91,8 +95,7 @@ def get_mean_period(pulse, tmin, tmax): def get_jitter_local(pulse, tmin, tmax): tmin, tmax = unidirectional_autowindow(pulse, tmin, tmax) - left, right = bisect.bisect(pulse['t'], tmin), bisect.bisect(pulse['t'], tmax) - first, last = left, right - bool(right) # decrease right if it is not zero + first, last = get_window_points(pulse, tmin, tmax) numberOfPeriods = max(0, last - first) if numberOfPeriods < 2: return None diff --git a/lingvodoc/views/v2/phonology.py b/lingvodoc/views/v2/phonology.py index c0a04a81..f9626ac6 100644 --- a/lingvodoc/views/v2/phonology.py +++ b/lingvodoc/views/v2/phonology.py @@ -97,6 +97,7 @@ from lingvodoc.utils import sanitize_worksheet_name from lingvodoc.views.v2.utils import anonymous_userid, as_storage_file, message, storage_file, unimplemented from lingvodoc.views.v2.jitter import pitch_to_point, get_jitter_local +from lingvodoc.views.v2.shimmer import get_shimmer_local from pdb import set_trace as A @@ -2179,6 +2180,7 @@ def __init__(self, max_length_str, max_length_r_length, max_length_jt_local, + max_length_sh_local, max_length_p_mean, max_length_i_list, max_length_f_list, @@ -2186,6 +2188,7 @@ def __init__(self, max_intensity_str, max_intensity_r_length, max_intensity_jt_local, + max_intensity_sh_local, max_intensity_p_mean, max_intensity_i_list, max_intensity_f_list, @@ -2202,6 +2205,7 @@ def __init__(self, self.max_length_str = max_length_str self.max_length_r_length = max_length_r_length self.max_length_jt_local = max_length_jt_local + self.max_length_sh_local = max_length_sh_local self.max_length_p_mean = max_length_p_mean self.max_length_i_list = max_length_i_list self.max_length_f_list = max_length_f_list @@ -2210,6 +2214,7 @@ def __init__(self, self.max_intensity_str = max_intensity_str self.max_intensity_r_length = max_intensity_r_length self.max_intensity_jt_local = max_intensity_jt_local + self.max_intensity_sh_local = max_intensity_sh_local self.max_intensity_p_mean = max_intensity_p_mean self.max_intensity_i_list = max_intensity_i_list self.max_intensity_f_list = max_intensity_f_list @@ -2231,11 +2236,12 @@ def format(self): '{0:.2f}%'.format(r_length * 100), is_max_length, is_max_intensity, source_index], j_local, + s_local, p_mean, i_list, f_list) - for interval_str, r_length, j_local, p_mean, i_list, f_list, is_max_length, is_max_intensity, source_index in + for interval_str, r_length, j_local, s_local, p_mean, i_list, f_list, is_max_length, is_max_intensity, source_index in self.interval_data_list] return pprint.pformat( @@ -2374,6 +2380,8 @@ def process_sound(tier_data_list, sound, vowel_selection = None): max_length_f_list = [] max_intensity_jt_local = 0 max_length_jt_local = 0 + max_intensity_sh_local = 0 + max_length_sh_local = 0 interval_data_list = [] if vowel_selection is None or vowel_selection == True: @@ -2397,6 +2405,12 @@ def process_sound(tier_data_list, sound, vowel_selection = None): max_length_jt_local = ( get_jitter_local(pulse, *max_length_interval[:2])) + max_intensity_sh_local = ( + get_shimmer_local(pulse, sound_dict, *max_intensity_interval[:2])) + + max_length_sh_local = ( + get_shimmer_local(pulse, sound_dict, *max_length_interval[:2])) + if vowel_selection is None or vowel_selection == False: print('Calculating lists for all intervals...') @@ -2411,9 +2425,13 @@ def process_sound(tier_data_list, sound, vowel_selection = None): jitter_list = [ get_jitter_local(pulse, begin_sec, end_sec) for begin_sec, end_sec, text in interval_list] - log.debug(f"jitter_list: {jitter_list}") + shimmer_list = [ + get_shimmer_local(pulse, sound_dict, begin_sec, end_sec) + for begin_sec, end_sec, text in interval_list] + log.debug(f"shimmer_list: {shimmer_list}") + # Preparing data of all other intervals. str_list = [ @@ -2438,6 +2456,7 @@ def process_sound(tier_data_list, sound, vowel_selection = None): (interval_str, (end - begin) / mean_interval_length, f'{j_local:.3f}', + f'{s_local:.3f}', f'{p_mean:.3f}', [f'{i_min:.3f}', f'{i_max:.3f}', f'{i_max - i_min:.3f}'], list(map('{0:.3f}'.format, f_list)), @@ -2449,6 +2468,7 @@ def process_sound(tier_data_list, sound, vowel_selection = None): index, (interval_str, j_local, + s_local, p_mean, (_, i_min, i_max), f_list, @@ -2459,6 +2479,7 @@ def process_sound(tier_data_list, sound, vowel_selection = None): zip( str_list, jitter_list, + shimmer_list, pitch_means, intensity_list, formant_list, @@ -2493,6 +2514,7 @@ def process_sound(tier_data_list, sound, vowel_selection = None): max_length_str, max_length / mean_interval_length, f'{max_length_jt_local:.3f}', + f'{max_length_sh_local:.3f}', f'{xl_p_mean:.3f}', list(map('{0:.3f}'.format, max_length_i_list)), list(map('{0:.3f}'.format, max_length_f_list)), @@ -2500,6 +2522,7 @@ def process_sound(tier_data_list, sound, vowel_selection = None): max_intensity_str, (max_intensity_interval[1] - max_intensity_interval[0]) / mean_interval_length, f'{max_intensity_jt_local:.3f}', + f'{max_intensity_sh_local:.3f}', f'{xi_p_mean:.3f}', list(map('{0:.3f}'.format, max_intensity_i_list)), list(map('{0:.3f}'.format, max_intensity_f_list)), @@ -3185,6 +3208,7 @@ def compile_workbook( 'Longest (seconds) interval', 'Relative length', 'Jitter local', + 'Shimmer local', 'Pitch mean (Hz)', 'Intensity minimum (dB)', 'Intensity maximum (dB)', 'Intensity range (dB)', 'F1 mean (Hz)', 'F2 mean (Hz)', 'F3 mean (Hz)', @@ -3193,6 +3217,7 @@ def compile_workbook( 'Highest intensity (dB) interval', 'Relative length', 'Jitter local', + 'Shimmer local', 'Pitch mean (Hz)', 'Intensity minimum (dB)', 'Intensity maximum (dB)', 'Intensity range (dB)', 'F1 mean (Hz)', 'F2 mean (Hz)', 'F3 mean (Hz)', @@ -3207,6 +3232,7 @@ def compile_workbook( 'Interval', 'Relative length', 'Jitter local', + 'Shimmer local', 'Pitch mean (Hz)', 'Intensity minimum (dB)', 'Intensity maximum (dB)', 'Intensity range (dB)', 'F1 mean (Hz)', 'F2 mean (Hz)', 'F3 mean (Hz)', @@ -3231,25 +3257,25 @@ def compile_workbook( if args.vowel_selection: worksheet_results.set_column(0, 2, 20) - worksheet_results.set_column(3, 4, 8, format_percent) - worksheet_results.set_column(5, 8, 8) - worksheet_results.set_column(9, 11, 10) - worksheet_results.set_column(12, 12, 4) - worksheet_results.set_column(13, 13, 20) - worksheet_results.set_column(14, 14, 8, format_percent) - worksheet_results.set_column(15, 18, 8) - worksheet_results.set_column(19, 21, 10) - worksheet_results.set_column(22, 23, 4) - worksheet_results.set_column(24, 24, 8, format_percent) + worksheet_results.set_column(3, 5, 8, format_percent) + worksheet_results.set_column(6, 9, 8) + worksheet_results.set_column(10, 12, 10) + worksheet_results.set_column(13, 13, 4) + worksheet_results.set_column(14, 14, 20) + worksheet_results.set_column(15, 17, 8, format_percent) + worksheet_results.set_column(18, 21, 8) + worksheet_results.set_column(22, 24, 10) + worksheet_results.set_column(25, 26, 4) + worksheet_results.set_column(27, 27, 8, format_percent) else: worksheet_results.set_column(0, 2, 20) - worksheet_results.set_column(3, 4, 8, format_percent) - worksheet_results.set_column(5, 8, 8) - worksheet_results.set_column(9, 11, 10) - worksheet_results.set_column(12, 14, 4) - worksheet_results.set_column(15, 15, 8, format_percent) + worksheet_results.set_column(3, 5, 8, format_percent) + worksheet_results.set_column(6, 9, 8) + worksheet_results.set_column(10, 12, 10) + worksheet_results.set_column(13, 15, 4) + worksheet_results.set_column(16, 16, 8, format_percent) worksheet_dict[group] = ( @@ -3378,6 +3404,7 @@ def next_text(): round(tier_result.max_length_r_length, 4)] + [ float(tier_result.max_length_jt_local)] + + [ float(tier_result.max_length_sh_local)] + [ float(tier_result.max_length_p_mean) ] + i_list_a + f_list_a + @@ -3387,6 +3414,7 @@ def next_text(): round(tier_result.max_intensity_r_length, 4)] + [ float(tier_result.max_intensity_jt_local)] + + [ float(tier_result.max_intensity_sh_local)] + [ float(tier_result.max_intensity_p_mean) ] + i_list_b + f_list_b + @@ -3437,7 +3465,7 @@ def next_text(): else: for index, (interval_str, interval_r_length, - j_local, p_mean, i_list, f_list, sign_longest, sign_highest, source_index) in ( + j_local, s_local, p_mean, i_list, f_list, sign_longest, sign_highest, source_index) in ( enumerate(tier_result.interval_data_list)): @@ -3460,6 +3488,7 @@ def next_text(): round(interval_r_length, 4)] + [float(j_local)] + + [float(s_local)] + [float(p_mean)] + i_list + f_list + diff --git a/lingvodoc/views/v2/shimmer.py b/lingvodoc/views/v2/shimmer.py new file mode 100644 index 00000000..2bd9790b --- /dev/null +++ b/lingvodoc/views/v2/shimmer.py @@ -0,0 +1,114 @@ +import numpy as np +from lingvodoc.views.v2.jitter import ( + unidirectional_autowindow, + get_window_points, + pmin, pmax, + maximumPeriodFactor, +) +from math import pi +from pdb import set_trace as A + +maximumAmplitudeFactor = 1.6 + +class AmplitudeTier: + def __init__(self, tmin, tmax): + self.xmin = tmin + self.xmax = tmax + self.points = [] + + def addPoint(self, t, value): + self.points.append({ + 'number': t, + 'value': value + }) + + +def get_hann_windowed_rms(sound, tmid, widthLeft, widthRight): + + if (edges := get_window_samples(sound, tmid - widthLeft, tmid + widthRight)) is None: + return None + + imin, imax = edges + + sumOfSquares = 0.0 + windowSumOfSquares = 0.0 + for i in range(imin, imax + 1): + t = sound['x1'] + sound['dx'] * i + width = widthLeft if t < tmid else widthRight + windowPhase = (t - tmid) / width # in [-1 .. 1] + window = 0.5 + 0.5 * np.cos(pi * windowPhase) # Hann + if sound['ny'] == 1: + windowedValue = sound['z'][0][i] * window + else: + windowedValue = 0.5 * (sound['z'][0][i] + sound['z'][1][i]) * window + sumOfSquares += windowedValue ** 2 + windowSumOfSquares += window ** 2 + + return np.sqrt(sumOfSquares / windowSumOfSquares) + + +def get_window_samples(sound, tmin, tmax): + imin = np.ceil((tmin - sound['x1']) / sound['dx']) + imax = np.floor((tmax - sound['x1']) / sound['dx']) + imin = int(max(0.0, imin)) + imax = int(min(sound['nx'], imax)) + if imax - imin < 2: + return None + return imin, imax + + +def point_to_amplitude_period(pulse, sound, tmin, tmax): + try: + tmin, tmax = unidirectional_autowindow(pulse, tmin, tmax) + first, last = get_window_points(pulse, tmin, tmax) + if last - first < 2: + raise ValueError(f"Too few pulses between {tmin} and {tmax} seconds.") + + amplitude = AmplitudeTier(tmin, tmax) + for i in range(first + 1, last): + p1 = pulse['t'][i] - pulse['t'][i - 1] + p2 = pulse['t'][i + 1] - pulse['t'][i] + intervalFactor = p1 / p2 if p1 > p2 else p2 / p1 + if pmin == pmax or (pmin <= p1 <= pmax and pmin <= p2 <= pmax and intervalFactor <= maximumPeriodFactor): + peak = get_hann_windowed_rms(sound, pulse['t'][i], 0.2 * p1, 0.2 * p2) + if peak is not None and peak > 0.0: + amplitude.addPoint(pulse['t'][i], peak) + + return amplitude + except Exception as e: + raise ValueError(f"{pulse} & {sound}: not converted to AmplitudeTier.") from e + + +def get_shimmer_local(pulse, sound, tmin, tmax): + numberOfPeaks = 0 + numerator = 0.0 + denominator = 0.0 + try: + tmin, tmax = unidirectional_autowindow(pulse, tmin, tmax) + peaks = point_to_amplitude_period(pulse, sound, tmin, tmax) + points = peaks.points + for i in range(1, len(points)): + p = points[i]['number'] - points[i-1]['number'] + if pmin == pmax or (p >= pmin and p <= pmax): + a1 = points[i-1]['value'] + a2 = points[i]['value'] + amplitudeFactor = a1 / a2 if a1 > a2 else a2 / a1 + if amplitudeFactor <= maximumAmplitudeFactor: + numerator += abs(a1 - a2) + numberOfPeaks += 1 + if numberOfPeaks < 1: + return None + numerator /= numberOfPeaks + numberOfPeaks = 0 + for point in points[:-1]: # why -1? + denominator += point['value'] + numberOfPeaks += 1 + denominator /= numberOfPeaks + if denominator == 0.0: + return None + return numerator / denominator + except Exception as e: + if "Too few pulses between" in str(e): + return None + else: + raise Exception(f"{pulse} & {sound}: shimmer (local) not computed.")