• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
1# Lint as: python2, python3
2# Copyright 2016 The Chromium OS Authors. All rights reserved.
3# Use of this source code is governed by a BSD-style license that can be
4# found in the LICENSE file.
5
6"""This module provides utilities to detect some artifacts and measure the
7    quality of audio."""
8
9from __future__ import absolute_import
10from __future__ import division
11from __future__ import print_function
12import logging
13import math
14import numpy
15from six.moves import range
16
17# Normal autotest environment.
18try:
19    import common
20    from autotest_lib.client.cros.audio import audio_analysis
21# Standalone execution without autotest environment.
22except ImportError:
23    import audio_analysis
24
25
26# The input signal should be one sine wave with fixed frequency which
27# can have silence before and/or after sine wave.
28# For example:
29#   silence      sine wave      silence
30#  -----------|VVVVVVVVVVVVV|-----------
31#     (a)           (b)           (c)
32# This module detects these artifacts:
33#   1. Detect noise in (a) and (c).
34#   2. Detect delay in (b).
35#   3. Detect burst in (b).
36# Assume the transitions between (a)(b) and (b)(c) are smooth and
37# amplitude increases/decreases linearly.
38# This module will detect artifacts in the sine wave.
39# This module also estimates the equivalent noise level by teager operator.
40# This module also detects volume changes in the sine wave. However, volume
41# changes may be affected by delay or burst.
42# Some artifacts may cause each other.
43
44# In this module, amplitude and frequency are derived from Hilbert transform.
45# Both amplitude and frequency are a function of time.
46
47# To detect each artifact, each point will be compared with
48# average amplitude of its block. The block size will be 1.5 ms.
49# Using average amplitude can mitigate the error caused by
50# Hilbert transform and noise.
51# In some case, for more accuracy, the block size may be modified
52# to other values.
53DEFAULT_BLOCK_SIZE_SECS = 0.0015
54
55# If the difference between average frequency of this block and
56# dominant frequency of full signal is less than 0.5 times of
57# dominant frequency, this block is considered to be within the
58# sine wave. In most cases, if there is no sine wave(only noise),
59# average frequency will be much greater than 5 times of
60# dominant frequency.
61# Also, for delay during playback, the frequency will be about 0
62# in perfect situation or much greater than 5 times of dominant
63# frequency if it's noised.
64DEFAULT_FREQUENCY_ERROR = 0.5
65
66# If the amplitude of some sample is less than 0.6 times of the
67# average amplitude of its left/right block, it will be considered
68# as a delay during playing.
69DEFAULT_DELAY_AMPLITUDE_THRESHOLD = 0.6
70
71# If the average amplitude of the block before or after playing
72# is more than 0.5 times to the average amplitude of the wave,
73# it will be considered as a noise artifact.
74DEFAULT_NOISE_AMPLITUDE_THRESHOLD = 0.5
75
76# In the sine wave, if the amplitude is more than 1.4 times of
77# its left side and its right side, it will be considered as
78# a burst.
79DEFAULT_BURST_AMPLITUDE_THRESHOLD = 1.4
80
81# When detecting burst, if the amplitude is lower than 0.5 times
82# average amplitude, we ignore it.
83DEFAULT_BURST_TOO_SMALL = 0.5
84
85# For a signal which is the combination of sine wave with fixed frequency f and
86# amplitude 1 and standard noise with amplitude k, the average teager value is
87# nearly linear to the noise level k.
88# Given frequency f, we simulate a sine wave with default noise level and
89# calculate its average teager value. Then, we can estimate the equivalent
90# noise level of input signal by the average teager value of input signal.
91DEFAULT_STANDARD_NOISE = 0.005
92
93# For delay, burst, volume increasing/decreasing, if two delay(
94# burst, volume increasing/decreasing) happen within
95# DEFAULT_SAME_EVENT_SECS seconds, we consider they are the
96# same event.
97DEFAULT_SAME_EVENT_SECS = 0.001
98
99# When detecting increasing/decreasing volume of signal, if the amplitude
100# is lower than 0.1 times average amplitude, we ignore it.
101DEFAULT_VOLUME_CHANGE_TOO_SMALL = 0.1
102
103# If average amplitude of right block is less/more than average
104# amplitude of left block times DEFAULT_VOLUME_CHANGE_AMPLITUDE, it will be
105# considered as decreasing/increasing on volume.
106DEFAULT_VOLUME_CHANGE_AMPLITUDE = 0.1
107
108# If the increasing/decreasing volume event is too close to the start or the end
109# of sine wave, we consider its volume change as part of rising/falling phase in
110# the start/end.
111NEAR_START_OR_END_SECS = 0.01
112
113# After applying Hilbert transform, the resulting amplitude and frequency may be
114# extremely large in the start and/or the end part. Thus, we will append zeros
115# before and after the whole wave for 0.1 secs.
116APPEND_ZEROS_SECS = 0.1
117
118# If the noise event is too close to the start or the end of the data, we
119# consider its noise as part of artifacts caused by edge effect of Hilbert
120# transform.
121# For example, originally, the data duration is 10 seconds.
122# We append 0.1 seconds of zeros in the beginning and the end of the data, so
123# the data becomes 10.2 seocnds long.
124# Then, we apply Hilbert transform to 10.2 seconds of data.
125# Near 0.1 seconds and 10.1 seconds, there will be edge effect of Hilbert
126# transform. We do not want these be treated as noise.
127# If NEAR_DATA_START_OR_END_SECS is set to 0.01, then the noise happened
128# at [0, 0.11] and [10.09, 10.1] will be ignored.
129NEAR_DATA_START_OR_END_SECS = 0.01
130
131# If the noise event is too close to the start or the end of the sine wave in
132# the data, we consider its noise as part of artifacts caused by edge effect of
133# Hilbert transform.
134# A |-------------|vvvvvvvvvvvvvvvvvvvvvvv|-------------|
135# B |ooooooooo| d |                       | d |ooooooooo|
136#
137# A is full signal. It contains a sine wave and silence before and after sine
138# wave.
139# In B, |oooo| shows the parts that we are going to check for noise before/after
140# sine wave. | d | is determined by NEAR_SINE_START_OR_END_SECS.
141NEAR_SINE_START_OR_END_SECS = 0.01
142
143
144class SineWaveNotFound(Exception):
145    """Error when there's no sine wave found in the signal"""
146    pass
147
148
149def hilbert(x):
150    """Hilbert transform copied from scipy.
151
152    More information can be found here:
153    http://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.hilbert.html
154
155    @param x: Real signal data to transform.
156
157    @returns: Analytic signal of x, we can further extract amplitude and
158              frequency from it.
159
160    """
161    x = numpy.asarray(x)
162    if numpy.iscomplexobj(x):
163        raise ValueError("x must be real.")
164    axis = -1
165    N = x.shape[axis]
166    if N <= 0:
167        raise ValueError("N must be positive.")
168
169    Xf = numpy.fft.fft(x, N, axis=axis)
170    h = numpy.zeros(N)
171    if N % 2 == 0:
172        h[0] = h[N // 2] = 1
173        h[1:N // 2] = 2
174    else:
175        h[0] = 1
176        h[1:(N + 1) // 2] = 2
177
178    if len(x.shape) > 1:
179        ind = [newaxis] * x.ndim
180        ind[axis] = slice(None)
181        h = h[ind]
182    x = numpy.fft.ifft(Xf * h, axis=axis)
183    return x
184
185
186def noised_sine_wave(frequency, rate, noise_level):
187    """Generates a sine wave of 2 second with specified noise level.
188
189    @param frequency: Frequency of sine wave.
190    @param rate: Sampling rate.
191    @param noise_level: Required noise level.
192
193    @returns: A sine wave with specified noise level.
194
195    """
196    wave = []
197    for index in range(0, rate * 2):
198        sample = 2.0 * math.pi * frequency * float(index) / float(rate)
199        sine_wave = math.sin(sample)
200        noise = noise_level * numpy.random.standard_normal()
201        wave.append(sine_wave + noise)
202    return wave
203
204
205def average_teager_value(wave, amplitude):
206    """Computes the normalized average teager value.
207
208    After averaging the teager value, we will normalize the value by
209    dividing square of amplitude.
210
211    @param wave: Wave to apply teager operator.
212    @param amplitude: Average amplitude of given wave.
213
214    @returns: Average teager value.
215
216    """
217    teager_value, length = 0, len(wave)
218    for i in range(1, length - 1):
219        ith_teager_value = abs(wave[i] * wave[i] - wave[i - 1] * wave[i + 1])
220        ith_teager_value *= max(1, abs(wave[i]))
221        teager_value += ith_teager_value
222    teager_value = (float(teager_value) / length) / (amplitude ** 2)
223    return teager_value
224
225
226def noise_level(amplitude, frequency, rate, teager_value_of_input):
227    """Computes the noise level compared with standard_noise.
228
229    For a signal which is the combination of sine wave with fixed frequency f
230    and amplitude 1 and standard noise with amplitude k, the average teager
231    value is nearly linear to the noise level k.
232    Thus, we can compute the average teager value of a sine wave with
233    standard_noise. Then, we can estimate the noise level of given input.
234
235    @param amplitude: Amplitude of input audio.
236    @param frequency: Dominant frequency of input audio.
237    @param rate: Sampling rate.
238    @param teager_value_of_input: Average teager value of input audio.
239
240    @returns: A float value denotes the audio is equivalent to have
241              how many times of noise compared with its amplitude.
242              For example, 0.02 denotes that the wave has a noise which
243              has standard distribution with standard deviation being
244              0.02 times the amplitude of the wave.
245
246    """
247    standard_noise = DEFAULT_STANDARD_NOISE
248
249    # Generates the standard sine wave with stdandard_noise level of noise.
250    standard_wave = noised_sine_wave(frequency, rate, standard_noise)
251
252    # Calculates the average teager value.
253    teager_value_of_std_wave = average_teager_value(standard_wave, amplitude)
254
255    return (teager_value_of_input / teager_value_of_std_wave) * standard_noise
256
257
258def error(f1, f2):
259    """Calculates the relative error between f1 and f2.
260
261    @param f1: Exact value.
262    @param f2: Test value.
263
264    @returns: Relative error between f1 and f2.
265
266    """
267    return abs(float(f1) - float(f2)) / float(f1)
268
269
270def hilbert_analysis(signal, rate, block_size):
271    """Finds amplitude and frequency of each time of signal by Hilbert transform.
272
273    @param signal: The wave to analyze.
274    @param rate: Sampling rate.
275    @param block_size: The size of block to transform.
276
277    @returns: A tuple of list: (amplitude, frequency) composed of
278              amplitude and frequency of each time.
279
280    """
281    # To apply Hilbert transform, the wave will be transformed
282    # segment by segment. For each segment, its size will be
283    # block_size and we will only take middle part of it.
284    # Thus, each segment looks like: |-----|=====|=====|-----|.
285    # "=...=" part will be taken while "-...-" part will be ignored.
286    #
287    # The whole size of taken part will be half of block_size
288    # which will be hilbert_block.
289    # The size of each ignored part will be half of hilbert_block
290    # which will be half_hilbert_block.
291    hilbert_block = block_size // 2
292    half_hilbert_block = hilbert_block // 2
293    # As mentioned above, for each block, we will only take middle
294    # part of it. Thus, the whole transformation will be completed as:
295    # |=====|=====|-----|           |-----|=====|=====|-----|
296    #       |-----|=====|=====|-----|           |-----|=====|=====|
297    #                   |-----|=====|=====|-----|
298    # Specially, beginning and ending part may not have ignored part.
299    length = len(signal)
300    result = []
301    for left_border in range(0, length, hilbert_block):
302        right_border = min(length, left_border + hilbert_block)
303        temp_left_border = max(0, left_border - half_hilbert_block)
304        temp_right_border = min(length, right_border + half_hilbert_block)
305        temp = hilbert(signal[temp_left_border:temp_right_border])
306        for index in range(left_border, right_border):
307            result.append(temp[index - temp_left_border])
308    result = numpy.asarray(result)
309    amplitude = numpy.abs(result)
310    phase = numpy.unwrap(numpy.angle(result))
311    frequency = numpy.diff(phase) / (2.0 * numpy.pi) * rate
312    #frequency.append(frequency[len(frequency)-1])
313    frequecny = numpy.append(frequency, frequency[len(frequency) - 1])
314    return (amplitude, frequency)
315
316
317def find_block_average_value(arr, side_block_size, block_size):
318    """For each index, finds average value of its block, left block, right block.
319
320    It will find average value for each index in the range.
321
322    For each index, the range of its block is
323        [max(0, index - block_size / 2), min(length - 1, index + block_size / 2)]
324    For each index, the range of its left block is
325        [max(0, index - size_block_size), index]
326    For each index, the range of its right block is
327        [index, min(length - 1, index + side_block_size)]
328
329    @param arr: The array to be computed.
330    @param side_block_size: the size of the left_block and right_block.
331    @param block_size: the size of the block.
332
333    @returns: A tuple of lists: (left_block_average_array,
334                                 right_block_average_array,
335                                 block_average_array)
336    """
337    length = len(arr)
338    left_border, right_border = 0, 1
339    left_block_sum = arr[0]
340    right_block_sum = arr[0]
341    left_average_array = numpy.zeros(length)
342    right_average_array = numpy.zeros(length)
343    block_average_array = numpy.zeros(length)
344    for index in range(0, length):
345        while left_border < index - side_block_size:
346            left_block_sum -= arr[left_border]
347            left_border += 1
348        while right_border < min(length, index + side_block_size):
349            right_block_sum += arr[right_border]
350            right_border += 1
351
352        left_average_value = float(left_block_sum) / (index - left_border + 1)
353        right_average_value = float(right_block_sum) / (right_border - index)
354        left_average_array[index] = left_average_value
355        right_average_array[index] = right_average_value
356
357        if index + 1 < length:
358            left_block_sum += arr[index + 1]
359        right_block_sum -= arr[index]
360    left_border, right_border = 0, 1
361    block_sum = 0
362    for index in range(0, length):
363        while left_border < index - block_size // 2:
364            block_sum -= arr[left_border]
365            left_border += 1
366        while right_border < min(length, index + block_size // 2):
367            block_sum += arr[right_border]
368            right_border += 1
369
370        average_value = float(block_sum) / (right_border - left_border)
371        block_average_array[index] = average_value
372    return (left_average_array, right_average_array, block_average_array)
373
374
375def find_start_end_index(dominant_frequency,
376                         block_frequency_delta,
377                         block_size,
378                         frequency_error_threshold):
379    """Finds start and end index of sine wave.
380
381    For each block with size of block_size, we check that whether its frequency
382    is close enough to the dominant_frequency. If yes, we will consider this
383    block to be within the sine wave.
384    Then, it will return the start and end index of sine wave indicating that
385    sine wave is between [start_index, end_index)
386    It's okay if the whole signal only contains sine wave.
387
388    @param dominant_frequency: Dominant frequency of signal.
389    @param block_frequency_delta: Average absolute difference between dominant
390                                  frequency and frequency of each block. For
391                                  each index, its block is
392                                  [max(0, index - block_size / 2),
393                                   min(length - 1, index + block_size / 2)]
394    @param block_size: Block size in samples.
395
396    @returns: A tuple composed of (start_index, end_index)
397
398    """
399    length = len(block_frequency_delta)
400
401    # Finds the start/end time index of playing based on dominant frequency
402    start_index, end_index = length - 1 , 0
403    for index in range(0, length):
404        left_border = max(0, index - block_size // 2)
405        right_border = min(length - 1, index + block_size // 2)
406        frequency_error = block_frequency_delta[index] / dominant_frequency
407        if frequency_error < frequency_error_threshold:
408            start_index = min(start_index, left_border)
409            end_index = max(end_index, right_border + 1)
410    return (start_index, end_index)
411
412
413def noise_detection(start_index, end_index,
414                    block_amplitude, average_amplitude,
415                    rate, noise_amplitude_threshold):
416    """Detects noise before/after sine wave.
417
418    If average amplitude of some sample's block before start of wave or after
419    end of wave is more than average_amplitude times noise_amplitude_threshold,
420    it will be considered as a noise.
421
422    @param start_index: Start index of sine wave.
423    @param end_index: End index of sine wave.
424    @param block_amplitude: An array for average amplitude of each block, where
425                            amplitude is computed from Hilbert transform.
426    @param average_amplitude: Average amplitude of sine wave.
427    @param rate: Sampling rate
428    @param noise_amplitude_threshold: If the average amplitude of a block is
429                        higher than average amplitude of the wave times
430                        noise_amplitude_threshold, it will be considered as
431                        noise before/after playback.
432
433    @returns: A tuple of lists indicating the time that noise happens:
434            (noise_before_playing, noise_after_playing).
435
436    """
437    length = len(block_amplitude)
438    amplitude_threshold = average_amplitude * noise_amplitude_threshold
439    same_event_samples = rate * DEFAULT_SAME_EVENT_SECS
440
441    # Detects noise before playing.
442    noise_time_point = []
443    last_noise_end_time_point = []
444    previous_noise_index = None
445    times = 0
446    for index in range(0, length):
447        # Ignore noise too close to the beginning or the end of sine wave.
448        # Check the docstring of NEAR_SINE_START_OR_END_SECS.
449        if ((start_index - rate * NEAR_SINE_START_OR_END_SECS) <= index and
450            (index < end_index + rate * NEAR_SINE_START_OR_END_SECS)):
451            continue
452
453        # Ignore noise too close to the beginning or the end of original data.
454        # Check the docstring of NEAR_DATA_START_OR_END_SECS.
455        if (float(index) / rate <=
456            NEAR_DATA_START_OR_END_SECS + APPEND_ZEROS_SECS):
457            continue
458        if (float(length - index) / rate <=
459            NEAR_DATA_START_OR_END_SECS + APPEND_ZEROS_SECS):
460            continue
461        if block_amplitude[index] > amplitude_threshold:
462            same_event = False
463            if previous_noise_index:
464                same_event = (index - previous_noise_index) < same_event_samples
465            if not same_event:
466                index_start_sec = float(index) / rate - APPEND_ZEROS_SECS
467                index_end_sec = float(index + 1) / rate - APPEND_ZEROS_SECS
468                noise_time_point.append(index_start_sec)
469                last_noise_end_time_point.append(index_end_sec)
470                times += 1
471            index_end_sec = float(index + 1) / rate - APPEND_ZEROS_SECS
472            last_noise_end_time_point[times - 1] = index_end_sec
473            previous_noise_index = index
474
475    noise_before_playing, noise_after_playing = [], []
476    for i in range(times):
477        duration = last_noise_end_time_point[i] - noise_time_point[i]
478        if noise_time_point[i] < float(start_index) / rate - APPEND_ZEROS_SECS:
479            noise_before_playing.append((noise_time_point[i], duration))
480        else:
481            noise_after_playing.append((noise_time_point[i], duration))
482
483    return (noise_before_playing, noise_after_playing)
484
485
486def delay_detection(start_index, end_index,
487                    block_amplitude, average_amplitude,
488                    dominant_frequency,
489                    rate,
490                    left_block_amplitude,
491                    right_block_amplitude,
492                    block_frequency_delta,
493                    delay_amplitude_threshold,
494                    frequency_error_threshold):
495    """Detects delay during playing.
496
497    For each sample, we will check whether the average amplitude of its block
498    is less than average amplitude of its left block and its right block times
499    delay_amplitude_threshold. Also, we will check whether the frequency of
500    its block is far from the dominant frequency.
501    If at least one constraint fulfilled, it will be considered as a delay.
502
503    @param start_index: Start index of sine wave.
504    @param end_index: End index of sine wave.
505    @param block_amplitude: An array for average amplitude of each block, where
506                            amplitude is computed from Hilbert transform.
507    @param average_amplitude: Average amplitude of sine wave.
508    @param dominant_frequency: Dominant frequency of signal.
509    @param rate: Sampling rate
510    @param left_block_amplitude: Average amplitude of left block of each index.
511                                Ref to find_block_average_value function.
512    @param right_block_amplitude: Average amplitude of right block of each index.
513                                Ref to find_block_average_value function.
514    @param block_frequency_delta: Average absolute difference frequency to
515                                dominant frequency of block of each index.
516                                Ref to find_block_average_value function.
517    @param delay_amplitude_threshold: If the average amplitude of a block is
518                        lower than average amplitude of the wave times
519                        delay_amplitude_threshold, it will be considered
520                        as delay.
521    @param frequency_error_threshold: Ref to DEFAULT_FREQUENCY_ERROR
522
523    @returns: List of delay occurrence:
524                [(time_1, duration_1), (time_2, duration_2), ...],
525              where time and duration are in seconds.
526
527    """
528    delay_time_points = []
529    last_delay_end_time_points = []
530    previous_delay_index = None
531    times = 0
532    same_event_samples = rate * DEFAULT_SAME_EVENT_SECS
533    start_time = float(start_index) / rate - APPEND_ZEROS_SECS
534    end_time = float(end_index) / rate - APPEND_ZEROS_SECS
535    for index in range(start_index, end_index):
536        if block_amplitude[index] > average_amplitude * delay_amplitude_threshold:
537            continue
538        now_time = float(index) / rate - APPEND_ZEROS_SECS
539        if abs(now_time - start_time) < NEAR_START_OR_END_SECS:
540            continue
541        if abs(now_time - end_time) < NEAR_START_OR_END_SECS:
542            continue
543        # If amplitude less than its left/right side and small enough,
544        # it will be considered as a delay.
545        amp_threshold = average_amplitude * delay_amplitude_threshold
546        left_threshold = delay_amplitude_threshold * left_block_amplitude[index]
547        amp_threshold = min(amp_threshold, left_threshold)
548        right_threshold = delay_amplitude_threshold * right_block_amplitude[index]
549        amp_threshold = min(amp_threshold, right_threshold)
550
551        frequency_error = block_frequency_delta[index] / dominant_frequency
552
553        amplitude_too_small = block_amplitude[index] < amp_threshold
554        frequency_not_match = frequency_error > frequency_error_threshold
555
556        if amplitude_too_small or frequency_not_match:
557            same_event = False
558            if previous_delay_index:
559                same_event = (index - previous_delay_index) < same_event_samples
560            if not same_event:
561                index_start_sec = float(index) / rate - APPEND_ZEROS_SECS
562                index_end_sec = float(index + 1) / rate - APPEND_ZEROS_SECS
563                delay_time_points.append(index_start_sec)
564                last_delay_end_time_points.append(index_end_sec)
565                times += 1
566            previous_delay_index = index
567            index_end_sec = float(index + 1) / rate - APPEND_ZEROS_SECS
568            last_delay_end_time_points[times - 1] = index_end_sec
569
570    delay_list = []
571    for i in range(len(delay_time_points)):
572        duration = last_delay_end_time_points[i] - delay_time_points[i]
573        delay_list.append( (delay_time_points[i], duration) )
574    return delay_list
575
576
577def burst_detection(start_index, end_index,
578                    block_amplitude, average_amplitude,
579                    dominant_frequency,
580                    rate,
581                    left_block_amplitude,
582                    right_block_amplitude,
583                    block_frequency_delta,
584                    burst_amplitude_threshold,
585                    frequency_error_threshold):
586    """Detects burst during playing.
587
588    For each sample, we will check whether the average amplitude of its block is
589    more than average amplitude of its left block and its right block times
590    burst_amplitude_threshold. Also, we will check whether the frequency of
591    its block is not compatible to the dominant frequency.
592    If at least one constraint fulfilled, it will be considered as a burst.
593
594    @param start_index: Start index of sine wave.
595    @param end_index: End index of sine wave.
596    @param block_amplitude: An array for average amplitude of each block, where
597                            amplitude is computed from Hilbert transform.
598    @param average_amplitude: Average amplitude of sine wave.
599    @param dominant_frequency: Dominant frequency of signal.
600    @param rate: Sampling rate
601    @param left_block_amplitude: Average amplitude of left block of each index.
602                                Ref to find_block_average_value function.
603    @param right_block_amplitude: Average amplitude of right block of each index.
604                                Ref to find_block_average_value function.
605    @param block_frequency_delta: Average absolute difference frequency to
606                                dominant frequency of block of each index.
607    @param burst_amplitude_threshold: If the amplitude is higher than average
608                            amplitude of its left block and its right block
609                            times burst_amplitude_threshold. It will be
610                            considered as a burst.
611    @param frequency_error_threshold: Ref to DEFAULT_FREQUENCY_ERROR
612
613    @returns: List of burst occurence: [time_1, time_2, ...],
614              where time is in seconds.
615
616    """
617    burst_time_points = []
618    previous_burst_index = None
619    same_event_samples = rate * DEFAULT_SAME_EVENT_SECS
620    for index in range(start_index, end_index):
621        # If amplitude higher than its left/right side and large enough,
622        # it will be considered as a burst.
623        if block_amplitude[index] <= average_amplitude * DEFAULT_BURST_TOO_SMALL:
624            continue
625        if abs(index - start_index) < rate * NEAR_START_OR_END_SECS:
626            continue
627        if abs(index - end_index) < rate * NEAR_START_OR_END_SECS:
628            continue
629        amp_threshold = average_amplitude * DEFAULT_BURST_TOO_SMALL
630        left_threshold = burst_amplitude_threshold * left_block_amplitude[index]
631        amp_threshold = max(amp_threshold, left_threshold)
632        right_threshold = burst_amplitude_threshold * right_block_amplitude[index]
633        amp_threshold = max(amp_threshold, right_threshold)
634
635        frequency_error = block_frequency_delta[index] / dominant_frequency
636
637        amplitude_too_large = block_amplitude[index] > amp_threshold
638        frequency_not_match = frequency_error > frequency_error_threshold
639
640        if amplitude_too_large or frequency_not_match:
641            same_event = False
642            if previous_burst_index:
643                same_event = index - previous_burst_index < same_event_samples
644            if not same_event:
645                burst_time_points.append(float(index) / rate - APPEND_ZEROS_SECS)
646            previous_burst_index = index
647
648    return burst_time_points
649
650
651def changing_volume_detection(start_index, end_index,
652                              average_amplitude,
653                              rate,
654                              left_block_amplitude,
655                              right_block_amplitude,
656                              volume_changing_amplitude_threshold):
657    """Finds volume changing during playback.
658
659    For each index, we will compare average amplitude of its left block and its
660    right block. If average amplitude of right block is more than average
661    amplitude of left block times (1 + DEFAULT_VOLUME_CHANGE_AMPLITUDE), it will
662    be considered as an increasing volume. If the one of right block is less
663    than that of left block times (1 - DEFAULT_VOLUME_CHANGE_AMPLITUDE), it will
664    be considered as a decreasing volume.
665
666    @param start_index: Start index of sine wave.
667    @param end_index: End index of sine wave.
668    @param average_amplitude: Average amplitude of sine wave.
669    @param rate: Sampling rate
670    @param left_block_amplitude: Average amplitude of left block of each index.
671                                Ref to find_block_average_value function.
672    @param right_block_amplitude: Average amplitude of right block of each index.
673                                Ref to find_block_average_value function.
674    @param volume_changing_amplitude_threshold: If the average amplitude of right
675                                                block is higher or lower than
676                                                that of left one times this
677                                                value, it will be considered as
678                                                a volume change.
679                                                Also refer to
680                                                DEFAULT_VOLUME_CHANGE_AMPLITUDE
681
682    @returns: List of volume changing composed of 1 for increasing and
683              -1 for decreasing.
684
685    """
686    length = len(left_block_amplitude)
687
688    # Detects rising and/or falling volume.
689    previous_rising_index, previous_falling_index = None, None
690    changing_time = []
691    changing_events = []
692    amplitude_threshold = average_amplitude * DEFAULT_VOLUME_CHANGE_TOO_SMALL
693    same_event_samples = rate * DEFAULT_SAME_EVENT_SECS
694    for index in range(start_index, end_index):
695        # Skips if amplitude is too small.
696        if left_block_amplitude[index] < amplitude_threshold:
697            continue
698        if right_block_amplitude[index] < amplitude_threshold:
699            continue
700        # Skips if changing is from start or end time
701        if float(abs(start_index - index)) / rate < NEAR_START_OR_END_SECS:
702            continue
703        if float(abs(end_index   - index)) / rate < NEAR_START_OR_END_SECS:
704            continue
705
706        delta_margin = volume_changing_amplitude_threshold
707        if left_block_amplitude[index] > 0:
708            delta_margin *= left_block_amplitude[index]
709
710        increasing_threshold = left_block_amplitude[index] + delta_margin
711        decreasing_threshold = left_block_amplitude[index] - delta_margin
712
713        if right_block_amplitude[index] > increasing_threshold:
714            same_event = False
715            if previous_rising_index:
716                same_event = index - previous_rising_index < same_event_samples
717            if not same_event:
718                changing_time.append(float(index) / rate - APPEND_ZEROS_SECS)
719                changing_events.append(+1)
720            previous_rising_index = index
721        if right_block_amplitude[index] < decreasing_threshold:
722            same_event = False
723            if previous_falling_index:
724                same_event = index - previous_falling_index < same_event_samples
725            if not same_event:
726                changing_time.append(float(index) / rate - APPEND_ZEROS_SECS)
727                changing_events.append(-1)
728            previous_falling_index = index
729
730    # Combines consecutive increasing/decreasing event.
731    combined_changing_events, prev = [], 0
732    for i in range(len(changing_events)):
733        if changing_events[i] == prev:
734            continue
735        combined_changing_events.append((changing_time[i], changing_events[i]))
736        prev = changing_events[i]
737    return combined_changing_events
738
739
740def quality_measurement(
741        signal, rate,
742        dominant_frequency=None,
743        block_size_secs=DEFAULT_BLOCK_SIZE_SECS,
744        frequency_error_threshold=DEFAULT_FREQUENCY_ERROR,
745        delay_amplitude_threshold=DEFAULT_DELAY_AMPLITUDE_THRESHOLD,
746        noise_amplitude_threshold=DEFAULT_NOISE_AMPLITUDE_THRESHOLD,
747        burst_amplitude_threshold=DEFAULT_BURST_AMPLITUDE_THRESHOLD,
748        volume_changing_amplitude_threshold=DEFAULT_VOLUME_CHANGE_AMPLITUDE):
749    """Detects several artifacts and estimates the noise level.
750
751    This method detects artifact before playing, after playing, and delay
752    during playing. Also, it estimates the noise level of the signal.
753    To avoid the influence of noise, it calculates amplitude and frequency
754    block by block.
755
756    @param signal: A list of numbers for one-channel PCM data. The data should
757                   be normalized to [-1,1].
758    @param rate: Sampling rate
759    @param dominant_frequency: Dominant frequency of signal. Set None to
760                               recalculate the frequency in this function.
761    @param block_size_secs: Block size in seconds. The measurement will be done
762                            block-by-block using average amplitude and frequency
763                            in each block to avoid noise.
764    @param frequency_error_threshold: Ref to DEFAULT_FREQUENCY_ERROR.
765    @param delay_amplitude_threshold: If the average amplitude of a block is
766                                      lower than average amplitude of the wave
767                                      times delay_amplitude_threshold, it will
768                                      be considered as delay.
769                                      Also refer to delay_detection and
770                                      DEFAULT_DELAY_AMPLITUDE_THRESHOLD.
771    @param noise_amplitude_threshold: If the average amplitude of a block is
772                                      higher than average amplitude of the wave
773                                      times noise_amplitude_threshold, it will
774                                      be considered as noise before/after
775                                      playback.
776                                      Also refer to noise_detection and
777                                      DEFAULT_NOISE_AMPLITUDE_THRESHOLD.
778    @param burst_amplitude_threshold: If the average amplitude of a block is
779                                      higher than average amplitude of its left
780                                      block and its right block times
781                                      burst_amplitude_threshold. It will be
782                                      considered as a burst.
783                                      Also refer to burst_detection and
784                                      DEFAULT_BURST_AMPLITUDE_THRESHOLD.
785    @param volume_changing_amplitude_threshold: If the average amplitude of right
786                                                block is higher or lower than
787                                                that of left one times this
788                                                value, it will be considered as
789                                                a volume change.
790                                                Also refer to
791                                                changing_volume_detection and
792                                                DEFAULT_VOLUME_CHANGE_AMPLITUDE
793
794    @returns: A dictoinary of detection/estimation:
795              {'artifacts':
796                {'noise_before_playback':
797                    [(time_1, duration_1), (time_2, duration_2), ...],
798                 'noise_after_playback':
799                    [(time_1, duration_1), (time_2, duration_2), ...],
800                 'delay_during_playback':
801                    [(time_1, duration_1), (time_2, duration_2), ...],
802                 'burst_during_playback':
803                    [time_1, time_2, ...]
804                },
805               'volume_changes':
806                 [(time_1, flag_1), (time_2, flag_2), ...],
807               'equivalent_noise_level': level
808              }
809              where durations and time points are in seconds. And,
810              equivalence_noise_level is the quotient of noise and
811              wave which refers to DEFAULT_STANDARD_NOISE.
812              volume_changes is a list of tuples containing time
813              stamps and decreasing/increasing flags for volume
814              change events.
815
816    """
817    # Calculates the block size, from seconds to samples.
818    block_size = int(block_size_secs * rate)
819
820    signal = numpy.concatenate((numpy.zeros(int(rate * APPEND_ZEROS_SECS)),
821                                signal,
822                                numpy.zeros(int(rate * APPEND_ZEROS_SECS))))
823    signal = numpy.array(signal, dtype=float)
824    length = len(signal)
825
826    # Calculates the amplitude and frequency.
827    amplitude, frequency = hilbert_analysis(signal, rate, block_size)
828
829    # Finds the dominant frequency.
830    if not dominant_frequency:
831        dominant_frequency = audio_analysis.spectral_analysis(signal, rate)[0][0]
832
833    # Finds the array which contains absolute difference between dominant
834    # frequency and frequency at each time point.
835    frequency_delta = abs(frequency - dominant_frequency)
836
837    # Computes average amplitude of each type of block
838    res = find_block_average_value(amplitude, block_size * 2, block_size)
839    left_block_amplitude, right_block_amplitude, block_amplitude = res
840
841    # Computes average absolute difference of frequency and dominant frequency
842    # of the block of each index
843    _, _, block_frequency_delta = find_block_average_value(frequency_delta,
844                                                           block_size * 2,
845                                                           block_size)
846
847    # Finds start and end index of sine wave.
848    start_index, end_index = find_start_end_index(dominant_frequency,
849                                                  block_frequency_delta,
850                                                  block_size,
851                                                  frequency_error_threshold)
852
853    if start_index > end_index:
854        raise SineWaveNotFound('No sine wave found in signal')
855
856    logging.debug('Found sine wave: start: %s, end: %s',
857                  float(start_index) / rate - APPEND_ZEROS_SECS,
858                  float(end_index) / rate - APPEND_ZEROS_SECS)
859
860    sum_of_amplitude = float(sum(amplitude[start_index:end_index]))
861    # Finds average amplitude of sine wave.
862    average_amplitude = sum_of_amplitude / (end_index - start_index)
863
864    # Finds noise before and/or after playback.
865    noise_before_playing, noise_after_playing = noise_detection(
866            start_index, end_index,
867            block_amplitude, average_amplitude,
868            rate,
869            noise_amplitude_threshold)
870
871    # Finds delay during playback.
872    delays = delay_detection(start_index, end_index,
873                             block_amplitude, average_amplitude,
874                             dominant_frequency,
875                             rate,
876                             left_block_amplitude,
877                             right_block_amplitude,
878                             block_frequency_delta,
879                             delay_amplitude_threshold,
880                             frequency_error_threshold)
881
882    # Finds burst during playback.
883    burst_time_points = burst_detection(start_index, end_index,
884                                        block_amplitude, average_amplitude,
885                                        dominant_frequency,
886                                        rate,
887                                        left_block_amplitude,
888                                        right_block_amplitude,
889                                        block_frequency_delta,
890                                        burst_amplitude_threshold,
891                                        frequency_error_threshold)
892
893    # Finds volume changing during playback.
894    volume_changes = changing_volume_detection(
895            start_index, end_index,
896            average_amplitude,
897            rate,
898            left_block_amplitude,
899            right_block_amplitude,
900            volume_changing_amplitude_threshold)
901
902    # Calculates the average teager value.
903    teager_value = average_teager_value(signal[start_index:end_index],
904                                        average_amplitude)
905
906    # Finds out the noise level.
907    noise = noise_level(average_amplitude, dominant_frequency,
908                        rate,
909                        teager_value)
910
911    return {'artifacts':
912            {'noise_before_playback': noise_before_playing,
913             'noise_after_playback': noise_after_playing,
914             'delay_during_playback': delays,
915             'burst_during_playback': burst_time_points
916            },
917            'volume_changes': volume_changes,
918            'equivalent_noise_level': noise
919           }
920