ajouter du bruit à un signal en python

Question

Je veux ajouter un bruit aléatoire à un signal de 100 cases que je simule en Python, pour le rendre plus réaliste.

Au niveau de base, ma première idée était de procéder bac par bac et de générer un nombre aléatoire compris dans une certaine fourchette, puis de l'ajouter ou de le soustraire au signal.

J'espérais (comme il s'agit de python) qu'il y aurait un moyen plus intelligent de le faire via numpy ou quelque chose comme ça. (Je suppose qu'idéalement, un nombre tiré d'une distribution gaussienne et ajouté à chaque bin serait également mieux).

Merci d'avance pour toute réponse.

---

Je suis juste au stade de la planification de mon code, donc je n'ai rien à montrer. Je pensais juste qu'il pourrait y avoir une façon plus sophistiquée de générer le bruit.

En termes de résultats, si j'avais 10 cases avec les valeurs suivantes :

Bin 1 : 1 Bac 2 : 4 Bin 3 : 9 Bin 4 : 16 Bin 5 : 25 Bin 6 : 25 Bin 7 : 16 Bin 8 : 9 Bin 9 : 4 Bac 10 : 1

Je me demandais juste s'il existait une fonction prédéfinie qui pourrait ajouter du bruit pour me donner quelque chose comme :

Bac 1 : 1,13 Bin 2 : 4,21 Bin 3 : 8.79 Bin 4 : 16,08 Bin 5 : 24,97 Bac 6 : 25,14 Bac 7 : 16,22 Case 8 : 8.90 Bac 9 : 4,02 Bac 10 : 0,91

Sinon, je vais simplement procéder bac par bac et ajouter à chacun d'eux un nombre choisi dans une distribution gaussienne.

Merci.

---

C'est en fait le signal d'un radiotélescope que je simule. Je veux pouvoir éventuellement choisir le rapport signal/bruit de ma simulation.

Answer 1

Dans la vie réelle, vous souhaitez simuler un signal avec un bruit blanc. Vous devez ajouter à votre signal des points aléatoires qui ont une distribution normale gaussienne. Si nous parlons d'un appareil dont la sensibilité est donnée en unité/SQRT(Hz), alors vous devez déterminer l'écart-type de vos points par rapport à celui-ci. Ici, je donne la fonction "white_noise" qui fait cela pour vous, et le reste du code est la démonstration et la vérification si elle fait ce qu'elle devrait.

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

"""
parameters: 
rhp - spectral noise density unit/SQRT(Hz)
sr  - sample rate
n   - no of points
mu  - mean value, optional

returns:
n points of noise signal with spectral noise density of rho
"""
def white_noise(rho, sr, n, mu=0):
    sigma = rho * np.sqrt(sr/2)
    noise = np.random.normal(mu, sigma, n)
    return noise

rho = 1 
sr = 1000
n = 1000
period = n/sr
time = np.linspace(0, period, n)
signal_pure = 100*np.sin(2*np.pi*13*time)
noise = white_noise(rho, sr, n)
signal_with_noise = signal_pure + noise

f, psd = signal.periodogram(signal_with_noise, sr)

print("Mean spectral noise density = ",np.sqrt(np.mean(psd[50:])), "arb.u/SQRT(Hz)")

plt.plot(time, signal_with_noise)
plt.plot(time, signal_pure)
plt.xlabel("time (s)")
plt.ylabel("signal (arb.u.)")
plt.show()

plt.semilogy(f[1:], np.sqrt(psd[1:]))
plt.xlabel("frequency (Hz)")
plt.ylabel("psd (arb.u./SQRT(Hz))")
#plt.axvline(13, ls="dashed", color="g")
plt.axhline(rho, ls="dashed", color="r")
plt.show()

Answer 2

Réponses impressionnantes d'Akavall et de Noel (c'est ce qui a fonctionné pour moi). J'ai également vu quelques commentaires sur les différentes distributions. Une solution que j'ai aussi essayée était de faire un test sur ma variable et de trouver de quelle distribution elle était la plus proche.

numpy.random

a différentes distributions qui peuvent être utilisées, cela peut être vu dans sa documentation :

documentation numpy.random

A titre d'exemple tiré d'une autre distribution (exemple référencé dans la réponse de Noel) :

import numpy as np
pure = np.linspace(-1, 1, 100)
noise = np.random.lognormal(0, 1, 100)
signal = pure + noise
print(pure[:10])
print(signal[:10])

J'espère que cela pourra aider quelqu'un qui recherche cette branche spécifique de la question originale.

Answer 3

Superbes réponses ci-dessus. J'ai récemment eu besoin de générer des données simulées et c'est ce que j'ai fini par utiliser. Je partage au cas où cela serait utile à d'autres,

import logging
__name__ = "DataSimulator"
logging.basicConfig(level=logging.INFO)
logger = logging.getLogger(__name__)

import numpy as np
import pandas as pd

def generate_simulated_data(add_anomalies:bool=True, random_state:int=42):
    rnd_state = np.random.RandomState(random_state)
    time = np.linspace(0, 200, num=2000)
    pure = 20*np.sin(time/(2*np.pi))

    # concatenate on the second axis; this will allow us to mix different data 
    # distribution
    data = np.c_[pure]
    mu = np.mean(data)
    sd = np.std(data)
    logger.info(f"Data shape : {data.shape}. mu: {mu} with sd: {sd}")
    data_df = pd.DataFrame(data, columns=['Value'])
    data_df['Index'] = data_df.index.values

    # Adding gaussian jitter
    jitter = 0.3*rnd_state.normal(mu, sd, size=data_df.shape[0])
    data_df['with_jitter'] = data_df['Value'] + jitter

    index_further_away = None
    if add_anomalies:
        # As per the 68-95-99.7 rule(also known as the empirical rule) mu+-2*sd 
        # covers 95.4% of the dataset.
        # Since, anomalies are considered to be rare and typically within the 
        # 5-10% of the data; this filtering
        # technique might work 
        #for us(https://en.wikipedia.org/wiki/68%E2%80%9395%E2%80%9399.7_rule)
        indexes_furhter_away = np.where(np.abs(data_df['with_jitter']) > (mu + 
         2*sd))[0]
        logger.info(f"Number of points further away : 
        {len(indexes_furhter_away)}. Indexes: {indexes_furhter_away}")
        # Generate a point uniformly and embed it into the dataset
        random = rnd_state.uniform(0, 5, 1)
        data_df.loc[indexes_furhter_away, 'with_jitter'] +=  
        random*data_df.loc[indexes_furhter_away, 'with_jitter']
    return data_df, indexes_furhter_away