Coverage for src/thunderfish/harmonics.py: 84%
919 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-19 20:05 +0000
1"""
2Extract and analyze harmonic frequencies from power spectra.
4## Harmonic group extraction
6- `harmonic_groups()`: detect peaks in power spectrum and group them
7 according to their harmonic structure.
8- `expand_group()`: add more harmonics to harmonic group.
9- `extract_fundamentals()`: collect harmonic groups from
10 lists of power spectrum peaks.
11- `threshold_estimate()`: estimates thresholds for peak detection
12 in a power spectrum.
14## Helper functions for harmonic group extraction
16- `build_harmonic_group()`: find all the harmonics belonging to the largest peak in the power spectrum.
17- `retrieve_harmonic_group()`: find all the harmonics belonging to a given fundamental.
18- `group_candidate()`: candidate harmonic frequencies belonging to a fundamental frequency.
19- `update_group()`: update frequency lists and harmonic group.
21## Handling of lists of harmonic groups
23- `fundamental_freqs()`: extract fundamental frequencies from
24 lists of harmonic groups.
25- `fundamental_freqs_and_power()`: extract fundamental frequencies and their
26 power in decibel from lists of harmonic groups.
28## Handling of lists of fundamental frequencies
30- `add_relative_power()`: add a column with relative power.
31- `add_power_ranks()`: add a column with power ranks.
32- `similar_indices()`: indices of similar frequencies.
33- `unique_indices()`: for each unique fundamental frequency a list of indices indicating where it is found.
34- `collapse()`: collapse harmonic groups with similar fundamental frequencies to a single one.
35- `consistent()`: harmonic groups whose fundamental frequencies consistently appear everywhere.
36- `closest()`: harmonic groups whose fundamental frequencies are closest to each other.
37- `unique_mask()`: mark similar frequencies from different recordings as dublicate.
38- `unique()`: remove similar frequencies from different recordings.
40## Visualization
42- `colors_markers()`: Generate a list of colors and markers for plotting.
43- `plot_harmonic_groups()`: Mark decibel power of fundamentals and their
44 harmonics.
45- `plot_psd_harmonic_groups()`: Plot decibel power-spectrum with detected peaks,
46 harmonic groups, and mains frequencies.
47- `annotate_harmonic_group()`: annotate peaks of a harmonic group in a spectrum.
48- `plot_unique_frequencies()`: plot all fundamental frequencies and mark unique frequencies.
49- `plot_selected_groups()`: plot all fundamental frequencies and mark selected frequencies with a bar.
51## Configuration
53- `add_psd_peak_detection_config()`: add parameters for the detection of
54 peaks in power spectra to configuration.
55- `psd_peak_detection_args()`: retrieve parameters for the detection of peaks
56 in power spectra from configuration.
57- `add_harmonic_groups_config()`: add parameters for the detection of
58 harmonic groups to configuration.
59- `harmonic_groups_args()`: retrieve parameters for the detection of
60 harmonic groups from configuration.
61"""
63import math as m
64import numpy as np
65import scipy.signal as sig
66try:
67 import matplotlib.pyplot as plt
68 import matplotlib.cm as cm
69 import matplotlib.colors as mc
70except ImportError:
71 pass
73from thunderlab.eventdetection import detect_peaks, trim, hist_threshold
74from thunderlab.powerspectrum import decibel, power, plot_decibel_psd
77def group_candidate(good_freqs, all_freqs, freq, divisor,
78 freq_tol, max_freq_tol, min_group_size, verbose):
79 """Candidate harmonic frequencies belonging to a fundamental frequency.
81 Parameters
82 ----------
83 good_freqs: 2-D array
84 Frequency, power, and use count (columns) of strong peaks detected
85 in a power spectrum.
86 all_freqs: 2-D array
87 Frequency, power, and use count (columns) of all peaks detected
88 in a power spectrum.
89 freq: float
90 Fundamental frequency for which a harmonic group should be assembled.
91 divisor: int
92 Fundamental frequency was obtained from a frequency divided by divisor.
93 0: Fundamental frequency is given as is.
94 freq_tol: float
95 Harmonics should fall within this frequency tolerance.
96 This should be in the range of the frequency resolution
97 and should not be smaller than half of the frequency resolution.
98 max_freq_tol: float
99 Maximum deviation of harmonics from their expected frequency.
100 Peaks with frequencies between `freq_tol` and `max_freq_tol`
101 get penalized the further away from their expected frequency.
102 min_group_size: int
103 Minimum number of harmonics of a harmonic group.
104 The harmonics from min_group_size/3 to max(min_group_size, divisor)
105 need to be in good_freqs.
106 verbose: int
107 Verbosity level.
109 Returns
110 -------
111 new_group: 1-D aray of indices
112 Frequencies in all_freqs belonging to the candidate group.
113 fzero: float
114 Adjusted fundamental frequency. If negative, no group was found.
115 fzero_harmonics: int
116 The highest harmonics that was used to recompute the fundamental frequency.
117 """
118 if verbose > 0:
119 print('')
121 dvs = True
122 if divisor <= 0:
123 divisor = 1
124 dvs = False
126 # 1. find harmonics in good_freqs and adjust fzero accordingly:
127 fzero = freq
128 fzero_harmonics = 1
129 if len(good_freqs[:,0]) > 0:
130 """
131 ff = good_freqs[:,0]
132 h = np.round(ff/fzero)
133 h_indices = np.unique(h, return_index=True)[1]
134 f_best = [hi+np.argmin(np.abs(fh/h - fzero)) for hi, fh in zip(h_indices, np.split(ff, h_indices))]
135 fe = np.abs(ff/h - fzero)
136 h = h[fe < freq_tol]
137 fe = fe[fe < freq_tol]
138 """
139 prev_freq = divisor * freq
140 for h in range(divisor+1, 2*min_group_size+1):
141 idx = np.argmin(np.abs(good_freqs[:,0]/h - fzero))
142 ff = good_freqs[idx,0]
143 if m.fabs(ff/h - fzero) > freq_tol:
144 continue
145 df = ff - prev_freq
146 dh = np.round(df/fzero)
147 fe = m.fabs(df/dh - fzero)
148 if fe > 2.0*freq_tol:
149 if h > min_group_size:
150 break
151 continue
152 # update fzero:
153 prev_freq = ff
154 fzero_harmonics = h
155 fzero = ff/fzero_harmonics
156 if verbose > 1:
157 print('adjusted fzero to %.2fHz from %d-th harmonic' % (fzero, h))
159 if verbose > 0:
160 ds = 'divisor: %d, ' % divisor if dvs else ''
161 print('# %sfzero=%7.2fHz adjusted from harmonics %d'
162 % (ds, fzero, fzero_harmonics))
164 # 2. check fzero:
165 # freq might not be in our group anymore, because fzero was adjusted:
166 if np.abs(freq - fzero) > freq_tol:
167 if verbose > 0:
168 print(' discarded: lost frequency')
169 return [], -1.0, fzero_harmonics
171 # 3. collect harmonics from all_freqs:
172 new_group = -np.ones(min_group_size, dtype=int)
173 new_penalties = np.ones(min_group_size)
174 freqs = []
175 prev_h = 0
176 prev_fe = 0.0
177 for h in range(1, min_group_size+1):
178 penalty = 0
179 i = np.argmin(np.abs(all_freqs[:,0]/h - fzero))
180 f = all_freqs[i,0]
181 if verbose > 2:
182 print(' check %7.2fHz as %d. harmonics' % (f, h))
183 fac = 1.0 if h >= divisor else 2.0
184 fe = m.fabs(f/h - fzero)
185 if fe > fac*max_freq_tol:
186 if verbose > 1 and fe < 2*fac*max_freq_tol:
187 print(' %d. harmonics at %7.2fHz is off by %7.2fHz (max between %5.2fHz and %5.2fHz) from %7.2fHz'
188 % (h, f, h*fe, h*fac*freq_tol, h*fac*max_freq_tol, h*fzero))
189 continue
190 if fe > fac*freq_tol:
191 penalty = np.interp(fe, [fac*freq_tol, fac*max_freq_tol], [0.0, 1.0])
192 if len(freqs) > 0:
193 pf = freqs[-1]
194 df = f - pf
195 if df < 0.5*fzero:
196 if len(freqs)>1:
197 pf = freqs[-2]
198 df = f - pf
199 else:
200 pf = 0.0
201 df = h*fzero
202 dh = m.floor(df/fzero + 0.5)
203 fe = m.fabs(df/dh - fzero)
204 if fe > 2*dh*fac*max_freq_tol:
205 if verbose > 1:
206 print(' %d. harmonics at %7.2fHz is off by %7.2fHz (max between %5.2fHz and %5.2fHz) from previous harmonics at %7.2fHz'
207 % (h, f, dh*fe, 2*fac*dh*freq_tol, 2*fac*dh*max_freq_tol, pf))
208 continue
209 if fe > 2*dh*fac*freq_tol:
210 penalty = np.interp(fe, [2*dh*fac*freq_tol, 2*dh*fac*max_freq_tol], [0.0, 1.0])
211 else:
212 fe = 0.0
213 if h > prev_h or fe < prev_fe:
214 if prev_h > 0 and h - prev_h > 1:
215 if verbose > 1:
216 print(' previous harmonics %d more than 1 away from %d. harmonics at %7.2fHz'
217 % (prev_h, h, f))
218 break
219 if h == prev_h and len(freqs) > 0:
220 freqs.pop()
221 freqs.append(f)
222 new_group[int(h)-1] = i
223 new_penalties[int(h)-1] = penalty
224 prev_h = h
225 prev_fe = fe
226 if verbose > 1:
227 print(' %d. harmonics at %7.2fHz has been taken (peak %2d) with penalty %3.1f' % (h, f, i, penalty))
229 # 4. check new group:
231 # almost all harmonics in min_group_size required:
232 max_penalties = min_group_size/3
233 if np.sum(new_penalties) > max_penalties:
234 if verbose > 0:
235 print(' discarded group because sum of penalties %3.1f is more than %3.1f: indices' %
236 (np.sum(new_penalties), max_penalties), new_group, ' penalties', new_penalties)
237 return [], -1.0, fzero_harmonics
238 new_group = new_group[new_group>=0]
240 # check use count of frequencies:
241 double_use = np.sum(all_freqs[new_group, 2]>0)
242 if double_use >= 2: # XXX make this a parameter?
243 if verbose > 0:
244 print(' discarded group because of use count = %2d >= 2' % double_use)
245 return [], -1.0, fzero_harmonics
247 # 5. return results:
248 return new_group, fzero, fzero_harmonics
251def update_group(good_freqs, all_freqs, new_group, fzero,
252 freq_tol, verbose, group_str):
253 """Update good frequencies and harmonic group.
255 Remove frequencies from good_freqs, add missing fundamental to group.
257 Parameters
258 ----------
259 good_freqs: 2-D array
260 Frequency, power, and use count (columns) of strong peaks detected
261 in a power spectrum.
262 all_freqs: 2-D array
263 Frequency, power, and use count (columns) of all peaks detected
264 in a power spectrum.
265 new_group: 1-D aray of indices
266 Frequencies in all_freqs of an harmonic group.
267 fzero: float
268 Fundamental frequency for which frequencies are collected in good_freqs.
269 freq_tol: float
270 Harmonics need to fall within this frequency tolerance.
271 This should be in the range of the frequency resolution
272 and not be smaller than half of the frequency resolution.
273 verbose: int
274 Verbosity level.
275 group_str: string
276 String for debug message.
278 Returns
279 -------
280 good_freqs: 2-D array
281 Frequency, power, and use count (columns) of strong peaks detected
282 in a power spectrum with frequencies for harmonic group
283 of fundamental frequency fzero removed.
284 group: 2-D array
285 Frequency, power, and use count (columns) of harmonic group
286 for fundamental frequency fzero.
287 """
288 # initialize group:
289 group = all_freqs[new_group,:]
291 # indices of group in good_freqs:
292 freq_tol *= 1.1
293 indices = []
294 for f in group[:,0]:
295 idx = np.argmin(np.abs(good_freqs[:,0]-f))
296 if np.abs(good_freqs[idx,0]-f) <= freq_tol:
297 indices.append(idx)
298 indices = np.asarray(indices, dtype=int)
300 # harmonics in good_freqs:
301 nharm = np.round(good_freqs[:,0]/fzero)
302 idxs = np.where(np.abs(good_freqs[:,0] - nharm*fzero) <= freq_tol)[0]
303 indices = np.unique(np.concatenate((indices, idxs)))
305 # report:
306 if verbose > 1:
307 print('# good freqs: ', indices,
308 '[', ', '.join(['%.2f' % f for f in good_freqs[indices,0]]), ']')
309 print('# all freqs : ', new_group,
310 '[', ', '.join(['%.2f' % f for f in all_freqs[new_group,0]]), ']')
311 if verbose > 0:
312 refi = np.argmax(group[:,1] > 0.0)
313 print('')
314 print(group_str)
315 for i in range(len(group)):
316 print('f=%8.2fHz n=%5.2f: power=%9.3g power/pmax=%6.4f=%5.1fdB'
317 % (group[i,0], group[i,0]/fzero,
318 group[i,1], group[i,1]/group[refi,1], decibel(group[i,1], group[refi,1])))
319 print('')
321 # erase group from good_freqs:
322 good_freqs = np.delete(good_freqs, indices, axis=0)
324 # adjust frequencies to fzero:
325 group[:,0] = np.round(group[:,0]/fzero)*fzero
327 # insert missing fzero:
328 if np.round(group[0,0]/fzero) != 1.0:
329 group = np.vstack(((fzero, group[0,1], -2.0), group))
331 return good_freqs, group
334def build_harmonic_group(good_freqs, all_freqs, freq_tol, max_freq_tol,
335 verbose=0, min_group_size=3, max_divisor=4):
336 """Find all the harmonics belonging to the largest peak in a list of frequency peaks.
338 Parameters
339 ----------
340 good_freqs: 2-D array
341 Frequency, power, and use count (columns) of strong peaks detected
342 in a power spectrum.
343 all_freqs: 2-D array
344 Frequency, power, and use count (columns) of all peaks detected
345 in a power spectrum.
346 freq_tol: float
347 Harmonics should fall within this frequency tolerance.
348 This should be in the range of the frequency resolution
349 and should not be smaller than half of the frequency resolution.
350 max_freq_tol: float
351 Maximum deviation of harmonics from their expected frequency.
352 Peaks with frequencies between `freq_tol` and `max_freq_tol`
353 get penalized the further away from their expected frequency.
354 verbose: int
355 Verbosity level.
356 min_group_size: int
357 Minimum number of harmonics of a harmonic group.
358 The harmonics from min_group_size/3 to max(min_group_size, divisor)
359 need to be in good_freqs.
360 max_divisor: int
361 Maximum divisor used for checking for sub-harmonics.
363 Returns
364 -------
365 good_freqs: 2-D array
366 Frequency, power, and use count (columns) of strong peaks detected
367 in a power spectrum with frequencies of the returned harmonic group removed.
368 group: 2-D array
369 The detected harmonic group. Might be empty.
370 indices: 1-D array of indices
371 Indices of the harmonic group in all_freqs.
372 best_fzero_harmonics: int
373 The highest harmonics that was used to recompute
374 the fundamental frequency.
375 fmax: float
376 The frequency of the largest peak in good_freqs
377 for which the harmonic group was detected.
378 """
379 # select strongest frequency for building the harmonic group:
380 fmaxinx = np.argmax(good_freqs[:,1])
381 fmax = good_freqs[fmaxinx,0]
382 if verbose > 0:
383 print('')
384 print('%s build_harmonic_group %s' % (10*'#', 38*'#'))
385 print('%s fmax=%7.2fHz, power=%9.3g %s'
386 % (10*'#', fmax, good_freqs[fmaxinx,1], 27*'#'))
387 print('good_freqs: ', '[',
388 ', '.join(['%.2f' % f for f in good_freqs[:,0]]), ']')
390 # container for harmonic groups:
391 best_group = []
392 best_value = -1e6
393 best_divisor = 0
394 best_fzero = 0.0
395 best_fzero_harmonics = 0
397 # check for integer fractions of the frequency:
398 for divisor in range(1, max_divisor + 1):
399 # 1. hypothesized fundamental:
400 freq = fmax / divisor
402 # 2. find harmonics in good_freqs and adjust fzero accordingly:
403 group_size = min_group_size if divisor <= min_group_size else divisor
404 new_group, fzero, fzero_harmonics = group_candidate(good_freqs, all_freqs,
405 freq, divisor,
406 freq_tol, max_freq_tol,
407 group_size, verbose)
408 # no group found:
409 if fzero < 0.0:
410 continue
412 # 3. compare new group to best group:
413 peaksum = decibel(np.sum(all_freqs[new_group, 1])*min_group_size/len(new_group))
414 diff = np.std(np.diff(decibel(all_freqs[new_group, 1])))
415 new_group_value = peaksum - diff
416 counts = int(np.sum(all_freqs[new_group, 2]))
417 if verbose > 0:
418 print(f' new group: fzero={fzero:7.2f}Hz, nharmonics={len(new_group)}, value={new_group_value:6.1f}dB, peaksum={peaksum:5.1f}dB, diff={diff:6.1f}dB, use count={counts:2d}, peaks:', new_group)
419 if verbose > 1 and len(best_group) > 0:
420 print(f' best group: divisor={best_divisor}, fzero={best_fzero:7.2f}Hz, nharmonics={len(best_group)}, value={best_value:6.1f}dB, peaks:', best_group)
421 # select new group if sum of peak power minus diff is larger:
422 if (len(new_group) >= min_group_size and len(best_group) < min_group_size) or \
423 (len(new_group) >= len(best_group) and new_group_value >= best_value):
424 best_value = new_group_value
425 best_group = new_group
426 best_divisor = divisor
427 best_fzero = fzero
428 best_fzero_harmonics = fzero_harmonics
429 if verbose > 1:
430 print(f' new best group: divisor={best_divisor}, fzero={best_fzero:7.2f}Hz, nharmonics={len(best_group)}, value={best_value:6.1f}dB, peaks:', best_group)
431 elif verbose > 0:
432 print(' took as new best group')
433 elif verbose > 0:
434 if not (len(new_group) >= min_group_size and len(best_group) < min_group_size):
435 s = f', number of harmonics {len(new_group)} not greater than minimum group size {min_group_size}'
436 if len(new_group) < len(best_group):
437 print(f' keep best group: number of harmonics {len(new_group)} is SMALLER than in best group {len(best_group)}{s}')
438 elif new_group_value < best_value:
439 print(f' keep best group: value {new_group_value:6.1f}dB is SMALLER than in best group {best_value:6.1f}dB{s}')
440 else:
441 print(f' keep best group: {s[2:]}')
443 # no group found:
444 if len(best_group) == 0:
445 # erase freq:
446 good_freqs = np.delete(good_freqs, fmaxinx, axis=0)
447 group = np.zeros((0, 3))
448 return good_freqs, group, best_group, 1, fmax
450 # update frequencies and group:
451 if verbose > 1:
452 print('')
453 print('# best group found for fmax=%.2fHz, fzero=%.2fHz, divisor=%d:'
454 % (fmax, best_fzero, best_divisor))
455 group_str = '%s resulting harmonic group for fmax=%.2fHz' % (10*'#', fmax)
456 good_freqs, group = update_group(good_freqs, all_freqs, best_group, best_fzero, freq_tol, verbose, group_str)
458 # good_freqs: removed all frequencies of bestgroup
459 return good_freqs, group, best_group, best_fzero_harmonics, fmax
462def retrieve_harmonic_group(freq, good_freqs, all_freqs,
463 freq_tol, max_freq_tol, verbose=0,
464 min_group_size=3):
465 """Find all the harmonics belonging to a given fundamental.
467 Parameters
468 ----------
469 freq: float
470 Fundamental frequency for which harmonics are to be retrieved.
471 good_freqs: 2-D array
472 Frequency, power, and use count (columns) of strong peaks detected
473 in a power spectrum. All harmonics of `freq` will be
474 removed from `good_freqs`.
475 all_freqs: 2-D array
476 Frequency, power, and use count (columns) of all peaks detected
477 in a power spectrum.
478 freq_tol: float
479 Harmonics should fall within this frequency tolerance.
480 This should be in the range of the frequency resolution
481 and should not be smaller than half of the frequency resolution.
482 max_freq_tol: float
483 Maximum deviation of harmonics from their expected frequency.
484 Peaks with frequencies between `freq_tol` and `max_freq_tol`
485 get penalized the further away from their expected frequency.
486 verbose: int
487 Verbosity level.
488 min_group_size: int
489 Minimum number of harmonics of a harmonic group.
490 The harmonics from min_group_size/3 to max(min_group_size, divisor)
491 need to be in good_freqs.
493 Returns
494 -------
495 good_freqs: 2-D array
496 Frequency, power, and use count (columns) of strong peaks detected
497 in a power spectrum with frequencies of the returned harmonic group removed.
498 group: 2-D array
499 The detected harmonic group. Might be empty.
500 indices: 1-D array of indices
501 Indices of the harmonic group in all_freqs.
502 fzero_harmonics: int
503 The highest harmonics that was used to recompute
504 the fundamental frequency.
505 """
506 if verbose > 0:
507 print('')
508 print('%s retrieve harmonic group %s' % (10*'#', 35*'#'))
509 print('%s freq=%7.2fHz %s' % (10*'#', freq, 44*'#'))
510 print('good_freqs: ', '[',
511 ', '.join(['%.2f' % f for f in good_freqs[:,0]]), ']')
513 # find harmonics in good_freqs and adjust fzero accordingly:
514 new_group, fzero, fzero_harmonics = group_candidate(good_freqs, all_freqs,
515 freq, 0,
516 freq_tol, max_freq_tol,
517 min_group_size, verbose)
519 # no group found:
520 if fzero < 0.0:
521 return good_freqs, np.zeros((0, 2)), np.zeros(0), fzero_harmonics
523 if verbose > 1:
524 print('')
525 print('# group found for freq=%.2fHz, fzero=%.2fHz:'
526 % (freq, fzero))
527 group_str = '#### resulting harmonic group for freq=%.2fHz' % freq
528 good_freqs, group = update_group(good_freqs, all_freqs, new_group,
529 fzero, freq_tol, verbose, group_str)
531 # good_freqs: removed all frequencies of bestgroup
532 return good_freqs, group, new_group, fzero_harmonics
535def expand_group(group, indices, freqs, freq_tol, max_harmonics=0):
536 """Add more harmonics to harmonic group.
538 Parameters
539 ----------
540 group: 2-D array
541 Group of fundamental frequency and harmonics
542 as returned by build_harmonic_group.
543 indices: 1-D array of indices
544 Indices of the harmonics in group in `freqs`.
545 freqs: 2-D array
546 Frequency, power, and use count (columns) of all peaks detected
547 in a power spectrum.
548 freq_tol: float
549 Harmonics need to fall within this frequency tolerance.
550 This should be in the range of the frequency resolution
551 and not be smaller than half of the frequency resolution.
552 max_harmonics: int
553 Maximum number of harmonics to be returned for each group.
555 Returns
556 -------
557 group: 2-D array
558 Expanded group of fundamental frequency and harmonics.
559 indices: 1-D array of indices
560 Indices of the harmonics in the expanded group in all_freqs.
561 """
562 if len(group) == 0:
563 return group, indices
564 fzero = group[0,0]
565 if max_harmonics <= 0:
566 max_harmonics = m.floor(freqs[-1,0]/fzero + 0.5) # round
567 if max_harmonics <= len(group):
568 return group, indices
569 group_freqs = list(group[:,0])
570 indices = list(indices)
571 last_h = m.floor(group_freqs[-1]/fzero + 0.5) # round
572 for h in range(last_h+1, max_harmonics+1):
573 i = np.argmin(np.abs(freqs[:,0]/h - fzero))
574 f = freqs[i,0]
575 if m.fabs(f/h - fzero) > freq_tol:
576 continue
577 df = f - group_freqs[-1]
578 dh = m.floor(df/fzero + 0.5) # round
579 if dh < 1:
580 continue
581 fe = m.fabs(df/dh - fzero)
582 if fe > 2*freq_tol:
583 continue
584 group_freqs.append(f)
585 indices.append(i)
586 # assemble group:
587 new_group = freqs[indices,:group.shape[1]]
588 # keep filled in fundamental:
589 if group[0,2] == -2:
590 new_group = np.vstack((group[0,:], new_group))
591 return new_group, np.array(indices, dtype=int)
594def extract_fundamentals(good_freqs, all_freqs, freq_tol, max_freq_tol,
595 verbose=0, check_freqs=[],
596 mains_freq=60.0, mains_freq_tol=1.0,
597 min_freq=0.0, max_freq=2000.0, max_db_diff=20.0,
598 max_divisor=4, min_group_size=3, max_harmonics_db=-5.0,
599 max_harmonics=0, max_groups=0, **kwargs):
600 """Extract fundamental frequencies from power-spectrum peaks.
602 Parameters
603 ----------
604 good_freqs: 2-D array
605 Frequency, power, and use count (columns) of strong peaks detected
606 in a power spectrum.
607 all_freqs: 2-D array
608 Frequency, power, and use count (columns) of all peaks detected
609 in a power spectrum.
610 freq_tol: float
611 Harmonics should fall within this frequency tolerance.
612 This should be in the range of the frequency resolution
613 and should not be smaller than half of the frequency resolution.
614 max_freq_tol: float
615 Maximum deviation of harmonics from their expected frequency.
616 Peaks with frequencies between `freq_tol` and `max_freq_tol`
617 get penalized the further away from their expected frequency.
618 verbose: int
619 Verbosity level.
620 check_freqs: list of float
621 List of fundamental frequencies that will be checked
622 first for being present and valid harmonic groups in the peak frequencies
623 of a power spectrum.
624 mains_freq: float
625 Frequency of the mains power supply.
626 mains_freq_tol: float
627 Tolarerance around harmonics of the mains frequency,
628 within which peaks are removed.
629 min_freq: float
630 Minimum frequency accepted as a fundamental frequency.
631 max_freq: float
632 Maximum frequency accepted as a fundamental frequency.
633 max_db_diff: float
634 If larger than zero, maximum standard deviation of differences between
635 logarithmic powers of harmonics in decibel.
636 Low values enforce smoother power spectra.
637 max_divisor: int
638 Maximum divisor used for checking for sub-harmonics.
639 min_group_size: int
640 Minimum number of harmonics of a harmonic group.
641 The harmonics from min_group_size/3 to max(min_group_size, divisor)
642 need to be in good_freqs.
643 max_harmonics_db: float
644 Maximum allowed power of the `min_group_size`-th and higher harmonics
645 after the peak (in decibel relative to peak power withn the first
646 `min_group_size` harmonics, i.e. if harmonics are required to be
647 smaller than fundamental then this is a negative number).
648 Make it a large positive number to effectively not check for relative power.
649 max_harmonics: int
650 Maximum number of harmonics to be returned for each group.
651 max_groups: int
652 If not zero the maximum number of most powerful harmonic groups.
654 Returns
655 -------
656 group_list: list of 2-D arrays
657 List of all harmonic groups found sorted by fundamental frequency.
658 Each harmonic group is a 2-D array with the first dimension the harmonics
659 and the second dimension containing frequency and power of each harmonic.
660 If the power is zero, there was no corresponding peak in the power spectrum.
661 fzero_harmonics_list: list of int
662 The harmonics from which the fundamental frequencies were computed.
663 mains_freqs: 2-d array
664 Array of mains peaks found in all_freqs (frequency, power).
665 """
666 if verbose > 0:
667 print('')
669 # set use count to zero:
670 all_freqs[:,2] = 0.0
672 # remove power line harmonics from good_freqs:
673 if mains_freq > 0.0:
674 indices = np.where(np.abs(good_freqs[:,0] - np.round(good_freqs[:,0]/mains_freq)*mains_freq) < mains_freq_tol)[0]
675 if len(indices)>0:
676 if verbose > 1:
677 print('remove power line frequencies',
678 ', '.join(['%.1f' % f for f in good_freqs[indices,0]]))
679 good_freqs = np.delete(good_freqs, indices, axis=0)
681 if verbose > 1:
682 print('all_freqs: ', '[',
683 ', '.join(['%.2f' % f for f in all_freqs[:,0] if f < 3000.0]), ']')
685 group_list = []
686 fzero_harmonics_list = []
687 first = True
688 # as long as there are frequencies left in good_freqs:
689 fi = 0
690 while len(good_freqs) > 0:
691 if fi < len(check_freqs):
692 # check for harmonic group of a given fundamental frequency:
693 fmax = check_freqs[fi]
694 f0s = 'freq'
695 good_freqs, harm_group, harm_indices, fzero_harmonics = \
696 retrieve_harmonic_group(fmax, good_freqs, all_freqs,
697 freq_tol, max_freq_tol,
698 verbose-1, min_group_size)
699 fi += 1
700 else:
701 # check for harmonic groups:
702 f0s = 'fmax'
703 good_freqs, harm_group, harm_indices, fzero_harmonics, fmax = \
704 build_harmonic_group(good_freqs, all_freqs,
705 freq_tol, max_freq_tol, verbose-1,
706 min_group_size, max_divisor)
708 # nothing found:
709 if len(harm_group) == 0:
710 if verbose > 1 - first:
711 s = ' (largest peak)' if first else ''
712 print('No harmonic group for %7.2fHz%s' % (fmax, s))
713 first = False
714 continue
715 first = False
717 # fill up harmonic group:
718 harm_group, harm_indices = expand_group(harm_group, harm_indices,
719 all_freqs, freq_tol,
720 max_harmonics)
722 # check whether fundamental was filled in:
723 first_h = 0 if harm_group[0,2] > -2 else 1
725 # increment use count:
726 all_freqs[harm_indices,2] += 1
727 if harm_group[0,2] == -1:
728 harm_group[0,2] = -2
730 # check frequency range of fundamental:
731 fundamental_ok = (harm_group[0, 0] >= min_freq and
732 harm_group[0, 0] <= max_freq)
733 # check power hum:
734 mains_ok = ((mains_freq <= 0.0) or
735 (m.fabs(harm_group[0,0] - mains_freq) > freq_tol))
736 # check smoothness:
737 db_powers = decibel(harm_group[first_h:,1])
738 diff = np.std(np.diff(db_powers))
739 smooth_ok = max_db_diff <= 0.0 or diff < max_db_diff
740 # check relative power of higher harmonics:
741 p_max = np.argmax(db_powers[:min_group_size])
742 db_powers -= db_powers[p_max]
743 amplitude_ok = len(db_powers[p_max+min_group_size:]) == 0
744 if amplitude_ok:
745 pi = len(db_powers) - 1
746 else:
747 amplitude_ok = np.all((db_powers[p_max+min_group_size:] < max_harmonics_db) |
748 (harm_group[first_h+p_max+min_group_size:,2] > 1))
749 if amplitude_ok:
750 pi = p_max+min_group_size-1
751 else:
752 pi = np.where((db_powers[p_max+min_group_size:] >= max_harmonics_db) &
753 (harm_group[first_h+p_max+min_group_size:,2] <= 1))[0][0]
755 # check:
756 if fundamental_ok and mains_ok and smooth_ok and amplitude_ok:
757 if verbose > 0:
758 print('accept harmonic group from %s=%7.2fHz: %7.2fHz power=%9.3g nharmonics=%2d, use count=%d'
759 % (f0s, fmax, harm_group[0,0], np.sum(harm_group[first_h:,1]),
760 len(harm_group), np.sum(harm_group[first_h:,2])))
761 group_list.append(harm_group[:,0:2])
762 fzero_harmonics_list.append(fzero_harmonics)
763 else:
764 if verbose > 0:
765 fs = 'is ' if fundamental_ok else 'NOT'
766 ms = 'not ' if mains_ok else 'IS'
767 ss = 'smaller' if smooth_ok else 'LARGER '
768 ps = 'smaller' if amplitude_ok else 'LARGER '
769 print('discard 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'
770 % (f0s, fmax, harm_group[0,0], np.sum(harm_group[first_h:,1]),
771 fs, ms, diff, ss, max_db_diff,
772 pi, db_powers[pi], ps, max_harmonics_db))
774 # select most powerful harmonic groups:
775 if max_groups > 0 and len(group_list) > max_groups:
776 n = len(group_list)
777 powers = [np.sum(group[:,1]) for group in group_list]
778 powers_inx = np.argsort(powers)
779 group_list = [group_list[pi] for pi in powers_inx[-max_groups:]]
780 fzero_harmonics_list = [fzero_harmonics_list[pi] for pi in powers_inx[-max_groups:]]
781 if verbose > 0:
782 print(f'selected from {n} groups the {max_groups} most powerful groups.')
784 # sort groups by fundamental frequency:
785 freqs = [group[0,0] for group in group_list]
786 freq_inx = np.argsort(freqs)
787 group_list = [group_list[fi] for fi in freq_inx]
788 fzero_harmonics_list = [fzero_harmonics_list[fi] for fi in freq_inx]
790 if verbose > 1:
791 print('')
792 if len(group_list) > 0:
793 print('## FUNDAMENTALS FOUND: ##')
794 for group, fzero_h in zip(group_list, fzero_harmonics_list):
795 print('%7.2fHz: power=%9.3g fzero-h=%2d'
796 % (group[0,0], np.sum(group[:,1]), fzero_h))
797 else:
798 print('## NO FUNDAMENTALS FOUND ##')
800 # assemble mains frequencies from all_freqs:
801 if mains_freq > 0.0:
802 mains_freqs = all_freqs[np.abs(all_freqs[:,0] - np.round(all_freqs[:,0]/mains_freq)*mains_freq) < mains_freq_tol,:2]
803 else:
804 mains_freqs = np.zeros((0, 2))
806 return group_list, fzero_harmonics_list, mains_freqs
809def threshold_estimate(psd_data, low_thresh_factor=6.0, high_thresh_factor=10.0,
810 nbins=100, hist_height=1.0/ np.sqrt(np.e)):
811 """Estimate thresholds for peak detection from histogram of power spectrum.
813 The standard deviation of the noise floor without peaks is estimated from
814 the width of the histogram of the power spectrum at `hist_height` relative height.
815 The histogram is computed in the third quarter of the linearly detrended power spectrum.
817 Parameters
818 ----------
819 psd_data: 1-D array
820 The power spectrum from which to estimate the thresholds.
821 low_thresh_factor: float
822 Factor by which the estimated standard deviation of the noise floor
823 is multiplied to set the `low_threshold`.
824 high_thresh_factor: float
825 Factor by which the estimated standard deviation of the noise floor
826 is multiplied to set the `high_threshold`.
827 nbins: int or list of floats
828 Number of bins or the bins for computing the histogram.
829 hist_height: float
830 Height between 0 and 1 at which the standard deviation of the histogram is estimated.
832 Returns
833 -------
834 low_threshold: float
835 The threshold for peaks just above the noise floor.
836 high_threshold: float
837 The threshold for distinct peaks.
838 center: float
839 The baseline level of the power spectrum.
840 """
841 n = len(psd_data)
842 psd_data_seg = psd_data[n//2:n*3//4]
843 psd_data_seg = psd_data_seg[~np.isinf(psd_data_seg)]
844 psd_data_seg = np.mean(psd_data_seg) + \
845 sig.detrend(psd_data_seg, type='linear')
846 noise_std, center = hist_threshold(psd_data_seg, thresh_fac=1.0, nbins=nbins)
847 low_threshold = noise_std * low_thresh_factor
848 high_threshold = noise_std * high_thresh_factor
849 return low_threshold, high_threshold, center
852def harmonic_groups(psd_freqs, psd, verbose=0, check_freqs=[],
853 low_threshold=0.0, high_threshold=0.0, thresh_bins=100,
854 low_thresh_factor=6.0, high_thresh_factor=10.0,
855 freq_tol_fac=1.0, max_freq_tol=1.0,
856 mains_freq=60.0, mains_freq_tol=1.0,
857 min_freq=0.0, max_freq=2000.0, max_db_diff=20.0, max_divisor=4,
858 min_group_size=3, max_harmonics_db=-5.0,
859 max_harmonics=0, max_groups=0, **kwargs):
860 """Detect peaks in power spectrum and group them according to their harmonic structure.
862 Parameters
863 ----------
864 psd_freqs: 1-D array
865 Frequencies of the power spectrum.
866 psd: 1-D array
867 Power spectrum (linear, not decible).
868 verbose: int
869 Verbosity level.
870 check_freqs: list of float
871 List of fundamental frequencies that will be checked first for being
872 present and valid harmonic groups in the peak frequencies
873 of the power spectrum.
874 low_threshold: float
875 The relative threshold for detecting all peaks in the decibel spectrum.
876 high_threshold: float
877 The relative threshold for detecting good peaks in the decibel spectrum.
878 thresh_bins: int or list of floats
879 Number of bins or the bins for computing the histogram from
880 which the standard deviation of the noise level in the `psd` is estimated.
881 low_thresh_factor: float
882 Factor by which the estimated standard deviation of the noise floor
883 is multiplied to set the `low_threshold`.
884 high_thresh_factor: float
885 Factor by which the estimated standard deviation of the noise floor
886 is multiplied to set the `high_threshold`.
887 freq_tol_fac: float
888 Harmonics should fall within `deltaf*freq_tol_fac`.
889 max_freq_tol: float
890 Maximum absolute frequency deviation of harmonics in Hertz..
891 mains_freq: float
892 Frequency of the mains power supply.
893 mains_freq_tol: float
894 Tolarerance around harmonics of the mains frequency,
895 within which peaks are removed.
896 min_freq: float
897 Minimum frequency accepted as a fundamental frequency.
898 max_freq: float
899 Maximum frequency accepted as a fundamental frequency.
900 max_db_diff: float
901 If larger than zero, maximum standard deviation of differences between
902 logarithmic powers of harmonics in decibel (larger than zero).
903 Low values enforce smoother power spectra.
904 max_divisor: int
905 Maximum divisor used for checking for sub-harmonics.
906 min_group_size: int
907 Minimum number of harmonics of a harmonic group.
908 The harmonics from min_group_size/3 to max(min_group_size, divisor)
909 need to be in good_freqs.
910 max_harmonics_db: float
911 Maximum allowed power of the `min_group_size`-th and higher harmonics
912 after the peak (in decibel relative to peak power withn the first
913 `min_group_size` harmonics, i.e. if harmonics are required to be
914 smaller than fundamental then this is a negative number).
915 Make it a large positive number to effectively not check for relative power.
916 max_harmonics: int
917 Maximum number of harmonics to be returned for each group.
918 max_groups: int
919 If not zero the maximum number of most powerful harmonic groups.
921 Returns
922 -------
923 groups: list of 2-D arrays
924 List of all extracted harmonic groups, sorted by fundamental frequency.
925 Each harmonic group is a 2-D array with the harmonics in the rows,
926 frequency of each harmonic in the first column and corresponding
927 power in the second column. If the power is zero, there was no
928 corresponding peak in the power spectrum.
929 fzero_harmonics: list of ints
930 The harmonics from which the fundamental frequencies were computed.
931 mains: 2-d array
932 Frequencies and power of multiples of the mains frequency found in the power spectrum.
933 all_freqs: 2-D array
934 Frequency, power, and use count (columns) of all peaks detected
935 in the power spectrum.
936 good_freqs: 1-D array
937 Frequencies of strong peaks detected in the power spectrum.
938 low_threshold: float
939 The relative threshold for detecting all peaks in the decibel spectrum.
940 high_threshold: float
941 The relative threshold for detecting good peaks in the decibel spectrum.
942 center: float
943 The baseline level of the power spectrum.
944 """
945 if verbose > 0:
946 print('')
947 if verbose > 1:
948 print(70*'#')
949 print('##### harmonic_groups', 48*'#')
951 # decibel power spectrum:
952 log_psd = decibel(psd)
953 max_idx = np.argmax(~np.isfinite(log_psd))
954 if max_idx > 0:
955 log_psd = log_psd[:max_idx]
956 delta_f = psd_freqs[1] - psd_freqs[0]
958 # thresholds:
959 center = np.nan
960 if low_threshold <= 0.0 or high_threshold <= 0.0:
961 low_th, high_th, center = threshold_estimate(log_psd, low_thresh_factor,
962 high_thresh_factor,
963 thresh_bins)
964 if low_threshold <= 0.0:
965 low_threshold = low_th
966 if high_threshold <= 0.0:
967 high_threshold = high_th
968 if verbose > 1:
969 print('')
970 print(f'low_threshold ={low_threshold:4.1f}dB, center+low_threshold ={center + low_threshold:6.1f}dB')
971 print(f'high_threshold={high_threshold:4.1f}dB, center+high_threshold={center + high_threshold:6.1f}dB')
972 print(f' center={center:6.1f}dB')
974 # detect peaks in decibel power spectrum:
975 peaks, troughs = detect_peaks(log_psd, low_threshold)
976 peaks, troughs = trim(peaks, troughs)
977 all_freqs = np.zeros((len(peaks), 3))
978 all_freqs[:,0] = psd_freqs[peaks]
979 all_freqs[:,1] = psd[peaks]
981 if len(all_freqs) == 0:
982 return [], [], [], np.zeros((0, 3)), [], low_threshold, high_threshold, center
984 # select good peaks:
985 good_freqs = all_freqs[(log_psd[peaks] - log_psd[troughs] > high_threshold) &
986 (all_freqs[:,0] >= min_freq) &
987 (all_freqs[:,0] < max_freq*max_divisor),:]
989 # detect harmonic groups:
990 freq_tol = delta_f*freq_tol_fac
991 if max_freq_tol < 1.1*freq_tol:
992 max_freq_tol = 1.1*freq_tol
993 groups, fzero_harmonics, mains = \
994 extract_fundamentals(good_freqs, all_freqs,
995 freq_tol, max_freq_tol,
996 verbose, check_freqs, mains_freq, mains_freq_tol,
997 min_freq, max_freq, max_db_diff, max_divisor, min_group_size,
998 max_harmonics_db, max_harmonics, max_groups)
1000 return (groups, fzero_harmonics, mains, all_freqs, good_freqs[:,0],
1001 low_threshold, high_threshold, center)
1004def fundamental_freqs(group_list):
1005 """
1006 Extract fundamental frequencies from lists of harmonic groups.
1008 The inner list of 2-D arrays of the input argument is transformed into
1009 a 1-D array containig the fundamental frequencies extracted from
1010 the 2-D arrays.
1012 Parameters
1013 ----------
1014 group_list: (list of (list of ...)) list of 2-D arrays
1015 Arbitrarily nested lists of harmonic groups as returned by
1016 extract_fundamentals() and harmonic_groups() with the element
1017 [0, 0] being the fundamental frequency.
1019 Returns
1020 -------
1021 fundamentals: (list of (list of ...)) 1-D array
1022 Nested list (corresponding to `group_list`) of 1-D arrays
1023 with the fundamental frequencies extracted from the harmonic groups.
1024 """
1025 if len(group_list) == 0:
1026 return np.zeros(0)
1028 # check whether group_list is list of harmonic groups:
1029 list_of_groups = True
1030 for group in group_list:
1031 if not ( hasattr(group, 'shape') and len(group.shape) == 2 ):
1032 list_of_groups = False
1033 break
1035 if list_of_groups:
1036 fundamentals = np.array([group[0, 0] for group in group_list if len(group) > 0])
1037 else:
1038 fundamentals = []
1039 for groups in group_list:
1040 f = fundamental_freqs(groups)
1041 fundamentals.append(f)
1042 return fundamentals
1045def fundamental_freqs_and_power(group_list, power=False,
1046 ref_power=1.0, min_power=1e-20):
1047 """Extract fundamental frequencies and their power in decibel from lists
1048 of harmonic groups.
1050 The inner list of 2-D arrays of the input argument is transformed
1051 into a 2-D array containig for each fish (1st dimension) the
1052 fundamental frequencies and powers (summed over all harmonics)
1053 extracted from the 2-D arrays.
1055 Parameters
1056 ----------
1057 group_list: (list of (list of ...)) list of 2-D arrays
1058 Arbitrarily nested lists of harmonic groups as returned by
1059 extract_fundamentals() and harmonic_groups() with the element
1060 [0, 0] being the fundamental frequency and the elements [:,1] being
1061 the powers of each harmonics.
1062 power: boolean
1063 If `False` convert the power into decibel using the
1064 `thunderlab.powerspectrum.decibel()` function.
1065 ref_power: float
1066 Reference power for computing decibel.
1067 If set to `None` the maximum power is used.
1068 min_power: float
1069 Power values smaller than `min_power` are set to `-np.inf`.
1071 Returns
1072 -------
1073 fundamentals: (list of (list of ...)) 2-D array
1074 Nested list (corresponding to `group_list`) of 2-D arrays
1075 with fundamental frequencies in first column and
1076 corresponding power in second column.
1077 """
1078 if len(group_list) == 0:
1079 return np.zeros((0, 2))
1081 # check whether group_list is list of harmonic groups:
1082 list_of_groups = True
1083 for group in group_list:
1084 if not ( hasattr(group, 'shape') and len(group.shape) == 2 ):
1085 list_of_groups = False
1086 break
1088 if list_of_groups:
1089 fundamentals = np.array([[group[0, 0], np.sum(group[:, 1])]
1090 for group in group_list if len(group) > 0])
1091 if not power:
1092 fundamentals[:, 1] = decibel(fundamentals[:, 1],
1093 ref_power, min_power)
1094 else:
1095 fundamentals = []
1096 for groups in group_list:
1097 f = fundamental_freqs_and_power(groups, power,
1098 ref_power, min_power)
1099 fundamentals.append(f)
1100 return fundamentals
1103def add_relative_power(freqs):
1104 """Add a column with relative power.
1106 For each element in `freqs`, its maximum power is subtracted
1107 from all powers.
1109 Parameters
1110 ----------
1111 freqs: list of 2-D arrays
1112 First column in the ndarrays is fundamental frequency and
1113 second column the corresponding power.
1114 Further columns are optional and kept in the returned list.
1115 fundamental_freqs_and_power() returns such a list.
1117 Returns
1118 -------
1119 power_freqs: list of 2-D arrays
1120 Same as freqs, but with an added column containing the relative power.
1121 """
1122 return [np.column_stack((f, f[:, 1] - np.max(f[:, 1]))) for f in freqs]
1125def add_power_ranks(freqs):
1126 """Add a column with power ranks.
1128 Parameters
1129 ----------
1130 freqs: list of 2-D arrays
1131 First column in the arrays is fundamental frequency and
1132 second column the corresponding power.
1133 Further columns are optional and kept in the returned list.
1134 fundamental_freqs_and_power() returns such a list.
1136 Returns
1137 -------
1138 rank_freqs: list of 2-D arrays
1139 Same as freqs, but with an added column containing the ranks.
1140 The highest power is assinged to zero,
1141 lower powers are assigned negative integers.
1142 """
1143 rank_freqs = []
1144 for f in freqs:
1145 i = np.argsort(f[:, 1])[::-1]
1146 ranks = np.empty_like(i)
1147 ranks[i] = -np.arange(len(i))
1148 rank_freqs.append(np.column_stack((f, ranks)))
1149 return rank_freqs
1152def similar_indices(freqs, df_thresh, nextfs=0):
1153 """Indices of similar frequencies.
1155 If two frequencies from different elements in the inner lists of `freqs` are
1156 reciprocally the closest to each other and closer than `df_thresh`,
1157 then two indices (element, frequency) of the respective other frequency
1158 are appended.
1160 Parameters
1161 ----------
1162 freqs: (list of (list of ...)) list of 1-D or 2-D arrays
1163 First column in the inner arrays is fundamental frequency.
1164 df_thresh: float
1165 Fundamental frequencies reciprocally being the closest to each other
1166 and closer than this threshold are considered equal.
1167 nextfs: int
1168 If zero, compare all elements in freqs with each other. Otherwise,
1169 only compare with the `nextfs` next elements in freqs.
1171 Returns
1172 -------
1173 indices: (list of (list of ...)) list of list of two-tuples of int
1174 For each frequency of each element in `freqs` a list of two tuples
1175 containing the indices of elements and frequencies that are similar.
1176 """
1177 if len(freqs) == 0:
1178 return []
1180 # check whether freqs is list of array of fundamental frequencies:
1181 list_of_freqs = True
1182 for group in freqs:
1183 if not (hasattr(group, 'shape') and 1 <= len(group.shape) <= 2):
1184 list_of_freqs = False
1185 break
1187 if list_of_freqs:
1188 indices = [ [[] for j in range(len(freqs[i]))] for i in range(len(freqs))]
1189 for j in range(len(freqs) - 1):
1190 freqsj = freqs[j]
1191 nn = len(freqs) if nextfs == 0 else j + nextfs + 1
1192 if nn > len(freqs):
1193 nn = len(freqs)
1194 for m in range(len(freqsj)):
1195 freq1 = freqsj[m]
1196 for k in range(j + 1, nn):
1197 freqsk = freqs[k]
1198 if len(freqsk) == 0:
1199 continue
1200 n = np.argmin(np.abs(freqsk[:, 0] - freq1[0]))
1201 freq2 = freqsk[n]
1202 if np.argmin(np.abs(freqsj[:, 0] - freq2[0])) != m:
1203 continue
1204 if np.abs(freq1[0] - freq2[0]) < df_thresh:
1205 indices[j][m].append((k, n))
1206 indices[k][n].append((j, m))
1207 else:
1208 indices = []
1209 for groups in freqs:
1210 indices.append(similar_indices(groups, df_thresh, nextfs))
1211 return indices
1214def unique_indices(freqs, df_thresh, nextfs=0):
1215 """For each unique fundamental frequency a list of indices indicating where it is found.
1217 Parameters
1218 ----------
1219 freqs: (list of (list of ...)) list of 1-D or 2-D arrays
1220 First column in the inner arrays is fundamental frequency.
1221 df_thresh: float
1222 Fundamental frequencies reciprocally being the closest to each other
1223 and closer than this threshold are considered equal.
1224 nextfs: int
1225 If zero, compare all elements in `freqs` with each other. Otherwise,
1226 only compare with the `nextfs` next elements in `freqs`.
1228 Returns
1229 -------
1230 freq_indices: (list of (list of ...)) list of list of 2-tuple of int
1231 For each unique frequency in `freqs`, a list of two-tuples of indices
1232 into `freqs`, indicating where in `freqs` these frequencies are.
1233 """
1234 if len(freqs) == 0:
1235 return []
1237 # check whether freqs is list of array of fundamental frequencies:
1238 list_of_freqs = True
1239 for group in freqs:
1240 if not (hasattr(group, 'shape') and 1 <= len(group.shape) <= 2):
1241 list_of_freqs = False
1242 break
1244 if list_of_freqs:
1245 freq_indices = []
1246 indices = similar_indices(freqs, df_thresh, 0)
1247 if len(indices) == 0:
1248 return freq_indices
1249 for j in range(len(indices)):
1250 for m in range(len(indices[j])):
1251 if indices[j][m] is not None:
1252 # store:
1253 freq_indices.append([(j, m)])
1254 freq_indices[-1].extend(indices[j][m])
1255 # remove from indices:
1256 for inx in indices[j][m]:
1257 indices[inx[0]][inx[1]] = None
1258 else:
1259 indices = []
1260 for groups in freqs:
1261 indices.append(unique_indices(groups, df_thresh, nextfs))
1262 return freq_indices
1265def collapse(group_list):
1266 """Collapse harmonic groups with similar fundamental frequencies to a single one.
1268 The group with the highest power (second column summed over all harmonics) is returned.
1270 Parameters
1271 ----------
1272 group_list: list of 2-D arrays of float
1273 List of harmonic groups that are supposed to have similar fundamental frequencies.
1275 Returns
1276 -------
1277 group: 2-D array of float
1278 Harmonic group with first column frequency of harmonics, second column corresponding powers.
1279 Further columns are optional.
1280 index: int
1281 Index of the returned harmonic group.
1282 """
1283 if len(group_list) == 0:
1284 return np.zeros((0, 2)), 0
1285 rgroup = group_list[0]
1286 rindex = 0
1287 for k, group in enumerate(group_list[1:]):
1288 if np.sum(group[:, 1]) > np.sum(rgroup[:, 1]):
1289 rgroup = group
1290 rindex = 1 + k
1291 return rgroup, rindex
1294def consistent(group_list, df_thresh):
1295 """Harmonic groups whose fundamental frequencies consistently appear everywhere.
1297 Parameters
1298 ----------
1299 group_list: (list of (list of ...)) list of list of 2-D arrays of float
1300 Harmonic groups as returned by several calls to extract_fundamentals()
1301 or harmonic_groups() with the element [0, 0] of the harmonic groups
1302 being the fundamental frequency, and element[0, 1] the corresponding
1303 power.
1304 df_thresh: float
1305 Fundamental frequencies reciprocally being the closest to each other
1306 and closer than this threshold are considered equal.
1308 Returns
1309 -------
1310 consistent_groups: (list of (list of ...)) list of 2-D arrays of float
1311 List of harmonic groups that consistently occur in all elements of
1312 `group_list`. Each element has frequency of the harmonics in the
1313 first column and corresponding power in the second column.
1314 """
1315 if len(group_list) == 0:
1316 return []
1318 # check whether group_list is list of harmonic groups:
1319 list_of_groups = True
1320 for groups in group_list:
1321 for group in groups:
1322 if not ( hasattr(group, 'shape') and len(group.shape) == 2 ):
1323 list_of_groups = False
1324 break
1325 if not list_of_groups:
1326 break
1328 if list_of_groups:
1329 consistent_groups = []
1330 freqs = fundamental_freqs_and_power(group_list)
1331 indices = unique_indices(freqs, df_thresh, 0)
1332 for idx in indices:
1333 # frequency appears everywhere:
1334 if len(idx) == len(freqs):
1335 cgroup, _ = collapse([group_list[j][k] for j, k in idx])
1336 consistent_groups.append(cgroup)
1337 else:
1338 consistent_groups = []
1339 for groups in group_list:
1340 consistent_groups.append(consistent(groups, df_thresh))
1341 return consistent_groups
1344def closest(group_list, df_thresh, close_thresh=0, min_closest=1):
1345 """Harmonic groups whose fundamental frequencies are closest to each other.
1347 Parameters
1348 ----------
1349 group_list: (list of (list of ...)) list of list of 2-D arrays of float
1350 Harmonic groups as returned by several calls to extract_fundamentals()
1351 or harmonic_groups() with the element [0, 0] of the harmonic groups
1352 being the fundamental frequency, and element[0, 1] the corresponding
1353 power.
1354 df_thresh: float
1355 Fundamental frequencies reciprocally being the closest to each other
1356 and closer than this threshold are considered equal and form a unique group.
1357 close_thresh: float
1358 Subsequent fundamental frequencies within a unique group that are closer
1359 than this threshold are considered close.
1360 min_closest: int
1361 Only return harmonic groups that have at minimum that many
1362 harmonic groups whose fundamental frequencies are closest to
1363 each other.
1365 Returns
1366 -------
1367 closest_groups: (list of (list of ...)) list of 2-D arrays of float
1368 For each unique fundamental frequency in `group_list` the
1369 harmonic group with the largest power among subsequent ones whose
1370 fundamental frequencies are closest to each other. Each element
1371 has frequency of the harmonics in the first column and
1372 corresponding power in the second column.
1373 closest_indices: (list of (list of ...)) 2-D array of int
1374 For each unique fundamental frequency in `group_list` (rows)
1375 the first and last index (two columns) into `group_list`
1376 where the closest frequencies have been found.
1378 """
1379 if len(group_list) == 0:
1380 return [], np.zeros((0, 2))
1382 # check whether group_list is list of harmonic groups:
1383 list_of_groups = True
1384 for groups in group_list:
1385 for group in groups:
1386 if not ( hasattr(group, 'shape') and len(group.shape) == 2 ):
1387 list_of_groups = False
1388 break
1389 if not list_of_groups:
1390 break
1392 if list_of_groups:
1393 if close_thresh < 1e-8:
1394 close_thresh = 1e-8
1395 closest_groups = []
1396 closest_indices = []
1397 freqs = fundamental_freqs_and_power(group_list)
1398 freq_indices = unique_indices(freqs, df_thresh, 0)
1399 for fidx in freq_indices:
1400 if len(fidx) == 0:
1401 continue
1402 if len(fidx) < 2:
1403 if min_closest < 2:
1404 cgroup = group_list[fidx[0][0]][fidx[0][1]]
1405 closest_groups.append(cgroup)
1406 closest_indices.append((fidx[0][0], fidx[0][0]))
1407 continue
1408 # find minimum difference in subsequent fundamental frequencies:
1409 fp = []
1410 fpi = []
1411 fi = 0
1412 for i in range(len(freqs)):
1413 if fi < len(fidx):
1414 j, k = fidx[fi]
1415 else:
1416 j = len(freqs)
1417 fpi. append(fi)
1418 if j == i:
1419 fp. append(freqs[j][k][0])
1420 fi += 1
1421 else:
1422 fp. append(np.nan)
1423 fp = np.array(fp)
1424 dfs = np.abs(np.diff(fp))
1425 if np.all(np.isnan(dfs)):
1426 if min_closest <= 1:
1427 cgroup, cinx = collapse([group_list[j][k] for j, k in fidx])
1428 closest_groups.append(cgroup)
1429 closest_indices.append((fidx[cinx][0], fidx[cinx][0]))
1430 continue
1431 i0 = fpi[np.nanargmin(dfs)]
1432 i1 = i0 + 1
1433 dfs = np.abs(np.diff(fp[np.isfinite(fp)]))
1434 # expand:
1435 for k in range(1, len(fidx)):
1436 go_left = i0 > 0 and fidx[i0 - 1][0] + 1 == fidx[i0][0] and abs(dfs[i0 - 1]) < close_thresh
1437 go_right = i1 < len(dfs) and fidx[i1][0] + 1 == fidx[i1 + 1][0] and abs(dfs[i1]) < close_thresh
1438 if go_left and go_right:
1439 if dfs[i0 - 1] <= dfs[i1]:
1440 i0 -= 1
1441 else:
1442 i1 += 1
1443 elif go_left:
1444 i0 -= 1
1445 elif go_right:
1446 i1 += 1
1447 else:
1448 break
1449 if i1 + 1 - i0 >= min_closest:
1450 cgroup, _ = collapse([group_list[j][k] for j, k in fidx[i0:i1 + 1]])
1451 closest_groups.append(cgroup)
1452 closest_indices.append((fidx[i0][0], fidx[i1][0]))
1453 closest_indices = np.array(closest_indices)
1454 else:
1455 closest_groups = []
1456 closest_indices = []
1457 for groups in group_list:
1458 cg, ci = closest(groups, df_thresh, min_closest)
1459 closest_groups.append(cg)
1460 closest_indices.append(ci)
1461 return closest_groups, closest_indices
1464def unique_mask(freqs, df_thresh, nextfs=0):
1465 """Mark similar frequencies from different recordings as dublicate.
1467 If two frequencies from different elements in `freqs` are
1468 reciprocally the closest to each other and closer than `df_thresh`,
1469 then the one with the smaller power is marked for removal.
1471 Parameters
1472 ----------
1473 freqs: list of 2-D arrays
1474 First column in the arrays is fundamental frequency and
1475 second column the corresponding power or equivalent.
1476 If values in the second column are equal (e.g. they are the same ranks),
1477 and there is a third column (e.g. power),
1478 the third column is used to decide, which element should be removed.
1479 df_thresh: float
1480 Fundamental frequencies closer than this threshold are considered
1481 equal.
1482 nextfs: int
1483 If zero, compare all elements in freqs with each other. Otherwise,
1484 only compare with the `nextfs` next elements in freqs.
1486 Returns
1487 -------
1488 mask: list of boolean arrays
1489 For each element in `freqs` True if that frequency should be kept.
1490 """
1491 mask = [np.ones(len(freqs[i]), dtype=bool) for i in range(len(freqs))]
1492 for j in range(len(freqs)-1):
1493 freqsj = np.asarray(freqs[j])
1494 for m in range(len(freqsj)):
1495 freq1 = freqsj[m]
1496 nn = len(freqs) if nextfs == 0 else j+1+nextfs
1497 if nn > len(freqs):
1498 nn = len(freqs)
1499 for k in range(j+1, nn):
1500 freqsk = np.asarray(freqs[k])
1501 if len(freqsk) == 0:
1502 continue
1503 n = np.argmin(np.abs(freqsk[:,0] - freq1[0]))
1504 freq2 = freqsk[n]
1505 if np.argmin(np.abs(freqsj[:,0] - freq2[0])) != m:
1506 continue
1507 if np.abs(freq1[0] - freq2[0]) < df_thresh:
1508 if freq1[1] > freq2[1]:
1509 mask[k][n] = False
1510 elif freq1[1] < freq2[1]:
1511 mask[j][m] = False
1512 elif len(freq1) > 2:
1513 if freq1[2] > freq2[2]:
1514 mask[k][n] = False
1515 else:
1516 mask[j][m] = False
1517 else:
1518 mask[j][m] = False
1519 return mask
1522def unique(freqs, df_thresh, mode='power', nextfs=0):
1523 """Remove similar frequencies from different recordings.
1525 If two frequencies from different elements in the inner lists of `freqs`
1526 are reciprocally the closest to each other and closer than `df_thresh`,
1527 then the one with the smaller power is removed. As power, either the
1528 absolute power as provided in the second column of the data elements
1529 in `freqs` is taken (mode=='power'), or the relative power
1530 (mode='relpower'), or the power rank (mode='rank').
1532 Parameters
1533 ----------
1534 freqs: (list of (list of ...)) list of 2-D arrays
1535 First column in the arrays is fundamental frequency and
1536 second column the corresponding power, as returned by
1537 fundamental_freqs_and_power().
1538 df_thresh: float
1539 Fundamental frequencies closer than this threshold are considered
1540 equal.
1541 mode: string
1542 - 'power': use second column of freqs elements as power.
1543 - 'relpower': use relative power computed from the second column
1544 of freqs elements for deciding which frequency to delete.
1545 - 'rank': use rank of second column of freqs elements
1546 for deciding which frequency to delete.
1547 nextfs: int
1548 If zero, compare all elements in freqs with each other. Otherwise,
1549 only compare with the `nextfs` next elements in freqs.
1551 Returns
1552 -------
1553 unique_freqs: (list of (list of ...)) list of 2-D arrays
1554 Same as `freqs` but elements with similar fundamental frequencies
1555 removed.
1557 Raises
1558 ------
1559 ValueError:
1560 Invalid mode string.
1561 """
1562 if len(freqs) == 0:
1563 return []
1565 # check whether freqs is array of fundamental frequencies and powers:
1566 list_of_freq_power = True
1567 for group in freqs:
1568 if not (hasattr(group, 'shape') and len(group.shape) == 2):
1569 list_of_freq_power = False
1570 break
1572 if list_of_freq_power:
1573 if mode == 'power':
1574 mask = unique_mask(freqs, df_thresh, nextfs)
1575 elif mode == 'relpower':
1576 power_freqs = [f[:, [0, 2, 1]] for f in add_relative_power(freqs)]
1577 mask = unique_mask(power_freqs, df_thresh, nextfs)
1578 elif mode == 'rank':
1579 rank_freqs = [f[:, [0, 2, 1]] for f in add_power_ranks(freqs)]
1580 mask = unique_mask(rank_freqs, df_thresh, nextfs)
1581 else:
1582 raise ValueError(f'{mode} is not a valid mode for unique(). Choose one of "power", "relpower", or "rank"')
1583 unique_freqs = []
1584 for f, m in zip(freqs, mask):
1585 unique_freqs.append(f[m])
1586 else:
1587 unique_freqs = []
1588 for groups in freqs:
1589 unique_freqs.append(unique(groups, df_thresh, mode, nextfs))
1590 return unique_freqs
1593def colors_markers():
1594 """Generate a list of colors and markers for plotting.
1596 Returns
1597 -------
1598 colors: list
1599 list of colors
1600 markers: list
1601 list of markers
1602 """
1603 # color and marker range:
1604 colors = []
1605 markers = []
1606 mr2 = []
1608 # first color range:
1609 cc0 = cm.gist_rainbow(np.linspace(0.0, 1.0, 8))
1611 # shuffle it:
1612 for k in range((len(cc0) + 1) // 2):
1613 colors.extend(cc0[k::(len(cc0) + 1) // 2])
1614 markers.extend(len(cc0) * 'o')
1615 mr2.extend(len(cc0) * 'v')
1616 # second darker color range:
1617 cc1 = cm.gist_rainbow(np.linspace(0.33 / 7.0, 1.0, 7))
1618 cc1 = mc.hsv_to_rgb(mc.rgb_to_hsv(np.array([cc1[:, :3]])) * np.array([1.0, 0.9, 0.7]))[0]
1619 cc1 = np.hstack((cc1, np.ones((len(cc1),1))))
1620 # shuffle it:
1621 for k in range((len(cc1) + 1) // 2):
1622 colors.extend(cc1[k::(len(cc1) + 1) // 2])
1623 markers.extend(len(cc1) * '^')
1624 mr2.extend(len(cc1) * '*')
1625 # third lighter color range:
1626 cc2 = cm.gist_rainbow(np.linspace(0.67 / 6.0, 1.0, 6))
1627 cc2 = mc.hsv_to_rgb(mc.rgb_to_hsv(np.array([cc1[:, :3]])) * np.array([1.0, 0.5, 1.0]))[0]
1628 cc2 = np.hstack((cc2, np.ones((len(cc2),1))))
1629 # shuffle it:
1630 for k in range((len(cc2) + 1) // 2):
1631 colors.extend(cc2[k::(len(cc2) + 1) // 2])
1632 markers.extend(len(cc2) * 'D')
1633 mr2.extend(len(cc2) * 'x')
1634 markers.extend(mr2)
1635 return colors, markers
1638def plot_harmonic_groups(ax, group_list, indices=None, max_groups=0,
1639 skip_bad=False, sort_by_freq=True, label_power=False,
1640 colors=None, markers=None, marker_size=18,
1641 legend_rows=8, **kwargs):
1642 """Mark decibel power of fundamentals and their harmonics in a plot.
1644 Parameters
1645 ----------
1646 ax: axis for plot
1647 Axis used for plotting.
1648 group_list: list of 2-D arrays
1649 Lists of harmonic groups as returned by extract_fundamentals() and
1650 harmonic_groups() with the element [0, 0] of the harmonic groups being
1651 the fundamental frequency, and element[0, 1] being the corresponding power.
1652 indices: list of int or None
1653 If smaller than zero then set the legend label of the corresponding group in brackets.
1654 max_groups: int
1655 If not zero plot only the max_groups most powerful groups.
1656 skip_bad: bool
1657 Skip harmonic groups without index (entry in indices is negative).
1658 sort_by_freq: boolean
1659 If True sort legend by frequency, otherwise by power.
1660 label_power: boolean
1661 If `True` put the power in decibel in addition to the frequency into the legend.
1662 colors: list of colors or None
1663 If not None list of colors for plotting each group
1664 markers: list of markers or None
1665 If not None list of markers for plotting each group
1666 marker_size: float
1667 Maximum size of markers used to mark the pulses.
1668 legend_rows: int
1669 Maximum number of rows to be used for the legend.
1670 kwargs:
1671 Key word arguments for the legend of the plot.
1673 Returns
1674 -------
1675 groups_dict: dict
1676 Dictionary that maps artist handles marking the harmonic groups
1677 with tuples containing the index and the 2D arrays from the
1678 `group_list`.
1679 groups_colors: list of matplotlib colors
1680 For each group the color that was used to mark it in
1681 the power spectrum.
1682 Color can be None if not plotted.
1683 groups_markers: list of string
1684 For each group the marker symbol that was used to mark it in
1685 the power spectrum.
1686 Marker symbol can be None if not plotted.
1687 groups_sizes: list of float
1688 For each group the marker size that was used to mark it in
1689 the power spectrum relative to the maximumpossible size.
1690 Marker size can be None if not plotted.
1691 """
1692 if len(group_list) == 0:
1693 return {}, [], [], []
1695 # sort by power:
1696 powers = np.array([np.sum(group[:, 1]) for group in group_list])
1697 max_power = np.max(powers)
1698 idx = np.argsort(-powers)
1699 if max_groups > 0 and len(idx > max_groups):
1700 idx = idx[:max_groups]
1702 # sort by frequency:
1703 if sort_by_freq:
1704 freqs = np.array([group_list[group][0, 0] for group in idx])
1705 idx = idx[np.argsort(freqs)]
1707 # plot:
1708 groups_dict = {}
1709 artists = []
1710 labels = []
1711 groups_colors = [None]*len(group_list)
1712 groups_markers = [None]*len(group_list)
1713 groups_sizes = [None]*len(group_list)
1714 k = 0
1715 for i in idx:
1716 if indices is not None and skip_bad and indices[i] < 0:
1717 continue
1718 group = group_list[i]
1719 x = group[:, 0]
1720 y = decibel(group[:, 1])
1721 msize = marker_size*(0.3 + 0.7*(powers[i]/max_power)**0.25)
1722 color_kwargs = {}
1723 if colors is not None:
1724 color_kwargs = {'color': colors[k%len(colors)]}
1725 label = f'{group[0, 0]:6.1f} Hz'
1726 if label_power:
1727 label += f' {decibel(np.array([np.sum(group[:,1])]))[0]:6.1f}dB'
1728 if indices is not None:
1729 if indices[i] < 0:
1730 label = '(' + label + ')'
1731 else:
1732 label = ' ' + label + ' '
1733 if markers is None:
1734 aa = ax.plot(x, y, 'o', label=label, ms=msize, **color_kwargs)
1735 else:
1736 if k >= len(markers):
1737 break
1738 aa = ax.plot(x, y, label=label, linestyle='none', marker=markers[k],
1739 mec=None, mew=0.0, ms=msize, **color_kwargs)
1740 groups_colors[i] = aa[0].get_markerfacecolor()
1741 groups_markers[i] = aa[0].get_marker()
1742 groups_sizes[i] = aa[0].get_markersize()/marker_size
1743 for a in aa:
1744 a.set_picker(8)
1745 groups_dict[a] = [i, group]
1746 k += 1
1748 # legend:
1749 ncol = (k - 1) // legend_rows + 1 if legend_rows > 0 else 1
1750 if ncol < 1:
1751 ncol = 1
1752 try:
1753 handles = ax.legend(numpoints=1, ncols=ncol, **kwargs)
1754 except TypeError:
1755 handles = ax.legend(numpoints=1, ncol=ncol, **kwargs)
1756 lines = handles.get_lines()
1757 for k in range(min(len(lines), len(idx))):
1758 a = lines[k]
1759 a.set_picker(8)
1760 i = idx[k]
1761 groups_dict[a] = [i, group_list[i]]
1762 return groups_dict, groups_colors, groups_markers, groups_sizes
1765def plot_psd_harmonic_groups(ax, psd_freqs, psd, group_list,
1766 mains=None, all_freqs=None,
1767 good_freqs=None, log_freq=False,
1768 min_freq=0.0, max_freq=2000.0, ymarg=0.0,
1769 sstyle=dict(color='tab:blue', lw=1)):
1770 """Plot decibel power-spectrum with detected peaks, harmonic groups,
1771 and mains frequencies.
1773 Parameters
1774 ----------
1775 ax: axis for plot
1776 Axis used for plotting.
1777 psd_freqs: array
1778 Frequencies of the power spectrum.
1779 psd: array
1780 Power spectrum (linear, not decible).
1781 group_list: list of 2-D arrays
1782 Lists of harmonic groups as returned by extract_fundamentals() and
1783 harmonic_groups() with the element [0, 0] of the harmonic groups being the fundamental frequency,
1784 and element[0, 1] being the corresponding power.
1785 mains: 2-D array
1786 Frequencies and power of multiples of the mains frequency found in the power spectrum.
1787 all_freqs: 2-D array
1788 Peaks in the power spectrum detected with low threshold.
1789 good_freqs: 1-D array
1790 Frequencies of peaks detected with high threshold.
1791 log_freq: boolean
1792 Logarithmic (True) or linear (False) frequency axis.
1793 min_freq: float
1794 Limits of frequency axis are set to `(min_freq, max_freq)`
1795 if `max_freq` is greater than zero
1796 max_freq: float
1797 Limits of frequency axis are set to `(min_freq, max_freq)`
1798 and limits of power axis are computed from powers below max_freq
1799 if `max_freq` is greater than zero
1800 ymarg: float
1801 Add this to the maximum decibel power for setting the ylim.
1802 sstyle: dict
1803 Arguments passed on to the plot command for the power spectrum.
1804 """
1806 # mark all and good psd peaks:
1807 pmin, pmax = ax.get_ylim()
1808 doty = pmax - 5.0
1809 if all_freqs is not None:
1810 ax.plot(all_freqs[:, 0], np.zeros(len(all_freqs[:, 0])) + doty, 'o', color='#ffffff')
1811 if good_freqs is not None:
1812 ax.plot(good_freqs, np.zeros(len(good_freqs)) + doty, 'o', color='#888888')
1813 # mark mains frequencies:
1814 if mains is not None and len(mains) > 0:
1815 fpeaks = mains[:, 0]
1816 fpeakinx = [int(np.round(fp/(psd_freqs[1]-psd_freqs[0]))) for fp in fpeaks if fp < psd_freqs[-1]]
1817 ax.plot(fpeaks[:len(fpeakinx)], decibel(psd[fpeakinx]), linestyle='None',
1818 marker='.', color='k', ms=10, mec=None, mew=0.0,
1819 label=f'{mains[0, 0]:3.0f} Hz mains')
1820 # mark harmonic groups:
1821 colors, markers = colors_markers()
1822 plot_harmonic_groups(ax, group_list, max_groups=0, sort_by_freq=True,
1823 colors=colors, markers=markers, legend_rows=8,
1824 loc='upper right')
1825 # plot power spectrum:
1826 plot_decibel_psd(ax, psd_freqs, psd, log_freq=log_freq, min_freq=min_freq,
1827 max_freq=max_freq, ymarg=ymarg, sstyle=sstyle)
1830def annotate_harmonic_group(ax, group, index=None, all_peaks=None,
1831 freq_thresh=None):
1832 """Annotate peaks of a harmonic group in a spectrum.
1834 Parameters
1835 ----------
1836 ax: axis for plot
1837 Axis used for plotting.
1838 group: 2D array of float
1839 All harmonics of a harmonic group first column their
1840 frequencies and second comun their power. Element [0, 0] is
1841 its fundamental frequency.
1842 index: None or ont
1843 Index of the group.
1844 all_freqs: None or 2D array of float
1845 All peaks in the power spectrum detected with low threshold.
1846 freq_thresh: None or float
1847 Threshold for matching peak frequencies. Something like
1848 `0.8*delta_f`, where `delta_f` is the frequency resolution
1849 of the power spectrum.
1851 Returns
1852 -------
1853 artists: list of matplotlib.Artist
1854 List of the annotations. Use it to remove the annotations:
1855 ```
1856 for a in artists:
1857 a.remove()
1858 artists = []
1859 ax.get_figure().canvas.draw()
1860 ```
1862 """
1863 fmin, fmax = ax.get_xlim()
1864 fwidth = fmax - fmin
1865 artists = []
1866 f1 = group[0, 0]
1867 for harmonic, (freq, power) in enumerate(group):
1868 count = None
1869 if all_peaks is not None and freq_thresh is not None:
1870 peak = all_peaks[np.abs(all_peaks[:, 0] - freq) < freq_thresh, :]
1871 if len(peak) == 0:
1872 continue
1873 count = peak[2]
1874 pl = []
1875 if index is not None:
1876 pl.append(f'{index}: {f1:.1f}Hz')
1877 pl.append(f'$f{harmonic + 1} = ${freq:.1f}Hz')
1878 pl.append(f'$p=${power:g}')
1879 if count is not None:
1880 pl.append(f'$c=${count}')
1881 a = ax.annotate('\n'.join(pl), xy=(freq, decibel(power)),
1882 xytext=(freq + 0.03*fwidth, decibel(power)),
1883 bbox=dict(boxstyle='round', facecolor='white'),
1884 arrowprops=dict(arrowstyle='-'))
1885 artists.append(a)
1886 return artists
1889def plot_selected_groups(ax, group_list, selected_groups, indices=None,
1890 label=None,
1891 freq_style=dict(ls='none', marker='o',
1892 color='k', ms=10),
1893 bar_style=dict(ls='-', color='c', lw=12)):
1894 """
1895 Plot all fundamental frequencies and mark selected frequencies with a bar.
1897 Parameters
1898 ----------
1899 ax: matplotlib.Axes
1900 Axes for plotting.
1901 group_list: list of list of 2-D arrays of float
1902 Harmonic groups as returned by several calls to extract_fundamentals()
1903 or harmonic_groups() with the element [0, 0] of the harmonic groups
1904 being the fundamental frequency, and element[0, 1] the corresponding
1905 power.
1906 selected_groups: list of 2-D arrays of float
1907 Frequencies and power of selected harmonic groups from `group_list`.
1908 indices: None or list of two-tuple of int
1909 Ranges of the selected groups.
1910 label: str or None
1911 Label for the selected groups that is added to the legend.
1912 freq_style: dict
1913 Plot style for marking fundamental frequencies of `group_list`.
1914 bar_style: dict
1915 Plot style for marking frequencies from `selected_groups`.
1916 """
1917 # mark selected frequencies:
1918 if indices is None:
1919 for index in range(len(selected_groups)):
1920 x = [-0.25, len(group_list) - 0.75]
1921 y = [selected_groups[index][0, 0]]*2
1922 if index == 0 and label is not None:
1923 ax.plot(x, y, label=label, **bar_style)
1924 else:
1925 ax.plot(x, y, **bar_style)
1926 else:
1927 for index in range(len(selected_groups)):
1928 x = [indices[index][0] - 0.25,
1929 indices[index][1] + 0.25]
1930 y = [selected_groups[index][0, 0]]*2
1931 if index == 0 and label is not None:
1932 ax.plot(x, y, label=label, **bar_style)
1933 else:
1934 ax.plot(x, y, **bar_style)
1935 # mark all frequencies:
1936 for index in range(len(group_list)):
1937 for group in group_list[index]:
1938 ax.plot(index, group[0, 0], **freq_style)
1939 ax.set_xlim([-0.5, len(group_list) - 0.5])
1940 ax.xaxis.set_major_locator(plt.MultipleLocator(1))
1941 ax.set_ylabel('frequency [Hz]')
1942 ax.set_xlabel('group index')
1945def plot_unique_frequencies(ax, freqs_list, unique_freqs,
1946 freq_style=dict(ls='none', marker='o',
1947 color='k', ms=10),
1948 unique_style=dict(ls='none', marker='o',
1949 color='r', ms=16)):
1950 """
1951 Plot all fundamental frequencies and mark unique frequencies.
1953 Parameters
1954 ----------
1955 ax: matplotlib.Axes
1956 Axes for plotting.
1957 freqs_list: list of 2-D arrays of float
1958 Harmonic groups as returned by severla calls to extract_fundamentals()
1959 or harmonic_groups() with the element [0, 0] of the harmonic groups
1960 being the fundamental frequency, and element[0, 1] the corresponding
1961 power.
1962 unique_freqs: list of 2-D arrays of float
1963 Frequencies and power of unique harmonic groups from `freqs_list`.
1964 freq_style: dict
1965 Plot style for marking fundamental frequencies of `freqs_list`.
1966 unique_style: dict
1967 Plot style for marking frequencies from `unique_freqs`.
1968 """
1969 for index in range(len(freqs_list)):
1970 freqs = freqs_list[index][:, 0]
1971 for freq in freqs:
1972 if len(unique_freqs[index][:, 0]) > 0 and \
1973 np.min(np.abs(unique_freqs[index][:, 0] - freq)) < 1e-8:
1974 ax.plot(index, freq, **unique_style)
1975 ax.plot(np.zeros(len(freqs)) + index, freqs, **freq_style)
1976 ax.set_xlim([-0.5, len(freqs_list) - 0.5])
1977 ax.xaxis.set_major_locator(plt.MultipleLocator(1))
1978 ax.set_ylabel('frequency [Hz]')
1979 ax.set_xlabel('group index')
1982def add_psd_peak_detection_config(cfg, low_threshold=0.0, high_threshold=0.0,
1983 thresh_bins=100,
1984 low_thresh_factor=6.0, high_thresh_factor=10.0):
1985 """Add parameter needed for detection of peaks in power spectrum used
1986 by harmonic_groups() as a new section to a configuration.
1988 Parameters
1989 ----------
1990 cfg: ConfigFile
1991 The configuration.
1992 """
1993 cfg.add_section('Thresholds for peak detection in power spectra:')
1994 cfg.add('lowThreshold', low_threshold, 'dB', 'Threshold for all peaks.\n If 0.0 estimate threshold from histogram.')
1995 cfg.add('highThreshold', high_threshold, 'dB', 'Threshold for good peaks. If 0.0 estimate threshold from histogram.')
1997 cfg.add_section('Threshold estimation:\nIf no thresholds are specified they are estimated from the histogram of the decibel power spectrum.')
1998 cfg.add('thresholdBins', thresh_bins, '', 'Number of bins used to compute the histogram used for threshold estimation.')
1999 cfg.add('lowThresholdFactor', low_thresh_factor, '', 'Factor for multiplying standard deviation of noise floor for lower threshold.')
2000 cfg.add('highThresholdFactor', high_thresh_factor, '', 'Factor for multiplying standard deviation of noise floor for higher threshold.')
2003def psd_peak_detection_args(cfg):
2004 """Translates a configuration to the respective parameter names for
2005 the detection of peaks in power spectrum used by
2006 harmonic_groups(). The return value can then be passed as
2007 key-word arguments to this function.
2009 Parameters
2010 ----------
2011 cfg: ConfigFile
2012 The configuration.
2014 Returns
2015 -------
2016 a: dict
2017 Dictionary with names of arguments of the `harmonic-group()` function
2018 and their values as supplied by `cfg`.
2019 """
2020 return cfg.map(low_threshold='lowThreshold',
2021 high_threshold='highThreshold',
2022 thresh_bins='thresholdBins',
2023 low_thresh_factor='lowThresholdFactor',
2024 high_thresh_factor='highThresholdFactor')
2027def add_harmonic_groups_config(cfg, mains_freq=60.0, mains_freq_tol=1.0,
2028 max_divisor=4, freq_tol_fac=1.0,
2029 max_freq_tol=1.0, min_group_size=3,
2030 min_freq=20.0, max_freq=2000.0, max_db_diff=20.0,
2031 max_harmonics_db=-5.0, max_harmonics=0, max_groups=0):
2032 """Add parameter needed for detection of harmonic groups as a new
2033 section to a configuration.
2035 Parameters
2036 ----------
2037 cfg: ConfigFile
2038 The configuration.
2039 """
2040 cfg.add_section('Harmonic groups:')
2041 cfg.add('mainsFreq', mains_freq, 'Hz', 'Mains frequency to be excluded.')
2042 cfg.add('mainsFreqTolerance', mains_freq_tol, 'Hz', 'Exclude peaks within this tolerance around multiples of the mains frequency.')
2043 cfg.add('minimumGroupSize', min_group_size, '',
2044'Minimum number of harmonics (inclusively fundamental) that make up a harmonic group.')
2045 cfg.add('maxDivisor', max_divisor, '', 'Maximum ratio between the frequency of the largest peak and its fundamental.')
2046 cfg.add('freqTolerance', freq_tol_fac, '',
2047 'Harmonics should be within this factor times the frequency resolution of the power spectrum. Needs to be higher than 0.5!')
2048 cfg.add('maximumFreqTolerance', max_freq_tol, 'Hz',
2049 'Maximum deviation of harmonics from their expected value.')
2051 cfg.add_section('Acceptance of best harmonic groups:')
2052 cfg.add('minimumFrequency', min_freq, 'Hz', 'Minimum frequency allowed for the fundamental.')
2053 cfg.add('maximumFrequency', max_freq, 'Hz', 'Maximum frequency allowed for the fundamental.')
2054 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.')
2055 cfg.add('maximumHarmonicsPower', max_harmonics_db, 'dB', 'Maximum allowed power of the minimumGroupSize-th and higher harmonics relative to peak power.')
2056 cfg.add('maximumHarmonics', max_harmonics, '', 'Only keep the first # harmonics. If 0 keep all.')
2057 cfg.add('maximumGroups', max_groups, '', 'Maximum number of harmonic groups. If 0 process all.')
2060def harmonic_groups_args(cfg):
2061 """Translates a configuration to the respective parameter names of
2062 the harmonic-group detection functions. The return value can then
2063 be passed as key-word arguments to this function.
2065 Parameters
2066 ----------
2067 cfg: ConfigFile
2068 The configuration.
2070 Returns
2071 -------
2072 a: dict
2073 Dictionary with names of arguments of the harmonic-group detection functions
2074 and their values as supplied by `cfg`.
2075 """
2076 return cfg.map(mains_freq='mainsFreq',
2077 mains_freq_tol='mainsFreqTolerance',
2078 freq_tol_fac='freqTolerance',
2079 max_freq_tol='maximumFreqTolerance',
2080 max_divisor='maxDivisor',
2081 min_group_size='minimumGroupSize',
2082 min_freq='minimumFrequency',
2083 max_freq='maximumFrequency',
2084 max_db_diff='maximumPowerDifference',
2085 max_harmonics_db='maximumHarmonicsPower',
2086 max_harmonics='maximumHarmonics',
2087 max_groups='maximumGroups')
2090def main(data_file=None):
2091 import matplotlib.pyplot as plt
2092 from thunderlab.powerspectrum import psd
2093 from .fakefish import wavefish_eods
2095 def mkg(f, p, n=5):
2096 return np.array([(h*f, p/h) for h in range(1, n)])
2098 if data_file is None:
2099 # generate data:
2100 title = 'simulation'
2101 rate = 44100.0
2102 d = 20.0
2103 noise = 0.01
2104 eodfs = [123.0, 333.0, 666.0, 666.5, 792.8]
2105 fish1 = 0.5*wavefish_eods('Eigenmannia', eodfs[0], rate,
2106 duration=d, noise_std=noise)
2107 fish2 = 1.0*wavefish_eods('Eigenmannia', eodfs[1], rate,
2108 duration=d, noise_std=noise)
2109 fish3 = 10.0*wavefish_eods('Alepto', eodfs[2], rate,
2110 duration=d, noise_std=noise)
2111 fish4 = 6.0*wavefish_eods('Alepto', eodfs[3], rate,
2112 duration=d, noise_std=noise)
2113 fish5 = 2.7*wavefish_eods('Alepto', eodfs[4], rate,
2114 duration=d, noise_std=noise)
2115 data = fish1 + fish2 + fish3 + fish4 + fish5
2116 else:
2117 from thunderlab.dataloader import load_data
2118 print(f"load {data_file} ...")
2119 data, rate, unit, amax = load_data(data_file)
2120 data = data[:, 0]
2121 title = data_file
2123 # retrieve fundamentals from power spectrum:
2124 psd_data = psd(data, rate, freq_resolution=0.1)
2125 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)
2126 fundamentals = fundamental_freqs(groups)
2127 np.set_printoptions(formatter={'float': lambda x: f'{x:5.1f}'})
2128 print('fundamental frequencies extracted from power spectrum:')
2129 print(fundamentals)
2130 fig, ax = plt.subplots()
2131 plot_psd_harmonic_groups(ax, psd_data[0], psd_data[1], groups, mains,
2132 all_freqs, good_freqs, max_freq=3000.0)
2133 annotate_harmonic_group(ax, groups[1])
2134 ax.set_title(title)
2135 plt.show()
2136 print()
2138 # add groups:
2139 group_list = [groups]
2140 group_list.append([mkg(44.0, 2.0), mkg(123.4, 17.2), mkg(320.5, 294.3), mkg(665.5, 486.1), mkg(666.2, 782.3), mkg(793.1, 74.9)])
2141 group_list.append([mkg(123.3, 163.8), mkg(320.2, 82.5), mkg(666.1, 184.9), mkg(668.4, 205.1)])
2142 group_list.append([mkg(123.1, 44.5), mkg(666.1, 86.4)])
2143 fundamentals = fundamental_freqs(group_list)
2144 print('fundamental frequencies of group_list:')
2145 print('\n'.join((str(f) for f in fundamentals)))
2146 print()
2148 # unify fundamental frequencies:
2149 freqs = fundamental_freqs_and_power(group_list)
2150 rank_freqs = add_relative_power(freqs)
2151 freqs = add_power_ranks(rank_freqs)
2152 print('all frequencies (frequency, power, relpower, rank):')
2153 print('\n'.join((str(f) for f in freqs)))
2154 print()
2155 indices = similar_indices(freqs, 1.0)
2156 print('similar indices:')
2157 for i in range(len(indices)):
2158 for inx in indices[i]:
2159 if len(inx) > 0:
2160 print(f'{i} {freqs[inx[0][0]][inx[0][1]][0]:6.1f}Hz:', inx)
2161 else:
2162 print(f'{i} :', inx)
2163 print()
2164 indices = unique_indices(freqs, 1.0)
2165 print('unique indices:')
2166 for inx in indices:
2167 print(f'{freqs[inx[0][0]][inx[0][1]][0]:6.1f}Hz:', inx)
2168 print()
2170 # consistent groups:
2171 cgroups = consistent(group_list, 1.0)
2172 print('fundamental frequencies of consistent groups:')
2173 print(fundamental_freqs(cgroups))
2174 print()
2175 fig, axs = plt.subplots(2, 2, layout='constrained', sharey=True)
2176 axs[0, 0].set_title('consistent groups')
2177 plot_selected_groups(axs[0, 0], group_list, cgroups)
2179 # check almost empty group_list:
2180 group2_list = [[], group_list[1], []]
2181 cgroups = consistent(group2_list, 1.0)
2182 axs[0, 1].set_title('consistent in mostly empty group list')
2183 plot_selected_groups(axs[0, 1], group2_list, cgroups)
2185 # check single group_list:
2186 group3_list = [group_list[1]]
2187 cgroups = consistent(group3_list, 1.0)
2188 axs[1, 0].set_title('consistent in single group list')
2189 plot_selected_groups(axs[1, 0], group3_list, cgroups)
2191 # closest groups:
2192 cgroups, cindices = closest(group_list, 1.0, 3, 2)
2193 print('closest frequencies:')
2194 print(fundamental_freqs(cgroups))
2195 print(cindices)
2196 axs[1, 1].set_title('closest groups')
2197 plot_selected_groups(axs[1, 1], group_list, cgroups, cindices)
2198 plt.show()
2199 print()
2201 fig, axs = plt.subplots(2, 2, layout='constrained')
2202 unique_freqs = unique(freqs, 1.0, 'power')
2203 axs[0, 0].set_title('unique power')
2204 plot_unique_frequencies(axs[0, 0], freqs, unique_freqs)
2205 print('unique power:')
2206 print('\n'.join((str(f) for f in unique_freqs)))
2207 print()
2209 unique_freqs = unique(freqs, 1.0, 'relpower')
2210 axs[0, 1].set_title('unique relative power')
2211 plot_unique_frequencies(axs[0, 1], freqs, unique_freqs)
2212 print('unique relative power:')
2213 print('\n'.join((str(f) for f in unique_freqs)))
2214 print()
2216 unique_freqs = unique(freqs, 1.0, 'rank')
2217 axs[1, 0].set_title('unique rank')
2218 plot_unique_frequencies(axs[1, 0], freqs, unique_freqs)
2219 print('unique rank:')
2220 print('\n'.join((str(f) for f in unique_freqs)))
2221 print()
2223 unique_freqs = unique(freqs, 1.0, 'rank', 1)
2224 axs[1, 1].set_title('unique rank for next neighbor')
2225 plot_unique_frequencies(axs[1, 1], freqs, unique_freqs)
2226 print('unique rank for next neighour only:')
2227 print('\n'.join((str(f) for f in unique_freqs)))
2228 print()
2229 plt.show()
2232if __name__ == "__main__":
2233 import sys
2234 data_file = sys.argv[1] if len(sys.argv) > 1 else None
2235 main(data_file)