Module thunderfish.harmonics
Extract and analyze harmonic frequencies from power spectra.
Harmonic group extraction
harmonic_groups()
: detect peaks in power spectrum and group them according to their harmonic structure.expand_group()
: add more harmonics to harmonic group.extract_fundamentals()
: collect harmonic groups from lists of power spectrum peaks.threshold_estimate()
: estimates thresholds for peak detection in a power spectrum.
Helper functions for harmonic group extraction
build_harmonic_group()
: find all the harmonics belonging to the largest peak in the power spectrum.retrieve_harmonic_group()
: find all the harmonics belonging to a given fundamental.group_candidate()
: candidate harmonic frequencies belonging to a fundamental frequency.update_group()
: update frequency lists and harmonic group.
Handling of lists of harmonic groups
fundamental_freqs()
: extract fundamental frequencies from lists of harmonic groups.fundamental_freqs_and_power()
: extract fundamental frequencies and their power in dB from lists of harmonic groups.
Handling of lists of fundamental frequencies
add_relative_power()
: add a column with relative power.add_power_ranks()
: add a column with power ranks.similar_indices()
: indices of similar frequencies.unique_mask()
: mark similar frequencies from different recordings as dublicate.unique()
: remove similar frequencies from different recordings.
Visualization
colors_markers()
: Generate a list of colors and markers for plotting.plot_harmonic_groups()
: Mark decibel power of fundamentals and their harmonics.plot_psd_harmonic_groups()
: Plot decibel powerspectrum with detected peaks, harmonic groups, and mains frequencies.
Configuration
add_psd_peak_detection_config()
: add parameters for the detection of peaks in power spectra to configuration.psd_peak_detection_args()
: retrieve parameters for the detection of peaks in power spectra from configuration.add_harmonic_groups_config()
: add parameters for the detection of harmonic groups to configuration.harmonic_groups_args()
: retrieve parameters for the detection of harmonic groups from configuration.
Expand source code
"""
Extract and analyze harmonic frequencies from power spectra.
## Harmonic group extraction
 `harmonic_groups()`: detect peaks in power spectrum and group them
according to their harmonic structure.
 `expand_group()`: add more harmonics to harmonic group.
 `extract_fundamentals()`: collect harmonic groups from
lists of power spectrum peaks.
 `threshold_estimate()`: estimates thresholds for peak detection
in a power spectrum.
## Helper functions for harmonic group extraction
 `build_harmonic_group()`: find all the harmonics belonging to the largest peak in the power spectrum.
 `retrieve_harmonic_group()`: find all the harmonics belonging to a given fundamental.
 `group_candidate()`: candidate harmonic frequencies belonging to a fundamental frequency.
 `update_group()`: update frequency lists and harmonic group.
## Handling of lists of harmonic groups
 `fundamental_freqs()`: extract fundamental frequencies from
lists of harmonic groups.
 `fundamental_freqs_and_power()`: extract fundamental frequencies and their
power in dB from lists of harmonic groups.
## Handling of lists of fundamental frequencies
 `add_relative_power()`: add a column with relative power.
 `add_power_ranks()`: add a column with power ranks.
 `similar_indices()`: indices of similar frequencies.
 `unique_mask()`: mark similar frequencies from different recordings as dublicate.
 `unique()`: remove similar frequencies from different recordings.
## Visualization
 `colors_markers()`: Generate a list of colors and markers for plotting.
 `plot_harmonic_groups()`: Mark decibel power of fundamentals and their
harmonics.
 `plot_psd_harmonic_groups()`: Plot decibel powerspectrum with detected peaks,
harmonic groups, and mains frequencies.
## Configuration
 `add_psd_peak_detection_config()`: add parameters for the detection of
peaks in power spectra to configuration.
 `psd_peak_detection_args()`: retrieve parameters for the detection of peaks
in power spectra from configuration.
 `add_harmonic_groups_config()`: add parameters for the detection of
harmonic groups to configuration.
 `harmonic_groups_args()`: retrieve parameters for the detection of
harmonic groups from configuration.
"""
import math as m
import numpy as np
import scipy.signal as sig
from .eventdetection import detect_peaks, trim, hist_threshold
from .powerspectrum import decibel, power, plot_decibel_psd
try:
import matplotlib.cm as cm
import matplotlib.colors as mc
except ImportError:
pass
def group_candidate(good_freqs, all_freqs, freq, divisor,
freq_tol, max_freq_tol, min_group_size, verbose):
"""Candidate harmonic frequencies belonging to a fundamental frequency.
Parameters

good_freqs: 2D array
Frequency, power, and use count (columns) of strong peaks detected
in a power spectrum.
all_freqs: 2D array
Frequency, power, and use count (columns) of all peaks detected
in a power spectrum.
freq: float
Fundamental frequency for which a harmonic group should be assembled.
divisor: int
Fundamental frequency was obtained from a frequency divided by divisor.
0: Fundamental frequency is given as is.
freq_tol: float
Harmonics should fall within this frequency tolerance.
This should be in the range of the frequency resolution
and should not be smaller than half of the frequency resolution.
max_freq_tol: float
Maximum deviation of harmonics from their expected frequency.
Peaks with frequencies between `freq_tol` and `max_freq_tol`
get penalized the further away from their expected frequency.
min_group_size: int
Minimum number of harmonics of a harmonic group.
The harmonics from min_group_size/3 to max(min_group_size, divisor)
need to be in good_freqs.
verbose: int
Verbosity level.
Returns

new_group: 1D aray of indices
Frequencies in all_freqs belonging to the candidate group.
fzero: float
Adjusted fundamental frequency. If negative, no group was found.
fzero_harmonics: int
The highest harmonics that was used to recompute the fundamental frequency.
"""
if verbose > 0:
print('')
dvs = True
if divisor <= 0:
divisor = 1
dvs = False
# 1. find harmonics in good_freqs and adjust fzero accordingly:
fzero = freq
fzero_harmonics = 1
if len(good_freqs[:,0]) > 0:
"""
ff = good_freqs[:,0]
h = np.round(ff/fzero)
h_indices = np.unique(h, return_index=True)[1]
f_best = [hi+np.argmin(np.abs(fh/h  fzero)) for hi, fh in zip(h_indices, np.split(ff, h_indices))]
fe = np.abs(ff/h  fzero)
h = h[fe < freq_tol]
fe = fe[fe < freq_tol]
"""
prev_freq = divisor * freq
for h in range(divisor+1, 2*min_group_size+1):
idx = np.argmin(np.abs(good_freqs[:,0]/h  fzero))
ff = good_freqs[idx,0]
if m.fabs(ff/h  fzero) > freq_tol:
continue
df = ff  prev_freq
dh = np.round(df/fzero)
fe = m.fabs(df/dh  fzero)
if fe > 2.0*freq_tol:
if h > min_group_size:
break
continue
# update fzero:
prev_freq = ff
fzero_harmonics = h
fzero = ff/fzero_harmonics
if verbose > 1:
print('adjusted fzero to %.2fHz from %dth harmonic' % (fzero, h))
if verbose > 0:
ds = 'divisor: %d, ' % divisor if dvs else ''
print('# %sfzero=%7.2fHz adjusted from harmonics %d'
% (ds, fzero, fzero_harmonics))
# 2. check fzero:
# freq might not be in our group anymore, because fzero was adjusted:
if np.abs(freq  fzero) > freq_tol:
if verbose > 0:
print(' discarded: lost frequency')
return [], 1.0, fzero_harmonics
# 3. collect harmonics from all_freqs:
new_group = np.ones(min_group_size, dtype=np.int)
new_penalties = np.ones(min_group_size)
freqs = []
prev_h = 0
prev_fe = 0.0
for h in range(1, min_group_size+1):
penalty = 0
i = np.argmin(np.abs(all_freqs[:,0]/h  fzero))
f = all_freqs[i,0]
if verbose > 2:
print(' check %7.2fHz as %d. harmonics' % (f, h))
fac = 1.0 if h >= divisor else 2.0
fe = m.fabs(f/h  fzero)
if fe > fac*max_freq_tol:
if verbose > 1 and fe < 2*fac*max_freq_tol:
print(' %d. harmonics at %7.2fHz is off by %7.2fHz (max between %5.2fHz and %5.2fHz) from %7.2fHz'
% (h, f, h*fe, h*fac*freq_tol, h*fac*max_freq_tol, h*fzero))
continue
if fe > fac*freq_tol:
penalty = np.interp(fe, [fac*freq_tol, fac*max_freq_tol], [0.0, 1.0])
if len(freqs) > 0:
pf = freqs[1]
df = f  pf
if df < 0.5*fzero:
if len(freqs)>1:
pf = freqs[2]
df = f  pf
else:
pf = 0.0
df = h*fzero
dh = m.floor(df/fzero + 0.5)
fe = m.fabs(df/dh  fzero)
if fe > 2*dh*fac*max_freq_tol:
if verbose > 1:
print(' %d. harmonics at %7.2fHz is off by %7.2fHz (max between %5.2fHz and %5.2fHz) from previous harmonics at %7.2fHz'
% (h, f, dh*fe, 2*fac*dh*freq_tol, 2*fac*dh*max_freq_tol, pf))
continue
if fe > 2*dh*fac*freq_tol:
penalty = np.interp(fe, [2*dh*fac*freq_tol, 2*dh*fac*max_freq_tol], [0.0, 1.0])
else:
fe = 0.0
if h > prev_h or fe < prev_fe:
if prev_h > 0 and h  prev_h > 1:
if verbose > 1:
print(' previous harmonics %d more than 1 away from %d. harmonics at %7.2fHz'
% (prev_h, h, f))
break
if h == prev_h and len(freqs) > 0:
freqs.pop()
freqs.append(f)
new_group[int(h)1] = i
new_penalties[int(h)1] = penalty
prev_h = h
prev_fe = fe
if verbose > 1:
print(' %d. harmonics at %7.2fHz has been taken (peak %2d) with penalty %3.1f' % (h, f, i, penalty))
# 4. check new group:
# almost all harmonics in min_group_size required:
max_penalties = min_group_size/3
if np.sum(new_penalties) > max_penalties:
if verbose > 0:
print(' discarded group because sum of penalties %3.1f is more than %3.1f: indices' %
(np.sum(new_penalties), max_penalties), new_group, ' penalties', new_penalties)
return [], 1.0, fzero_harmonics
new_group = new_group[new_group>=0]
# check use count of frequencies:
double_use = np.sum(all_freqs[new_group, 2]>0)
if double_use >= 2: # XXX make this a parameter?
if verbose > 0:
print(' discarded group because of use count = %2d >= 2' % double_use)
return [], 1.0, fzero_harmonics
# 5. return results:
return new_group, fzero, fzero_harmonics
def update_group(good_freqs, all_freqs, new_group, fzero,
freq_tol, verbose, group_str):
"""Update good frequencies and harmonic group.
Remove frequencies from good_freqs, add missing fundamental to group.
Parameters

good_freqs: 2D array
Frequency, power, and use count (columns) of strong peaks detected
in a power spectrum.
all_freqs: 2D array
Frequency, power, and use count (columns) of all peaks detected
in a power spectrum.
new_group: 1D aray of indices
Frequencies in all_freqs of an harmonic group.
fzero: float
Fundamental frequency for which frequencies are collected in good_freqs.
freq_tol: float
Harmonics need to fall within this frequency tolerance.
This should be in the range of the frequency resolution
and not be smaller than half of the frequency resolution.
verbose: int
Verbosity level.
group_str: string
String for debug message.
Returns

good_freqs: 2D array
Frequency, power, and use count (columns) of strong peaks detected
in a power spectrum with frequencies for harmonic group
of fundamental frequency fzero removed.
group: 2D array
Frequency, power, and use count (columns) of harmonic group
for fundamental frequency fzero.
"""
# initialize group:
group = all_freqs[new_group,:]
# indices of group in good_freqs:
freq_tol *= 1.1
indices = []
for f in group[:,0]:
idx = np.argmin(np.abs(good_freqs[:,0]f))
if np.abs(good_freqs[idx,0]f) <= freq_tol:
indices.append(idx)
indices = np.asarray(indices, dtype=np.int)
# harmonics in good_freqs:
nharm = np.round(good_freqs[:,0]/fzero)
idxs = np.where(np.abs(good_freqs[:,0]  nharm*fzero) <= freq_tol)[0]
indices = np.unique(np.concatenate((indices, idxs)))
# report:
if verbose > 1:
print('# good freqs: ', indices,
'[', ', '.join(['%.2f' % f for f in good_freqs[indices,0]]), ']')
print('# all freqs : ', new_group,
'[', ', '.join(['%.2f' % f for f in all_freqs[new_group,0]]), ']')
if verbose > 0:
refi = np.argmax(group[:,1] > 0.0)
print('')
print(group_str)
for i in range(len(group)):
print('f=%8.2fHz n=%5.2f: power=%9.3g power/pmax=%6.4f=%5.1fdB'
% (group[i,0], group[i,0]/fzero,
group[i,1], group[i,1]/group[refi,1], decibel(group[i,1], group[refi,1])))
print('')
# erase group from good_freqs:
good_freqs = np.delete(good_freqs, indices, axis=0)
# adjust frequencies to fzero:
group[:,0] = np.round(group[:,0]/fzero)*fzero
# insert missing fzero:
if np.round(group[0,0]/fzero) != 1.0:
group = np.vstack(((fzero, group[0,1], 2.0), group))
return good_freqs, group
def build_harmonic_group(good_freqs, all_freqs, freq_tol, max_freq_tol,
verbose=0, min_group_size=3, max_divisor=4):
"""Find all the harmonics belonging to the largest peak in a list of frequency peaks.
Parameters

good_freqs: 2D array
Frequency, power, and use count (columns) of strong peaks detected
in a power spectrum.
all_freqs: 2D array
Frequency, power, and use count (columns) of all peaks detected
in a power spectrum.
freq_tol: float
Harmonics should fall within this frequency tolerance.
This should be in the range of the frequency resolution
and should not be smaller than half of the frequency resolution.
max_freq_tol: float
Maximum deviation of harmonics from their expected frequency.
Peaks with frequencies between `freq_tol` and `max_freq_tol`
get penalized the further away from their expected frequency.
verbose: int
Verbosity level.
min_group_size: int
Minimum number of harmonics of a harmonic group.
The harmonics from min_group_size/3 to max(min_group_size, divisor)
need to be in good_freqs.
max_divisor: int
Maximum divisor used for checking for subharmonics.
Returns

good_freqs: 2D array
Frequency, power, and use count (columns) of strong peaks detected
in a power spectrum with frequencies of the returned harmonic group removed.
group: 2D array
The detected harmonic group. Might be empty.
indices: 1D array of indices
Indices of the harmonic group in all_freqs.
best_fzero_harmonics: int
The highest harmonics that was used to recompute
the fundamental frequency.
fmax: float
The frequency of the largest peak in good_freqs
for which the harmonic group was detected.
"""
# select strongest frequency for building the harmonic group:
fmaxinx = np.argmax(good_freqs[:,1])
fmax = good_freqs[fmaxinx,0]
if verbose > 0:
print('')
print('%s build_harmonic_group %s' % (10*'#', 38*'#'))
print('%s fmax=%7.2fHz, power=%9.3g %s'
% (10*'#', fmax, good_freqs[fmaxinx,1], 27*'#'))
print('good_freqs: ', '[',
', '.join(['%.2f' % f for f in good_freqs[:,0]]), ']')
# container for harmonic groups:
best_group = []
best_value = 1e6
best_divisor = 0
best_fzero = 0.0
best_fzero_harmonics = 0
# check for integer fractions of the frequency:
for divisor in range(1, max_divisor + 1):
# 1. hypothesized fundamental:
freq = fmax / divisor
# 2. find harmonics in good_freqs and adjust fzero accordingly:
group_size = min_group_size if divisor <= min_group_size else divisor
new_group, fzero, fzero_harmonics = group_candidate(good_freqs, all_freqs,
freq, divisor,
freq_tol, max_freq_tol,
group_size, verbose)
# no group found:
if fzero < 0.0:
continue
# 3. compare new group to best group:
peaksum = decibel(np.sum(all_freqs[new_group, 1])*min_group_size/len(new_group))
diff = np.std(np.diff(decibel(all_freqs[new_group, 1])))
new_group_value = peaksum  diff
counts = np.sum(all_freqs[new_group, 2])
if verbose > 0:
print(' new group: fzero=%7.2fHz, nharmonics=%d, value=%6.1fdB, peaksum=%5.1fdB, diff=%6.1fdB, use count=%2d, peaks:'
% (fzero, len(new_group), new_group_value, peaksum, diff, counts), new_group)
if verbose > 1:
print(' best group: divisor=%d, fzero=%7.2fHz, nharmonics=%d, value=%6.1fdB, peaks:'
% (best_divisor, best_fzero, len(best_group), best_value), best_group)
# select new group if sum of peak power minus diff is larger:
if len(new_group) >= len(best_group) and new_group_value >= best_value:
best_value = new_group_value
best_group = new_group
best_divisor = divisor
best_fzero = fzero
best_fzero_harmonics = fzero_harmonics
if verbose > 1:
print(' new best group: divisor=%d, fzero=%7.2fHz, value=%6.1fdB, peaks:'
% (best_divisor, best_fzero, best_value), best_group)
elif verbose > 0:
print(' took as new best group')
# no group found:
if len(best_group) == 0:
# erase freq:
good_freqs = np.delete(good_freqs, fmaxinx, axis=0)
group = np.zeros((0, 3))
return good_freqs, group, best_group, 1, fmax
# update frequencies and group:
if verbose > 1:
print('')
print('# best group found for fmax=%.2fHz, fzero=%.2fHz, divisor=%d:'
% (fmax, best_fzero, best_divisor))
group_str = '%s resulting harmonic group for fmax=%.2fHz' % (10*'#', fmax)
good_freqs, group = update_group(good_freqs, all_freqs, best_group, best_fzero, freq_tol, verbose, group_str)
# good_freqs: removed all frequencies of bestgroup
return good_freqs, group, best_group, best_fzero_harmonics, fmax
def retrieve_harmonic_group(freq, good_freqs, all_freqs,
freq_tol, max_freq_tol, verbose=0,
min_group_size=3):
"""Find all the harmonics belonging to a given fundamental.
Parameters

freq: float
Fundamental frequency for which harmonics are to be retrieved.
good_freqs: 2D array
Frequency, power, and use count (columns) of strong peaks detected
in a power spectrum. All harmonics of `freq` will be
removed from `good_freqs`.
all_freqs: 2D array
Frequency, power, and use count (columns) of all peaks detected
in a power spectrum.
freq_tol: float
Harmonics should fall within this frequency tolerance.
This should be in the range of the frequency resolution
and should not be smaller than half of the frequency resolution.
max_freq_tol: float
Maximum deviation of harmonics from their expected frequency.
Peaks with frequencies between `freq_tol` and `max_freq_tol`
get penalized the further away from their expected frequency.
verbose: int
Verbosity level.
min_group_size: int
Minimum number of harmonics of a harmonic group.
The harmonics from min_group_size/3 to max(min_group_size, divisor)
need to be in good_freqs.
Returns

good_freqs: 2D array
Frequency, power, and use count (columns) of strong peaks detected
in a power spectrum with frequencies of the returned harmonic group removed.
group: 2D array
The detected harmonic group. Might be empty.
indices: 1D array of indices
Indices of the harmonic group in all_freqs.
fzero_harmonics: int
The highest harmonics that was used to recompute
the fundamental frequency.
"""
if verbose > 0:
print('')
print('%s retrieve harmonic group %s' % (10*'#', 35*'#'))
print('%s freq=%7.2fHz %s' % (10*'#', freq, 44*'#'))
print('good_freqs: ', '[',
', '.join(['%.2f' % f for f in good_freqs[:,0]]), ']')
# find harmonics in good_freqs and adjust fzero accordingly:
new_group, fzero, fzero_harmonics = group_candidate(good_freqs, all_freqs,
freq, 0,
freq_tol, max_freq_tol,
min_group_size, verbose)
# no group found:
if fzero < 0.0:
return good_freqs, np.zeros((0, 2)), np.zeros(0), fzero_harmonics
if verbose > 1:
print('')
print('# group found for freq=%.2fHz, fzero=%.2fHz:'
% (freq, fzero))
group_str = '#### resulting harmonic group for freq=%.2fHz' % freq
good_freqs, group = update_group(good_freqs, all_freqs, new_group,
fzero, freq_tol, verbose, group_str)
# good_freqs: removed all frequencies of bestgroup
return good_freqs, group, new_group, fzero_harmonics
def expand_group(group, indices, freqs, freq_tol, max_harmonics=0):
"""Add more harmonics to harmonic group.
Parameters

group: 2D array
Group of fundamental frequency and harmonics
as returned by build_harmonic_group.
indices: 1D array of indices
Indices of the harmonics in group in `freqs`.
freqs: 2D array
Frequency, power, and use count (columns) of all peaks detected
in a power spectrum.
freq_tol: float
Harmonics need to fall within this frequency tolerance.
This should be in the range of the frequency resolution
and not be smaller than half of the frequency resolution.
max_harmonics: int
Maximum number of harmonics to be returned for each group.
Returns

group: 2D array
Expanded group of fundamental frequency and harmonics.
indices: 1D array of indices
Indices of the harmonics in the expanded group in all_freqs.
"""
if len(group) == 0:
return group, indices
fzero = group[0,0]
if max_harmonics <= 0:
max_harmonics = m.floor(freqs[1,0]/fzero + 0.5) # round
if max_harmonics <= len(group):
return group, indices
group_freqs = list(group[:,0])
indices = list(indices)
last_h = m.floor(group_freqs[1]/fzero + 0.5) # round
for h in range(last_h+1, max_harmonics+1):
i = np.argmin(np.abs(freqs[:,0]/h  fzero))
f = freqs[i,0]
if m.fabs(f/h  fzero) > freq_tol:
continue
df = f  group_freqs[1]
dh = m.floor(df/fzero + 0.5) # round
if dh < 1:
continue
fe = m.fabs(df/dh  fzero)
if fe > 2*freq_tol:
continue
group_freqs.append(f)
indices.append(i)
# assemble group:
new_group = freqs[indices,:group.shape[1]]
# keep filled in fundamental:
if group[0,2] == 2:
new_group = np.vstack((group[0,:], new_group))
return new_group, np.array(indices, dtype=np.int)
def extract_fundamentals(good_freqs, all_freqs, freq_tol, max_freq_tol,
verbose=0, check_freqs=[],
mains_freq=60.0, mains_freq_tol=1.0,
min_freq=0.0, max_freq=2000.0, max_db_diff=20.0,
max_divisor=4, min_group_size=3, max_harmonics_db=5.0,
max_harmonics=0, max_groups=0, **kwargs):
"""Extract fundamental frequencies from powerspectrum peaks.
Parameters

good_freqs: 2D array
Frequency, power, and use count (columns) of strong peaks detected
in a power spectrum.
all_freqs: 2D array
Frequency, power, and use count (columns) of all peaks detected
in a power spectrum.
freq_tol: float
Harmonics should fall within this frequency tolerance.
This should be in the range of the frequency resolution
and should not be smaller than half of the frequency resolution.
max_freq_tol: float
Maximum deviation of harmonics from their expected frequency.
Peaks with frequencies between `freq_tol` and `max_freq_tol`
get penalized the further away from their expected frequency.
verbose: int
Verbosity level.
check_freqs: list of float
List of fundamental frequencies that will be checked
first for being present and valid harmonic groups in the peak frequencies
of a power spectrum.
mains_freq: float
Frequency of the mains power supply.
mains_freq_tol: float
Tolarerance around harmonics of the mains frequency,
within which peaks are removed.
min_freq: float
Minimum frequency accepted as a fundamental frequency.
max_freq: float
Maximum frequency accepted as a fundamental frequency.
max_db_diff: float
If larger than zero, maximum standard deviation of differences between
logarithmic powers of harmonics in decibel.
Low values enforce smoother power spectra.
max_divisor: int
Maximum divisor used for checking for subharmonics.
min_group_size: int
Minimum number of harmonics of a harmonic group.
The harmonics from min_group_size/3 to max(min_group_size, divisor)
need to be in good_freqs.
max_harmonics_db: float
Maximum allowed power of the `min_group_size`th and higher harmonics
after the peak (in decibel relative to peak power withn the first
`min_group_size` harmonics, i.e. if harmonics are required to be
smaller than fundamental then this is a negative number).
Make it a large positive number to effectively not check for relative power.
max_harmonics: int
Maximum number of harmonics to be returned for each group.
max_groups: int
If not zero the maximum number of most powerful harmonic groups.
Returns

group_list: list of 2D arrays
List of all harmonic groups found sorted by fundamental frequency.
Each harmonic group is a 2D array with the first dimension the harmonics
and the second dimension containing frequency and power of each harmonic.
If the power is zero, there was no corresponding peak in the power spectrum.
fzero_harmonics_list: list of int
The harmonics from which the fundamental frequencies were computed.
mains_freqs: 2d array
Array of mains peaks found in all_freqs (frequency, power).
"""
if verbose > 0:
print('')
# set use count to zero:
all_freqs[:,2] = 0.0
# remove power line harmonics from good_freqs:
if mains_freq > 0.0:
indices = np.where(np.abs(good_freqs[:,0]  np.round(good_freqs[:,0]/mains_freq)*mains_freq) < mains_freq_tol)[0]
if len(indices)>0:
if verbose > 1:
print('remove power line frequencies',
', '.join(['%.1f' % f for f in good_freqs[indices,0]]))
good_freqs = np.delete(good_freqs, indices, axis=0)
if verbose > 1:
print('all_freqs: ', '[',
', '.join(['%.2f' % f for f in all_freqs[:,0] if f < 3000.0]), ']')
group_list = []
fzero_harmonics_list = []
first = True
# as long as there are frequencies left in good_freqs:
fi = 0
while len(good_freqs) > 0:
if fi < len(check_freqs):
# check for harmonic group of a given fundamental frequency:
fmax = check_freqs[fi]
f0s = 'freq'
good_freqs, harm_group, harm_indices, fzero_harmonics = \
retrieve_harmonic_group(fmax, good_freqs, all_freqs,
freq_tol, max_freq_tol,
verbose1, min_group_size)
fi += 1
else:
# check for harmonic groups:
f0s = 'fmax'
good_freqs, harm_group, harm_indices, fzero_harmonics, fmax = \
build_harmonic_group(good_freqs, all_freqs,
freq_tol, max_freq_tol, verbose1,
min_group_size, max_divisor)
# nothing found:
if len(harm_group) == 0:
if verbose > 1  first:
s = ' (largest peak)' if first else ''
print('No harmonic group for %7.2fHz%s' % (fmax, s))
first = False
continue
first = False
# fill up harmonic group:
harm_group, harm_indices = expand_group(harm_group, harm_indices,
all_freqs, freq_tol,
max_harmonics)
# check whether fundamental was filled in:
first_h = 0 if harm_group[0,2] > 2 else 1
# increment use count:
all_freqs[harm_indices,2] += 1
if harm_group[0,2] == 1:
harm_group[0,2] = 2
# check frequency range of fundamental:
fundamental_ok = (harm_group[0, 0] >= min_freq and
harm_group[0, 0] <= max_freq)
# check power hum:
mains_ok = ((mains_freq <= 0.0) or
(m.fabs(harm_group[0,0]  mains_freq) > freq_tol))
# check smoothness:
db_powers = decibel(harm_group[first_h:,1])
diff = np.std(np.diff(db_powers))
smooth_ok = max_db_diff <= 0.0 or diff < max_db_diff
# check relative power of higher harmonics:
p_max = np.argmax(db_powers[:min_group_size])
db_powers = db_powers[p_max]
amplitude_ok = len(db_powers[p_max+min_group_size:]) == 0
if amplitude_ok:
pi = len(db_powers)  1
else:
amplitude_ok = np.all((db_powers[p_max+min_group_size:] < max_harmonics_db) 
(harm_group[first_h+p_max+min_group_size:,2] > 1))
if amplitude_ok:
pi = p_max+min_group_size1
else:
pi = np.where((db_powers[p_max+min_group_size:] >= max_harmonics_db) &
(harm_group[first_h+p_max+min_group_size:,2] <= 1))[0][0]
# check:
if fundamental_ok and mains_ok and smooth_ok and amplitude_ok:
if verbose > 0:
print('Accepting harmonic group from %s=%7.2fHz: %7.2fHz power=%9.3g nharmonics=%2d, use count=%d'
% (f0s, fmax, harm_group[0,0], np.sum(harm_group[first_h:,1]),
len(harm_group), np.sum(harm_group[first_h:,2])))
group_list.append(harm_group[:,0:2])
fzero_harmonics_list.append(fzero_harmonics)
else:
if verbose > 0:
fs = 'is ' if fundamental_ok else 'NOT'
ms = 'not ' if mains_ok else 'IS'
ss = 'smaller' if smooth_ok else 'LARGER '
ps = 'smaller' if amplitude_ok else 'LARGER '
print('Discarded harmonic group from %s=%7.2fHz: %7.2fHz power=%9.3g: %s in frequency range, %s mains frequency, smooth=%4.1fdB %s than %4.1fdB, relpower[%d]=%5.1fdB %s than %5.1fdB'
% (f0s, fmax, harm_group[0,0], np.sum(harm_group[first_h:,1]),
fs, ms, diff, ss, max_db_diff,
pi, db_powers[pi], ps, max_harmonics_db))
# select most powerful harmonic groups:
if max_groups > 0 and len(group_list) > max_groups:
n = len(group_list)
powers = [np.sum(group[:,1]) for group in group_list]
powers_inx = np.argsort(powers)
group_list = [group_list[pi] for pi in powers_inx[max_groups:]]
fzero_harmonics_list = [fzero_harmonics_list[pi] for pi in powers_inx[max_groups:]]
if verbose > 0:
print('Selected from %d groups the %d most powerful groups.' % (n, max_groups))
# sort groups by fundamental frequency:
freqs = [group[0,0] for group in group_list]
freq_inx = np.argsort(freqs)
group_list = [group_list[fi] for fi in freq_inx]
fzero_harmonics_list = [fzero_harmonics_list[fi] for fi in freq_inx]
if verbose > 1:
print('')
if len(group_list) > 0:
print('## FUNDAMENTALS FOUND: ##')
for group, fzero_h in zip(group_list, fzero_harmonics_list):
print('%7.2fHz: power=%9.3g fzeroh=%2d'
% (group[0,0], np.sum(group[:,1]), fzero_h))
else:
print('## NO FUNDAMENTALS FOUND ##')
# assemble mains frequencies from all_freqs:
if mains_freq > 0.0:
mains_freqs = all_freqs[np.abs(all_freqs[:,0]  np.round(all_freqs[:,0]/mains_freq)*mains_freq) < mains_freq_tol,:2]
else:
mains_freqs = np.zeros((0, 2))
return group_list, fzero_harmonics_list, mains_freqs
def threshold_estimate(psd_data, low_thresh_factor=6.0, high_thresh_factor=10.0,
nbins=100, hist_height=1.0/ np.sqrt(np.e)):
"""Estimate thresholds for peak detection from histogram of power spectrum.
The standard deviation of the noise floor without peaks is estimated from
the width of the histogram of the power spectrum at `hist_height` relative height.
The histogram is computed in the third quarter of the linearly detrended power spectrum.
Parameters

psd_data: 1D array
The power spectrum from which to estimate the thresholds.
low_thresh_factor: float
Factor by which the estimated standard deviation of the noise floor
is multiplied to set the `low_threshold`.
high_thresh_factor: float
Factor by which the estimated standard deviation of the noise floor
is multiplied to set the `high_threshold`.
nbins: int or list of floats
Number of bins or the bins for computing the histogram.
hist_height: float
Height between 0 and 1 at which the standard deviation of the histogram is estimated.
Returns

low_threshold: float
The threshold for peaks just above the noise floor.
high_threshold: float
The threshold for distinct peaks.
center: float
The baseline level of the power spectrum.
"""
n = len(psd_data)
psd_data_seg = psd_data[n//2:n*3//4]
psd_data_seg = psd_data_seg[~np.isinf(psd_data_seg)]
psd_data_seg = np.mean(psd_data_seg) + \
sig.detrend(psd_data_seg, type='linear')
noise_std, center = hist_threshold(psd_data_seg, thresh_fac=1.0, nbins=nbins)
low_threshold = noise_std * low_thresh_factor
high_threshold = noise_std * high_thresh_factor
return low_threshold, high_threshold, center
def harmonic_groups(psd_freqs, psd, verbose=0, check_freqs=[],
low_threshold=0.0, high_threshold=0.0, thresh_bins=100,
low_thresh_factor=6.0, high_thresh_factor=10.0,
freq_tol_fac=1.0, max_freq_tol=1.0,
mains_freq=60.0, mains_freq_tol=1.0,
min_freq=0.0, max_freq=2000.0, max_db_diff=20.0, max_divisor=4,
min_group_size=3, max_harmonics_db=5.0,
max_harmonics=0, max_groups=0, **kwargs):
"""Detect peaks in power spectrum and group them according to their harmonic structure.
Parameters

psd_freqs: 1D array
Frequencies of the power spectrum.
psd: 1D array
Power spectrum (linear, not decible).
verbose: int
Verbosity level.
check_freqs: list of float
List of fundamental frequencies that will be checked first for being
present and valid harmonic groups in the peak frequencies
of the power spectrum.
low_threshold: float
The relative threshold for detecting all peaks in the decibel spectrum.
high_threshold: float
The relative threshold for detecting good peaks in the decibel spectrum.
thresh_bins: int or list of floats
Number of bins or the bins for computing the histogram from
which the standard deviation of the noise level in the `psd` is estimated.
low_thresh_factor: float
Factor by which the estimated standard deviation of the noise floor
is multiplied to set the `low_threshold`.
high_thresh_factor: float
Factor by which the estimated standard deviation of the noise floor
is multiplied to set the `high_threshold`.
freq_tol_fac: float
Harmonics should fall within `deltaf*freq_tol_fac`.
max_freq_tol: float
Maximum absolute frequency deviation of harmonics in Hertz..
mains_freq: float
Frequency of the mains power supply.
mains_freq_tol: float
Tolarerance around harmonics of the mains frequency,
within which peaks are removed.
min_freq: float
Minimum frequency accepted as a fundamental frequency.
max_freq: float
Maximum frequency accepted as a fundamental frequency.
max_db_diff: float
If larger than zero, maximum standard deviation of differences between
logarithmic powers of harmonics in decibel (larger than zero).
Low values enforce smoother power spectra.
max_divisor: int
Maximum divisor used for checking for subharmonics.
min_group_size: int
Minimum number of harmonics of a harmonic group.
The harmonics from min_group_size/3 to max(min_group_size, divisor)
need to be in good_freqs.
max_harmonics_db: float
Maximum allowed power of the `min_group_size`th and higher harmonics
after the peak (in decibel relative to peak power withn the first
`min_group_size` harmonics, i.e. if harmonics are required to be
smaller than fundamental then this is a negative number).
Make it a large positive number to effectively not check for relative power.
max_harmonics: int
Maximum number of harmonics to be returned for each group.
max_groups: int
If not zero the maximum number of most powerful harmonic groups.
Returns

groups: list of 2D arrays
List of all extracted harmonic groups, sorted by fundamental frequency.
Each harmonic group is a 2D array with the first dimension the harmonics
and the second dimension containing frequency and power of each harmonic.
If the power is zero, there was no corresponding peak in the power spectrum.
fzero_harmonics: list of ints
The harmonics from which the fundamental frequencies were computed.
mains: 2d array
Frequencies and power of multiples of the mains frequency found in the power spectrum.
all_freqs: 2D array
Frequency, power, and use count (columns) of all peaks detected
in the power spectrum.
good_freqs: 1D array
Frequencies of strong peaks detected in the power spectrum.
low_threshold: float
The relative threshold for detecting all peaks in the decibel spectrum.
high_threshold: float
The relative threshold for detecting good peaks in the decibel spectrum.
center: float
The baseline level of the power spectrum.
"""
if verbose > 0:
print('')
if verbose > 1:
print(70*'#')
print('##### harmonic_groups', 48*'#')
# decibel power spectrum:
log_psd = decibel(psd)
max_idx = np.argmax(~np.isfinite(log_psd))
if max_idx > 0:
log_psd = log_psd[:max_idx]
delta_f = psd_freqs[1]  psd_freqs[0]
# thresholds:
center = np.NaN
if low_threshold <= 0.0 or high_threshold <= 0.0:
low_th, high_th, center = threshold_estimate(log_psd, low_thresh_factor,
high_thresh_factor,
thresh_bins)
if low_threshold <= 0.0:
low_threshold = low_th
if high_threshold <= 0.0:
high_threshold = high_th
if verbose > 1:
print('')
print('low_threshold =%4.1fdB, center+low_threshold =%6.1fdB' % (low_threshold, center + low_threshold))
print('high_threshold=%4.1fdB, center+high_threshold=%6.1fdB' % (high_threshold, center + high_threshold))
print(' center=%6.1fdB' % center)
# detect peaks in decibel power spectrum:
peaks, troughs = detect_peaks(log_psd, low_threshold)
peaks, troughs = trim(peaks, troughs)
all_freqs = np.zeros((len(peaks), 3))
all_freqs[:,0] = psd_freqs[peaks]
all_freqs[:,1] = psd[peaks]
if len(all_freqs) == 0:
return [], [], [], np.zeros((0, 3)), [], low_threshold, high_threshold, center
# select good peaks:
good_freqs = all_freqs[(log_psd[peaks]  log_psd[troughs] > high_threshold) &
(all_freqs[:,0] >= min_freq) &
(all_freqs[:,0] < max_freq*max_divisor),:]
# detect harmonic groups:
freq_tol = delta_f*freq_tol_fac
if max_freq_tol < 1.1*freq_tol:
max_freq_tol = 1.1*freq_tol
groups, fzero_harmonics, mains = \
extract_fundamentals(good_freqs, all_freqs,
freq_tol, max_freq_tol,
verbose, check_freqs, mains_freq, mains_freq_tol,
min_freq, max_freq, max_db_diff, max_divisor, min_group_size,
max_harmonics_db, max_harmonics, max_groups)
return (groups, fzero_harmonics, mains, all_freqs, good_freqs[:,0],
low_threshold, high_threshold, center)
def fundamental_freqs(group_list):
"""
Extract fundamental frequencies from lists of harmonic groups.
The inner list of 2D arrays of the input argument is transformed into
a 1D array containig the fundamental frequencies extracted from
the 2D arrays.
Parameters

group_list: (list of (list of ...)) list of 2D arrays
Arbitrarily nested lists of harmonic groups as returned by
extract_fundamentals() and harmonic_groups() with the element
[0, 0] being the fundamental frequency.
Returns

fundamentals: (list of (list of ...)) 1D array
Nested list (corresponding to `group_list`) of 1D arrays
with the fundamental frequencies extracted from the harmonic groups.
"""
if len(group_list) == 0:
return np.array([])
# check whether group_list is list of harmonic groups:
list_of_groups = True
for group in group_list:
if not ( hasattr(group, 'shape') and len(group.shape) == 2 ):
list_of_groups = False
break
if list_of_groups:
fundamentals = np.array([group[0, 0] for group in group_list if len(group) > 0])
else:
fundamentals = []
for groups in group_list:
f = fundamental_freqs(groups)
fundamentals.append(f)
return fundamentals
def fundamental_freqs_and_power(group_list, power=False,
ref_power=1.0, min_power=1e20):
"""Extract fundamental frequencies and their power in dB from lists
of harmonic groups.
The inner list of 2D arrays of the input argument is transformed
into a 2D array containig for each fish (1st dimension) the
fundamental frequencies and powers (summed over all harmonics)
extracted from the 2D arrays.
Parameters

group_list: (list of (list of ...)) list of 2D arrays
Arbitrarily nested lists of harmonic groups as returned by
extract_fundamentals() and harmonic_groups() with the element
[0, 0] being the fundamental frequency and the elements [:,1] being
the powers of each harmonics.
power: boolean
If `False` convert the power into decibel using the
powerspectrum.decibel() function.
ref_power: float
Reference power for computing decibel.
If set to `None` the maximum power is used.
min_power: float
Power values smaller than `min_power` are set to `np.inf`.
Returns

fundamentals: (list of (list of ...)) 2D array
Nested list (corresponding to `group_list`) of 2D arrays
with fundamental frequencies in first column and
corresponding power in second column.
"""
if len(group_list) == 0:
return np.array([])
# check whether group_list is list of harmonic groups:
list_of_groups = True
for group in group_list:
if not ( hasattr(group, 'shape') and len(group.shape) == 2 ):
list_of_groups = False
break
if list_of_groups:
fundamentals = np.array([[group[0, 0], np.sum(group[:, 1])]
for group in group_list if len(group) > 0])
if not power:
fundamentals[:, 1] = decibel(fundamentals[:, 1],
ref_power, min_power)
else:
fundamentals = []
for groups in group_list:
f = fundamental_freqs_and_power(groups, power,
ref_power, min_power)
fundamentals.append(f)
return fundamentals
def add_relative_power(freqs):
"""Add a column with relative power.
For each element in `freqs`, its maximum power is subtracted
from all powers.
Parameters

freqs: list of 2D arrays
First column in the ndarrays is fundamental frequency and
second column the corresponding power.
Further columns are optional and kept in the returned list.
fundamental_freqs_and_power() returns such a list.
Returns

power_freqs: list of 2D arrays
Same as freqs, but with an added column containing the relative power.
"""
return [np.column_stack((f, f[:,1]  np.max(f[:,1]))) for f in freqs]
def add_power_ranks(freqs):
"""Add a column with power ranks.
Parameters

freqs: list of 2D arrays
First column in the arrays is fundamental frequency and
second column the corresponding power.
Further columns are optional and kept in the returned list.
fundamental_freqs_and_power() returns such a list.
Returns

rank_freqs: list of 2D arrays
Same as freqs, but with an added column containing the ranks.
The highest power is assinged to zero,
lower powers are assigned negative integers.
"""
rank_freqs = []
for f in freqs:
i = np.argsort(f[:,1])[::1]
ranks = np.empty_like(i)
ranks[i] = np.arange(len(i))
rank_freqs.append(np.column_stack((f, ranks)))
return rank_freqs
def similar_indices(freqs, df_thresh, nextfs=0):
"""Indices of similar frequencies.
If two frequencies from different elements in the inner lists of `freqs` are
reciprocally the closest to each other and closer than `df_thresh`,
then two indices (element, frequency) of the respective other frequency
are appended.
Parameters

freqs: (list of (list of ...)) list of 2D arrays
First column in the arrays is fundamental frequency.
df_thresh: float
Fundamental frequencies closer than this threshold are considered
equal.
nextfs: int
If zero, compare all elements in freqs with each other. Otherwise,
only compare with the `nextfs` next elements in freqs.
Returns

indices: (list of (list of ...)) list of list of twotuples of int
For each frequency of each element in `freqs` a list of two tuples containing
the indices of elements and frequencies that are similar.
"""
if len(freqs) == 0:
return []
# check whether freqs is list of fundamental frequencies and powers:
list_of_freq_power = True
for group in freqs:
if not (hasattr(group, 'shape') and len(group.shape) == 2):
list_of_freq_power = False
break
if list_of_freq_power:
indices = [ [[] for j in range(len(freqs[i]))] for i in range(len(freqs))]
for j in range(len(freqs)1):
freqsj = np.asarray(freqs[j])
for m in range(len(freqsj)):
freq1 = freqsj[m]
nn = len(freqs) if nextfs == 0 else j+1+nextfs
if nn > len(freqs):
nn = len(freqs)
for k in range(j+1, nn):
freqsk = np.asarray(freqs[k])
if len(freqsk) == 0:
continue
n = np.argmin(np.abs(freqsk[:,0]  freq1[0]))
freq2 = freqsk[n]
if np.argmin(np.abs(freqsj[:,0]  freq2[0])) != m:
continue
if np.abs(freq1[0]  freq2[0]) < df_thresh:
indices[k][n].append((j, m))
indices[j][m].append((k, n))
else:
indices = []
for groups in freqs:
indices.append(similar_indices(groups, df_thresh, nextfs))
return indices
def unique_mask(freqs, df_thresh, nextfs=0):
"""Mark similar frequencies from different recordings as dublicate.
If two frequencies from different elements in `freqs` are
reciprocally the closest to each other and closer than `df_thresh`,
then the one with the smaller power is marked for removal.
Parameters

freqs: list of 2D arrays
First column in the arrays is fundamental frequency and
second column the corresponding power or equivalent.
If values in the second column are equal (e.g. they are the same ranks),
and there is a third column (e.g. power),
the third column is used to decide, which element should be removed.
df_thresh: float
Fundamental frequencies closer than this threshold are considered
equal.
nextfs: int
If zero, compare all elements in freqs with each other. Otherwise,
only compare with the `nextfs` next elements in freqs.
Returns

mask: list of boolean arrays
For each element in `freqs` True if that frequency should be kept.
"""
mask = [np.ones(len(freqs[i]), dtype=bool) for i in range(len(freqs))]
for j in range(len(freqs)1):
freqsj = np.asarray(freqs[j])
for m in range(len(freqsj)):
freq1 = freqsj[m]
nn = len(freqs) if nextfs == 0 else j+1+nextfs
if nn > len(freqs):
nn = len(freqs)
for k in range(j+1, nn):
freqsk = np.asarray(freqs[k])
if len(freqsk) == 0:
continue
n = np.argmin(np.abs(freqsk[:,0]  freq1[0]))
freq2 = freqsk[n]
if np.argmin(np.abs(freqsj[:,0]  freq2[0])) != m:
continue
if np.abs(freq1[0]  freq2[0]) < df_thresh:
if freq1[1] > freq2[1]:
mask[k][n] = False
elif freq1[1] < freq2[1]:
mask[j][m] = False
elif len(freq1) > 2:
if freq1[2] > freq2[2]:
mask[k][n] = False
else:
mask[j][m] = False
else:
mask[j][m] = False
return mask
def unique(freqs, df_thresh, mode='power', nextfs=0):
"""Remove similar frequencies from different recordings.
If two frequencies from different elements in the inner lists of `freqs`
are reciprocally the closest to each other and closer than `df_thresh`,
then the one with the smaller power is removed. As power, either the
absolute power as provided in the second column of the data elements
in `freqs` is taken (mode=='power'), or the relative power
(mode='relpower'), or the power rank (mode='rank').
Parameters

freqs: (list of (list of ...)) list of 2D arrays
First column in the arrays is fundamental frequency and
second column the corresponding power, as returned by
fundamental_freqs_and_power().
df_thresh: float
Fundamental frequencies closer than this threshold are considered
equal.
mode: string
 'power': use second column of freqs elements as power.
 'relpower': use relative power computed from the second column
of freqs elements for deciding which frequency to delete.
 'rank': use rank of second column of freqs elements
for deciding which frequency to delete.
nextfs: int
If zero, compare all elements in freqs with each other. Otherwise,
only compare with the `nextfs` next elements in freqs.
Returns

uniqe_freqs: (list of (list of ...)) list of 2D arrays
Same as `freqs` but elements with similar fundamental frequencies
removed.
"""
if len(freqs) == 0:
return []
# check whether freqs is list of fundamental frequencies and powers:
list_of_freq_power = True
for group in freqs:
if not (hasattr(group, 'shape') and len(group.shape) == 2):
list_of_freq_power = False
break
if list_of_freq_power:
if mode == 'power':
mask = unique_mask(freqs, df_thresh, nextfs)
elif mode == 'relpower':
power_freqs = [f[:,[0, 2, 1]] for f in add_relative_power(freqs)]
mask = unique_mask(power_freqs, df_thresh, nextfs)
elif mode == 'rank':
rank_freqs = [f[:,[0, 2, 1]] for f in add_power_ranks(freqs)]
mask = unique_mask(rank_freqs, df_thresh, nextfs)
else:
raise ValueError('%s is not a valid mode for unique(). Choose one of "power" or "rank"')
unique_freqs = []
for f, m in zip(freqs, mask):
unique_freqs.append(f[m])
else:
unique_freqs = []
for groups in freqs:
unique_freqs.append(unique(groups, df_thresh, mode, nextfs))
return unique_freqs
def colors_markers():
"""Generate a list of colors and markers for plotting.
Returns

colors: list
list of colors
markers: list
list of markers
"""
# color and marker range:
colors = []
markers = []
mr2 = []
# first color range:
cc0 = cm.gist_rainbow(np.linspace(0.0, 1.0, 8))
# shuffle it:
for k in range((len(cc0) + 1) // 2):
colors.extend(cc0[k::(len(cc0) + 1) // 2])
markers.extend(len(cc0) * 'o')
mr2.extend(len(cc0) * 'v')
# second darker color range:
cc1 = cm.gist_rainbow(np.linspace(0.33 / 7.0, 1.0, 7))
cc1 = mc.hsv_to_rgb(mc.rgb_to_hsv(np.array([cc1[:, :3]])) * np.array([1.0, 0.9, 0.7]))[0]
cc1 = np.hstack((cc1, np.ones((len(cc1),1))))
# shuffle it:
for k in range((len(cc1) + 1) // 2):
colors.extend(cc1[k::(len(cc1) + 1) // 2])
markers.extend(len(cc1) * '^')
mr2.extend(len(cc1) * '*')
# third lighter color range:
cc2 = cm.gist_rainbow(np.linspace(0.67 / 6.0, 1.0, 6))
cc2 = mc.hsv_to_rgb(mc.rgb_to_hsv(np.array([cc1[:, :3]])) * np.array([1.0, 0.5, 1.0]))[0]
cc2 = np.hstack((cc2, np.ones((len(cc2),1))))
# shuffle it:
for k in range((len(cc2) + 1) // 2):
colors.extend(cc2[k::(len(cc2) + 1) // 2])
markers.extend(len(cc2) * 'D')
mr2.extend(len(cc2) * 'x')
markers.extend(mr2)
return colors, markers
def plot_harmonic_groups(ax, group_list, indices=None, max_groups=0,
skip_bad=False, sort_by_freq=True, label_power=False,
colors=None, markers=None, legend_rows=8, **kwargs):
"""Mark decibel power of fundamentals and their harmonics in a plot.
Parameters

ax: axis for plot
Axis used for plotting.
group_list: list of 2D arrays
Lists of harmonic groups as returned by extract_fundamentals() and
harmonic_groups() with the element [0, 0] of the harmonic groups being
the fundamental frequency, and element[0, 1] being the corresponding power.
indices: list of int or None
If smaller than zero then set the legend label of the corresponding group in brackets.
max_groups: int
If not zero plot only the max_groups most powerful groups.
skip_bad: bool
Skip harmonic groups without index (entry in indices is negative).
sort_by_freq: boolean
If True sort legend by frequency, otherwise by power.
label_power: boolean
If `True` put the power in decibel in addition to the frequency into the legend.
colors: list of colors or None
If not None list of colors for plotting each group
markers: list of markers or None
If not None list of markers for plotting each group
legend_rows: int
Maximum number of rows to be used for the legend.
kwargs:
Key word arguments for the legend of the plot.
"""
if len(group_list) == 0:
return
# sort by power:
powers = np.array([np.sum(group[:,1]) for group in group_list])
max_power = np.max(powers)
idx = np.argsort(powers)
if max_groups > 0 and len(idx > max_groups):
idx = idx[:max_groups]
# sort by frequency:
if sort_by_freq:
freqs = np.array([group_list[group][0,0] for group in idx])
if legend_rows > 0 and legend_rows < len(freqs):
idx = idx[np.argsort(freqs[:legend_rows])]
else:
idx = idx[np.argsort(freqs)]
# plot:
k = 0
for i in idx:
if indices is not None and skip_bad and indices[i] < 0:
continue
group = group_list[i]
x = group[:,0]
y = decibel(group[:,1])
msize = 7.0 + 10.0*(powers[i]/max_power)**0.25
color_kwargs = {}
if colors is not None:
color_kwargs = {'color': colors[k%len(colors)]}
label = '%6.1f Hz' % group[0, 0]
if label_power:
label += ' %6.1f dB' % decibel(np.array([np.sum(group[:,1])]))[0]
if indices is not None:
if indices[i] < 0:
label = '(' + label + ')'
else:
label = ' ' + label + ' '
if legend_rows > 5 and k >= legend_rows:
label = None
if markers is None:
ax.plot(x, y, 'o', ms=msize, label=label, **color_kwargs)
else:
if k >= len(markers):
break
ax.plot(x, y, linestyle='None', marker=markers[k],
mec=None, mew=0.0, ms=msize, label=label, **color_kwargs)
k += 1
# legend:
if legend_rows > 0:
if legend_rows > 5:
ncol = 1
else:
ncol = (len(idx)1) // legend_rows + 1
ax.legend(numpoints=1, ncol=ncol, **kwargs)
else:
ax.legend(numpoints=1, **kwargs)
def plot_psd_harmonic_groups(ax, psd_freqs, psd, group_list,
mains=None, all_freqs=None, good_freqs=None,
log_freq=False, min_freq=0.0, max_freq=2000.0, ymarg=0.0):
"""Plot decibel powerspectrum with detected peaks, harmonic groups,
and mains frequencies.
Parameters

psd_freqs: array
Frequencies of the power spectrum.
psd: array
Power spectrum (linear, not decible).
group_list: list of 2D arrays
Lists of harmonic groups as returned by extract_fundamentals() and
harmonic_groups() with the element [0, 0] of the harmonic groups being the fundamental frequency,
and element[0, 1] being the corresponding power.
mains: 2D array
Frequencies and power of multiples of the mains frequency found in the power spectrum.
all_freqs: 2D array
Peaks in the power spectrum detected with low threshold.
good_freqs: 1D array
Frequencies of peaks detected with high threshold.
log_freq: boolean
Logarithmic (True) or linear (False) frequency axis.
min_freq: float
Limits of frequency axis are set to `(min_freq, max_freq)`
if `max_freq` is greater than zero
max_freq: float
Limits of frequency axis are set to `(min_freq, max_freq)`
and limits of power axis are computed from powers below max_freq
if `max_freq` is greater than zero
ymarg: float
Add this to the maximum decibel power for setting the ylim.
"""
# mark all and good psd peaks:
pmin, pmax = ax.get_ylim()
doty = pmax  5.0
if all_freqs is not None:
ax.plot(all_freqs[:, 0], np.zeros(len(all_freqs[:, 0])) + doty, 'o', color='#ffffff')
if good_freqs is not None:
ax.plot(good_freqs, np.zeros(len(good_freqs)) + doty, 'o', color='#888888')
# mark mains frequencies:
if mains is not None and len(mains) > 0:
fpeaks = mains[:, 0]
fpeakinx = [int(np.round(fp/(psd_freqs[1]psd_freqs[0]))) for fp in fpeaks if fp < psd_freqs[1]]
ax.plot(fpeaks[:len(fpeakinx)], decibel(psd[fpeakinx]), linestyle='None',
marker='.', color='k', ms=10, mec=None, mew=0.0,
label='%3.0f Hz mains' % mains[0, 0])
# mark harmonic groups:
colors, markers = colors_markers()
plot_harmonic_groups(ax, group_list, max_groups=0, sort_by_freq=True,
colors=colors, markers=markers, legend_rows=8,
loc='upper right')
# plot power spectrum:
plot_decibel_psd(ax, psd_freqs, psd, log_freq=log_freq, min_freq=min_freq,
max_freq=max_freq, ymarg=ymarg, color='blue')
def add_psd_peak_detection_config(cfg, low_threshold=0.0, high_threshold=0.0,
thresh_bins=100,
low_thresh_factor=6.0, high_thresh_factor=10.0):
"""Add parameter needed for detection of peaks in power spectrum used
by harmonic_groups() as a new section to a configuration.
Parameters

cfg: ConfigFile
The configuration.
"""
cfg.add_section('Thresholds for peak detection in power spectra:')
cfg.add('lowThreshold', low_threshold, 'dB', 'Threshold for all peaks.\n If 0.0 estimate threshold from histogram.')
cfg.add('highThreshold', high_threshold, 'dB', 'Threshold for good peaks. If 0.0 estimate threshold from histogram.')
cfg.add_section('Threshold estimation:\nIf no thresholds are specified they are estimated from the histogram of the decibel power spectrum.')
cfg.add('thresholdBins', thresh_bins, '', 'Number of bins used to compute the histogram used for threshold estimation.')
cfg.add('lowThresholdFactor', low_thresh_factor, '', 'Factor for multiplying standard deviation of noise floor for lower threshold.')
cfg.add('highThresholdFactor', high_thresh_factor, '', 'Factor for multiplying standard deviation of noise floor for higher threshold.')
def psd_peak_detection_args(cfg):
"""Translates a configuration to the respective parameter names for
the detection of peaks in power spectrum used by
harmonic_groups(). The return value can then be passed as
keyword arguments to this function.
Parameters

cfg: ConfigFile
The configuration.
Returns

a: dict
Dictionary with names of arguments of the `harmonicgroup()` function
and their values as supplied by `cfg`.
"""
return cfg.map({'low_threshold': 'lowThreshold',
'high_threshold': 'highThreshold',
'thresh_bins': 'thresholdBins',
'low_thresh_factor': 'lowThresholdFactor',
'high_thresh_factor': 'highThresholdFactor'})
def add_harmonic_groups_config(cfg, mains_freq=60.0, mains_freq_tol=1.0,
max_divisor=4, freq_tol_fac=1.0,
max_freq_tol=1.0, min_group_size=3,
min_freq=20.0, max_freq=2000.0, max_db_diff=20.0,
max_harmonics_db=5.0, max_harmonics=0, max_groups=0):
"""Add parameter needed for detection of harmonic groups as a new
section to a configuration.
Parameters

cfg: ConfigFile
The configuration.
"""
cfg.add_section('Harmonic groups:')
cfg.add('mainsFreq', mains_freq, 'Hz', 'Mains frequency to be excluded.')
cfg.add('mainsFreqTolerance', mains_freq_tol, 'Hz', 'Exclude peaks within this tolerance around multiples of the mains frequency.')
cfg.add('minimumGroupSize', min_group_size, '',
'Minimum number of harmonics (inclusively fundamental) that make up a harmonic group.')
cfg.add('maxDivisor', max_divisor, '', 'Maximum ratio between the frequency of the largest peak and its fundamental.')
cfg.add('freqTolerance', freq_tol_fac, '',
'Harmonics should be within this factor times the frequency resolution of the power spectrum. Needs to be higher than 0.5!')
cfg.add('maximumFreqTolerance', max_freq_tol, 'Hz',
'Maximum deviation of harmonics from their expected value.')
cfg.add_section('Acceptance of best harmonic groups:')
cfg.add('minimumFrequency', min_freq, 'Hz', 'Minimum frequency allowed for the fundamental.')
cfg.add('maximumFrequency', max_freq, 'Hz', 'Maximum frequency allowed for the fundamental.')
cfg.add('maximumPowerDifference', max_db_diff, 'dB', 'If larger than zero, maximum standard deviation allowed for difference in logarithmic power between successive harmonics. Smaller values enforce smoother spectra.')
cfg.add('maximumHarmonicsPower', max_harmonics_db, 'dB', 'Maximum allowed power of the minimumGroupSizeth and higher harmonics relative to peak power.')
cfg.add('maximumHarmonics', max_harmonics, '', '0: keep all, >0 only keep the first # harmonics.')
cfg.add('maximumGroups', max_groups, '', 'Maximum number of harmonic groups. If 0 process all.')
def harmonic_groups_args(cfg):
"""Translates a configuration to the respective parameter names of
the harmonicgroup detection functions. The return value can then
be passed as keyword arguments to this function.
Parameters

cfg: ConfigFile
The configuration.
Returns

a: dict
Dictionary with names of arguments of the harmonicgroup detection functions
and their values as supplied by `cfg`.
"""
return cfg.map({'mains_freq': 'mainsFreq',
'mains_freq_tol': 'mainsFreqTolerance',
'freq_tol_fac': 'freqTolerance',
'max_freq_tol': 'maximumFreqTolerance',
'max_divisor': 'maxDivisor',
'min_group_size': 'minimumGroupSize',
'min_freq': 'minimumFrequency',
'max_freq': 'maximumFrequency',
'max_db_diff': 'maximumPowerDifference',
'max_harmonics_db': 'maximumHarmonicsPower',
'max_harmonics': 'maximumHarmonics',
'max_groups': 'maximumGroups'})
def main(data_file=None):
import matplotlib.pyplot as plt
from .fakefish import wavefish_eods
from .powerspectrum import psd
if data_file is None:
# generate data:
title = 'simulation'
samplerate = 44100.0
d = 20.0
noise = 0.01
eodfs = [123.0, 333.0, 666.0, 666.5]
fish1 = 0.5*wavefish_eods('Eigenmannia', eodfs[0], samplerate,
duration=d, noise_std=noise)
fish2 = 1.0*wavefish_eods('Eigenmannia', eodfs[1], samplerate,
duration=d, noise_std=noise)
fish3 = 10.0*wavefish_eods('Alepto', eodfs[2], samplerate,
duration=d, noise_std=noise)
fish4 = 6.0*wavefish_eods('Alepto', eodfs[3], samplerate,
duration=d, noise_std=noise)
data = fish1 + fish2 + fish3 + fish4
else:
from .dataloader import load_data
print("load %s ..." % data_file)
data, samplerate, unit = load_data(data_file, 0)
title = data_file
# retrieve fundamentals from power spectrum:
psd_data = psd(data, samplerate, freq_resolution=0.1)
groups, _, mains, all_freqs, good_freqs, _, _, _ = harmonic_groups(psd_data[0], psd_data[1], check_freqs=[123.0, 666.0], max_db_diff=30.0, verbose=0)
fig, ax = plt.subplots()
plot_psd_harmonic_groups(ax, psd_data[0], psd_data[1], groups, mains,
all_freqs, good_freqs, max_freq=3000.0)
ax.set_title(title)
plt.show()
# unify fundamental frequencies:
fundamentals = fundamental_freqs(groups)
np.set_printoptions(formatter={'float': lambda x: '%5.1f' % x})
print('fundamental frequencies extracted from power spectrum:')
print(fundamentals)
print('')
freqs = fundamental_freqs_and_power([groups])
freqs.append(np.array([[44.0, 20.0], [44.2, 10.0], [320.5, 2.5], [665.5, 5.0], [666.2, 10.0]]))
freqs.append(np.array([[123.3, 1.0], [320.2, 2.0], [668.4, 2.0]]))
rank_freqs = add_relative_power(freqs)
rank_freqs = add_power_ranks(rank_freqs)
print('all frequencies (frequency, power, relpower, rank):')
print('\n'.join(( str(f) for f in rank_freqs)))
print('')
indices = similar_indices(freqs, 1.0)
print('similar indices:')
print('\n'.join(( ('\n '.join((str(f) for f in g)) for g in indices))))
print('')
unique_freqs = unique(freqs, 1.0, 'power')
print('unique power:')
print('\n'.join(( str(f) for f in unique_freqs)))
print('')
unique_freqs = unique(freqs, 1.0, 'relpower')
print('unique relative power:')
print('\n'.join(( str(f) for f in unique_freqs)))
print('')
unique_freqs = unique(freqs, 1.0, 'rank')
print('unique rank:')
print('\n'.join(( str(f) for f in unique_freqs)))
print('')
unique_freqs = unique(freqs, 1.0, 'rank', 1)
print('unique rank for next neighour only:')
print('\n'.join(( str(f) for f in unique_freqs)))
print('')
if __name__ == "__main__":
import sys
data_file = sys.argv[1] if len(sys.argv) > 1 else None
main(data_file)
Functions
def add_harmonic_groups_config(cfg, mains_freq=60.0, mains_freq_tol=1.0, max_divisor=4, freq_tol_fac=1.0, max_freq_tol=1.0, min_group_size=3, min_freq=20.0, max_freq=2000.0, max_db_diff=20.0, max_harmonics_db=5.0, max_harmonics=0, max_groups=0)

Add parameter needed for detection of harmonic groups as a new section to a configuration.
Parameters
cfg
:ConfigFile
 The configuration.
Expand source code
def add_harmonic_groups_config(cfg, mains_freq=60.0, mains_freq_tol=1.0, max_divisor=4, freq_tol_fac=1.0, max_freq_tol=1.0, min_group_size=3, min_freq=20.0, max_freq=2000.0, max_db_diff=20.0, max_harmonics_db=5.0, max_harmonics=0, max_groups=0): """Add parameter needed for detection of harmonic groups as a new section to a configuration. Parameters  cfg: ConfigFile The configuration. """ cfg.add_section('Harmonic groups:') cfg.add('mainsFreq', mains_freq, 'Hz', 'Mains frequency to be excluded.') cfg.add('mainsFreqTolerance', mains_freq_tol, 'Hz', 'Exclude peaks within this tolerance around multiples of the mains frequency.') cfg.add('minimumGroupSize', min_group_size, '', 'Minimum number of harmonics (inclusively fundamental) that make up a harmonic group.') cfg.add('maxDivisor', max_divisor, '', 'Maximum ratio between the frequency of the largest peak and its fundamental.') cfg.add('freqTolerance', freq_tol_fac, '', 'Harmonics should be within this factor times the frequency resolution of the power spectrum. Needs to be higher than 0.5!') cfg.add('maximumFreqTolerance', max_freq_tol, 'Hz', 'Maximum deviation of harmonics from their expected value.') cfg.add_section('Acceptance of best harmonic groups:') cfg.add('minimumFrequency', min_freq, 'Hz', 'Minimum frequency allowed for the fundamental.') cfg.add('maximumFrequency', max_freq, 'Hz', 'Maximum frequency allowed for the fundamental.') cfg.add('maximumPowerDifference', max_db_diff, 'dB', 'If larger than zero, maximum standard deviation allowed for difference in logarithmic power between successive harmonics. Smaller values enforce smoother spectra.') cfg.add('maximumHarmonicsPower', max_harmonics_db, 'dB', 'Maximum allowed power of the minimumGroupSizeth and higher harmonics relative to peak power.') cfg.add('maximumHarmonics', max_harmonics, '', '0: keep all, >0 only keep the first # harmonics.') cfg.add('maximumGroups', max_groups, '', 'Maximum number of harmonic groups. If 0 process all.')
def add_power_ranks(freqs)

Add a column with power ranks.
Parameters
freqs
:list
of2D arrays
 First column in the arrays is fundamental frequency and second column the corresponding power. Further columns are optional and kept in the returned list. fundamental_freqs_and_power() returns such a list.
Returns
rank_freqs
:list
of2D arrays
 Same as freqs, but with an added column containing the ranks. The highest power is assinged to zero, lower powers are assigned negative integers.
Expand source code
def add_power_ranks(freqs): """Add a column with power ranks. Parameters  freqs: list of 2D arrays First column in the arrays is fundamental frequency and second column the corresponding power. Further columns are optional and kept in the returned list. fundamental_freqs_and_power() returns such a list. Returns  rank_freqs: list of 2D arrays Same as freqs, but with an added column containing the ranks. The highest power is assinged to zero, lower powers are assigned negative integers. """ rank_freqs = [] for f in freqs: i = np.argsort(f[:,1])[::1] ranks = np.empty_like(i) ranks[i] = np.arange(len(i)) rank_freqs.append(np.column_stack((f, ranks))) return rank_freqs
def add_psd_peak_detection_config(cfg, low_threshold=0.0, high_threshold=0.0, thresh_bins=100, low_thresh_factor=6.0, high_thresh_factor=10.0)

Add parameter needed for detection of peaks in power spectrum used by harmonic_groups() as a new section to a configuration.
Parameters
cfg
:ConfigFile
 The configuration.
Expand source code
def add_psd_peak_detection_config(cfg, low_threshold=0.0, high_threshold=0.0, thresh_bins=100, low_thresh_factor=6.0, high_thresh_factor=10.0): """Add parameter needed for detection of peaks in power spectrum used by harmonic_groups() as a new section to a configuration. Parameters  cfg: ConfigFile The configuration. """ cfg.add_section('Thresholds for peak detection in power spectra:') cfg.add('lowThreshold', low_threshold, 'dB', 'Threshold for all peaks.\n If 0.0 estimate threshold from histogram.') cfg.add('highThreshold', high_threshold, 'dB', 'Threshold for good peaks. If 0.0 estimate threshold from histogram.') cfg.add_section('Threshold estimation:\nIf no thresholds are specified they are estimated from the histogram of the decibel power spectrum.') cfg.add('thresholdBins', thresh_bins, '', 'Number of bins used to compute the histogram used for threshold estimation.') cfg.add('lowThresholdFactor', low_thresh_factor, '', 'Factor for multiplying standard deviation of noise floor for lower threshold.') cfg.add('highThresholdFactor', high_thresh_factor, '', 'Factor for multiplying standard deviation of noise floor for higher threshold.')
def add_relative_power(freqs)

Add a column with relative power.
For each element in
freqs
, its maximum power is subtracted from all powers.Parameters
freqs
:list
of2D arrays
 First column in the ndarrays is fundamental frequency and second column the corresponding power. Further columns are optional and kept in the returned list. fundamental_freqs_and_power() returns such a list.
Returns
power_freqs
:list
of2D arrays
 Same as freqs, but with an added column containing the relative power.
Expand source code
def add_relative_power(freqs): """Add a column with relative power. For each element in `freqs`, its maximum power is subtracted from all powers. Parameters  freqs: list of 2D arrays First column in the ndarrays is fundamental frequency and second column the corresponding power. Further columns are optional and kept in the returned list. fundamental_freqs_and_power() returns such a list. Returns  power_freqs: list of 2D arrays Same as freqs, but with an added column containing the relative power. """ return [np.column_stack((f, f[:,1]  np.max(f[:,1]))) for f in freqs]
def build_harmonic_group(good_freqs, all_freqs, freq_tol, max_freq_tol, verbose=0, min_group_size=3, max_divisor=4)

Find all the harmonics belonging to the largest peak in a list of frequency peaks.
Parameters
good_freqs
:2D array
 Frequency, power, and use count (columns) of strong peaks detected in a power spectrum.
all_freqs
:2D array
 Frequency, power, and use count (columns) of all peaks detected in a power spectrum.
freq_tol
:float
 Harmonics should fall within this frequency tolerance. This should be in the range of the frequency resolution and should not be smaller than half of the frequency resolution.
max_freq_tol
:float
 Maximum deviation of harmonics from their expected frequency.
Peaks with frequencies between
freq_tol
andmax_freq_tol
get penalized the further away from their expected frequency. verbose
:int
 Verbosity level.
min_group_size
:int
 Minimum number of harmonics of a harmonic group. The harmonics from min_group_size/3 to max(min_group_size, divisor) need to be in good_freqs.
max_divisor
:int
 Maximum divisor used for checking for subharmonics.
Returns
good_freqs
:2D array
 Frequency, power, and use count (columns) of strong peaks detected in a power spectrum with frequencies of the returned harmonic group removed.
group
:2D array
 The detected harmonic group. Might be empty.
indices
:1D array
ofindices
 Indices of the harmonic group in all_freqs.
best_fzero_harmonics
:int
 The highest harmonics that was used to recompute the fundamental frequency.
fmax
:float
 The frequency of the largest peak in good_freqs for which the harmonic group was detected.
Expand source code
def build_harmonic_group(good_freqs, all_freqs, freq_tol, max_freq_tol, verbose=0, min_group_size=3, max_divisor=4): """Find all the harmonics belonging to the largest peak in a list of frequency peaks. Parameters  good_freqs: 2D array Frequency, power, and use count (columns) of strong peaks detected in a power spectrum. all_freqs: 2D array Frequency, power, and use count (columns) of all peaks detected in a power spectrum. freq_tol: float Harmonics should fall within this frequency tolerance. This should be in the range of the frequency resolution and should not be smaller than half of the frequency resolution. max_freq_tol: float Maximum deviation of harmonics from their expected frequency. Peaks with frequencies between `freq_tol` and `max_freq_tol` get penalized the further away from their expected frequency. verbose: int Verbosity level. min_group_size: int Minimum number of harmonics of a harmonic group. The harmonics from min_group_size/3 to max(min_group_size, divisor) need to be in good_freqs. max_divisor: int Maximum divisor used for checking for subharmonics. Returns  good_freqs: 2D array Frequency, power, and use count (columns) of strong peaks detected in a power spectrum with frequencies of the returned harmonic group removed. group: 2D array The detected harmonic group. Might be empty. indices: 1D array of indices Indices of the harmonic group in all_freqs. best_fzero_harmonics: int The highest harmonics that was used to recompute the fundamental frequency. fmax: float The frequency of the largest peak in good_freqs for which the harmonic group was detected. """ # select strongest frequency for building the harmonic group: fmaxinx = np.argmax(good_freqs[:,1]) fmax = good_freqs[fmaxinx,0] if verbose > 0: print('') print('%s build_harmonic_group %s' % (10*'#', 38*'#')) print('%s fmax=%7.2fHz, power=%9.3g %s' % (10*'#', fmax, good_freqs[fmaxinx,1], 27*'#')) print('good_freqs: ', '[', ', '.join(['%.2f' % f for f in good_freqs[:,0]]), ']') # container for harmonic groups: best_group = [] best_value = 1e6 best_divisor = 0 best_fzero = 0.0 best_fzero_harmonics = 0 # check for integer fractions of the frequency: for divisor in range(1, max_divisor + 1): # 1. hypothesized fundamental: freq = fmax / divisor # 2. find harmonics in good_freqs and adjust fzero accordingly: group_size = min_group_size if divisor <= min_group_size else divisor new_group, fzero, fzero_harmonics = group_candidate(good_freqs, all_freqs, freq, divisor, freq_tol, max_freq_tol, group_size, verbose) # no group found: if fzero < 0.0: continue # 3. compare new group to best group: peaksum = decibel(np.sum(all_freqs[new_group, 1])*min_group_size/len(new_group)) diff = np.std(np.diff(decibel(all_freqs[new_group, 1]))) new_group_value = peaksum  diff counts = np.sum(all_freqs[new_group, 2]) if verbose > 0: print(' new group: fzero=%7.2fHz, nharmonics=%d, value=%6.1fdB, peaksum=%5.1fdB, diff=%6.1fdB, use count=%2d, peaks:' % (fzero, len(new_group), new_group_value, peaksum, diff, counts), new_group) if verbose > 1: print(' best group: divisor=%d, fzero=%7.2fHz, nharmonics=%d, value=%6.1fdB, peaks:' % (best_divisor, best_fzero, len(best_group), best_value), best_group) # select new group if sum of peak power minus diff is larger: if len(new_group) >= len(best_group) and new_group_value >= best_value: best_value = new_group_value best_group = new_group best_divisor = divisor best_fzero = fzero best_fzero_harmonics = fzero_harmonics if verbose > 1: print(' new best group: divisor=%d, fzero=%7.2fHz, value=%6.1fdB, peaks:' % (best_divisor, best_fzero, best_value), best_group) elif verbose > 0: print(' took as new best group') # no group found: if len(best_group) == 0: # erase freq: good_freqs = np.delete(good_freqs, fmaxinx, axis=0) group = np.zeros((0, 3)) return good_freqs, group, best_group, 1, fmax # update frequencies and group: if verbose > 1: print('') print('# best group found for fmax=%.2fHz, fzero=%.2fHz, divisor=%d:' % (fmax, best_fzero, best_divisor)) group_str = '%s resulting harmonic group for fmax=%.2fHz' % (10*'#', fmax) good_freqs, group = update_group(good_freqs, all_freqs, best_group, best_fzero, freq_tol, verbose, group_str) # good_freqs: removed all frequencies of bestgroup return good_freqs, group, best_group, best_fzero_harmonics, fmax
def colors_markers()

Generate a list of colors and markers for plotting.
Returns
colors
:list
 list of colors
markers
:list
 list of markers
Expand source code
def colors_markers(): """Generate a list of colors and markers for plotting. Returns  colors: list list of colors markers: list list of markers """ # color and marker range: colors = [] markers = [] mr2 = [] # first color range: cc0 = cm.gist_rainbow(np.linspace(0.0, 1.0, 8)) # shuffle it: for k in range((len(cc0) + 1) // 2): colors.extend(cc0[k::(len(cc0) + 1) // 2]) markers.extend(len(cc0) * 'o') mr2.extend(len(cc0) * 'v') # second darker color range: cc1 = cm.gist_rainbow(np.linspace(0.33 / 7.0, 1.0, 7)) cc1 = mc.hsv_to_rgb(mc.rgb_to_hsv(np.array([cc1[:, :3]])) * np.array([1.0, 0.9, 0.7]))[0] cc1 = np.hstack((cc1, np.ones((len(cc1),1)))) # shuffle it: for k in range((len(cc1) + 1) // 2): colors.extend(cc1[k::(len(cc1) + 1) // 2]) markers.extend(len(cc1) * '^') mr2.extend(len(cc1) * '*') # third lighter color range: cc2 = cm.gist_rainbow(np.linspace(0.67 / 6.0, 1.0, 6)) cc2 = mc.hsv_to_rgb(mc.rgb_to_hsv(np.array([cc1[:, :3]])) * np.array([1.0, 0.5, 1.0]))[0] cc2 = np.hstack((cc2, np.ones((len(cc2),1)))) # shuffle it: for k in range((len(cc2) + 1) // 2): colors.extend(cc2[k::(len(cc2) + 1) // 2]) markers.extend(len(cc2) * 'D') mr2.extend(len(cc2) * 'x') markers.extend(mr2) return colors, markers
def expand_group(group, indices, freqs, freq_tol, max_harmonics=0)

Add more harmonics to harmonic group.
Parameters
group
:2D array
 Group of fundamental frequency and harmonics as returned by build_harmonic_group.
indices
:1D array
ofindices
 Indices of the harmonics in group in
freqs
. freqs
:2D array
 Frequency, power, and use count (columns) of all peaks detected in a power spectrum.
freq_tol
:float
 Harmonics need to fall within this frequency tolerance. This should be in the range of the frequency resolution and not be smaller than half of the frequency resolution.
max_harmonics
:int
 Maximum number of harmonics to be returned for each group.
Returns
group
:2D array
 Expanded group of fundamental frequency and harmonics.
indices
:1D array
ofindices
 Indices of the harmonics in the expanded group in all_freqs.
Expand source code
def expand_group(group, indices, freqs, freq_tol, max_harmonics=0): """Add more harmonics to harmonic group. Parameters  group: 2D array Group of fundamental frequency and harmonics as returned by build_harmonic_group. indices: 1D array of indices Indices of the harmonics in group in `freqs`. freqs: 2D array Frequency, power, and use count (columns) of all peaks detected in a power spectrum. freq_tol: float Harmonics need to fall within this frequency tolerance. This should be in the range of the frequency resolution and not be smaller than half of the frequency resolution. max_harmonics: int Maximum number of harmonics to be returned for each group. Returns  group: 2D array Expanded group of fundamental frequency and harmonics. indices: 1D array of indices Indices of the harmonics in the expanded group in all_freqs. """ if len(group) == 0: return group, indices fzero = group[0,0] if max_harmonics <= 0: max_harmonics = m.floor(freqs[1,0]/fzero + 0.5) # round if max_harmonics <= len(group): return group, indices group_freqs = list(group[:,0]) indices = list(indices) last_h = m.floor(group_freqs[1]/fzero + 0.5) # round for h in range(last_h+1, max_harmonics+1): i = np.argmin(np.abs(freqs[:,0]/h  fzero)) f = freqs[i,0] if m.fabs(f/h  fzero) > freq_tol: continue df = f  group_freqs[1] dh = m.floor(df/fzero + 0.5) # round if dh < 1: continue fe = m.fabs(df/dh  fzero) if fe > 2*freq_tol: continue group_freqs.append(f) indices.append(i) # assemble group: new_group = freqs[indices,:group.shape[1]] # keep filled in fundamental: if group[0,2] == 2: new_group = np.vstack((group[0,:], new_group)) return new_group, np.array(indices, dtype=np.int)
def extract_fundamentals(good_freqs, all_freqs, freq_tol, max_freq_tol, verbose=0, check_freqs=[], mains_freq=60.0, mains_freq_tol=1.0, min_freq=0.0, max_freq=2000.0, max_db_diff=20.0, max_divisor=4, min_group_size=3, max_harmonics_db=5.0, max_harmonics=0, max_groups=0, **kwargs)

Extract fundamental frequencies from powerspectrum peaks.
Parameters
good_freqs
:2D array
 Frequency, power, and use count (columns) of strong peaks detected in a power spectrum.
all_freqs
:2D array
 Frequency, power, and use count (columns) of all peaks detected in a power spectrum.
freq_tol
:float
 Harmonics should fall within this frequency tolerance. This should be in the range of the frequency resolution and should not be smaller than half of the frequency resolution.
max_freq_tol
:float
 Maximum deviation of harmonics from their expected frequency.
Peaks with frequencies between
freq_tol
andmax_freq_tol
get penalized the further away from their expected frequency. verbose
:int
 Verbosity level.
check_freqs
:list
offloat
 List of fundamental frequencies that will be checked first for being present and valid harmonic groups in the peak frequencies of a power spectrum.
mains_freq
:float
 Frequency of the mains power supply.
mains_freq_tol
:float
 Tolarerance around harmonics of the mains frequency, within which peaks are removed.
min_freq
:float
 Minimum frequency accepted as a fundamental frequency.
max_freq
:float
 Maximum frequency accepted as a fundamental frequency.
max_db_diff
:float
 If larger than zero, maximum standard deviation of differences between logarithmic powers of harmonics in decibel. Low values enforce smoother power spectra.
max_divisor
:int
 Maximum divisor used for checking for subharmonics.
min_group_size
:int
 Minimum number of harmonics of a harmonic group. The harmonics from min_group_size/3 to max(min_group_size, divisor) need to be in good_freqs.
max_harmonics_db
:float
 Maximum allowed power of the
min_group_size
th and higher harmonics after the peak (in decibel relative to peak power withn the firstmin_group_size
harmonics, i.e. if harmonics are required to be smaller than fundamental then this is a negative number). Make it a large positive number to effectively not check for relative power. max_harmonics
:int
 Maximum number of harmonics to be returned for each group.
max_groups
:int
 If not zero the maximum number of most powerful harmonic groups.
Returns
group_list
:list
of2D arrays
 List of all harmonic groups found sorted by fundamental frequency. Each harmonic group is a 2D array with the first dimension the harmonics and the second dimension containing frequency and power of each harmonic. If the power is zero, there was no corresponding peak in the power spectrum.
fzero_harmonics_list
:list
ofint
 The harmonics from which the fundamental frequencies were computed.
mains_freqs
:2d array
 Array of mains peaks found in all_freqs (frequency, power).
Expand source code
def extract_fundamentals(good_freqs, all_freqs, freq_tol, max_freq_tol, verbose=0, check_freqs=[], mains_freq=60.0, mains_freq_tol=1.0, min_freq=0.0, max_freq=2000.0, max_db_diff=20.0, max_divisor=4, min_group_size=3, max_harmonics_db=5.0, max_harmonics=0, max_groups=0, **kwargs): """Extract fundamental frequencies from powerspectrum peaks. Parameters  good_freqs: 2D array Frequency, power, and use count (columns) of strong peaks detected in a power spectrum. all_freqs: 2D array Frequency, power, and use count (columns) of all peaks detected in a power spectrum. freq_tol: float Harmonics should fall within this frequency tolerance. This should be in the range of the frequency resolution and should not be smaller than half of the frequency resolution. max_freq_tol: float Maximum deviation of harmonics from their expected frequency. Peaks with frequencies between `freq_tol` and `max_freq_tol` get penalized the further away from their expected frequency. verbose: int Verbosity level. check_freqs: list of float List of fundamental frequencies that will be checked first for being present and valid harmonic groups in the peak frequencies of a power spectrum. mains_freq: float Frequency of the mains power supply. mains_freq_tol: float Tolarerance around harmonics of the mains frequency, within which peaks are removed. min_freq: float Minimum frequency accepted as a fundamental frequency. max_freq: float Maximum frequency accepted as a fundamental frequency. max_db_diff: float If larger than zero, maximum standard deviation of differences between logarithmic powers of harmonics in decibel. Low values enforce smoother power spectra. max_divisor: int Maximum divisor used for checking for subharmonics. min_group_size: int Minimum number of harmonics of a harmonic group. The harmonics from min_group_size/3 to max(min_group_size, divisor) need to be in good_freqs. max_harmonics_db: float Maximum allowed power of the `min_group_size`th and higher harmonics after the peak (in decibel relative to peak power withn the first `min_group_size` harmonics, i.e. if harmonics are required to be smaller than fundamental then this is a negative number). Make it a large positive number to effectively not check for relative power. max_harmonics: int Maximum number of harmonics to be returned for each group. max_groups: int If not zero the maximum number of most powerful harmonic groups. Returns  group_list: list of 2D arrays List of all harmonic groups found sorted by fundamental frequency. Each harmonic group is a 2D array with the first dimension the harmonics and the second dimension containing frequency and power of each harmonic. If the power is zero, there was no corresponding peak in the power spectrum. fzero_harmonics_list: list of int The harmonics from which the fundamental frequencies were computed. mains_freqs: 2d array Array of mains peaks found in all_freqs (frequency, power). """ if verbose > 0: print('') # set use count to zero: all_freqs[:,2] = 0.0 # remove power line harmonics from good_freqs: if mains_freq > 0.0: indices = np.where(np.abs(good_freqs[:,0]  np.round(good_freqs[:,0]/mains_freq)*mains_freq) < mains_freq_tol)[0] if len(indices)>0: if verbose > 1: print('remove power line frequencies', ', '.join(['%.1f' % f for f in good_freqs[indices,0]])) good_freqs = np.delete(good_freqs, indices, axis=0) if verbose > 1: print('all_freqs: ', '[', ', '.join(['%.2f' % f for f in all_freqs[:,0] if f < 3000.0]), ']') group_list = [] fzero_harmonics_list = [] first = True # as long as there are frequencies left in good_freqs: fi = 0 while len(good_freqs) > 0: if fi < len(check_freqs): # check for harmonic group of a given fundamental frequency: fmax = check_freqs[fi] f0s = 'freq' good_freqs, harm_group, harm_indices, fzero_harmonics = \ retrieve_harmonic_group(fmax, good_freqs, all_freqs, freq_tol, max_freq_tol, verbose1, min_group_size) fi += 1 else: # check for harmonic groups: f0s = 'fmax' good_freqs, harm_group, harm_indices, fzero_harmonics, fmax = \ build_harmonic_group(good_freqs, all_freqs, freq_tol, max_freq_tol, verbose1, min_group_size, max_divisor) # nothing found: if len(harm_group) == 0: if verbose > 1  first: s = ' (largest peak)' if first else '' print('No harmonic group for %7.2fHz%s' % (fmax, s)) first = False continue first = False # fill up harmonic group: harm_group, harm_indices = expand_group(harm_group, harm_indices, all_freqs, freq_tol, max_harmonics) # check whether fundamental was filled in: first_h = 0 if harm_group[0,2] > 2 else 1 # increment use count: all_freqs[harm_indices,2] += 1 if harm_group[0,2] == 1: harm_group[0,2] = 2 # check frequency range of fundamental: fundamental_ok = (harm_group[0, 0] >= min_freq and harm_group[0, 0] <= max_freq) # check power hum: mains_ok = ((mains_freq <= 0.0) or (m.fabs(harm_group[0,0]  mains_freq) > freq_tol)) # check smoothness: db_powers = decibel(harm_group[first_h:,1]) diff = np.std(np.diff(db_powers)) smooth_ok = max_db_diff <= 0.0 or diff < max_db_diff # check relative power of higher harmonics: p_max = np.argmax(db_powers[:min_group_size]) db_powers = db_powers[p_max] amplitude_ok = len(db_powers[p_max+min_group_size:]) == 0 if amplitude_ok: pi = len(db_powers)  1 else: amplitude_ok = np.all((db_powers[p_max+min_group_size:] < max_harmonics_db)  (harm_group[first_h+p_max+min_group_size:,2] > 1)) if amplitude_ok: pi = p_max+min_group_size1 else: pi = np.where((db_powers[p_max+min_group_size:] >= max_harmonics_db) & (harm_group[first_h+p_max+min_group_size:,2] <= 1))[0][0] # check: if fundamental_ok and mains_ok and smooth_ok and amplitude_ok: if verbose > 0: print('Accepting harmonic group from %s=%7.2fHz: %7.2fHz power=%9.3g nharmonics=%2d, use count=%d' % (f0s, fmax, harm_group[0,0], np.sum(harm_group[first_h:,1]), len(harm_group), np.sum(harm_group[first_h:,2]))) group_list.append(harm_group[:,0:2]) fzero_harmonics_list.append(fzero_harmonics) else: if verbose > 0: fs = 'is ' if fundamental_ok else 'NOT' ms = 'not ' if mains_ok else 'IS' ss = 'smaller' if smooth_ok else 'LARGER ' ps = 'smaller' if amplitude_ok else 'LARGER ' print('Discarded harmonic group from %s=%7.2fHz: %7.2fHz power=%9.3g: %s in frequency range, %s mains frequency, smooth=%4.1fdB %s than %4.1fdB, relpower[%d]=%5.1fdB %s than %5.1fdB' % (f0s, fmax, harm_group[0,0], np.sum(harm_group[first_h:,1]), fs, ms, diff, ss, max_db_diff, pi, db_powers[pi], ps, max_harmonics_db)) # select most powerful harmonic groups: if max_groups > 0 and len(group_list) > max_groups: n = len(group_list) powers = [np.sum(group[:,1]) for group in group_list] powers_inx = np.argsort(powers) group_list = [group_list[pi] for pi in powers_inx[max_groups:]] fzero_harmonics_list = [fzero_harmonics_list[pi] for pi in powers_inx[max_groups:]] if verbose > 0: print('Selected from %d groups the %d most powerful groups.' % (n, max_groups)) # sort groups by fundamental frequency: freqs = [group[0,0] for group in group_list] freq_inx = np.argsort(freqs) group_list = [group_list[fi] for fi in freq_inx] fzero_harmonics_list = [fzero_harmonics_list[fi] for fi in freq_inx] if verbose > 1: print('') if len(group_list) > 0: print('## FUNDAMENTALS FOUND: ##') for group, fzero_h in zip(group_list, fzero_harmonics_list): print('%7.2fHz: power=%9.3g fzeroh=%2d' % (group[0,0], np.sum(group[:,1]), fzero_h)) else: print('## NO FUNDAMENTALS FOUND ##') # assemble mains frequencies from all_freqs: if mains_freq > 0.0: mains_freqs = all_freqs[np.abs(all_freqs[:,0]  np.round(all_freqs[:,0]/mains_freq)*mains_freq) < mains_freq_tol,:2] else: mains_freqs = np.zeros((0, 2)) return group_list, fzero_harmonics_list, mains_freqs
def fundamental_freqs(group_list)

Extract fundamental frequencies from lists of harmonic groups.
The inner list of 2D arrays of the input argument is transformed into a 1D array containig the fundamental frequencies extracted from the 2D arrays.
Parameters
group_list
:(list
of(list
of…)) list
of2D arrays
 Arbitrarily nested lists of harmonic groups as returned by extract_fundamentals() and harmonic_groups() with the element [0, 0] being the fundamental frequency.
Returns
fundamentals
:(list
of(list
of...)) 1D array
 Nested list (corresponding to
group_list
) of 1D arrays with the fundamental frequencies extracted from the harmonic groups.
Expand source code
def fundamental_freqs(group_list): """ Extract fundamental frequencies from lists of harmonic groups. The inner list of 2D arrays of the input argument is transformed into a 1D array containig the fundamental frequencies extracted from the 2D arrays. Parameters  group_list: (list of (list of ...)) list of 2D arrays Arbitrarily nested lists of harmonic groups as returned by extract_fundamentals() and harmonic_groups() with the element [0, 0] being the fundamental frequency. Returns  fundamentals: (list of (list of ...)) 1D array Nested list (corresponding to `group_list`) of 1D arrays with the fundamental frequencies extracted from the harmonic groups. """ if len(group_list) == 0: return np.array([]) # check whether group_list is list of harmonic groups: list_of_groups = True for group in group_list: if not ( hasattr(group, 'shape') and len(group.shape) == 2 ): list_of_groups = False break if list_of_groups: fundamentals = np.array([group[0, 0] for group in group_list if len(group) > 0]) else: fundamentals = [] for groups in group_list: f = fundamental_freqs(groups) fundamentals.append(f) return fundamentals
def fundamental_freqs_and_power(group_list, power=False, ref_power=1.0, min_power=1e20)

Extract fundamental frequencies and their power in dB from lists of harmonic groups.
The inner list of 2D arrays of the input argument is transformed into a 2D array containig for each fish (1st dimension) the fundamental frequencies and powers (summed over all harmonics) extracted from the 2D arrays.
Parameters
group_list
:(list
of(list
of…)) list
of2D arrays
 Arbitrarily nested lists of harmonic groups as returned by extract_fundamentals() and harmonic_groups() with the element [0, 0] being the fundamental frequency and the elements [:,1] being the powers of each harmonics.
power
:boolean
 If
False
convert the power into decibel using the powerspectrum.decibel() function. ref_power
:float
 Reference power for computing decibel.
If set to
None
the maximum power is used. min_power
:float
 Power values smaller than
min_power
are set tonp.inf
.
Returns
fundamentals
:(list
of(list
of...)) 2D array
 Nested list (corresponding to
group_list
) of 2D arrays with fundamental frequencies in first column and corresponding power in second column.
Expand source code
def fundamental_freqs_and_power(group_list, power=False, ref_power=1.0, min_power=1e20): """Extract fundamental frequencies and their power in dB from lists of harmonic groups. The inner list of 2D arrays of the input argument is transformed into a 2D array containig for each fish (1st dimension) the fundamental frequencies and powers (summed over all harmonics) extracted from the 2D arrays. Parameters  group_list: (list of (list of ...)) list of 2D arrays Arbitrarily nested lists of harmonic groups as returned by extract_fundamentals() and harmonic_groups() with the element [0, 0] being the fundamental frequency and the elements [:,1] being the powers of each harmonics. power: boolean If `False` convert the power into decibel using the powerspectrum.decibel() function. ref_power: float Reference power for computing decibel. If set to `None` the maximum power is used. min_power: float Power values smaller than `min_power` are set to `np.inf`. Returns  fundamentals: (list of (list of ...)) 2D array Nested list (corresponding to `group_list`) of 2D arrays with fundamental frequencies in first column and corresponding power in second column. """ if len(group_list) == 0: return np.array([]) # check whether group_list is list of harmonic groups: list_of_groups = True for group in group_list: if not ( hasattr(group, 'shape') and len(group.shape) == 2 ): list_of_groups = False break if list_of_groups: fundamentals = np.array([[group[0, 0], np.sum(group[:, 1])] for group in group_list if len(group) > 0]) if not power: fundamentals[:, 1] = decibel(fundamentals[:, 1], ref_power, min_power) else: fundamentals = [] for groups in group_list: f = fundamental_freqs_and_power(groups, power, ref_power, min_power) fundamentals.append(f) return fundamentals
def group_candidate(good_freqs, all_freqs, freq, divisor, freq_tol, max_freq_tol, min_group_size, verbose)

Candidate harmonic frequencies belonging to a fundamental frequency.
Parameters
good_freqs
:2D array
 Frequency, power, and use count (columns) of strong peaks detected in a power spectrum.
all_freqs
:2D array
 Frequency, power, and use count (columns) of all peaks detected in a power spectrum.
freq
:float
 Fundamental frequency for which a harmonic group should be assembled.
divisor
:int
 Fundamental frequency was obtained from a frequency divided by divisor. 0: Fundamental frequency is given as is.
freq_tol
:float
 Harmonics should fall within this frequency tolerance. This should be in the range of the frequency resolution and should not be smaller than half of the frequency resolution.
max_freq_tol
:float
 Maximum deviation of harmonics from their expected frequency.
Peaks with frequencies between
freq_tol
andmax_freq_tol
get penalized the further away from their expected frequency. min_group_size
:int
 Minimum number of harmonics of a harmonic group. The harmonics from min_group_size/3 to max(min_group_size, divisor) need to be in good_freqs.
verbose
:int
 Verbosity level.
Returns
new_group
:1D aray
ofindices
 Frequencies in all_freqs belonging to the candidate group.
fzero
:float
 Adjusted fundamental frequency. If negative, no group was found.
fzero_harmonics
:int
 The highest harmonics that was used to recompute the fundamental frequency.
Expand source code
def group_candidate(good_freqs, all_freqs, freq, divisor, freq_tol, max_freq_tol, min_group_size, verbose): """Candidate harmonic frequencies belonging to a fundamental frequency. Parameters  good_freqs: 2D array Frequency, power, and use count (columns) of strong peaks detected in a power spectrum. all_freqs: 2D array Frequency, power, and use count (columns) of all peaks detected in a power spectrum. freq: float Fundamental frequency for which a harmonic group should be assembled. divisor: int Fundamental frequency was obtained from a frequency divided by divisor. 0: Fundamental frequency is given as is. freq_tol: float Harmonics should fall within this frequency tolerance. This should be in the range of the frequency resolution and should not be smaller than half of the frequency resolution. max_freq_tol: float Maximum deviation of harmonics from their expected frequency. Peaks with frequencies between `freq_tol` and `max_freq_tol` get penalized the further away from their expected frequency. min_group_size: int Minimum number of harmonics of a harmonic group. The harmonics from min_group_size/3 to max(min_group_size, divisor) need to be in good_freqs. verbose: int Verbosity level. Returns  new_group: 1D aray of indices Frequencies in all_freqs belonging to the candidate group. fzero: float Adjusted fundamental frequency. If negative, no group was found. fzero_harmonics: int The highest harmonics that was used to recompute the fundamental frequency. """ if verbose > 0: print('') dvs = True if divisor <= 0: divisor = 1 dvs = False # 1. find harmonics in good_freqs and adjust fzero accordingly: fzero = freq fzero_harmonics = 1 if len(good_freqs[:,0]) > 0: """ ff = good_freqs[:,0] h = np.round(ff/fzero) h_indices = np.unique(h, return_index=True)[1] f_best = [hi+np.argmin(np.abs(fh/h  fzero)) for hi, fh in zip(h_indices, np.split(ff, h_indices))] fe = np.abs(ff/h  fzero) h = h[fe < freq_tol] fe = fe[fe < freq_tol] """ prev_freq = divisor * freq for h in range(divisor+1, 2*min_group_size+1): idx = np.argmin(np.abs(good_freqs[:,0]/h  fzero)) ff = good_freqs[idx,0] if m.fabs(ff/h  fzero) > freq_tol: continue df = ff  prev_freq dh = np.round(df/fzero) fe = m.fabs(df/dh  fzero) if fe > 2.0*freq_tol: if h > min_group_size: break continue # update fzero: prev_freq = ff fzero_harmonics = h fzero = ff/fzero_harmonics if verbose > 1: print('adjusted fzero to %.2fHz from %dth harmonic' % (fzero, h)) if verbose > 0: ds = 'divisor: %d, ' % divisor if dvs else '' print('# %sfzero=%7.2fHz adjusted from harmonics %d' % (ds, fzero, fzero_harmonics)) # 2. check fzero: # freq might not be in our group anymore, because fzero was adjusted: if np.abs(freq  fzero) > freq_tol: if verbose > 0: print(' discarded: lost frequency') return [], 1.0, fzero_harmonics # 3. collect harmonics from all_freqs: new_group = np.ones(min_group_size, dtype=np.int) new_penalties = np.ones(min_group_size) freqs = [] prev_h = 0 prev_fe = 0.0 for h in range(1, min_group_size+1): penalty = 0 i = np.argmin(np.abs(all_freqs[:,0]/h  fzero)) f = all_freqs[i,0] if verbose > 2: print(' check %7.2fHz as %d. harmonics' % (f, h)) fac = 1.0 if h >= divisor else 2.0 fe = m.fabs(f/h  fzero) if fe > fac*max_freq_tol: if verbose > 1 and fe < 2*fac*max_freq_tol: print(' %d. harmonics at %7.2fHz is off by %7.2fHz (max between %5.2fHz and %5.2fHz) from %7.2fHz' % (h, f, h*fe, h*fac*freq_tol, h*fac*max_freq_tol, h*fzero)) continue if fe > fac*freq_tol: penalty = np.interp(fe, [fac*freq_tol, fac*max_freq_tol], [0.0, 1.0]) if len(freqs) > 0: pf = freqs[1] df = f  pf if df < 0.5*fzero: if len(freqs)>1: pf = freqs[2] df = f  pf else: pf = 0.0 df = h*fzero dh = m.floor(df/fzero + 0.5) fe = m.fabs(df/dh  fzero) if fe > 2*dh*fac*max_freq_tol: if verbose > 1: print(' %d. harmonics at %7.2fHz is off by %7.2fHz (max between %5.2fHz and %5.2fHz) from previous harmonics at %7.2fHz' % (h, f, dh*fe, 2*fac*dh*freq_tol, 2*fac*dh*max_freq_tol, pf)) continue if fe > 2*dh*fac*freq_tol: penalty = np.interp(fe, [2*dh*fac*freq_tol, 2*dh*fac*max_freq_tol], [0.0, 1.0]) else: fe = 0.0 if h > prev_h or fe < prev_fe: if prev_h > 0 and h  prev_h > 1: if verbose > 1: print(' previous harmonics %d more than 1 away from %d. harmonics at %7.2fHz' % (prev_h, h, f)) break if h == prev_h and len(freqs) > 0: freqs.pop() freqs.append(f) new_group[int(h)1] = i new_penalties[int(h)1] = penalty prev_h = h prev_fe = fe if verbose > 1: print(' %d. harmonics at %7.2fHz has been taken (peak %2d) with penalty %3.1f' % (h, f, i, penalty)) # 4. check new group: # almost all harmonics in min_group_size required: max_penalties = min_group_size/3 if np.sum(new_penalties) > max_penalties: if verbose > 0: print(' discarded group because sum of penalties %3.1f is more than %3.1f: indices' % (np.sum(new_penalties), max_penalties), new_group, ' penalties', new_penalties) return [], 1.0, fzero_harmonics new_group = new_group[new_group>=0] # check use count of frequencies: double_use = np.sum(all_freqs[new_group, 2]>0) if double_use >= 2: # XXX make this a parameter? if verbose > 0: print(' discarded group because of use count = %2d >= 2' % double_use) return [], 1.0, fzero_harmonics # 5. return results: return new_group, fzero, fzero_harmonics
def harmonic_groups(psd_freqs, psd, verbose=0, check_freqs=[], low_threshold=0.0, high_threshold=0.0, thresh_bins=100, low_thresh_factor=6.0, high_thresh_factor=10.0, freq_tol_fac=1.0, max_freq_tol=1.0, mains_freq=60.0, mains_freq_tol=1.0, min_freq=0.0, max_freq=2000.0, max_db_diff=20.0, max_divisor=4, min_group_size=3, max_harmonics_db=5.0, max_harmonics=0, max_groups=0, **kwargs)

Detect peaks in power spectrum and group them according to their harmonic structure.
Parameters
psd_freqs
:1D array
 Frequencies of the power spectrum.
psd
:1D array
 Power spectrum (linear, not decible).
verbose
:int
 Verbosity level.
check_freqs
:list
offloat
 List of fundamental frequencies that will be checked first for being present and valid harmonic groups in the peak frequencies of the power spectrum.
low_threshold
:float
 The relative threshold for detecting all peaks in the decibel spectrum.
high_threshold
:float
 The relative threshold for detecting good peaks in the decibel spectrum.
thresh_bins
:int
orlist
offloats
 Number of bins or the bins for computing the histogram from
which the standard deviation of the noise level in the
psd
is estimated. low_thresh_factor
:float
 Factor by which the estimated standard deviation of the noise floor
is multiplied to set the
low_threshold
. high_thresh_factor
:float
 Factor by which the estimated standard deviation of the noise floor
is multiplied to set the
high_threshold
. freq_tol_fac
:float
 Harmonics should fall within
deltaf*freq_tol_fac
. max_freq_tol
:float
 Maximum absolute frequency deviation of harmonics in Hertz..
mains_freq
:float
 Frequency of the mains power supply.
mains_freq_tol
:float
 Tolarerance around harmonics of the mains frequency, within which peaks are removed.
min_freq
:float
 Minimum frequency accepted as a fundamental frequency.
max_freq
:float
 Maximum frequency accepted as a fundamental frequency.
max_db_diff
:float
 If larger than zero, maximum standard deviation of differences between logarithmic powers of harmonics in decibel (larger than zero). Low values enforce smoother power spectra.
max_divisor
:int
 Maximum divisor used for checking for subharmonics.
min_group_size
:int
 Minimum number of harmonics of a harmonic group. The harmonics from min_group_size/3 to max(min_group_size, divisor) need to be in good_freqs.
max_harmonics_db
:float
 Maximum allowed power of the
min_group_size
th and higher harmonics after the peak (in decibel relative to peak power withn the firstmin_group_size
harmonics, i.e. if harmonics are required to be smaller than fundamental then this is a negative number). Make it a large positive number to effectively not check for relative power. max_harmonics
:int
 Maximum number of harmonics to be returned for each group.
max_groups
:int
 If not zero the maximum number of most powerful harmonic groups.
Returns
groups
:list
of2D arrays
 List of all extracted harmonic groups, sorted by fundamental frequency. Each harmonic group is a 2D array with the first dimension the harmonics and the second dimension containing frequency and power of each harmonic. If the power is zero, there was no corresponding peak in the power spectrum.
fzero_harmonics
:list
ofints
 The harmonics from which the fundamental frequencies were computed.
mains
:2d array
 Frequencies and power of multiples of the mains frequency found in the power spectrum.
all_freqs
:2D array
 Frequency, power, and use count (columns) of all peaks detected in the power spectrum.
good_freqs
:1D array
 Frequencies of strong peaks detected in the power spectrum.
low_threshold
:float
 The relative threshold for detecting all peaks in the decibel spectrum.
high_threshold
:float
 The relative threshold for detecting good peaks in the decibel spectrum.
center
:float
 The baseline level of the power spectrum.
Expand source code
def harmonic_groups(psd_freqs, psd, verbose=0, check_freqs=[], low_threshold=0.0, high_threshold=0.0, thresh_bins=100, low_thresh_factor=6.0, high_thresh_factor=10.0, freq_tol_fac=1.0, max_freq_tol=1.0, mains_freq=60.0, mains_freq_tol=1.0, min_freq=0.0, max_freq=2000.0, max_db_diff=20.0, max_divisor=4, min_group_size=3, max_harmonics_db=5.0, max_harmonics=0, max_groups=0, **kwargs): """Detect peaks in power spectrum and group them according to their harmonic structure. Parameters  psd_freqs: 1D array Frequencies of the power spectrum. psd: 1D array Power spectrum (linear, not decible). verbose: int Verbosity level. check_freqs: list of float List of fundamental frequencies that will be checked first for being present and valid harmonic groups in the peak frequencies of the power spectrum. low_threshold: float The relative threshold for detecting all peaks in the decibel spectrum. high_threshold: float The relative threshold for detecting good peaks in the decibel spectrum. thresh_bins: int or list of floats Number of bins or the bins for computing the histogram from which the standard deviation of the noise level in the `psd` is estimated. low_thresh_factor: float Factor by which the estimated standard deviation of the noise floor is multiplied to set the `low_threshold`. high_thresh_factor: float Factor by which the estimated standard deviation of the noise floor is multiplied to set the `high_threshold`. freq_tol_fac: float Harmonics should fall within `deltaf*freq_tol_fac`. max_freq_tol: float Maximum absolute frequency deviation of harmonics in Hertz.. mains_freq: float Frequency of the mains power supply. mains_freq_tol: float Tolarerance around harmonics of the mains frequency, within which peaks are removed. min_freq: float Minimum frequency accepted as a fundamental frequency. max_freq: float Maximum frequency accepted as a fundamental frequency. max_db_diff: float If larger than zero, maximum standard deviation of differences between logarithmic powers of harmonics in decibel (larger than zero). Low values enforce smoother power spectra. max_divisor: int Maximum divisor used for checking for subharmonics. min_group_size: int Minimum number of harmonics of a harmonic group. The harmonics from min_group_size/3 to max(min_group_size, divisor) need to be in good_freqs. max_harmonics_db: float Maximum allowed power of the `min_group_size`th and higher harmonics after the peak (in decibel relative to peak power withn the first `min_group_size` harmonics, i.e. if harmonics are required to be smaller than fundamental then this is a negative number). Make it a large positive number to effectively not check for relative power. max_harmonics: int Maximum number of harmonics to be returned for each group. max_groups: int If not zero the maximum number of most powerful harmonic groups. Returns  groups: list of 2D arrays List of all extracted harmonic groups, sorted by fundamental frequency. Each harmonic group is a 2D array with the first dimension the harmonics and the second dimension containing frequency and power of each harmonic. If the power is zero, there was no corresponding peak in the power spectrum. fzero_harmonics: list of ints The harmonics from which the fundamental frequencies were computed. mains: 2d array Frequencies and power of multiples of the mains frequency found in the power spectrum. all_freqs: 2D array Frequency, power, and use count (columns) of all peaks detected in the power spectrum. good_freqs: 1D array Frequencies of strong peaks detected in the power spectrum. low_threshold: float The relative threshold for detecting all peaks in the decibel spectrum. high_threshold: float The relative threshold for detecting good peaks in the decibel spectrum. center: float The baseline level of the power spectrum. """ if verbose > 0: print('') if verbose > 1: print(70*'#') print('##### harmonic_groups', 48*'#') # decibel power spectrum: log_psd = decibel(psd) max_idx = np.argmax(~np.isfinite(log_psd)) if max_idx > 0: log_psd = log_psd[:max_idx] delta_f = psd_freqs[1]  psd_freqs[0] # thresholds: center = np.NaN if low_threshold <= 0.0 or high_threshold <= 0.0: low_th, high_th, center = threshold_estimate(log_psd, low_thresh_factor, high_thresh_factor, thresh_bins) if low_threshold <= 0.0: low_threshold = low_th if high_threshold <= 0.0: high_threshold = high_th if verbose > 1: print('') print('low_threshold =%4.1fdB, center+low_threshold =%6.1fdB' % (low_threshold, center + low_threshold)) print('high_threshold=%4.1fdB, center+high_threshold=%6.1fdB' % (high_threshold, center + high_threshold)) print(' center=%6.1fdB' % center) # detect peaks in decibel power spectrum: peaks, troughs = detect_peaks(log_psd, low_threshold) peaks, troughs = trim(peaks, troughs) all_freqs = np.zeros((len(peaks), 3)) all_freqs[:,0] = psd_freqs[peaks] all_freqs[:,1] = psd[peaks] if len(all_freqs) == 0: return [], [], [], np.zeros((0, 3)), [], low_threshold, high_threshold, center # select good peaks: good_freqs = all_freqs[(log_psd[peaks]  log_psd[troughs] > high_threshold) & (all_freqs[:,0] >= min_freq) & (all_freqs[:,0] < max_freq*max_divisor),:] # detect harmonic groups: freq_tol = delta_f*freq_tol_fac if max_freq_tol < 1.1*freq_tol: max_freq_tol = 1.1*freq_tol groups, fzero_harmonics, mains = \ extract_fundamentals(good_freqs, all_freqs, freq_tol, max_freq_tol, verbose, check_freqs, mains_freq, mains_freq_tol, min_freq, max_freq, max_db_diff, max_divisor, min_group_size, max_harmonics_db, max_harmonics, max_groups) return (groups, fzero_harmonics, mains, all_freqs, good_freqs[:,0], low_threshold, high_threshold, center)
def harmonic_groups_args(cfg)

Translates a configuration to the respective parameter names of the harmonicgroup detection functions. The return value can then be passed as keyword arguments to this function.
Parameters
cfg
:ConfigFile
 The configuration.
Returns
a
:dict
 Dictionary with names of arguments of the harmonicgroup detection functions
and their values as supplied by
cfg
.
Expand source code
def harmonic_groups_args(cfg): """Translates a configuration to the respective parameter names of the harmonicgroup detection functions. The return value can then be passed as keyword arguments to this function. Parameters  cfg: ConfigFile The configuration. Returns  a: dict Dictionary with names of arguments of the harmonicgroup detection functions and their values as supplied by `cfg`. """ return cfg.map({'mains_freq': 'mainsFreq', 'mains_freq_tol': 'mainsFreqTolerance', 'freq_tol_fac': 'freqTolerance', 'max_freq_tol': 'maximumFreqTolerance', 'max_divisor': 'maxDivisor', 'min_group_size': 'minimumGroupSize', 'min_freq': 'minimumFrequency', 'max_freq': 'maximumFrequency', 'max_db_diff': 'maximumPowerDifference', 'max_harmonics_db': 'maximumHarmonicsPower', 'max_harmonics': 'maximumHarmonics', 'max_groups': 'maximumGroups'})
def main(data_file=None)

Expand source code
def main(data_file=None): import matplotlib.pyplot as plt from .fakefish import wavefish_eods from .powerspectrum import psd if data_file is None: # generate data: title = 'simulation' samplerate = 44100.0 d = 20.0 noise = 0.01 eodfs = [123.0, 333.0, 666.0, 666.5] fish1 = 0.5*wavefish_eods('Eigenmannia', eodfs[0], samplerate, duration=d, noise_std=noise) fish2 = 1.0*wavefish_eods('Eigenmannia', eodfs[1], samplerate, duration=d, noise_std=noise) fish3 = 10.0*wavefish_eods('Alepto', eodfs[2], samplerate, duration=d, noise_std=noise) fish4 = 6.0*wavefish_eods('Alepto', eodfs[3], samplerate, duration=d, noise_std=noise) data = fish1 + fish2 + fish3 + fish4 else: from .dataloader import load_data print("load %s ..." % data_file) data, samplerate, unit = load_data(data_file, 0) title = data_file # retrieve fundamentals from power spectrum: psd_data = psd(data, samplerate, freq_resolution=0.1) groups, _, mains, all_freqs, good_freqs, _, _, _ = harmonic_groups(psd_data[0], psd_data[1], check_freqs=[123.0, 666.0], max_db_diff=30.0, verbose=0) fig, ax = plt.subplots() plot_psd_harmonic_groups(ax, psd_data[0], psd_data[1], groups, mains, all_freqs, good_freqs, max_freq=3000.0) ax.set_title(title) plt.show() # unify fundamental frequencies: fundamentals = fundamental_freqs(groups) np.set_printoptions(formatter={'float': lambda x: '%5.1f' % x}) print('fundamental frequencies extracted from power spectrum:') print(fundamentals) print('') freqs = fundamental_freqs_and_power([groups]) freqs.append(np.array([[44.0, 20.0], [44.2, 10.0], [320.5, 2.5], [665.5, 5.0], [666.2, 10.0]])) freqs.append(np.array([[123.3, 1.0], [320.2, 2.0], [668.4, 2.0]])) rank_freqs = add_relative_power(freqs) rank_freqs = add_power_ranks(rank_freqs) print('all frequencies (frequency, power, relpower, rank):') print('\n'.join(( str(f) for f in rank_freqs))) print('') indices = similar_indices(freqs, 1.0) print('similar indices:') print('\n'.join(( ('\n '.join((str(f) for f in g)) for g in indices)))) print('') unique_freqs = unique(freqs, 1.0, 'power') print('unique power:') print('\n'.join(( str(f) for f in unique_freqs))) print('') unique_freqs = unique(freqs, 1.0, 'relpower') print('unique relative power:') print('\n'.join(( str(f) for f in unique_freqs))) print('') unique_freqs = unique(freqs, 1.0, 'rank') print('unique rank:') print('\n'.join(( str(f) for f in unique_freqs))) print('') unique_freqs = unique(freqs, 1.0, 'rank', 1) print('unique rank for next neighour only:') print('\n'.join(( str(f) for f in unique_freqs))) print('')
def plot_harmonic_groups(ax, group_list, indices=None, max_groups=0, skip_bad=False, sort_by_freq=True, label_power=False, colors=None, markers=None, legend_rows=8, **kwargs)

Mark decibel power of fundamentals and their harmonics in a plot.
Parameters
ax
:axis for plot
 Axis used for plotting.
group_list
:list
of2D arrays
 Lists of harmonic groups as returned by extract_fundamentals() and harmonic_groups() with the element [0, 0] of the harmonic groups being the fundamental frequency, and element[0, 1] being the corresponding power.
indices
:list
ofint
orNone
 If smaller than zero then set the legend label of the corresponding group in brackets.
max_groups
:int
 If not zero plot only the max_groups most powerful groups.
skip_bad
:bool
 Skip harmonic groups without index (entry in indices is negative).
sort_by_freq
:boolean
 If True sort legend by frequency, otherwise by power.
label_power
:boolean
 If
True
put the power in decibel in addition to the frequency into the legend. colors
:list
ofcolors
orNone
 If not None list of colors for plotting each group
markers
:list
ofmarkers
orNone
 If not None list of markers for plotting each group
legend_rows
:int
 Maximum number of rows to be used for the legend.
kwargs
 Key word arguments for the legend of the plot.
Expand source code
def plot_harmonic_groups(ax, group_list, indices=None, max_groups=0, skip_bad=False, sort_by_freq=True, label_power=False, colors=None, markers=None, legend_rows=8, **kwargs): """Mark decibel power of fundamentals and their harmonics in a plot. Parameters  ax: axis for plot Axis used for plotting. group_list: list of 2D arrays Lists of harmonic groups as returned by extract_fundamentals() and harmonic_groups() with the element [0, 0] of the harmonic groups being the fundamental frequency, and element[0, 1] being the corresponding power. indices: list of int or None If smaller than zero then set the legend label of the corresponding group in brackets. max_groups: int If not zero plot only the max_groups most powerful groups. skip_bad: bool Skip harmonic groups without index (entry in indices is negative). sort_by_freq: boolean If True sort legend by frequency, otherwise by power. label_power: boolean If `True` put the power in decibel in addition to the frequency into the legend. colors: list of colors or None If not None list of colors for plotting each group markers: list of markers or None If not None list of markers for plotting each group legend_rows: int Maximum number of rows to be used for the legend. kwargs: Key word arguments for the legend of the plot. """ if len(group_list) == 0: return # sort by power: powers = np.array([np.sum(group[:,1]) for group in group_list]) max_power = np.max(powers) idx = np.argsort(powers) if max_groups > 0 and len(idx > max_groups): idx = idx[:max_groups] # sort by frequency: if sort_by_freq: freqs = np.array([group_list[group][0,0] for group in idx]) if legend_rows > 0 and legend_rows < len(freqs): idx = idx[np.argsort(freqs[:legend_rows])] else: idx = idx[np.argsort(freqs)] # plot: k = 0 for i in idx: if indices is not None and skip_bad and indices[i] < 0: continue group = group_list[i] x = group[:,0] y = decibel(group[:,1]) msize = 7.0 + 10.0*(powers[i]/max_power)**0.25 color_kwargs = {} if colors is not None: color_kwargs = {'color': colors[k%len(colors)]} label = '%6.1f Hz' % group[0, 0] if label_power: label += ' %6.1f dB' % decibel(np.array([np.sum(group[:,1])]))[0] if indices is not None: if indices[i] < 0: label = '(' + label + ')' else: label = ' ' + label + ' ' if legend_rows > 5 and k >= legend_rows: label = None if markers is None: ax.plot(x, y, 'o', ms=msize, label=label, **color_kwargs) else: if k >= len(markers): break ax.plot(x, y, linestyle='None', marker=markers[k], mec=None, mew=0.0, ms=msize, label=label, **color_kwargs) k += 1 # legend: if legend_rows > 0: if legend_rows > 5: ncol = 1 else: ncol = (len(idx)1) // legend_rows + 1 ax.legend(numpoints=1, ncol=ncol, **kwargs) else: ax.legend(numpoints=1, **kwargs)
def plot_psd_harmonic_groups(ax, psd_freqs, psd, group_list, mains=None, all_freqs=None, good_freqs=None, log_freq=False, min_freq=0.0, max_freq=2000.0, ymarg=0.0)

Plot decibel powerspectrum with detected peaks, harmonic groups, and mains frequencies.
Parameters
psd_freqs
:array
 Frequencies of the power spectrum.
psd
:array
 Power spectrum (linear, not decible).
group_list
:list
of2D arrays
 Lists of harmonic groups as returned by extract_fundamentals() and harmonic_groups() with the element [0, 0] of the harmonic groups being the fundamental frequency, and element[0, 1] being the corresponding power.
mains
:2D array
 Frequencies and power of multiples of the mains frequency found in the power spectrum.
all_freqs
:2D array
 Peaks in the power spectrum detected with low threshold.
good_freqs
:1D array
 Frequencies of peaks detected with high threshold.
log_freq
:boolean
 Logarithmic (True) or linear (False) frequency axis.
min_freq
:float
 Limits of frequency axis are set to
(min_freq, max_freq)
ifmax_freq
is greater than zero max_freq
:float
 Limits of frequency axis are set to
(min_freq, max_freq)
and limits of power axis are computed from powers below max_freq ifmax_freq
is greater than zero ymarg
:float
 Add this to the maximum decibel power for setting the ylim.
Expand source code
def plot_psd_harmonic_groups(ax, psd_freqs, psd, group_list, mains=None, all_freqs=None, good_freqs=None, log_freq=False, min_freq=0.0, max_freq=2000.0, ymarg=0.0): """Plot decibel powerspectrum with detected peaks, harmonic groups, and mains frequencies. Parameters  psd_freqs: array Frequencies of the power spectrum. psd: array Power spectrum (linear, not decible). group_list: list of 2D arrays Lists of harmonic groups as returned by extract_fundamentals() and harmonic_groups() with the element [0, 0] of the harmonic groups being the fundamental frequency, and element[0, 1] being the corresponding power. mains: 2D array Frequencies and power of multiples of the mains frequency found in the power spectrum. all_freqs: 2D array Peaks in the power spectrum detected with low threshold. good_freqs: 1D array Frequencies of peaks detected with high threshold. log_freq: boolean Logarithmic (True) or linear (False) frequency axis. min_freq: float Limits of frequency axis are set to `(min_freq, max_freq)` if `max_freq` is greater than zero max_freq: float Limits of frequency axis are set to `(min_freq, max_freq)` and limits of power axis are computed from powers below max_freq if `max_freq` is greater than zero ymarg: float Add this to the maximum decibel power for setting the ylim. """ # mark all and good psd peaks: pmin, pmax = ax.get_ylim() doty = pmax  5.0 if all_freqs is not None: ax.plot(all_freqs[:, 0], np.zeros(len(all_freqs[:, 0])) + doty, 'o', color='#ffffff') if good_freqs is not None: ax.plot(good_freqs, np.zeros(len(good_freqs)) + doty, 'o', color='#888888') # mark mains frequencies: if mains is not None and len(mains) > 0: fpeaks = mains[:, 0] fpeakinx = [int(np.round(fp/(psd_freqs[1]psd_freqs[0]))) for fp in fpeaks if fp < psd_freqs[1]] ax.plot(fpeaks[:len(fpeakinx)], decibel(psd[fpeakinx]), linestyle='None', marker='.', color='k', ms=10, mec=None, mew=0.0, label='%3.0f Hz mains' % mains[0, 0]) # mark harmonic groups: colors, markers = colors_markers() plot_harmonic_groups(ax, group_list, max_groups=0, sort_by_freq=True, colors=colors, markers=markers, legend_rows=8, loc='upper right') # plot power spectrum: plot_decibel_psd(ax, psd_freqs, psd, log_freq=log_freq, min_freq=min_freq, max_freq=max_freq, ymarg=ymarg, color='blue')
def psd_peak_detection_args(cfg)

Translates a configuration to the respective parameter names for the detection of peaks in power spectrum used by harmonic_groups(). The return value can then be passed as keyword arguments to this function.
Parameters
cfg
:ConfigFile
 The configuration.
Returns
a
:dict
 Dictionary with names of arguments of the
harmonicgroup()
function and their values as supplied bycfg
.
Expand source code
def psd_peak_detection_args(cfg): """Translates a configuration to the respective parameter names for the detection of peaks in power spectrum used by harmonic_groups(). The return value can then be passed as keyword arguments to this function. Parameters  cfg: ConfigFile The configuration. Returns  a: dict Dictionary with names of arguments of the `harmonicgroup()` function and their values as supplied by `cfg`. """ return cfg.map({'low_threshold': 'lowThreshold', 'high_threshold': 'highThreshold', 'thresh_bins': 'thresholdBins', 'low_thresh_factor': 'lowThresholdFactor', 'high_thresh_factor': 'highThresholdFactor'})
def retrieve_harmonic_group(freq, good_freqs, all_freqs, freq_tol, max_freq_tol, verbose=0, min_group_size=3)

Find all the harmonics belonging to a given fundamental.
Parameters
freq
:float
 Fundamental frequency for which harmonics are to be retrieved.
good_freqs
:2D array
 Frequency, power, and use count (columns) of strong peaks detected
in a power spectrum. All harmonics of
freq
will be removed fromgood_freqs
. all_freqs
:2D array
 Frequency, power, and use count (columns) of all peaks detected in a power spectrum.
freq_tol
:float
 Harmonics should fall within this frequency tolerance. This should be in the range of the frequency resolution and should not be smaller than half of the frequency resolution.
max_freq_tol
:float
 Maximum deviation of harmonics from their expected frequency.
Peaks with frequencies between
freq_tol
andmax_freq_tol
get penalized the further away from their expected frequency. verbose
:int
 Verbosity level.
min_group_size
:int
 Minimum number of harmonics of a harmonic group. The harmonics from min_group_size/3 to max(min_group_size, divisor) need to be in good_freqs.
Returns
good_freqs
:2D array
 Frequency, power, and use count (columns) of strong peaks detected in a power spectrum with frequencies of the returned harmonic group removed.
group
:2D array
 The detected harmonic group. Might be empty.
indices
:1D array
ofindices
 Indices of the harmonic group in all_freqs.
fzero_harmonics
:int
 The highest harmonics that was used to recompute the fundamental frequency.
Expand source code
def retrieve_harmonic_group(freq, good_freqs, all_freqs, freq_tol, max_freq_tol, verbose=0, min_group_size=3): """Find all the harmonics belonging to a given fundamental. Parameters  freq: float Fundamental frequency for which harmonics are to be retrieved. good_freqs: 2D array Frequency, power, and use count (columns) of strong peaks detected in a power spectrum. All harmonics of `freq` will be removed from `good_freqs`. all_freqs: 2D array Frequency, power, and use count (columns) of all peaks detected in a power spectrum. freq_tol: float Harmonics should fall within this frequency tolerance. This should be in the range of the frequency resolution and should not be smaller than half of the frequency resolution. max_freq_tol: float Maximum deviation of harmonics from their expected frequency. Peaks with frequencies between `freq_tol` and `max_freq_tol` get penalized the further away from their expected frequency. verbose: int Verbosity level. min_group_size: int Minimum number of harmonics of a harmonic group. The harmonics from min_group_size/3 to max(min_group_size, divisor) need to be in good_freqs. Returns  good_freqs: 2D array Frequency, power, and use count (columns) of strong peaks detected in a power spectrum with frequencies of the returned harmonic group removed. group: 2D array The detected harmonic group. Might be empty. indices: 1D array of indices Indices of the harmonic group in all_freqs. fzero_harmonics: int The highest harmonics that was used to recompute the fundamental frequency. """ if verbose > 0: print('') print('%s retrieve harmonic group %s' % (10*'#', 35*'#')) print('%s freq=%7.2fHz %s' % (10*'#', freq, 44*'#')) print('good_freqs: ', '[', ', '.join(['%.2f' % f for f in good_freqs[:,0]]), ']') # find harmonics in good_freqs and adjust fzero accordingly: new_group, fzero, fzero_harmonics = group_candidate(good_freqs, all_freqs, freq, 0, freq_tol, max_freq_tol, min_group_size, verbose) # no group found: if fzero < 0.0: return good_freqs, np.zeros((0, 2)), np.zeros(0), fzero_harmonics if verbose > 1: print('') print('# group found for freq=%.2fHz, fzero=%.2fHz:' % (freq, fzero)) group_str = '#### resulting harmonic group for freq=%.2fHz' % freq good_freqs, group = update_group(good_freqs, all_freqs, new_group, fzero, freq_tol, verbose, group_str) # good_freqs: removed all frequencies of bestgroup return good_freqs, group, new_group, fzero_harmonics
def similar_indices(freqs, df_thresh, nextfs=0)

Indices of similar frequencies.
If two frequencies from different elements in the inner lists of
freqs
are reciprocally the closest to each other and closer thandf_thresh
, then two indices (element, frequency) of the respective other frequency are appended.Parameters
freqs
:(list
of(list
of…)) list
of2D arrays
 First column in the arrays is fundamental frequency.
df_thresh
:float
 Fundamental frequencies closer than this threshold are considered equal.
nextfs
:int
 If zero, compare all elements in freqs with each other. Otherwise,
only compare with the
nextfs
next elements in freqs.
Returns
indices
:(list
of(list
of…)) list
oflist
oftwotuples
ofint
 For each frequency of each element in
freqs
a list of two tuples containing the indices of elements and frequencies that are similar.
Expand source code
def similar_indices(freqs, df_thresh, nextfs=0): """Indices of similar frequencies. If two frequencies from different elements in the inner lists of `freqs` are reciprocally the closest to each other and closer than `df_thresh`, then two indices (element, frequency) of the respective other frequency are appended. Parameters  freqs: (list of (list of ...)) list of 2D arrays First column in the arrays is fundamental frequency. df_thresh: float Fundamental frequencies closer than this threshold are considered equal. nextfs: int If zero, compare all elements in freqs with each other. Otherwise, only compare with the `nextfs` next elements in freqs. Returns  indices: (list of (list of ...)) list of list of twotuples of int For each frequency of each element in `freqs` a list of two tuples containing the indices of elements and frequencies that are similar. """ if len(freqs) == 0: return [] # check whether freqs is list of fundamental frequencies and powers: list_of_freq_power = True for group in freqs: if not (hasattr(group, 'shape') and len(group.shape) == 2): list_of_freq_power = False break if list_of_freq_power: indices = [ [[] for j in range(len(freqs[i]))] for i in range(len(freqs))] for j in range(len(freqs)1): freqsj = np.asarray(freqs[j]) for m in range(len(freqsj)): freq1 = freqsj[m] nn = len(freqs) if nextfs == 0 else j+1+nextfs if nn > len(freqs): nn = len(freqs) for k in range(j+1, nn): freqsk = np.asarray(freqs[k]) if len(freqsk) == 0: continue n = np.argmin(np.abs(freqsk[:,0]  freq1[0])) freq2 = freqsk[n] if np.argmin(np.abs(freqsj[:,0]  freq2[0])) != m: continue if np.abs(freq1[0]  freq2[0]) < df_thresh: indices[k][n].append((j, m)) indices[j][m].append((k, n)) else: indices = [] for groups in freqs: indices.append(similar_indices(groups, df_thresh, nextfs)) return indices
def threshold_estimate(psd_data, low_thresh_factor=6.0, high_thresh_factor=10.0, nbins=100, hist_height=0.6065306597126334)

Estimate thresholds for peak detection from histogram of power spectrum.
The standard deviation of the noise floor without peaks is estimated from the width of the histogram of the power spectrum at
hist_height
relative height. The histogram is computed in the third quarter of the linearly detrended power spectrum.Parameters
psd_data
:1D array
 The power spectrum from which to estimate the thresholds.
low_thresh_factor
:float
 Factor by which the estimated standard deviation of the noise floor
is multiplied to set the
low_threshold
. high_thresh_factor
:float
 Factor by which the estimated standard deviation of the noise floor
is multiplied to set the
high_threshold
. nbins
:int
orlist
offloats
 Number of bins or the bins for computing the histogram.
hist_height
:float
 Height between 0 and 1 at which the standard deviation of the histogram is estimated.
Returns
low_threshold
:float
 The threshold for peaks just above the noise floor.
high_threshold
:float
 The threshold for distinct peaks.
center
:float
 The baseline level of the power spectrum.
Expand source code
def threshold_estimate(psd_data, low_thresh_factor=6.0, high_thresh_factor=10.0, nbins=100, hist_height=1.0/ np.sqrt(np.e)): """Estimate thresholds for peak detection from histogram of power spectrum. The standard deviation of the noise floor without peaks is estimated from the width of the histogram of the power spectrum at `hist_height` relative height. The histogram is computed in the third quarter of the linearly detrended power spectrum. Parameters  psd_data: 1D array The power spectrum from which to estimate the thresholds. low_thresh_factor: float Factor by which the estimated standard deviation of the noise floor is multiplied to set the `low_threshold`. high_thresh_factor: float Factor by which the estimated standard deviation of the noise floor is multiplied to set the `high_threshold`. nbins: int or list of floats Number of bins or the bins for computing the histogram. hist_height: float Height between 0 and 1 at which the standard deviation of the histogram is estimated. Returns  low_threshold: float The threshold for peaks just above the noise floor. high_threshold: float The threshold for distinct peaks. center: float The baseline level of the power spectrum. """ n = len(psd_data) psd_data_seg = psd_data[n//2:n*3//4] psd_data_seg = psd_data_seg[~np.isinf(psd_data_seg)] psd_data_seg = np.mean(psd_data_seg) + \ sig.detrend(psd_data_seg, type='linear') noise_std, center = hist_threshold(psd_data_seg, thresh_fac=1.0, nbins=nbins) low_threshold = noise_std * low_thresh_factor high_threshold = noise_std * high_thresh_factor return low_threshold, high_threshold, center
def unique(freqs, df_thresh, mode='power', nextfs=0)

Remove similar frequencies from different recordings.
If two frequencies from different elements in the inner lists of
freqs
are reciprocally the closest to each other and closer thandf_thresh
, then the one with the smaller power is removed. As power, either the absolute power as provided in the second column of the data elements infreqs
is taken (mode=='power'), or the relative power (mode='relpower'), or the power rank (mode='rank').Parameters
freqs
:(list
of(list
of…)) list
of2D arrays
 First column in the arrays is fundamental frequency and second column the corresponding power, as returned by fundamental_freqs_and_power().
df_thresh
:float
 Fundamental frequencies closer than this threshold are considered equal.
mode
:string

 'power': use second column of freqs elements as power.
 'relpower': use relative power computed from the second column of freqs elements for deciding which frequency to delete.
 'rank': use rank of second column of freqs elements for deciding which frequency to delete.
nextfs
:int
 If zero, compare all elements in freqs with each other. Otherwise,
only compare with the
nextfs
next elements in freqs.
Returns
uniqe_freqs
:(list
of(list
of…)) list
of2D arrays
 Same as
freqs
but elements with similar fundamental frequencies removed.
Expand source code
def unique(freqs, df_thresh, mode='power', nextfs=0): """Remove similar frequencies from different recordings. If two frequencies from different elements in the inner lists of `freqs` are reciprocally the closest to each other and closer than `df_thresh`, then the one with the smaller power is removed. As power, either the absolute power as provided in the second column of the data elements in `freqs` is taken (mode=='power'), or the relative power (mode='relpower'), or the power rank (mode='rank'). Parameters  freqs: (list of (list of ...)) list of 2D arrays First column in the arrays is fundamental frequency and second column the corresponding power, as returned by fundamental_freqs_and_power(). df_thresh: float Fundamental frequencies closer than this threshold are considered equal. mode: string  'power': use second column of freqs elements as power.  'relpower': use relative power computed from the second column of freqs elements for deciding which frequency to delete.  'rank': use rank of second column of freqs elements for deciding which frequency to delete. nextfs: int If zero, compare all elements in freqs with each other. Otherwise, only compare with the `nextfs` next elements in freqs. Returns  uniqe_freqs: (list of (list of ...)) list of 2D arrays Same as `freqs` but elements with similar fundamental frequencies removed. """ if len(freqs) == 0: return [] # check whether freqs is list of fundamental frequencies and powers: list_of_freq_power = True for group in freqs: if not (hasattr(group, 'shape') and len(group.shape) == 2): list_of_freq_power = False break if list_of_freq_power: if mode == 'power': mask = unique_mask(freqs, df_thresh, nextfs) elif mode == 'relpower': power_freqs = [f[:,[0, 2, 1]] for f in add_relative_power(freqs)] mask = unique_mask(power_freqs, df_thresh, nextfs) elif mode == 'rank': rank_freqs = [f[:,[0, 2, 1]] for f in add_power_ranks(freqs)] mask = unique_mask(rank_freqs, df_thresh, nextfs) else: raise ValueError('%s is not a valid mode for unique(). Choose one of "power" or "rank"') unique_freqs = [] for f, m in zip(freqs, mask): unique_freqs.append(f[m]) else: unique_freqs = [] for groups in freqs: unique_freqs.append(unique(groups, df_thresh, mode, nextfs)) return unique_freqs
def unique_mask(freqs, df_thresh, nextfs=0)

Mark similar frequencies from different recordings as dublicate.
If two frequencies from different elements in
freqs
are reciprocally the closest to each other and closer thandf_thresh
, then the one with the smaller power is marked for removal.Parameters
freqs
:list
of2D arrays
 First column in the arrays is fundamental frequency and second column the corresponding power or equivalent. If values in the second column are equal (e.g. they are the same ranks), and there is a third column (e.g. power), the third column is used to decide, which element should be removed.
df_thresh
:float
 Fundamental frequencies closer than this threshold are considered equal.
nextfs
:int
 If zero, compare all elements in freqs with each other. Otherwise,
only compare with the
nextfs
next elements in freqs.
Returns
mask
:list
ofboolean arrays
 For each element in
freqs
True if that frequency should be kept.
Expand source code
def unique_mask(freqs, df_thresh, nextfs=0): """Mark similar frequencies from different recordings as dublicate. If two frequencies from different elements in `freqs` are reciprocally the closest to each other and closer than `df_thresh`, then the one with the smaller power is marked for removal. Parameters  freqs: list of 2D arrays First column in the arrays is fundamental frequency and second column the corresponding power or equivalent. If values in the second column are equal (e.g. they are the same ranks), and there is a third column (e.g. power), the third column is used to decide, which element should be removed. df_thresh: float Fundamental frequencies closer than this threshold are considered equal. nextfs: int If zero, compare all elements in freqs with each other. Otherwise, only compare with the `nextfs` next elements in freqs. Returns  mask: list of boolean arrays For each element in `freqs` True if that frequency should be kept. """ mask = [np.ones(len(freqs[i]), dtype=bool) for i in range(len(freqs))] for j in range(len(freqs)1): freqsj = np.asarray(freqs[j]) for m in range(len(freqsj)): freq1 = freqsj[m] nn = len(freqs) if nextfs == 0 else j+1+nextfs if nn > len(freqs): nn = len(freqs) for k in range(j+1, nn): freqsk = np.asarray(freqs[k]) if len(freqsk) == 0: continue n = np.argmin(np.abs(freqsk[:,0]  freq1[0])) freq2 = freqsk[n] if np.argmin(np.abs(freqsj[:,0]  freq2[0])) != m: continue if np.abs(freq1[0]  freq2[0]) < df_thresh: if freq1[1] > freq2[1]: mask[k][n] = False elif freq1[1] < freq2[1]: mask[j][m] = False elif len(freq1) > 2: if freq1[2] > freq2[2]: mask[k][n] = False else: mask[j][m] = False else: mask[j][m] = False return mask
def update_group(good_freqs, all_freqs, new_group, fzero, freq_tol, verbose, group_str)

Update good frequencies and harmonic group.
Remove frequencies from good_freqs, add missing fundamental to group.
Parameters
good_freqs
:2D array
 Frequency, power, and use count (columns) of strong peaks detected in a power spectrum.
all_freqs
:2D array
 Frequency, power, and use count (columns) of all peaks detected in a power spectrum.
new_group
:1D aray
ofindices
 Frequencies in all_freqs of an harmonic group.
fzero
:float
 Fundamental frequency for which frequencies are collected in good_freqs.
freq_tol
:float
 Harmonics need to fall within this frequency tolerance. This should be in the range of the frequency resolution and not be smaller than half of the frequency resolution.
verbose
:int
 Verbosity level.
group_str
:string
 String for debug message.
Returns
good_freqs
:2D array
 Frequency, power, and use count (columns) of strong peaks detected in a power spectrum with frequencies for harmonic group of fundamental frequency fzero removed.
group
:2D array
 Frequency, power, and use count (columns) of harmonic group for fundamental frequency fzero.
Expand source code
def update_group(good_freqs, all_freqs, new_group, fzero, freq_tol, verbose, group_str): """Update good frequencies and harmonic group. Remove frequencies from good_freqs, add missing fundamental to group. Parameters  good_freqs: 2D array Frequency, power, and use count (columns) of strong peaks detected in a power spectrum. all_freqs: 2D array Frequency, power, and use count (columns) of all peaks detected in a power spectrum. new_group: 1D aray of indices Frequencies in all_freqs of an harmonic group. fzero: float Fundamental frequency for which frequencies are collected in good_freqs. freq_tol: float Harmonics need to fall within this frequency tolerance. This should be in the range of the frequency resolution and not be smaller than half of the frequency resolution. verbose: int Verbosity level. group_str: string String for debug message. Returns  good_freqs: 2D array Frequency, power, and use count (columns) of strong peaks detected in a power spectrum with frequencies for harmonic group of fundamental frequency fzero removed. group: 2D array Frequency, power, and use count (columns) of harmonic group for fundamental frequency fzero. """ # initialize group: group = all_freqs[new_group,:] # indices of group in good_freqs: freq_tol *= 1.1 indices = [] for f in group[:,0]: idx = np.argmin(np.abs(good_freqs[:,0]f)) if np.abs(good_freqs[idx,0]f) <= freq_tol: indices.append(idx) indices = np.asarray(indices, dtype=np.int) # harmonics in good_freqs: nharm = np.round(good_freqs[:,0]/fzero) idxs = np.where(np.abs(good_freqs[:,0]  nharm*fzero) <= freq_tol)[0] indices = np.unique(np.concatenate((indices, idxs))) # report: if verbose > 1: print('# good freqs: ', indices, '[', ', '.join(['%.2f' % f for f in good_freqs[indices,0]]), ']') print('# all freqs : ', new_group, '[', ', '.join(['%.2f' % f for f in all_freqs[new_group,0]]), ']') if verbose > 0: refi = np.argmax(group[:,1] > 0.0) print('') print(group_str) for i in range(len(group)): print('f=%8.2fHz n=%5.2f: power=%9.3g power/pmax=%6.4f=%5.1fdB' % (group[i,0], group[i,0]/fzero, group[i,1], group[i,1]/group[refi,1], decibel(group[i,1], group[refi,1]))) print('') # erase group from good_freqs: good_freqs = np.delete(good_freqs, indices, axis=0) # adjust frequencies to fzero: group[:,0] = np.round(group[:,0]/fzero)*fzero # insert missing fzero: if np.round(group[0,0]/fzero) != 1.0: group = np.vstack(((fzero, group[0,1], 2.0), group)) return good_freqs, group