Coverage for src / thunderfish / harmonics.py: 85%
657 statements
« prev ^ index » next coverage.py v7.13.1, created at 2026-01-15 17:50 +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 dB 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_mask()`: mark similar frequencies from different recordings as dublicate.
34- `unique()`: remove similar frequencies from different recordings.
36## Visualization
38- `colors_markers()`: Generate a list of colors and markers for plotting.
39- `plot_harmonic_groups()`: Mark decibel power of fundamentals and their
40 harmonics.
41- `plot_psd_harmonic_groups()`: Plot decibel power-spectrum with detected peaks,
42 harmonic groups, and mains frequencies.
44## Configuration
46- `add_psd_peak_detection_config()`: add parameters for the detection of
47 peaks in power spectra to configuration.
48- `psd_peak_detection_args()`: retrieve parameters for the detection of peaks
49 in power spectra from configuration.
50- `add_harmonic_groups_config()`: add parameters for the detection of
51 harmonic groups to configuration.
52- `harmonic_groups_args()`: retrieve parameters for the detection of
53 harmonic groups from configuration.
54"""
56import math as m
57import numpy as np
58import scipy.signal as sig
59try:
60 import matplotlib.cm as cm
61 import matplotlib.colors as mc
62except ImportError:
63 pass
65from thunderlab.eventdetection import detect_peaks, trim, hist_threshold
66from thunderlab.powerspectrum import decibel, power, plot_decibel_psd
69def group_candidate(good_freqs, all_freqs, freq, divisor,
70 freq_tol, max_freq_tol, min_group_size, verbose):
71 """Candidate harmonic frequencies belonging to a fundamental frequency.
73 Parameters
74 ----------
75 good_freqs: 2-D array
76 Frequency, power, and use count (columns) of strong peaks detected
77 in a power spectrum.
78 all_freqs: 2-D array
79 Frequency, power, and use count (columns) of all peaks detected
80 in a power spectrum.
81 freq: float
82 Fundamental frequency for which a harmonic group should be assembled.
83 divisor: int
84 Fundamental frequency was obtained from a frequency divided by divisor.
85 0: Fundamental frequency is given as is.
86 freq_tol: float
87 Harmonics should fall within this frequency tolerance.
88 This should be in the range of the frequency resolution
89 and should not be smaller than half of the frequency resolution.
90 max_freq_tol: float
91 Maximum deviation of harmonics from their expected frequency.
92 Peaks with frequencies between `freq_tol` and `max_freq_tol`
93 get penalized the further away from their expected frequency.
94 min_group_size: int
95 Minimum number of harmonics of a harmonic group.
96 The harmonics from min_group_size/3 to max(min_group_size, divisor)
97 need to be in good_freqs.
98 verbose: int
99 Verbosity level.
101 Returns
102 -------
103 new_group: 1-D aray of indices
104 Frequencies in all_freqs belonging to the candidate group.
105 fzero: float
106 Adjusted fundamental frequency. If negative, no group was found.
107 fzero_harmonics: int
108 The highest harmonics that was used to recompute the fundamental frequency.
109 """
110 if verbose > 0:
111 print('')
113 dvs = True
114 if divisor <= 0:
115 divisor = 1
116 dvs = False
118 # 1. find harmonics in good_freqs and adjust fzero accordingly:
119 fzero = freq
120 fzero_harmonics = 1
121 if len(good_freqs[:,0]) > 0:
122 """
123 ff = good_freqs[:,0]
124 h = np.round(ff/fzero)
125 h_indices = np.unique(h, return_index=True)[1]
126 f_best = [hi+np.argmin(np.abs(fh/h - fzero)) for hi, fh in zip(h_indices, np.split(ff, h_indices))]
127 fe = np.abs(ff/h - fzero)
128 h = h[fe < freq_tol]
129 fe = fe[fe < freq_tol]
130 """
131 prev_freq = divisor * freq
132 for h in range(divisor+1, 2*min_group_size+1):
133 idx = np.argmin(np.abs(good_freqs[:,0]/h - fzero))
134 ff = good_freqs[idx,0]
135 if m.fabs(ff/h - fzero) > freq_tol:
136 continue
137 df = ff - prev_freq
138 dh = np.round(df/fzero)
139 fe = m.fabs(df/dh - fzero)
140 if fe > 2.0*freq_tol:
141 if h > min_group_size:
142 break
143 continue
144 # update fzero:
145 prev_freq = ff
146 fzero_harmonics = h
147 fzero = ff/fzero_harmonics
148 if verbose > 1:
149 print('adjusted fzero to %.2fHz from %d-th harmonic' % (fzero, h))
151 if verbose > 0:
152 ds = 'divisor: %d, ' % divisor if dvs else ''
153 print('# %sfzero=%7.2fHz adjusted from harmonics %d'
154 % (ds, fzero, fzero_harmonics))
156 # 2. check fzero:
157 # freq might not be in our group anymore, because fzero was adjusted:
158 if np.abs(freq - fzero) > freq_tol:
159 if verbose > 0:
160 print(' discarded: lost frequency')
161 return [], -1.0, fzero_harmonics
163 # 3. collect harmonics from all_freqs:
164 new_group = -np.ones(min_group_size, dtype=int)
165 new_penalties = np.ones(min_group_size)
166 freqs = []
167 prev_h = 0
168 prev_fe = 0.0
169 for h in range(1, min_group_size+1):
170 penalty = 0
171 i = np.argmin(np.abs(all_freqs[:,0]/h - fzero))
172 f = all_freqs[i,0]
173 if verbose > 2:
174 print(' check %7.2fHz as %d. harmonics' % (f, h))
175 fac = 1.0 if h >= divisor else 2.0
176 fe = m.fabs(f/h - fzero)
177 if fe > fac*max_freq_tol:
178 if verbose > 1 and fe < 2*fac*max_freq_tol:
179 print(' %d. harmonics at %7.2fHz is off by %7.2fHz (max between %5.2fHz and %5.2fHz) from %7.2fHz'
180 % (h, f, h*fe, h*fac*freq_tol, h*fac*max_freq_tol, h*fzero))
181 continue
182 if fe > fac*freq_tol:
183 penalty = np.interp(fe, [fac*freq_tol, fac*max_freq_tol], [0.0, 1.0])
184 if len(freqs) > 0:
185 pf = freqs[-1]
186 df = f - pf
187 if df < 0.5*fzero:
188 if len(freqs)>1:
189 pf = freqs[-2]
190 df = f - pf
191 else:
192 pf = 0.0
193 df = h*fzero
194 dh = m.floor(df/fzero + 0.5)
195 fe = m.fabs(df/dh - fzero)
196 if fe > 2*dh*fac*max_freq_tol:
197 if verbose > 1:
198 print(' %d. harmonics at %7.2fHz is off by %7.2fHz (max between %5.2fHz and %5.2fHz) from previous harmonics at %7.2fHz'
199 % (h, f, dh*fe, 2*fac*dh*freq_tol, 2*fac*dh*max_freq_tol, pf))
200 continue
201 if fe > 2*dh*fac*freq_tol:
202 penalty = np.interp(fe, [2*dh*fac*freq_tol, 2*dh*fac*max_freq_tol], [0.0, 1.0])
203 else:
204 fe = 0.0
205 if h > prev_h or fe < prev_fe:
206 if prev_h > 0 and h - prev_h > 1:
207 if verbose > 1:
208 print(' previous harmonics %d more than 1 away from %d. harmonics at %7.2fHz'
209 % (prev_h, h, f))
210 break
211 if h == prev_h and len(freqs) > 0:
212 freqs.pop()
213 freqs.append(f)
214 new_group[int(h)-1] = i
215 new_penalties[int(h)-1] = penalty
216 prev_h = h
217 prev_fe = fe
218 if verbose > 1:
219 print(' %d. harmonics at %7.2fHz has been taken (peak %2d) with penalty %3.1f' % (h, f, i, penalty))
221 # 4. check new group:
223 # almost all harmonics in min_group_size required:
224 max_penalties = min_group_size/3
225 if np.sum(new_penalties) > max_penalties:
226 if verbose > 0:
227 print(' discarded group because sum of penalties %3.1f is more than %3.1f: indices' %
228 (np.sum(new_penalties), max_penalties), new_group, ' penalties', new_penalties)
229 return [], -1.0, fzero_harmonics
230 new_group = new_group[new_group>=0]
232 # check use count of frequencies:
233 double_use = np.sum(all_freqs[new_group, 2]>0)
234 if double_use >= 2: # XXX make this a parameter?
235 if verbose > 0:
236 print(' discarded group because of use count = %2d >= 2' % double_use)
237 return [], -1.0, fzero_harmonics
239 # 5. return results:
240 return new_group, fzero, fzero_harmonics
243def update_group(good_freqs, all_freqs, new_group, fzero,
244 freq_tol, verbose, group_str):
245 """Update good frequencies and harmonic group.
247 Remove frequencies from good_freqs, add missing fundamental to group.
249 Parameters
250 ----------
251 good_freqs: 2-D array
252 Frequency, power, and use count (columns) of strong peaks detected
253 in a power spectrum.
254 all_freqs: 2-D array
255 Frequency, power, and use count (columns) of all peaks detected
256 in a power spectrum.
257 new_group: 1-D aray of indices
258 Frequencies in all_freqs of an harmonic group.
259 fzero: float
260 Fundamental frequency for which frequencies are collected in good_freqs.
261 freq_tol: float
262 Harmonics need to fall within this frequency tolerance.
263 This should be in the range of the frequency resolution
264 and not be smaller than half of the frequency resolution.
265 verbose: int
266 Verbosity level.
267 group_str: string
268 String for debug message.
270 Returns
271 -------
272 good_freqs: 2-D array
273 Frequency, power, and use count (columns) of strong peaks detected
274 in a power spectrum with frequencies for harmonic group
275 of fundamental frequency fzero removed.
276 group: 2-D array
277 Frequency, power, and use count (columns) of harmonic group
278 for fundamental frequency fzero.
279 """
280 # initialize group:
281 group = all_freqs[new_group,:]
283 # indices of group in good_freqs:
284 freq_tol *= 1.1
285 indices = []
286 for f in group[:,0]:
287 idx = np.argmin(np.abs(good_freqs[:,0]-f))
288 if np.abs(good_freqs[idx,0]-f) <= freq_tol:
289 indices.append(idx)
290 indices = np.asarray(indices, dtype=int)
292 # harmonics in good_freqs:
293 nharm = np.round(good_freqs[:,0]/fzero)
294 idxs = np.where(np.abs(good_freqs[:,0] - nharm*fzero) <= freq_tol)[0]
295 indices = np.unique(np.concatenate((indices, idxs)))
297 # report:
298 if verbose > 1:
299 print('# good freqs: ', indices,
300 '[', ', '.join(['%.2f' % f for f in good_freqs[indices,0]]), ']')
301 print('# all freqs : ', new_group,
302 '[', ', '.join(['%.2f' % f for f in all_freqs[new_group,0]]), ']')
303 if verbose > 0:
304 refi = np.argmax(group[:,1] > 0.0)
305 print('')
306 print(group_str)
307 for i in range(len(group)):
308 print('f=%8.2fHz n=%5.2f: power=%9.3g power/pmax=%6.4f=%5.1fdB'
309 % (group[i,0], group[i,0]/fzero,
310 group[i,1], group[i,1]/group[refi,1], decibel(group[i,1], group[refi,1])))
311 print('')
313 # erase group from good_freqs:
314 good_freqs = np.delete(good_freqs, indices, axis=0)
316 # adjust frequencies to fzero:
317 group[:,0] = np.round(group[:,0]/fzero)*fzero
319 # insert missing fzero:
320 if np.round(group[0,0]/fzero) != 1.0:
321 group = np.vstack(((fzero, group[0,1], -2.0), group))
323 return good_freqs, group
326def build_harmonic_group(good_freqs, all_freqs, freq_tol, max_freq_tol,
327 verbose=0, min_group_size=3, max_divisor=4):
328 """Find all the harmonics belonging to the largest peak in a list of frequency peaks.
330 Parameters
331 ----------
332 good_freqs: 2-D array
333 Frequency, power, and use count (columns) of strong peaks detected
334 in a power spectrum.
335 all_freqs: 2-D array
336 Frequency, power, and use count (columns) of all peaks detected
337 in a power spectrum.
338 freq_tol: float
339 Harmonics should fall within this frequency tolerance.
340 This should be in the range of the frequency resolution
341 and should not be smaller than half of the frequency resolution.
342 max_freq_tol: float
343 Maximum deviation of harmonics from their expected frequency.
344 Peaks with frequencies between `freq_tol` and `max_freq_tol`
345 get penalized the further away from their expected frequency.
346 verbose: int
347 Verbosity level.
348 min_group_size: int
349 Minimum number of harmonics of a harmonic group.
350 The harmonics from min_group_size/3 to max(min_group_size, divisor)
351 need to be in good_freqs.
352 max_divisor: int
353 Maximum divisor used for checking for sub-harmonics.
355 Returns
356 -------
357 good_freqs: 2-D array
358 Frequency, power, and use count (columns) of strong peaks detected
359 in a power spectrum with frequencies of the returned harmonic group removed.
360 group: 2-D array
361 The detected harmonic group. Might be empty.
362 indices: 1-D array of indices
363 Indices of the harmonic group in all_freqs.
364 best_fzero_harmonics: int
365 The highest harmonics that was used to recompute
366 the fundamental frequency.
367 fmax: float
368 The frequency of the largest peak in good_freqs
369 for which the harmonic group was detected.
370 """
371 # select strongest frequency for building the harmonic group:
372 fmaxinx = np.argmax(good_freqs[:,1])
373 fmax = good_freqs[fmaxinx,0]
374 if verbose > 0:
375 print('')
376 print('%s build_harmonic_group %s' % (10*'#', 38*'#'))
377 print('%s fmax=%7.2fHz, power=%9.3g %s'
378 % (10*'#', fmax, good_freqs[fmaxinx,1], 27*'#'))
379 print('good_freqs: ', '[',
380 ', '.join(['%.2f' % f for f in good_freqs[:,0]]), ']')
382 # container for harmonic groups:
383 best_group = []
384 best_value = -1e6
385 best_divisor = 0
386 best_fzero = 0.0
387 best_fzero_harmonics = 0
389 # check for integer fractions of the frequency:
390 for divisor in range(1, max_divisor + 1):
391 # 1. hypothesized fundamental:
392 freq = fmax / divisor
394 # 2. find harmonics in good_freqs and adjust fzero accordingly:
395 group_size = min_group_size if divisor <= min_group_size else divisor
396 new_group, fzero, fzero_harmonics = group_candidate(good_freqs, all_freqs,
397 freq, divisor,
398 freq_tol, max_freq_tol,
399 group_size, verbose)
400 # no group found:
401 if fzero < 0.0:
402 continue
404 # 3. compare new group to best group:
405 peaksum = decibel(np.sum(all_freqs[new_group, 1])*min_group_size/len(new_group))
406 diff = np.std(np.diff(decibel(all_freqs[new_group, 1])))
407 new_group_value = peaksum - diff
408 counts = np.sum(all_freqs[new_group, 2])
409 if verbose > 0:
410 print(' new group: fzero=%7.2fHz, nharmonics=%d, value=%6.1fdB, peaksum=%5.1fdB, diff=%6.1fdB, use count=%2d, peaks:'
411 % (fzero, len(new_group), new_group_value, peaksum, diff, counts), new_group)
412 if verbose > 1:
413 print(' best group: divisor=%d, fzero=%7.2fHz, nharmonics=%d, value=%6.1fdB, peaks:'
414 % (best_divisor, best_fzero, len(best_group), best_value), best_group)
415 # select new group if sum of peak power minus diff is larger:
416 if len(new_group) >= len(best_group) and new_group_value >= best_value:
417 best_value = new_group_value
418 best_group = new_group
419 best_divisor = divisor
420 best_fzero = fzero
421 best_fzero_harmonics = fzero_harmonics
422 if verbose > 1:
423 print(' new best group: divisor=%d, fzero=%7.2fHz, value=%6.1fdB, peaks:'
424 % (best_divisor, best_fzero, best_value), best_group)
425 elif verbose > 0:
426 print(' took as new best group')
428 # no group found:
429 if len(best_group) == 0:
430 # erase freq:
431 good_freqs = np.delete(good_freqs, fmaxinx, axis=0)
432 group = np.zeros((0, 3))
433 return good_freqs, group, best_group, 1, fmax
435 # update frequencies and group:
436 if verbose > 1:
437 print('')
438 print('# best group found for fmax=%.2fHz, fzero=%.2fHz, divisor=%d:'
439 % (fmax, best_fzero, best_divisor))
440 group_str = '%s resulting harmonic group for fmax=%.2fHz' % (10*'#', fmax)
441 good_freqs, group = update_group(good_freqs, all_freqs, best_group, best_fzero, freq_tol, verbose, group_str)
443 # good_freqs: removed all frequencies of bestgroup
444 return good_freqs, group, best_group, best_fzero_harmonics, fmax
447def retrieve_harmonic_group(freq, good_freqs, all_freqs,
448 freq_tol, max_freq_tol, verbose=0,
449 min_group_size=3):
450 """Find all the harmonics belonging to a given fundamental.
452 Parameters
453 ----------
454 freq: float
455 Fundamental frequency for which harmonics are to be retrieved.
456 good_freqs: 2-D array
457 Frequency, power, and use count (columns) of strong peaks detected
458 in a power spectrum. All harmonics of `freq` will be
459 removed from `good_freqs`.
460 all_freqs: 2-D array
461 Frequency, power, and use count (columns) of all peaks detected
462 in a power spectrum.
463 freq_tol: float
464 Harmonics should fall within this frequency tolerance.
465 This should be in the range of the frequency resolution
466 and should not be smaller than half of the frequency resolution.
467 max_freq_tol: float
468 Maximum deviation of harmonics from their expected frequency.
469 Peaks with frequencies between `freq_tol` and `max_freq_tol`
470 get penalized the further away from their expected frequency.
471 verbose: int
472 Verbosity level.
473 min_group_size: int
474 Minimum number of harmonics of a harmonic group.
475 The harmonics from min_group_size/3 to max(min_group_size, divisor)
476 need to be in good_freqs.
478 Returns
479 -------
480 good_freqs: 2-D array
481 Frequency, power, and use count (columns) of strong peaks detected
482 in a power spectrum with frequencies of the returned harmonic group removed.
483 group: 2-D array
484 The detected harmonic group. Might be empty.
485 indices: 1-D array of indices
486 Indices of the harmonic group in all_freqs.
487 fzero_harmonics: int
488 The highest harmonics that was used to recompute
489 the fundamental frequency.
490 """
491 if verbose > 0:
492 print('')
493 print('%s retrieve harmonic group %s' % (10*'#', 35*'#'))
494 print('%s freq=%7.2fHz %s' % (10*'#', freq, 44*'#'))
495 print('good_freqs: ', '[',
496 ', '.join(['%.2f' % f for f in good_freqs[:,0]]), ']')
498 # find harmonics in good_freqs and adjust fzero accordingly:
499 new_group, fzero, fzero_harmonics = group_candidate(good_freqs, all_freqs,
500 freq, 0,
501 freq_tol, max_freq_tol,
502 min_group_size, verbose)
504 # no group found:
505 if fzero < 0.0:
506 return good_freqs, np.zeros((0, 2)), np.zeros(0), fzero_harmonics
508 if verbose > 1:
509 print('')
510 print('# group found for freq=%.2fHz, fzero=%.2fHz:'
511 % (freq, fzero))
512 group_str = '#### resulting harmonic group for freq=%.2fHz' % freq
513 good_freqs, group = update_group(good_freqs, all_freqs, new_group,
514 fzero, freq_tol, verbose, group_str)
516 # good_freqs: removed all frequencies of bestgroup
517 return good_freqs, group, new_group, fzero_harmonics
520def expand_group(group, indices, freqs, freq_tol, max_harmonics=0):
521 """Add more harmonics to harmonic group.
523 Parameters
524 ----------
525 group: 2-D array
526 Group of fundamental frequency and harmonics
527 as returned by build_harmonic_group.
528 indices: 1-D array of indices
529 Indices of the harmonics in group in `freqs`.
530 freqs: 2-D array
531 Frequency, power, and use count (columns) of all peaks detected
532 in a power spectrum.
533 freq_tol: float
534 Harmonics need to fall within this frequency tolerance.
535 This should be in the range of the frequency resolution
536 and not be smaller than half of the frequency resolution.
537 max_harmonics: int
538 Maximum number of harmonics to be returned for each group.
540 Returns
541 -------
542 group: 2-D array
543 Expanded group of fundamental frequency and harmonics.
544 indices: 1-D array of indices
545 Indices of the harmonics in the expanded group in all_freqs.
546 """
547 if len(group) == 0:
548 return group, indices
549 fzero = group[0,0]
550 if max_harmonics <= 0:
551 max_harmonics = m.floor(freqs[-1,0]/fzero + 0.5) # round
552 if max_harmonics <= len(group):
553 return group, indices
554 group_freqs = list(group[:,0])
555 indices = list(indices)
556 last_h = m.floor(group_freqs[-1]/fzero + 0.5) # round
557 for h in range(last_h+1, max_harmonics+1):
558 i = np.argmin(np.abs(freqs[:,0]/h - fzero))
559 f = freqs[i,0]
560 if m.fabs(f/h - fzero) > freq_tol:
561 continue
562 df = f - group_freqs[-1]
563 dh = m.floor(df/fzero + 0.5) # round
564 if dh < 1:
565 continue
566 fe = m.fabs(df/dh - fzero)
567 if fe > 2*freq_tol:
568 continue
569 group_freqs.append(f)
570 indices.append(i)
571 # assemble group:
572 new_group = freqs[indices,:group.shape[1]]
573 # keep filled in fundamental:
574 if group[0,2] == -2:
575 new_group = np.vstack((group[0,:], new_group))
576 return new_group, np.array(indices, dtype=int)
579def extract_fundamentals(good_freqs, all_freqs, freq_tol, max_freq_tol,
580 verbose=0, check_freqs=[],
581 mains_freq=60.0, mains_freq_tol=1.0,
582 min_freq=0.0, max_freq=2000.0, max_db_diff=20.0,
583 max_divisor=4, min_group_size=3, max_harmonics_db=-5.0,
584 max_harmonics=0, max_groups=0, **kwargs):
585 """Extract fundamental frequencies from power-spectrum peaks.
587 Parameters
588 ----------
589 good_freqs: 2-D array
590 Frequency, power, and use count (columns) of strong peaks detected
591 in a power spectrum.
592 all_freqs: 2-D array
593 Frequency, power, and use count (columns) of all peaks detected
594 in a power spectrum.
595 freq_tol: float
596 Harmonics should fall within this frequency tolerance.
597 This should be in the range of the frequency resolution
598 and should not be smaller than half of the frequency resolution.
599 max_freq_tol: float
600 Maximum deviation of harmonics from their expected frequency.
601 Peaks with frequencies between `freq_tol` and `max_freq_tol`
602 get penalized the further away from their expected frequency.
603 verbose: int
604 Verbosity level.
605 check_freqs: list of float
606 List of fundamental frequencies that will be checked
607 first for being present and valid harmonic groups in the peak frequencies
608 of a power spectrum.
609 mains_freq: float
610 Frequency of the mains power supply.
611 mains_freq_tol: float
612 Tolarerance around harmonics of the mains frequency,
613 within which peaks are removed.
614 min_freq: float
615 Minimum frequency accepted as a fundamental frequency.
616 max_freq: float
617 Maximum frequency accepted as a fundamental frequency.
618 max_db_diff: float
619 If larger than zero, maximum standard deviation of differences between
620 logarithmic powers of harmonics in decibel.
621 Low values enforce smoother power spectra.
622 max_divisor: int
623 Maximum divisor used for checking for sub-harmonics.
624 min_group_size: int
625 Minimum number of harmonics of a harmonic group.
626 The harmonics from min_group_size/3 to max(min_group_size, divisor)
627 need to be in good_freqs.
628 max_harmonics_db: float
629 Maximum allowed power of the `min_group_size`-th and higher harmonics
630 after the peak (in decibel relative to peak power withn the first
631 `min_group_size` harmonics, i.e. if harmonics are required to be
632 smaller than fundamental then this is a negative number).
633 Make it a large positive number to effectively not check for relative power.
634 max_harmonics: int
635 Maximum number of harmonics to be returned for each group.
636 max_groups: int
637 If not zero the maximum number of most powerful harmonic groups.
639 Returns
640 -------
641 group_list: list of 2-D arrays
642 List of all harmonic groups found sorted by fundamental frequency.
643 Each harmonic group is a 2-D array with the first dimension the harmonics
644 and the second dimension containing frequency and power of each harmonic.
645 If the power is zero, there was no corresponding peak in the power spectrum.
646 fzero_harmonics_list: list of int
647 The harmonics from which the fundamental frequencies were computed.
648 mains_freqs: 2-d array
649 Array of mains peaks found in all_freqs (frequency, power).
650 """
651 if verbose > 0:
652 print('')
654 # set use count to zero:
655 all_freqs[:,2] = 0.0
657 # remove power line harmonics from good_freqs:
658 if mains_freq > 0.0:
659 indices = np.where(np.abs(good_freqs[:,0] - np.round(good_freqs[:,0]/mains_freq)*mains_freq) < mains_freq_tol)[0]
660 if len(indices)>0:
661 if verbose > 1:
662 print('remove power line frequencies',
663 ', '.join(['%.1f' % f for f in good_freqs[indices,0]]))
664 good_freqs = np.delete(good_freqs, indices, axis=0)
666 if verbose > 1:
667 print('all_freqs: ', '[',
668 ', '.join(['%.2f' % f for f in all_freqs[:,0] if f < 3000.0]), ']')
670 group_list = []
671 fzero_harmonics_list = []
672 first = True
673 # as long as there are frequencies left in good_freqs:
674 fi = 0
675 while len(good_freqs) > 0:
676 if fi < len(check_freqs):
677 # check for harmonic group of a given fundamental frequency:
678 fmax = check_freqs[fi]
679 f0s = 'freq'
680 good_freqs, harm_group, harm_indices, fzero_harmonics = \
681 retrieve_harmonic_group(fmax, good_freqs, all_freqs,
682 freq_tol, max_freq_tol,
683 verbose-1, min_group_size)
684 fi += 1
685 else:
686 # check for harmonic groups:
687 f0s = 'fmax'
688 good_freqs, harm_group, harm_indices, fzero_harmonics, fmax = \
689 build_harmonic_group(good_freqs, all_freqs,
690 freq_tol, max_freq_tol, verbose-1,
691 min_group_size, max_divisor)
693 # nothing found:
694 if len(harm_group) == 0:
695 if verbose > 1 - first:
696 s = ' (largest peak)' if first else ''
697 print('No harmonic group for %7.2fHz%s' % (fmax, s))
698 first = False
699 continue
700 first = False
702 # fill up harmonic group:
703 harm_group, harm_indices = expand_group(harm_group, harm_indices,
704 all_freqs, freq_tol,
705 max_harmonics)
707 # check whether fundamental was filled in:
708 first_h = 0 if harm_group[0,2] > -2 else 1
710 # increment use count:
711 all_freqs[harm_indices,2] += 1
712 if harm_group[0,2] == -1:
713 harm_group[0,2] = -2
715 # check frequency range of fundamental:
716 fundamental_ok = (harm_group[0, 0] >= min_freq and
717 harm_group[0, 0] <= max_freq)
718 # check power hum:
719 mains_ok = ((mains_freq <= 0.0) or
720 (m.fabs(harm_group[0,0] - mains_freq) > freq_tol))
721 # check smoothness:
722 db_powers = decibel(harm_group[first_h:,1])
723 diff = np.std(np.diff(db_powers))
724 smooth_ok = max_db_diff <= 0.0 or diff < max_db_diff
725 # check relative power of higher harmonics:
726 p_max = np.argmax(db_powers[:min_group_size])
727 db_powers -= db_powers[p_max]
728 amplitude_ok = len(db_powers[p_max+min_group_size:]) == 0
729 if amplitude_ok:
730 pi = len(db_powers) - 1
731 else:
732 amplitude_ok = np.all((db_powers[p_max+min_group_size:] < max_harmonics_db) |
733 (harm_group[first_h+p_max+min_group_size:,2] > 1))
734 if amplitude_ok:
735 pi = p_max+min_group_size-1
736 else:
737 pi = np.where((db_powers[p_max+min_group_size:] >= max_harmonics_db) &
738 (harm_group[first_h+p_max+min_group_size:,2] <= 1))[0][0]
740 # check:
741 if fundamental_ok and mains_ok and smooth_ok and amplitude_ok:
742 if verbose > 0:
743 print('Accepting harmonic group from %s=%7.2fHz: %7.2fHz power=%9.3g nharmonics=%2d, use count=%d'
744 % (f0s, fmax, harm_group[0,0], np.sum(harm_group[first_h:,1]),
745 len(harm_group), np.sum(harm_group[first_h:,2])))
746 group_list.append(harm_group[:,0:2])
747 fzero_harmonics_list.append(fzero_harmonics)
748 else:
749 if verbose > 0:
750 fs = 'is ' if fundamental_ok else 'NOT'
751 ms = 'not ' if mains_ok else 'IS'
752 ss = 'smaller' if smooth_ok else 'LARGER '
753 ps = 'smaller' if amplitude_ok else 'LARGER '
754 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'
755 % (f0s, fmax, harm_group[0,0], np.sum(harm_group[first_h:,1]),
756 fs, ms, diff, ss, max_db_diff,
757 pi, db_powers[pi], ps, max_harmonics_db))
759 # select most powerful harmonic groups:
760 if max_groups > 0 and len(group_list) > max_groups:
761 n = len(group_list)
762 powers = [np.sum(group[:,1]) for group in group_list]
763 powers_inx = np.argsort(powers)
764 group_list = [group_list[pi] for pi in powers_inx[-max_groups:]]
765 fzero_harmonics_list = [fzero_harmonics_list[pi] for pi in powers_inx[-max_groups:]]
766 if verbose > 0:
767 print('Selected from %d groups the %d most powerful groups.' % (n, max_groups))
769 # sort groups by fundamental frequency:
770 freqs = [group[0,0] for group in group_list]
771 freq_inx = np.argsort(freqs)
772 group_list = [group_list[fi] for fi in freq_inx]
773 fzero_harmonics_list = [fzero_harmonics_list[fi] for fi in freq_inx]
775 if verbose > 1:
776 print('')
777 if len(group_list) > 0:
778 print('## FUNDAMENTALS FOUND: ##')
779 for group, fzero_h in zip(group_list, fzero_harmonics_list):
780 print('%7.2fHz: power=%9.3g fzero-h=%2d'
781 % (group[0,0], np.sum(group[:,1]), fzero_h))
782 else:
783 print('## NO FUNDAMENTALS FOUND ##')
785 # assemble mains frequencies from all_freqs:
786 if mains_freq > 0.0:
787 mains_freqs = all_freqs[np.abs(all_freqs[:,0] - np.round(all_freqs[:,0]/mains_freq)*mains_freq) < mains_freq_tol,:2]
788 else:
789 mains_freqs = np.zeros((0, 2))
791 return group_list, fzero_harmonics_list, mains_freqs
794def threshold_estimate(psd_data, low_thresh_factor=6.0, high_thresh_factor=10.0,
795 nbins=100, hist_height=1.0/ np.sqrt(np.e)):
796 """Estimate thresholds for peak detection from histogram of power spectrum.
798 The standard deviation of the noise floor without peaks is estimated from
799 the width of the histogram of the power spectrum at `hist_height` relative height.
800 The histogram is computed in the third quarter of the linearly detrended power spectrum.
802 Parameters
803 ----------
804 psd_data: 1-D array
805 The power spectrum from which to estimate the thresholds.
806 low_thresh_factor: float
807 Factor by which the estimated standard deviation of the noise floor
808 is multiplied to set the `low_threshold`.
809 high_thresh_factor: float
810 Factor by which the estimated standard deviation of the noise floor
811 is multiplied to set the `high_threshold`.
812 nbins: int or list of floats
813 Number of bins or the bins for computing the histogram.
814 hist_height: float
815 Height between 0 and 1 at which the standard deviation of the histogram is estimated.
817 Returns
818 -------
819 low_threshold: float
820 The threshold for peaks just above the noise floor.
821 high_threshold: float
822 The threshold for distinct peaks.
823 center: float
824 The baseline level of the power spectrum.
825 """
826 n = len(psd_data)
827 psd_data_seg = psd_data[n//2:n*3//4]
828 psd_data_seg = psd_data_seg[~np.isinf(psd_data_seg)]
829 psd_data_seg = np.mean(psd_data_seg) + \
830 sig.detrend(psd_data_seg, type='linear')
831 noise_std, center = hist_threshold(psd_data_seg, thresh_fac=1.0, nbins=nbins)
832 low_threshold = noise_std * low_thresh_factor
833 high_threshold = noise_std * high_thresh_factor
834 return low_threshold, high_threshold, center
837def harmonic_groups(psd_freqs, psd, verbose=0, check_freqs=[],
838 low_threshold=0.0, high_threshold=0.0, thresh_bins=100,
839 low_thresh_factor=6.0, high_thresh_factor=10.0,
840 freq_tol_fac=1.0, max_freq_tol=1.0,
841 mains_freq=60.0, mains_freq_tol=1.0,
842 min_freq=0.0, max_freq=2000.0, max_db_diff=20.0, max_divisor=4,
843 min_group_size=3, max_harmonics_db=-5.0,
844 max_harmonics=0, max_groups=0, **kwargs):
845 """Detect peaks in power spectrum and group them according to their harmonic structure.
847 Parameters
848 ----------
849 psd_freqs: 1-D array
850 Frequencies of the power spectrum.
851 psd: 1-D array
852 Power spectrum (linear, not decible).
853 verbose: int
854 Verbosity level.
855 check_freqs: list of float
856 List of fundamental frequencies that will be checked first for being
857 present and valid harmonic groups in the peak frequencies
858 of the power spectrum.
859 low_threshold: float
860 The relative threshold for detecting all peaks in the decibel spectrum.
861 high_threshold: float
862 The relative threshold for detecting good peaks in the decibel spectrum.
863 thresh_bins: int or list of floats
864 Number of bins or the bins for computing the histogram from
865 which the standard deviation of the noise level in the `psd` is estimated.
866 low_thresh_factor: float
867 Factor by which the estimated standard deviation of the noise floor
868 is multiplied to set the `low_threshold`.
869 high_thresh_factor: float
870 Factor by which the estimated standard deviation of the noise floor
871 is multiplied to set the `high_threshold`.
872 freq_tol_fac: float
873 Harmonics should fall within `deltaf*freq_tol_fac`.
874 max_freq_tol: float
875 Maximum absolute frequency deviation of harmonics in Hertz..
876 mains_freq: float
877 Frequency of the mains power supply.
878 mains_freq_tol: float
879 Tolarerance around harmonics of the mains frequency,
880 within which peaks are removed.
881 min_freq: float
882 Minimum frequency accepted as a fundamental frequency.
883 max_freq: float
884 Maximum frequency accepted as a fundamental frequency.
885 max_db_diff: float
886 If larger than zero, maximum standard deviation of differences between
887 logarithmic powers of harmonics in decibel (larger than zero).
888 Low values enforce smoother power spectra.
889 max_divisor: int
890 Maximum divisor used for checking for sub-harmonics.
891 min_group_size: int
892 Minimum number of harmonics of a harmonic group.
893 The harmonics from min_group_size/3 to max(min_group_size, divisor)
894 need to be in good_freqs.
895 max_harmonics_db: float
896 Maximum allowed power of the `min_group_size`-th and higher harmonics
897 after the peak (in decibel relative to peak power withn the first
898 `min_group_size` harmonics, i.e. if harmonics are required to be
899 smaller than fundamental then this is a negative number).
900 Make it a large positive number to effectively not check for relative power.
901 max_harmonics: int
902 Maximum number of harmonics to be returned for each group.
903 max_groups: int
904 If not zero the maximum number of most powerful harmonic groups.
906 Returns
907 -------
908 groups: list of 2-D arrays
909 List of all extracted harmonic groups, sorted by fundamental frequency.
910 Each harmonic group is a 2-D array with the first dimension the harmonics
911 and the second dimension containing frequency and power of each harmonic.
912 If the power is zero, there was no corresponding peak in the power spectrum.
913 fzero_harmonics: list of ints
914 The harmonics from which the fundamental frequencies were computed.
915 mains: 2-d array
916 Frequencies and power of multiples of the mains frequency found in the power spectrum.
917 all_freqs: 2-D array
918 Frequency, power, and use count (columns) of all peaks detected
919 in the power spectrum.
920 good_freqs: 1-D array
921 Frequencies of strong peaks detected in the power spectrum.
922 low_threshold: float
923 The relative threshold for detecting all peaks in the decibel spectrum.
924 high_threshold: float
925 The relative threshold for detecting good peaks in the decibel spectrum.
926 center: float
927 The baseline level of the power spectrum.
928 """
929 if verbose > 0:
930 print('')
931 if verbose > 1:
932 print(70*'#')
933 print('##### harmonic_groups', 48*'#')
935 # decibel power spectrum:
936 log_psd = decibel(psd)
937 max_idx = np.argmax(~np.isfinite(log_psd))
938 if max_idx > 0:
939 log_psd = log_psd[:max_idx]
940 delta_f = psd_freqs[1] - psd_freqs[0]
942 # thresholds:
943 center = np.nan
944 if low_threshold <= 0.0 or high_threshold <= 0.0:
945 low_th, high_th, center = threshold_estimate(log_psd, low_thresh_factor,
946 high_thresh_factor,
947 thresh_bins)
948 if low_threshold <= 0.0:
949 low_threshold = low_th
950 if high_threshold <= 0.0:
951 high_threshold = high_th
952 if verbose > 1:
953 print('')
954 print('low_threshold =%4.1fdB, center+low_threshold =%6.1fdB' % (low_threshold, center + low_threshold))
955 print('high_threshold=%4.1fdB, center+high_threshold=%6.1fdB' % (high_threshold, center + high_threshold))
956 print(' center=%6.1fdB' % center)
958 # detect peaks in decibel power spectrum:
959 peaks, troughs = detect_peaks(log_psd, low_threshold)
960 peaks, troughs = trim(peaks, troughs)
961 all_freqs = np.zeros((len(peaks), 3))
962 all_freqs[:,0] = psd_freqs[peaks]
963 all_freqs[:,1] = psd[peaks]
965 if len(all_freqs) == 0:
966 return [], [], [], np.zeros((0, 3)), [], low_threshold, high_threshold, center
968 # select good peaks:
969 good_freqs = all_freqs[(log_psd[peaks] - log_psd[troughs] > high_threshold) &
970 (all_freqs[:,0] >= min_freq) &
971 (all_freqs[:,0] < max_freq*max_divisor),:]
973 # detect harmonic groups:
974 freq_tol = delta_f*freq_tol_fac
975 if max_freq_tol < 1.1*freq_tol:
976 max_freq_tol = 1.1*freq_tol
977 groups, fzero_harmonics, mains = \
978 extract_fundamentals(good_freqs, all_freqs,
979 freq_tol, max_freq_tol,
980 verbose, check_freqs, mains_freq, mains_freq_tol,
981 min_freq, max_freq, max_db_diff, max_divisor, min_group_size,
982 max_harmonics_db, max_harmonics, max_groups)
984 return (groups, fzero_harmonics, mains, all_freqs, good_freqs[:,0],
985 low_threshold, high_threshold, center)
988def fundamental_freqs(group_list):
989 """
990 Extract fundamental frequencies from lists of harmonic groups.
992 The inner list of 2-D arrays of the input argument is transformed into
993 a 1-D array containig the fundamental frequencies extracted from
994 the 2-D arrays.
996 Parameters
997 ----------
998 group_list: (list of (list of ...)) list of 2-D arrays
999 Arbitrarily nested lists of harmonic groups as returned by
1000 extract_fundamentals() and harmonic_groups() with the element
1001 [0, 0] being the fundamental frequency.
1003 Returns
1004 -------
1005 fundamentals: (list of (list of ...)) 1-D array
1006 Nested list (corresponding to `group_list`) of 1-D arrays
1007 with the fundamental frequencies extracted from the harmonic groups.
1008 """
1009 if len(group_list) == 0:
1010 return np.zeros(0)
1012 # check whether group_list is list of harmonic groups:
1013 list_of_groups = True
1014 for group in group_list:
1015 if not ( hasattr(group, 'shape') and len(group.shape) == 2 ):
1016 list_of_groups = False
1017 break
1019 if list_of_groups:
1020 fundamentals = np.array([group[0, 0] for group in group_list if len(group) > 0])
1021 else:
1022 fundamentals = []
1023 for groups in group_list:
1024 f = fundamental_freqs(groups)
1025 fundamentals.append(f)
1026 return fundamentals
1029def fundamental_freqs_and_power(group_list, power=False,
1030 ref_power=1.0, min_power=1e-20):
1031 """Extract fundamental frequencies and their power in dB from lists
1032 of harmonic groups.
1034 The inner list of 2-D arrays of the input argument is transformed
1035 into a 2-D array containig for each fish (1st dimension) the
1036 fundamental frequencies and powers (summed over all harmonics)
1037 extracted from the 2-D arrays.
1039 Parameters
1040 ----------
1041 group_list: (list of (list of ...)) list of 2-D arrays
1042 Arbitrarily nested lists of harmonic groups as returned by
1043 extract_fundamentals() and harmonic_groups() with the element
1044 [0, 0] being the fundamental frequency and the elements [:,1] being
1045 the powers of each harmonics.
1046 power: boolean
1047 If `False` convert the power into decibel using the
1048 `thunderlab.powerspectrum.decibel()` function.
1049 ref_power: float
1050 Reference power for computing decibel.
1051 If set to `None` the maximum power is used.
1052 min_power: float
1053 Power values smaller than `min_power` are set to `-np.inf`.
1055 Returns
1056 -------
1057 fundamentals: (list of (list of ...)) 2-D array
1058 Nested list (corresponding to `group_list`) of 2-D arrays
1059 with fundamental frequencies in first column and
1060 corresponding power in second column.
1061 """
1062 if len(group_list) == 0:
1063 return np.zeros((0, 2))
1065 # check whether group_list is list of harmonic groups:
1066 list_of_groups = True
1067 for group in group_list:
1068 if not ( hasattr(group, 'shape') and len(group.shape) == 2 ):
1069 list_of_groups = False
1070 break
1072 if list_of_groups:
1073 fundamentals = np.array([[group[0, 0], np.sum(group[:, 1])]
1074 for group in group_list if len(group) > 0])
1075 if not power:
1076 fundamentals[:, 1] = decibel(fundamentals[:, 1],
1077 ref_power, min_power)
1078 else:
1079 fundamentals = []
1080 for groups in group_list:
1081 f = fundamental_freqs_and_power(groups, power,
1082 ref_power, min_power)
1083 fundamentals.append(f)
1084 return fundamentals
1087def add_relative_power(freqs):
1088 """Add a column with relative power.
1090 For each element in `freqs`, its maximum power is subtracted
1091 from all powers.
1093 Parameters
1094 ----------
1095 freqs: list of 2-D arrays
1096 First column in the ndarrays is fundamental frequency and
1097 second column the corresponding power.
1098 Further columns are optional and kept in the returned list.
1099 fundamental_freqs_and_power() returns such a list.
1101 Returns
1102 -------
1103 power_freqs: list of 2-D arrays
1104 Same as freqs, but with an added column containing the relative power.
1105 """
1106 return [np.column_stack((f, f[:,1] - np.max(f[:,1]))) for f in freqs]
1109def add_power_ranks(freqs):
1110 """Add a column with power ranks.
1112 Parameters
1113 ----------
1114 freqs: list of 2-D arrays
1115 First column in the arrays is fundamental frequency and
1116 second column the corresponding power.
1117 Further columns are optional and kept in the returned list.
1118 fundamental_freqs_and_power() returns such a list.
1120 Returns
1121 -------
1122 rank_freqs: list of 2-D arrays
1123 Same as freqs, but with an added column containing the ranks.
1124 The highest power is assinged to zero,
1125 lower powers are assigned negative integers.
1126 """
1127 rank_freqs = []
1128 for f in freqs:
1129 i = np.argsort(f[:,1])[::-1]
1130 ranks = np.empty_like(i)
1131 ranks[i] = -np.arange(len(i))
1132 rank_freqs.append(np.column_stack((f, ranks)))
1133 return rank_freqs
1136def similar_indices(freqs, df_thresh, nextfs=0):
1137 """Indices of similar frequencies.
1139 If two frequencies from different elements in the inner lists of `freqs` are
1140 reciprocally the closest to each other and closer than `df_thresh`,
1141 then two indices (element, frequency) of the respective other frequency
1142 are appended.
1144 Parameters
1145 ----------
1146 freqs: (list of (list of ...)) list of 2-D arrays
1147 First column in the arrays is fundamental frequency.
1148 df_thresh: float
1149 Fundamental frequencies closer than this threshold are considered
1150 equal.
1151 nextfs: int
1152 If zero, compare all elements in freqs with each other. Otherwise,
1153 only compare with the `nextfs` next elements in freqs.
1155 Returns
1156 -------
1157 indices: (list of (list of ...)) list of list of two-tuples of int
1158 For each frequency of each element in `freqs` a list of two tuples containing
1159 the indices of elements and frequencies that are similar.
1160 """
1161 if len(freqs) == 0:
1162 return []
1164 # check whether freqs is list of fundamental frequencies and powers:
1165 list_of_freq_power = True
1166 for group in freqs:
1167 if not (hasattr(group, 'shape') and len(group.shape) == 2):
1168 list_of_freq_power = False
1169 break
1171 if list_of_freq_power:
1172 indices = [ [[] for j in range(len(freqs[i]))] for i in range(len(freqs))]
1173 for j in range(len(freqs)-1):
1174 freqsj = np.asarray(freqs[j])
1175 for m in range(len(freqsj)):
1176 freq1 = freqsj[m]
1177 nn = len(freqs) if nextfs == 0 else j+1+nextfs
1178 if nn > len(freqs):
1179 nn = len(freqs)
1180 for k in range(j+1, nn):
1181 freqsk = np.asarray(freqs[k])
1182 if len(freqsk) == 0:
1183 continue
1184 n = np.argmin(np.abs(freqsk[:,0] - freq1[0]))
1185 freq2 = freqsk[n]
1186 if np.argmin(np.abs(freqsj[:,0] - freq2[0])) != m:
1187 continue
1188 if np.abs(freq1[0] - freq2[0]) < df_thresh:
1189 indices[k][n].append((j, m))
1190 indices[j][m].append((k, n))
1191 else:
1192 indices = []
1193 for groups in freqs:
1194 indices.append(similar_indices(groups, df_thresh, nextfs))
1195 return indices
1198def unique_mask(freqs, df_thresh, nextfs=0):
1199 """Mark similar frequencies from different recordings as dublicate.
1201 If two frequencies from different elements in `freqs` are
1202 reciprocally the closest to each other and closer than `df_thresh`,
1203 then the one with the smaller power is marked for removal.
1205 Parameters
1206 ----------
1207 freqs: list of 2-D arrays
1208 First column in the arrays is fundamental frequency and
1209 second column the corresponding power or equivalent.
1210 If values in the second column are equal (e.g. they are the same ranks),
1211 and there is a third column (e.g. power),
1212 the third column is used to decide, which element should be removed.
1213 df_thresh: float
1214 Fundamental frequencies closer than this threshold are considered
1215 equal.
1216 nextfs: int
1217 If zero, compare all elements in freqs with each other. Otherwise,
1218 only compare with the `nextfs` next elements in freqs.
1220 Returns
1221 -------
1222 mask: list of boolean arrays
1223 For each element in `freqs` True if that frequency should be kept.
1224 """
1225 mask = [np.ones(len(freqs[i]), dtype=bool) for i in range(len(freqs))]
1226 for j in range(len(freqs)-1):
1227 freqsj = np.asarray(freqs[j])
1228 for m in range(len(freqsj)):
1229 freq1 = freqsj[m]
1230 nn = len(freqs) if nextfs == 0 else j+1+nextfs
1231 if nn > len(freqs):
1232 nn = len(freqs)
1233 for k in range(j+1, nn):
1234 freqsk = np.asarray(freqs[k])
1235 if len(freqsk) == 0:
1236 continue
1237 n = np.argmin(np.abs(freqsk[:,0] - freq1[0]))
1238 freq2 = freqsk[n]
1239 if np.argmin(np.abs(freqsj[:,0] - freq2[0])) != m:
1240 continue
1241 if np.abs(freq1[0] - freq2[0]) < df_thresh:
1242 if freq1[1] > freq2[1]:
1243 mask[k][n] = False
1244 elif freq1[1] < freq2[1]:
1245 mask[j][m] = False
1246 elif len(freq1) > 2:
1247 if freq1[2] > freq2[2]:
1248 mask[k][n] = False
1249 else:
1250 mask[j][m] = False
1251 else:
1252 mask[j][m] = False
1253 return mask
1256def unique(freqs, df_thresh, mode='power', nextfs=0):
1257 """Remove similar frequencies from different recordings.
1259 If two frequencies from different elements in the inner lists of `freqs`
1260 are reciprocally the closest to each other and closer than `df_thresh`,
1261 then the one with the smaller power is removed. As power, either the
1262 absolute power as provided in the second column of the data elements
1263 in `freqs` is taken (mode=='power'), or the relative power
1264 (mode='relpower'), or the power rank (mode='rank').
1266 Parameters
1267 ----------
1268 freqs: (list of (list of ...)) list of 2-D arrays
1269 First column in the arrays is fundamental frequency and
1270 second column the corresponding power, as returned by
1271 fundamental_freqs_and_power().
1272 df_thresh: float
1273 Fundamental frequencies closer than this threshold are considered
1274 equal.
1275 mode: string
1276 - 'power': use second column of freqs elements as power.
1277 - 'relpower': use relative power computed from the second column
1278 of freqs elements for deciding which frequency to delete.
1279 - 'rank': use rank of second column of freqs elements
1280 for deciding which frequency to delete.
1281 nextfs: int
1282 If zero, compare all elements in freqs with each other. Otherwise,
1283 only compare with the `nextfs` next elements in freqs.
1285 Returns
1286 -------
1287 uniqe_freqs: (list of (list of ...)) list of 2-D arrays
1288 Same as `freqs` but elements with similar fundamental frequencies
1289 removed.
1290 """
1291 if len(freqs) == 0:
1292 return []
1294 # check whether freqs is list of fundamental frequencies and powers:
1295 list_of_freq_power = True
1296 for group in freqs:
1297 if not (hasattr(group, 'shape') and len(group.shape) == 2):
1298 list_of_freq_power = False
1299 break
1301 if list_of_freq_power:
1302 if mode == 'power':
1303 mask = unique_mask(freqs, df_thresh, nextfs)
1304 elif mode == 'relpower':
1305 power_freqs = [f[:,[0, 2, 1]] for f in add_relative_power(freqs)]
1306 mask = unique_mask(power_freqs, df_thresh, nextfs)
1307 elif mode == 'rank':
1308 rank_freqs = [f[:,[0, 2, 1]] for f in add_power_ranks(freqs)]
1309 mask = unique_mask(rank_freqs, df_thresh, nextfs)
1310 else:
1311 raise ValueError('%s is not a valid mode for unique(). Choose one of "power" or "rank"')
1312 unique_freqs = []
1313 for f, m in zip(freqs, mask):
1314 unique_freqs.append(f[m])
1315 else:
1316 unique_freqs = []
1317 for groups in freqs:
1318 unique_freqs.append(unique(groups, df_thresh, mode, nextfs))
1319 return unique_freqs
1322def colors_markers():
1323 """Generate a list of colors and markers for plotting.
1325 Returns
1326 -------
1327 colors: list
1328 list of colors
1329 markers: list
1330 list of markers
1331 """
1332 # color and marker range:
1333 colors = []
1334 markers = []
1335 mr2 = []
1337 # first color range:
1338 cc0 = cm.gist_rainbow(np.linspace(0.0, 1.0, 8))
1340 # shuffle it:
1341 for k in range((len(cc0) + 1) // 2):
1342 colors.extend(cc0[k::(len(cc0) + 1) // 2])
1343 markers.extend(len(cc0) * 'o')
1344 mr2.extend(len(cc0) * 'v')
1345 # second darker color range:
1346 cc1 = cm.gist_rainbow(np.linspace(0.33 / 7.0, 1.0, 7))
1347 cc1 = mc.hsv_to_rgb(mc.rgb_to_hsv(np.array([cc1[:, :3]])) * np.array([1.0, 0.9, 0.7]))[0]
1348 cc1 = np.hstack((cc1, np.ones((len(cc1),1))))
1349 # shuffle it:
1350 for k in range((len(cc1) + 1) // 2):
1351 colors.extend(cc1[k::(len(cc1) + 1) // 2])
1352 markers.extend(len(cc1) * '^')
1353 mr2.extend(len(cc1) * '*')
1354 # third lighter color range:
1355 cc2 = cm.gist_rainbow(np.linspace(0.67 / 6.0, 1.0, 6))
1356 cc2 = mc.hsv_to_rgb(mc.rgb_to_hsv(np.array([cc1[:, :3]])) * np.array([1.0, 0.5, 1.0]))[0]
1357 cc2 = np.hstack((cc2, np.ones((len(cc2),1))))
1358 # shuffle it:
1359 for k in range((len(cc2) + 1) // 2):
1360 colors.extend(cc2[k::(len(cc2) + 1) // 2])
1361 markers.extend(len(cc2) * 'D')
1362 mr2.extend(len(cc2) * 'x')
1363 markers.extend(mr2)
1364 return colors, markers
1367def plot_harmonic_groups(ax, group_list, indices=None, max_groups=0,
1368 skip_bad=False, sort_by_freq=True, label_power=False,
1369 colors=None, markers=None, legend_rows=8, **kwargs):
1370 """Mark decibel power of fundamentals and their harmonics in a plot.
1372 Parameters
1373 ----------
1374 ax: axis for plot
1375 Axis used for plotting.
1376 group_list: list of 2-D arrays
1377 Lists of harmonic groups as returned by extract_fundamentals() and
1378 harmonic_groups() with the element [0, 0] of the harmonic groups being
1379 the fundamental frequency, and element[0, 1] being the corresponding power.
1380 indices: list of int or None
1381 If smaller than zero then set the legend label of the corresponding group in brackets.
1382 max_groups: int
1383 If not zero plot only the max_groups most powerful groups.
1384 skip_bad: bool
1385 Skip harmonic groups without index (entry in indices is negative).
1386 sort_by_freq: boolean
1387 If True sort legend by frequency, otherwise by power.
1388 label_power: boolean
1389 If `True` put the power in decibel in addition to the frequency into the legend.
1390 colors: list of colors or None
1391 If not None list of colors for plotting each group
1392 markers: list of markers or None
1393 If not None list of markers for plotting each group
1394 legend_rows: int
1395 Maximum number of rows to be used for the legend.
1396 kwargs:
1397 Key word arguments for the legend of the plot.
1398 """
1399 if len(group_list) == 0:
1400 return
1402 # sort by power:
1403 powers = np.array([np.sum(group[:,1]) for group in group_list])
1404 max_power = np.max(powers)
1405 idx = np.argsort(-powers)
1406 if max_groups > 0 and len(idx > max_groups):
1407 idx = idx[:max_groups]
1409 # sort by frequency:
1410 if sort_by_freq:
1411 freqs = np.array([group_list[group][0,0] for group in idx])
1412 if legend_rows > 0 and legend_rows < len(freqs):
1413 idx = idx[np.argsort(freqs[:legend_rows])]
1414 else:
1415 idx = idx[np.argsort(freqs)]
1417 # plot:
1418 k = 0
1419 for i in idx:
1420 if indices is not None and skip_bad and indices[i] < 0:
1421 continue
1422 group = group_list[i]
1423 x = group[:,0]
1424 y = decibel(group[:,1])
1425 msize = 7.0 + 10.0*(powers[i]/max_power)**0.25
1426 color_kwargs = {}
1427 if colors is not None:
1428 color_kwargs = {'color': colors[k%len(colors)]}
1429 label = f'{group[0, 0]:6.1f} Hz'
1430 if label_power:
1431 label += f' {decibel(np.array([np.sum(group[:,1])]))[0]:6.1f} dB'
1432 if indices is not None:
1433 if indices[i] < 0:
1434 label = '(' + label + ')'
1435 else:
1436 label = ' ' + label + ' '
1437 if legend_rows > 5 and k >= legend_rows:
1438 label = None
1439 if markers is None:
1440 ax.plot(x, y, 'o', ms=msize, label=label, **color_kwargs)
1441 else:
1442 if k >= len(markers):
1443 break
1444 ax.plot(x, y, linestyle='None', marker=markers[k],
1445 mec=None, mew=0.0, ms=msize, label=label, **color_kwargs)
1446 k += 1
1448 # legend:
1449 if legend_rows > 0:
1450 if legend_rows > 5:
1451 ncol = 1
1452 else:
1453 ncol = (len(idx)-1) // legend_rows + 1
1454 ax.legend(numpoints=1, ncol=ncol, **kwargs)
1455 else:
1456 ax.legend(numpoints=1, **kwargs)
1459def plot_psd_harmonic_groups(ax, psd_freqs, psd, group_list,
1460 mains=None, all_freqs=None,
1461 good_freqs=None, log_freq=False,
1462 min_freq=0.0, max_freq=2000.0, ymarg=0.0,
1463 sstyle=dict(color='tab:blue', lw=1)):
1464 """Plot decibel power-spectrum with detected peaks, harmonic groups,
1465 and mains frequencies.
1467 Parameters
1468 ----------
1469 psd_freqs: array
1470 Frequencies of the power spectrum.
1471 psd: array
1472 Power spectrum (linear, not decible).
1473 group_list: list of 2-D arrays
1474 Lists of harmonic groups as returned by extract_fundamentals() and
1475 harmonic_groups() with the element [0, 0] of the harmonic groups being the fundamental frequency,
1476 and element[0, 1] being the corresponding power.
1477 mains: 2-D array
1478 Frequencies and power of multiples of the mains frequency found in the power spectrum.
1479 all_freqs: 2-D array
1480 Peaks in the power spectrum detected with low threshold.
1481 good_freqs: 1-D array
1482 Frequencies of peaks detected with high threshold.
1483 log_freq: boolean
1484 Logarithmic (True) or linear (False) frequency axis.
1485 min_freq: float
1486 Limits of frequency axis are set to `(min_freq, max_freq)`
1487 if `max_freq` is greater than zero
1488 max_freq: float
1489 Limits of frequency axis are set to `(min_freq, max_freq)`
1490 and limits of power axis are computed from powers below max_freq
1491 if `max_freq` is greater than zero
1492 ymarg: float
1493 Add this to the maximum decibel power for setting the ylim.
1494 sstyle: dict
1495 Arguments passed on to the plot command for the power spectrum.
1496 """
1498 # mark all and good psd peaks:
1499 pmin, pmax = ax.get_ylim()
1500 doty = pmax - 5.0
1501 if all_freqs is not None:
1502 ax.plot(all_freqs[:, 0], np.zeros(len(all_freqs[:, 0])) + doty, 'o', color='#ffffff')
1503 if good_freqs is not None:
1504 ax.plot(good_freqs, np.zeros(len(good_freqs)) + doty, 'o', color='#888888')
1505 # mark mains frequencies:
1506 if mains is not None and len(mains) > 0:
1507 fpeaks = mains[:, 0]
1508 fpeakinx = [int(np.round(fp/(psd_freqs[1]-psd_freqs[0]))) for fp in fpeaks if fp < psd_freqs[-1]]
1509 ax.plot(fpeaks[:len(fpeakinx)], decibel(psd[fpeakinx]), linestyle='None',
1510 marker='.', color='k', ms=10, mec=None, mew=0.0,
1511 label=f'{mains[0, 0]:3.0f} Hz mains')
1512 # mark harmonic groups:
1513 colors, markers = colors_markers()
1514 plot_harmonic_groups(ax, group_list, max_groups=0, sort_by_freq=True,
1515 colors=colors, markers=markers, legend_rows=8,
1516 loc='upper right')
1517 # plot power spectrum:
1518 plot_decibel_psd(ax, psd_freqs, psd, log_freq=log_freq, min_freq=min_freq,
1519 max_freq=max_freq, ymarg=ymarg, sstyle=sstyle)
1522def add_psd_peak_detection_config(cfg, low_threshold=0.0, high_threshold=0.0,
1523 thresh_bins=100,
1524 low_thresh_factor=6.0, high_thresh_factor=10.0):
1525 """Add parameter needed for detection of peaks in power spectrum used
1526 by harmonic_groups() as a new section to a configuration.
1528 Parameters
1529 ----------
1530 cfg: ConfigFile
1531 The configuration.
1532 """
1533 cfg.add_section('Thresholds for peak detection in power spectra:')
1534 cfg.add('lowThreshold', low_threshold, 'dB', 'Threshold for all peaks.\n If 0.0 estimate threshold from histogram.')
1535 cfg.add('highThreshold', high_threshold, 'dB', 'Threshold for good peaks. If 0.0 estimate threshold from histogram.')
1537 cfg.add_section('Threshold estimation:\nIf no thresholds are specified they are estimated from the histogram of the decibel power spectrum.')
1538 cfg.add('thresholdBins', thresh_bins, '', 'Number of bins used to compute the histogram used for threshold estimation.')
1539 cfg.add('lowThresholdFactor', low_thresh_factor, '', 'Factor for multiplying standard deviation of noise floor for lower threshold.')
1540 cfg.add('highThresholdFactor', high_thresh_factor, '', 'Factor for multiplying standard deviation of noise floor for higher threshold.')
1543def psd_peak_detection_args(cfg):
1544 """Translates a configuration to the respective parameter names for
1545 the detection of peaks in power spectrum used by
1546 harmonic_groups(). The return value can then be passed as
1547 key-word arguments to this function.
1549 Parameters
1550 ----------
1551 cfg: ConfigFile
1552 The configuration.
1554 Returns
1555 -------
1556 a: dict
1557 Dictionary with names of arguments of the `harmonic-group()` function
1558 and their values as supplied by `cfg`.
1559 """
1560 return cfg.map({'low_threshold': 'lowThreshold',
1561 'high_threshold': 'highThreshold',
1562 'thresh_bins': 'thresholdBins',
1563 'low_thresh_factor': 'lowThresholdFactor',
1564 'high_thresh_factor': 'highThresholdFactor'})
1567def add_harmonic_groups_config(cfg, mains_freq=60.0, mains_freq_tol=1.0,
1568 max_divisor=4, freq_tol_fac=1.0,
1569 max_freq_tol=1.0, min_group_size=3,
1570 min_freq=20.0, max_freq=2000.0, max_db_diff=20.0,
1571 max_harmonics_db=-5.0, max_harmonics=0, max_groups=0):
1572 """Add parameter needed for detection of harmonic groups as a new
1573 section to a configuration.
1575 Parameters
1576 ----------
1577 cfg: ConfigFile
1578 The configuration.
1579 """
1580 cfg.add_section('Harmonic groups:')
1581 cfg.add('mainsFreq', mains_freq, 'Hz', 'Mains frequency to be excluded.')
1582 cfg.add('mainsFreqTolerance', mains_freq_tol, 'Hz', 'Exclude peaks within this tolerance around multiples of the mains frequency.')
1583 cfg.add('minimumGroupSize', min_group_size, '',
1584'Minimum number of harmonics (inclusively fundamental) that make up a harmonic group.')
1585 cfg.add('maxDivisor', max_divisor, '', 'Maximum ratio between the frequency of the largest peak and its fundamental.')
1586 cfg.add('freqTolerance', freq_tol_fac, '',
1587 'Harmonics should be within this factor times the frequency resolution of the power spectrum. Needs to be higher than 0.5!')
1588 cfg.add('maximumFreqTolerance', max_freq_tol, 'Hz',
1589 'Maximum deviation of harmonics from their expected value.')
1591 cfg.add_section('Acceptance of best harmonic groups:')
1592 cfg.add('minimumFrequency', min_freq, 'Hz', 'Minimum frequency allowed for the fundamental.')
1593 cfg.add('maximumFrequency', max_freq, 'Hz', 'Maximum frequency allowed for the fundamental.')
1594 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.')
1595 cfg.add('maximumHarmonicsPower', max_harmonics_db, 'dB', 'Maximum allowed power of the minimumGroupSize-th and higher harmonics relative to peak power.')
1596 cfg.add('maximumHarmonics', max_harmonics, '', '0: keep all, >0 only keep the first # harmonics.')
1597 cfg.add('maximumGroups', max_groups, '', 'Maximum number of harmonic groups. If 0 process all.')
1600def harmonic_groups_args(cfg):
1601 """Translates a configuration to the respective parameter names of
1602 the harmonic-group detection functions. The return value can then
1603 be passed as key-word arguments to this function.
1605 Parameters
1606 ----------
1607 cfg: ConfigFile
1608 The configuration.
1610 Returns
1611 -------
1612 a: dict
1613 Dictionary with names of arguments of the harmonic-group detection functions
1614 and their values as supplied by `cfg`.
1615 """
1616 return cfg.map({'mains_freq': 'mainsFreq',
1617 'mains_freq_tol': 'mainsFreqTolerance',
1618 'freq_tol_fac': 'freqTolerance',
1619 'max_freq_tol': 'maximumFreqTolerance',
1620 'max_divisor': 'maxDivisor',
1621 'min_group_size': 'minimumGroupSize',
1622 'min_freq': 'minimumFrequency',
1623 'max_freq': 'maximumFrequency',
1624 'max_db_diff': 'maximumPowerDifference',
1625 'max_harmonics_db': 'maximumHarmonicsPower',
1626 'max_harmonics': 'maximumHarmonics',
1627 'max_groups': 'maximumGroups'})
1630def main(data_file=None):
1631 import matplotlib.pyplot as plt
1632 from thunderlab.powerspectrum import psd
1633 from .fakefish import wavefish_eods
1635 if data_file is None:
1636 # generate data:
1637 title = 'simulation'
1638 rate = 44100.0
1639 d = 20.0
1640 noise = 0.01
1641 eodfs = [123.0, 333.0, 666.0, 666.5]
1642 fish1 = 0.5*wavefish_eods('Eigenmannia', eodfs[0], rate,
1643 duration=d, noise_std=noise)
1644 fish2 = 1.0*wavefish_eods('Eigenmannia', eodfs[1], rate,
1645 duration=d, noise_std=noise)
1646 fish3 = 10.0*wavefish_eods('Alepto', eodfs[2], rate,
1647 duration=d, noise_std=noise)
1648 fish4 = 6.0*wavefish_eods('Alepto', eodfs[3], rate,
1649 duration=d, noise_std=noise)
1650 data = fish1 + fish2 + fish3 + fish4
1651 else:
1652 from thunderlab.dataloader import load_data
1653 print(f"load {data_file} ...")
1654 data, rate, unit, amax = load_data(data_file)
1655 data = data[:,0]
1656 title = data_file
1658 # retrieve fundamentals from power spectrum:
1659 psd_data = psd(data, rate, freq_resolution=0.1)
1660 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)
1661 fig, ax = plt.subplots()
1662 plot_psd_harmonic_groups(ax, psd_data[0], psd_data[1], groups, mains,
1663 all_freqs, good_freqs, max_freq=3000.0)
1664 ax.set_title(title)
1665 plt.show()
1667 # unify fundamental frequencies:
1668 fundamentals = fundamental_freqs(groups)
1669 np.set_printoptions(formatter={'float': lambda x: f'{x:5.1f}'})
1670 print('fundamental frequencies extracted from power spectrum:')
1671 print(fundamentals)
1672 print('')
1673 freqs = fundamental_freqs_and_power([groups])
1674 freqs.append(np.array([[44.0, -20.0], [44.2, -10.0], [320.5, 2.5], [665.5, 5.0], [666.2, 10.0]]))
1675 freqs.append(np.array([[123.3, 1.0], [320.2, -2.0], [668.4, 2.0]]))
1676 rank_freqs = add_relative_power(freqs)
1677 rank_freqs = add_power_ranks(rank_freqs)
1678 print('all frequencies (frequency, power, relpower, rank):')
1679 print('\n'.join(( str(f) for f in rank_freqs)))
1680 print('')
1681 indices = similar_indices(freqs, 1.0)
1682 print('similar indices:')
1683 print('\n'.join(( ('\n '.join((str(f) for f in g)) for g in indices))))
1684 print('')
1685 unique_freqs = unique(freqs, 1.0, 'power')
1686 print('unique power:')
1687 print('\n'.join(( str(f) for f in unique_freqs)))
1688 print('')
1689 unique_freqs = unique(freqs, 1.0, 'relpower')
1690 print('unique relative power:')
1691 print('\n'.join(( str(f) for f in unique_freqs)))
1692 print('')
1693 unique_freqs = unique(freqs, 1.0, 'rank')
1694 print('unique rank:')
1695 print('\n'.join(( str(f) for f in unique_freqs)))
1696 print('')
1697 unique_freqs = unique(freqs, 1.0, 'rank', 1)
1698 print('unique rank for next neighour only:')
1699 print('\n'.join(( str(f) for f in unique_freqs)))
1700 print('')
1703if __name__ == "__main__":
1704 import sys
1705 data_file = sys.argv[1] if len(sys.argv) > 1 else None
1706 main(data_file)