Peaks extraction methods
Contents
Peaks extraction methods#
This notebook illustrates the use of the biotuner toolbox with the aims to show that we can retrieve harmonic information from a generated signal using different peaks extraction methods
This notebook is one step towards the exploration of biological harmonies and their relevance for the carving of new musical systems, as compared to pseudo-random generative engines.
Retrieving frequency information from simulated signals using the biotuner#
First, we will generate two signals composed of two sets of predetermined frequencies. The first one will contain harmonic frequencies, while the second will contain non-harmonic frequencies
import numpy as np
import matplotlib.pyplot as plt
import warnings
from biotuner.biotuner_utils import generate_signal
from neurodsp.sim import sim_combined
warnings.filterwarnings('ignore', category=DeprecationWarning)
warnings.filterwarnings('ignore', category=RuntimeWarning)
warnings.filterwarnings('ignore', category=FutureWarning)
warnings.filterwarnings('ignore', category=UserWarning)
sf = 1000
length =4
theta = 0
freqs_harm = [2, 6, 12, 18, 24, 36, 48]
freqs_noharm = [2, 5, 11, 17, 23, 33, 49]
# set amplitudes to random between 0.5 and 1
amps = np.random.uniform(0.5, 1, len(freqs_harm))
signal_harm = generate_signal(sf, length, freqs_harm, amps, show=False)
signal_noharm = generate_signal(sf, length, freqs_noharm, amps, show=False)
# Set up simulation parameters for neurodsp
n_seconds = 4
fs = 1000
components_harm = {
'sim_powerlaw': {'exponent': -0.8},
'sim_oscillation': [{'freq': freq} for freq in freqs_harm] # list of dictionaries
}
components_noharm = {
'sim_powerlaw': {'exponent': -0.8},
'sim_oscillation': [{'freq': freq} for freq in freqs_noharm] # list of dictionaries
}
# Each component (aperiodic & oscillatory) needs a variance
component_variances = [10 for _ in range(len(freqs_noharm) + 1)] # +1 for the 'sim_powerlaw' component
# Simulate the signal
brain_signal_harm = sim_combined(n_seconds, fs, components_harm, component_variances)
brain_signal_noharm = sim_combined(n_seconds, fs, components_noharm, component_variances)
# plot the signals with 2 rows of subplots
fig, axs = plt.subplots(2, 2, figsize=(10, 4)) # Create 2x2 subplots
# First row: Simulated oscillations
axs[0, 0].plot(signal_harm, color='darkblue')
axs[0, 0].set_title('Harmonic signal')
axs[0, 0].set_xlabel('Time')
axs[0, 0].set_ylabel('Amplitude')
axs[0, 0].set_xlim(0, 1000)
axs[0, 1].plot(signal_noharm, color='darkred')
axs[0, 1].set_title('Non-harmonic signal')
axs[0, 1].set_xlabel('Time')
axs[0, 1].set_ylabel('Amplitude')
axs[0, 1].set_xlim(0, 1000)
# Second row: Simulated brain signals
axs[1, 0].plot(brain_signal_harm, color='darkblue')
axs[1, 0].set_title('Simulated Harmonic Brain Signal')
axs[1, 0].set_xlabel('Time')
axs[1, 0].set_ylabel('Amplitude')
axs[1, 0].set_xlim(0, 1000)
axs[1, 1].plot(brain_signal_noharm, color='darkred')
axs[1, 1].set_title('Simulated Non-Harmonic Brain Signal')
axs[1, 1].set_xlabel('Time')
axs[1, 1].set_ylabel('Amplitude')
axs[1, 1].set_xlim(0, 1000)
plt.tight_layout()
plt.show()
Then, we want to retrieve the generator frequencies using different peaks_functions offered by the biotuner.
Peaks functions (see the doc)#
fixed : ranges of frequency bands are fixed
FOOOF : peaks rising above the aperiodic component
EMD: Intrinsic Mode Functions (IMFs) are derived from Empirical Mode Decomposition (EMD)
FFT is computed on each IMF
harmonic_recurrence: keeps peaks for which a maximum of other peaks are harmonics
EIMC: keeps peaks for which a maximum of other peaks are intermodulation components
Fixed#
And now with the ‘fixed’ bands method. With the specified bands, we would expect to find 2, 6 and 24Hz for the harmonic signal and 2, 5 and 23Hz for the non-harmonic signal
from biotuner.biotuner_object import compute_biotuner
FREQ_BANDS = [[1, 3], [4, 7], [10, 20], [20, 30]]
# Initialize biotuner object for harmonic signal
bt_harm = compute_biotuner(sf=sf, peaks_function = 'fixed', precision = 0.5)
# Extract spectral peaks
bt_harm.peaks_extraction(signal_harm, FREQ_BANDS = FREQ_BANDS, min_freq = 1, max_freq = 50, min_harms = 2, harm_limit = 1000, graph=True)
# Initialize biotuner object for non-harmonic signal
bt_noharm = compute_biotuner(sf=sf, peaks_function = 'fixed', precision = 0.5)
# Extract spectral peaks
bt_noharm.peaks_extraction(signal_noharm, FREQ_BANDS = FREQ_BANDS, min_freq = 1, max_freq = 50, min_harms = 2, harm_limit = 1000, graph=True)
# Compare results
print('HARMONIC PEAKS', '\nRETRIEVED', bt_harm.peaks, 'ORIGINAL', freqs_harm)
print('NON-HARMONIC PEAKS', '\nRETRIEVED', bt_noharm.peaks, 'ORIGINAL', freqs_noharm)
HARMONIC PEAKS
RETRIEVED [ 2. 6. 12. 24.] ORIGINAL [2, 6, 12, 18, 24, 36, 48]
NON-HARMONIC PEAKS
RETRIEVED [ 2. 5. 11. 23.] ORIGINAL [2, 5, 11, 17, 23, 33, 49]
FOOOF#
# Initialize biotuner object for harmonic signal
bt_harm = compute_biotuner(sf=sf, peaks_function = 'FOOOF', precision = 0.5)
# Extract spectral peaks
bt_harm.peaks_extraction(signal_harm, min_freq = 1, max_freq = 50, min_harms = 2, harm_limit = 1000, graph=True)
# Initialize biotuner object for non-harmonic signal
bt_noharm = compute_biotuner(sf=sf, peaks_function = 'FOOOF', precision = 0.5)
# Extract spectral peaks
bt_noharm.peaks_extraction(signal_noharm, min_freq = 1, max_freq = 50, min_harms = 2, harm_limit = 1000, graph=True)
# Compare results
print('HARMONIC PEAKS', '\nRETRIEVED', bt_harm.peaks, 'ORIGINAL', freqs_harm)
print('NON-HARMONIC PEAKS', '\nRETRIEVED', bt_noharm.peaks, 'ORIGINAL', freqs_noharm)
==================================================================================================
FOOOF - POWER SPECTRUM MODEL
The model was run on the frequency range 0 - 50 Hz
Frequency Resolution is 0.50 Hz
Aperiodic Parameters (offset, exponent):
-31.7637, -1.0123
7 peaks were found:
CF: 2.00, PW: 36.313, BW: 1.06
CF: 6.00, PW: 35.957, BW: 1.05
CF: 12.00, PW: 36.079, BW: 1.05
CF: 18.00, PW: 35.839, BW: 1.06
CF: 24.00, PW: 35.479, BW: 1.04
CF: 36.00, PW: 35.369, BW: 1.06
CF: 48.00, PW: 35.189, BW: 1.03
Goodness of fit metrics:
R^2 of model fit is 0.9062
Error of the fit is 2.4452
==================================================================================================
==================================================================================================
FOOOF - POWER SPECTRUM MODEL
The model was run on the frequency range 0 - 50 Hz
Frequency Resolution is 0.50 Hz
Aperiodic Parameters (offset, exponent):
-31.5186, -0.5956
7 peaks were found:
CF: 2.00, PW: 36.117, BW: 1.05
CF: 5.00, PW: 35.849, BW: 1.06
CF: 11.00, PW: 36.161, BW: 1.06
CF: 17.00, PW: 36.145, BW: 1.05
CF: 23.00, PW: 35.713, BW: 1.06
CF: 33.00, PW: 35.769, BW: 1.07
CF: 49.00, PW: 35.445, BW: 1.06
Goodness of fit metrics:
R^2 of model fit is 0.9128
Error of the fit is 2.3753
==================================================================================================
HARMONIC PEAKS
RETRIEVED [ 2. 12. 6. 18. 24.] ORIGINAL [2, 6, 12, 18, 24, 36, 48]
NON-HARMONIC PEAKS
RETRIEVED [11. 17. 2. 5. 33.] ORIGINAL [2, 5, 11, 17, 23, 33, 49]
Empirical Mode Decomposition#
# Initialize biotuner object for harmonic signal
bt_harm = compute_biotuner(sf=sf, peaks_function='EMD', precision=0.5)
# Extract spectral peaks
bt_harm.peaks_extraction(signal_harm, min_freq=1, max_freq=50, nIMFs=3, prominence=10, rel_height=1, keep_first_IMF=True,
)
# Initialize biotuner object for non-harmonic signal
bt_noharm = compute_biotuner(sf=sf, peaks_function='EMD', precision=0.5)
# Extract spectral peaks
bt_noharm.peaks_extraction(signal_noharm, min_freq=1, max_freq=50, nIMFs=3, prominence=10, rel_height=1, keep_first_IMF=True,
)
# Compare results
print('HARMONIC PEAKS', '\nRETRIEVED', bt_harm.peaks, 'ORIGINAL', freqs_harm)
print('NON-HARMONIC PEAKS', '\nRETRIEVED', bt_noharm.peaks, 'ORIGINAL', freqs_noharm)
HARMONIC PEAKS
RETRIEVED [ 2. 6. 18. 36.] ORIGINAL [2, 6, 12, 18, 24, 36, 48]
NON-HARMONIC PEAKS
RETRIEVED [ 5. 11. 23. 49.] ORIGINAL [2, 5, 11, 17, 23, 33, 49]
Harmonic recurrence#
This method identifies spectral peaks that have the highest recurrence in the spectrum based on their harmonic series. In the graph, each colored line reprensent a spectral peak, while the associated dashed lines represent harmonics of this peak.
# Initialize biotuner object for harmonic signal
bt_harm = compute_biotuner(sf=sf, peaks_function = 'harmonic_recurrence', precision = 0.5)
#set figzise
plt.figure(figsize=(5, 2))
# Extract spectral peaks
bt_harm.peaks_extraction(signal_harm, min_freq=1, max_freq=50, min_harms=2, harm_limit=1000, graph=True,
prominence=15, rel_height=1, n_peaks=4)
# Initialize biotuner object for non-harmonic signal
bt_noharm = compute_biotuner(sf=sf, peaks_function='harmonic_recurrence', precision=0.5)
# Extract spectral peaks
bt_noharm.peaks_extraction(signal_noharm, min_freq=1, max_freq=50, min_harms=2, harm_limit=1000, graph=True,
prominence=15, rel_height=1, n_peaks=4)
# Compare results
print('HARMONIC PEAKS', '\nRETRIEVED', bt_harm.peaks, 'ORIGINAL', freqs_harm)
print('NON-HARMONIC PEAKS', '\nRETRIEVED', bt_noharm.peaks, 'ORIGINAL', freqs_noharm)
<Figure size 500x200 with 0 Axes>
HARMONIC PEAKS
RETRIEVED [ 2. 6. 12. 48.] ORIGINAL [2, 6, 12, 18, 24, 36, 48]
NON-HARMONIC PEAKS
RETRIEVED [ 5. 2. 33. 11.] ORIGINAL [2, 5, 11, 17, 23, 33, 49]
Endogeneous Inter-Modulation Components#
# Initialize biotuner object for harmonic signal
bt_harm = compute_biotuner(sf=sf, peaks_function='EIMC', precision=0.5)
# Extract spectral peaks
bt_harm.peaks_extraction(signal_harm, min_freq=1, max_freq=50, graph=True,
n_peaks=4, prominence=10, rel_height=10)
# Initialize biotuner object for non-harmonic signal
bt_noharm = compute_biotuner(sf=sf, peaks_function='EIMC', precision=0.5)
# Extract spectral peaks
bt_noharm.peaks_extraction(signal_noharm, min_freq=1, max_freq=50, graph=True,
n_peaks=4, prominence=10, rel_height=10)
# Round the peaks according to the precision, which is a float
bt_harm.peaks = [round(x * 2) / 2 for x in bt_harm.peaks]
bt_noharm.peaks = [round(x * 2) / 2 for x in bt_noharm.peaks]
# Compare results
print('HARMONIC PEAKS', '\nRETRIEVED', bt_harm.peaks, 'ORIGINAL', freqs_harm)
print('NON-HARMONIC PEAKS', '\nRETRIEVED', bt_noharm.peaks, 'ORIGINAL', freqs_noharm)
HARMONIC PEAKS
RETRIEVED [12.0, 18.0, 6.0, 48.0] ORIGINAL [2, 6, 12, 18, 24, 36, 48]
NON-HARMONIC PEAKS
RETRIEVED [11.0, 33.0, 49.0, 40.5] ORIGINAL [2, 5, 11, 17, 23, 33, 49]
Effect of noise on peaks detection#
We will add noise to the signal to see how the algorithms perform
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import welch
# Set a target channel noise power to something very noisy
noise_db_list = range(0, 25, 5)
# Create subplots
fig, axs = plt.subplots(1, 4, figsize=(16, 4))
colors = plt.cm.inferno(np.linspace(0, 0.6, len(noise_db_list)))
for i, db in enumerate(noise_db_list):
# Convert to linear Watt units
target_noise_watts = 10 ** (db / 10)
# Generate noise samples
mean_noise = 0
noise_signal_harm = np.random.normal(mean_noise, np.sqrt(target_noise_watts), len(signal_harm))
noise_signal_brain_harm = np.random.normal(mean_noise, np.sqrt(target_noise_watts), len(brain_signal_harm))
# Noise up the original signal (oscillations and brain signals) and plot
noised_signal_harm = signal_harm + noise_signal_harm
noised_brain_signal_harm = brain_signal_harm + noise_signal_brain_harm
# For harmonic oscillation
axs[0].plot(noised_signal_harm[0:100], color = colors[i], label='{}dB noise'.format(db), alpha=0.8)
freqs, psd = welch(noised_signal_harm, sf, nperseg=1024)
axs[1].semilogy(freqs, psd, color = colors[i])
# For simulated brain signal
axs[2].plot(noised_brain_signal_harm[0:100], color = colors[i], alpha=1)
freqs, psd = welch(noised_brain_signal_harm, sf, nperseg=1024)
axs[3].semilogy(freqs, psd, color = colors[i])
# Set labels and titles
axs[0].set_title('Harmonic signal')
axs[1].set_title('Spectrum')
axs[2].set_title('Simulated brain signal')
axs[3].set_title('Spectrum')
for ax in axs[1:4]:
ax.set_xlim([2, 50])
ax.set_xlabel('Frequency (Hz)')
ax.set_ylabel('PSD [V**2/Hz]')
# Add legends
axs[0].legend(loc='upper right')
axs[1].legend(loc='upper right')
# Display the plot
plt.tight_layout()
plt.show()
No artists with labels found to put in legend. Note that artists whose label start with an underscore are ignored when legend() is called with no argument.
We will now quantify the effect of noise on peak detection by plotting the performance of the algorithms for different intensity of noise.#
The first plot is for a signal made from a set of harmonic frequencies.
import seaborn as sns
import pandas as pd
import matplotlib.pyplot as plt
from neurodsp.sim import sim_combined
from scipy.interpolate import make_interp_spline, BSpline
# Setting up parameters
sf = 1000
length = 4
n_times = 30 #number of times the bootstraping is executed
freqs = [2, 6, 12, 18, 24, 36]
amps = [1, 1, 1, 1, 1, 1]
FREQ_BANDS = [[1, 5], [5, 10], [10, 20], [20, 30]]
# Simulate the signals
components = {
'sim_powerlaw': {'exponent': -1},
'sim_oscillation': [{'freq': freq} for freq in freqs] # list of dictionaries
}
component_variances = [1 for _ in range(len(freqs) + 1)] # +1 for the 'sim_powerlaw' component
# Simulate original signal
simulated_EEG = sim_combined(length, sf, components, component_variances)
# Simulate EEG signal
simulated_oscillations = generate_signal(sf, length, freqs, amps, show=False, color = 'turquoise') # Assuming generate_signal is your custom function
signals = {'Simulated Oscillations': simulated_oscillations, 'Simulated EEG Signal': simulated_EEG}
for title, signal in signals.items():
peaks_methods = ['fixed', 'harmonic_recurrence', 'EIMC', 'EMD', 'FOOOF']
noise_db_list = range(0, 40, 1)
colors = ['darkred', 'darkblue', 'darkorange', 'black', 'cyan', 'deeppink']
data = []
for method, color in zip(peaks_methods, colors):
print(method)
median_distance_tot = []
for db in noise_db_list:
t = 0
median_distance = []
while t < n_times:
# Convert to linear Watt units
target_noise_watts = 10 ** (db / 10)
# Generate noise samples
mean_noise = 0
noise_signal = np.random.normal(mean_noise, np.sqrt(target_noise_watts), len(signal))
# Noise up the original signal
noised_signal = signal + noise_signal
# set dictionary of prominence and rel_height for each method
prominences = {'fixed': 10, 'harmonic_recurrence': 8, 'EMD': 1, 'FOOOF': 10, 'HH1D_max': 1, 'EIMC': 7}
rel_heights = {'fixed': 1, 'harmonic_recurrence': 1, 'EMD': 0.7, 'FOOOF': 1, 'HH1D_max': 1, 'EIMC': 1}
# Initialize biotuner object
biotuning = compute_biotuner(sf = sf, peaks_function=method, precision=0.5)
if method == 'EMD':
biotuning.peaks_extraction(noised_signal, min_freq=1, max_freq=50, min_harms=2, harm_limit=128,
FREQ_BANDS = FREQ_BANDS, verbose=False, graph=False, prominence=prominences[method],
rel_height=rel_heights[method], nIMFs=6, keep_first_IMF=False, n_peaks=6)
biotuning.peaks = biotuning.peaks[0:4] # select only the first 4 peaks
biotuning.peak = [x for x in biotuning.peaks if x <= freqs[-1]+1]
else:
try:
biotuning.peaks_extraction(noised_signal, min_freq=1, max_freq=50, min_harms=2, harm_limit=128,
FREQ_BANDS = FREQ_BANDS, verbose=False, graph=False, prominence=prominences[method],
rel_height=rel_heights[method], nIMFs=3, keep_first_IMF=True, n_peaks=4)
biotuning.peaks = [x for x in biotuning.peaks if x <= freqs[-1]+1]
except:
biotuning.peaks = []
if db == 0:
no_noise_peaks = biotuning.peaks
if len(biotuning.peaks) == 0:
median_distance.append(np.nan)
else:
list_distance = []
for p in biotuning.peaks:
takeClosest = lambda num,collection:min(collection,key=lambda x:abs(x-num))
closest = takeClosest(p,freqs)
list_distance.append(abs(p-closest))
median_distance.append(np.nanmedian(list_distance))
t += 1
median_distance_tot.append(np.nanmedian(median_distance))
# calculate standard deviation
std_dev = np.std(median_distance)
# Continue with your computations...
data.append([method, db, np.nanmedian(median_distance), np.std(median_distance)])
df = pd.DataFrame(data, columns=["Method", "Noise decibels", "Median distance", "Std deviation"])
# Add the option for std or 95% CI
error_option = 'CI' # change to 'std' for standard deviation
plt.figure(figsize=(10,6))
for method in peaks_methods:
df_method = df[df["Method"] == method]
df_method = df_method[df_method['Std deviation'].notna()]
# To create smooth line plot
xnew = np.linspace(df_method["Noise decibels"].min(), df_method["Noise decibels"].max(), 300)
spl = make_interp_spline(df_method["Noise decibels"], df_method["Median distance"], k=3) # type: BSpline
y_smooth = spl(xnew)
if error_option == 'std':
spl_std = make_interp_spline(df_method["Noise decibels"], df_method["Std deviation"], k=3) # type: BSpline
error_smooth = spl_std(xnew)
elif error_option == 'CI':
confidence_interval = df_method["Std deviation"] / np.sqrt(n_times) * 1.96
spl_ci = make_interp_spline(df_method["Noise decibels"], confidence_interval, k=3) # type: BSpline
error_smooth = spl_ci(xnew)
line, = plt.plot(xnew, y_smooth, label=method) # use same color for the line
color = line.get_color() # get color of the line
plt.fill_between(xnew, y_smooth-error_smooth, y_smooth+error_smooth, color=color, alpha=0.2)
plt.title(f'Impact of noise on peak detection: {title}', size=16)
plt.grid(True)
plt.legend(loc="upper left", title="Peak Extraction Method")
plt.xlabel("Noise decibels")
plt.ylabel("Median distance between generated and retrieved peaks")
plt.show()
plt.close()
fixed
harmonic_recurrence
No peaks were detected. Consider increasing precision or number of harmonics
No peak detected
No peaks were detected. Consider increasing precision or number of harmonics
No peak detected
No peaks were detected. Consider increasing precision or number of harmonics
No peak detected
No peaks were detected. Consider increasing precision or number of harmonics
No peak detected
EIMC
EMD
FOOOF
fixed
harmonic_recurrence
No peaks were detected. Consider increasing precision or number of harmonics
No peak detected
No peaks were detected. Consider increasing precision or number of harmonics
No peak detected
No peaks were detected. Consider increasing precision or number of harmonics
No peak detected
No peaks were detected. Consider increasing precision or number of harmonics
No peak detected
EIMC
EMD
FOOOF
import seaborn as sns
import pandas as pd
import matplotlib.pyplot as plt
from neurodsp.sim import sim_combined
from scipy.interpolate import make_interp_spline, BSpline
# Setting up parameters
sf = 1000
length = 4
n_times = 30 #number of times the bootstraping is executed
freqs = [2, 6.5, 15.5, 18.3, 24.9, 39.1]
amps = [1, 1, 1, 1, 1, 1]
FREQ_BANDS = [[1, 5], [5, 10], [10, 20], [20, 30]]
# Simulate the signals
components = {
'sim_powerlaw': {'exponent': -1},
'sim_oscillation': [{'freq': freq} for freq in freqs] # list of dictionaries
}
component_variances = [1 for _ in range(len(freqs) + 1)] # +1 for the 'sim_powerlaw' component
# Simulate original signal
simulated_EEG = sim_combined(length, sf, components, component_variances)
# Simulate EEG signal
simulated_oscillations = generate_signal(sf, length, freqs, amps, show=False, color = 'turquoise') # Assuming generate_signal is your custom function
signals = {'Simulated Oscillations': simulated_oscillations, 'Simulated EEG Signal': simulated_EEG}
for title, signal in signals.items():
peaks_methods = ['fixed', 'harmonic_recurrence', 'EIMC', 'EMD', 'FOOOF']
noise_db_list = range(0, 40, 1)
colors = ['darkred', 'darkblue', 'darkorange', 'black', 'cyan', 'deeppink']
data = []
for method, color in zip(peaks_methods, colors):
print(method)
median_distance_tot = []
for db in noise_db_list:
t = 0
median_distance = []
while t < n_times:
# Convert to linear Watt units
target_noise_watts = 10 ** (db / 10)
# Generate noise samples
mean_noise = 0
noise_signal = np.random.normal(mean_noise, np.sqrt(target_noise_watts), len(signal))
# Noise up the original signal
noised_signal = signal + noise_signal
# set dictionary of prominence and rel_height for each method
prominences = {'fixed': 10, 'harmonic_recurrence': 8, 'EMD': 1, 'FOOOF': 10, 'HH1D_max': 1, 'EIMC': 7}
rel_heights = {'fixed': 1, 'harmonic_recurrence': 1, 'EMD': 0.7, 'FOOOF': 1, 'HH1D_max': 1, 'EIMC': 1}
# Initialize biotuner object
biotuning = compute_biotuner(sf = sf, peaks_function=method, precision=0.5)
if method == 'EMD':
biotuning.peaks_extraction(noised_signal, min_freq=1, max_freq=50, min_harms=2, harm_limit=128,
FREQ_BANDS = FREQ_BANDS, verbose=False, graph=False, prominence=prominences[method],
rel_height=rel_heights[method], nIMFs=6, keep_first_IMF=False, n_peaks=6)
biotuning.peaks = biotuning.peaks[0:4] # select only the first 4 peaks
biotuning.peak = [x for x in biotuning.peaks if x <= freqs[-1]+1]
else:
try:
biotuning.peaks_extraction(noised_signal, min_freq=1, max_freq=50, min_harms=2, harm_limit=128,
FREQ_BANDS = FREQ_BANDS, verbose=False, graph=False, prominence=prominences[method],
rel_height=rel_heights[method], nIMFs=3, keep_first_IMF=True, n_peaks=4)
biotuning.peaks = [x for x in biotuning.peaks if x <= freqs[-1]+1]
except:
biotuning.peaks = []
if db == 0:
no_noise_peaks = biotuning.peaks
if len(biotuning.peaks) == 0:
median_distance.append(np.nan)
else:
list_distance = []
for p in biotuning.peaks:
takeClosest = lambda num,collection:min(collection,key=lambda x:abs(x-num))
closest = takeClosest(p,freqs)
list_distance.append(abs(p-closest))
median_distance.append(np.nanmedian(list_distance))
t += 1
median_distance_tot.append(np.nanmedian(median_distance))
# calculate standard deviation
std_dev = np.std(median_distance)
# Continue with your computations...
data.append([method, db, np.nanmedian(median_distance), np.std(median_distance)])
df = pd.DataFrame(data, columns=["Method", "Noise decibels", "Median distance", "Std deviation"])
# Add the option for std or 95% CI
error_option = 'CI' # change to 'std' for standard deviation
plt.figure(figsize=(10,6))
for method in peaks_methods:
df_method = df[df["Method"] == method]
df_method = df_method[df_method['Std deviation'].notna()]
# To create smooth line plot
xnew = np.linspace(df_method["Noise decibels"].min(), df_method["Noise decibels"].max(), 300)
spl = make_interp_spline(df_method["Noise decibels"], df_method["Median distance"], k=3) # type: BSpline
y_smooth = spl(xnew)
if error_option == 'std':
spl_std = make_interp_spline(df_method["Noise decibels"], df_method["Std deviation"], k=3) # type: BSpline
error_smooth = spl_std(xnew)
elif error_option == 'CI':
confidence_interval = df_method["Std deviation"] / np.sqrt(n_times) * 1.96
spl_ci = make_interp_spline(df_method["Noise decibels"], confidence_interval, k=3) # type: BSpline
error_smooth = spl_ci(xnew)
line, = plt.plot(xnew, y_smooth, label=method) # use same color for the line
color = line.get_color() # get color of the line
plt.fill_between(xnew, y_smooth-error_smooth, y_smooth+error_smooth, color=color, alpha=0.2)
plt.title(f'Impact of noise on peak detection: {title}')
plt.grid(True)
plt.legend(loc="upper left", title="Peak Extraction Method")
plt.xlabel("Noise decibels")
plt.ylabel("Median distance")
plt.show()
plt.close()
fixed
harmonic_recurrence
No peaks were detected. Consider increasing precision or number of harmonics
No peak detected
No peaks were detected. Consider increasing precision or number of harmonics
No peak detected
No peaks were detected. Consider increasing precision or number of harmonics
No peak detected
EIMC
EMD
FOOOF
fixed
harmonic_recurrence
No peaks were detected. Consider increasing precision or number of harmonics
No peak detected
No peaks were detected. Consider increasing precision or number of harmonics
No peak detected
EIMC
EMD
FOOOF
The second plot is for a signal made from a set of inharmonic frequencies.