Je veux ajouter du bruit aléatoire à quelque 100 signaux bin que je simule en Python - pour le rendre plus réaliste.
À un niveau de base, ma première pensée a été d'aller bac par bac et de générer simplement un nombre aléatoire entre une certaine plage et d'ajouter ou de soustraire cela du signal.
J'espérais (car c'est python) qu'il pourrait y avoir un moyen plus intelligent de le faire via numpy ou quelque chose. (Je suppose qu'idéalement, un nombre tiré d'une distribution gaussienne et ajouté à chaque casier serait également préférable.)
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 manière plus sophistiquée de générer le bruit.
En termes de sortie, si j'avais 10 bacs avec les valeurs suivantes:
Bac 1: 1 Bac 2: 4 Bac 3: 9 Bac 4: 16 Bac 5:25 Bac 6:25 Bac 7: 16 Bac 8: 9 Bac 9: 4 Bac 10: 1
Je me demandais simplement s'il y avait une fonction prédéfinie qui pourrait ajouter du bruit pour me donner quelque chose comme:
Bac 1: 1,13 Bac 2: 4,21 Bac 3: 8,79 Bac 4: 16,08 Bac 5: 24,97 Bac 6: 25,14 Bac 7: 16,22 Bac 8: 8,90 Bac 9: 4,02 Bac 10: 0,91
Sinon, je vais simplement aller bin par bin et ajouter un nombre sélectionné à partir d'une distribution gaussienne à chacun.
Merci.
C'est en fait un signal d'un radiotélescope que je simule. Je veux pouvoir éventuellement choisir le rapport signal sur bruit de ma simulation.
Réponses:
Vous pouvez générer un tableau de bruit et l'ajouter à votre signal
import numpy as np noise = np.random.normal(0,1,100) # 0 is the mean of the normal distribution you are choosing from # 1 is the standard deviation of the normal distribution # 100 is the number of elements you get in array noisela source
... Et pour ceux qui - comme moi - sont très tôt dans leur courbe d'apprentissage engourdie,
import numpy as np pure = np.linspace(-1, 1, 100) noise = np.random.normal(0, 1, 100) signal = pure + noisela source
Pour ceux qui essaient de faire la connexion entre SNR et une variable aléatoire normale générée par numpy:
[1]
, où il est important de garder à l'esprit que P est la puissance moyenne .
Ou en dB:
[2]
Dans ce cas, nous avons déjà un signal et nous voulons générer du bruit pour nous donner un SNR souhaité.
Bien que le bruit puisse prendre différentes saveurs en fonction de ce que vous modélisez, un bon début (en particulier pour cet exemple de radiotélescope) est le bruit gaussien blanc additif (AWGN) . Comme indiqué dans les réponses précédentes, pour modéliser l'AWGN, vous devez ajouter une variable aléatoire gaussienne de moyenne zéro à votre signal d'origine. La variance de cette variable aléatoire affectera la puissance moyenne du bruit.
Pour une variable aléatoire gaussienne X, la puissance moyenne
, également appelée second moment , est 
[3]
Donc pour le bruit blanc,
et la puissance moyenne est alors égale à la variance 
.
Lors de la modélisation en python, vous pouvez soit
1. Calculer la variance en fonction d'un SNR souhaité et d'un ensemble de mesures existantes, ce qui fonctionnerait si vous vous attendez à ce que vos mesures aient des valeurs d'amplitude assez cohérentes.
2. Vous pouvez également régler la puissance du bruit à un niveau connu pour correspondre à quelque chose comme le bruit du récepteur. Le bruit du récepteur peut être mesuré en pointant le télescope dans l'espace libre et en calculant la puissance moyenne.
Dans tous les cas, il est important de vous assurer d'ajouter du bruit à votre signal et de prendre des moyennes dans l'espace linéaire et non en unités dB.
Voici un code pour générer un signal et tracer la tension, la puissance en watts et la puissance en dB:
# Signal Generation # matplotlib inline import numpy as np import matplotlib.pyplot as plt t = np.linspace(1, 100, 1000) x_volts = 10*np.sin(t/(2*np.pi)) plt.subplot(3,1,1) plt.plot(t, x_volts) plt.title('Signal') plt.ylabel('Voltage (V)') plt.xlabel('Time (s)') plt.show() x_watts = x_volts ** 2 plt.subplot(3,1,2) plt.plot(t, x_watts) plt.title('Signal Power') plt.ylabel('Power (W)') plt.xlabel('Time (s)') plt.show() x_db = 10 * np.log10(x_watts) plt.subplot(3,1,3) plt.plot(t, x_db) plt.title('Signal Power in dB') plt.ylabel('Power (dB)') plt.xlabel('Time (s)') plt.show()Voici un exemple d'ajout d'AWGN basé sur un SNR souhaité:
# Adding noise using target SNR # Set a target SNR target_snr_db = 20 # Calculate signal power and convert to dB sig_avg_watts = np.mean(x_watts) sig_avg_db = 10 * np.log10(sig_avg_watts) # Calculate noise according to [2] then convert to watts noise_avg_db = sig_avg_db - target_snr_db noise_avg_watts = 10 ** (noise_avg_db / 10) # Generate an sample of white noise mean_noise = 0 noise_volts = np.random.normal(mean_noise, np.sqrt(noise_avg_watts), len(x_watts)) # Noise up the original signal y_volts = x_volts + noise_volts # Plot signal with noise plt.subplot(2,1,1) plt.plot(t, y_volts) plt.title('Signal with noise') plt.ylabel('Voltage (V)') plt.xlabel('Time (s)') plt.show() # Plot in dB y_watts = y_volts ** 2 y_db = 10 * np.log10(y_watts) plt.subplot(2,1,2) plt.plot(t, 10* np.log10(y_volts**2)) plt.title('Signal with noise (dB)') plt.ylabel('Power (dB)') plt.xlabel('Time (s)') plt.show()Et voici un exemple d'ajout d'AWGN basé sur une puissance de bruit connue:
# Adding noise using a target noise power # Set a target channel noise power to something very noisy target_noise_db = 10 # Convert to linear Watt units target_noise_watts = 10 ** (target_noise_db / 10) # Generate noise samples mean_noise = 0 noise_volts = np.random.normal(mean_noise, np.sqrt(target_noise_watts), len(x_watts)) # Noise up the original signal (again) and plot y_volts = x_volts + noise_volts # Plot signal with noise plt.subplot(2,1,1) plt.plot(t, y_volts) plt.title('Signal with noise') plt.ylabel('Voltage (V)') plt.xlabel('Time (s)') plt.show() # Plot in dB y_watts = y_volts ** 2 y_db = 10 * np.log10(y_watts) plt.subplot(2,1,2) plt.plot(t, 10* np.log10(y_volts**2)) plt.title('Signal with noise') plt.ylabel('Power (dB)') plt.xlabel('Time (s)') plt.show()la source
Pour ceux qui souhaitent ajouter du bruit à un ensemble de données multidimensionnel chargé dans une trame de données pandas ou même un ndarray numpy, voici un exemple:
import pandas as pd # create a sample dataset with dimension (2,2) # in your case you need to replace this with # clean_signal = pd.read_csv("your_data.csv") clean_signal = pd.DataFrame([[1,2],[3,4]], columns=list('AB'), dtype=float) print(clean_signal) """ print output: A B 0 1.0 2.0 1 3.0 4.0 """ import numpy as np mu, sigma = 0, 0.1 # creating a noise with the same dimension as the dataset (2,2) noise = np.random.normal(mu, sigma, [2,2]) print(noise) """ print output: array([[-0.11114313, 0.25927152], [ 0.06701506, -0.09364186]]) """ signal = clean_signal + noise print(signal) """ print output: A B 0 0.888857 2.259272 1 3.067015 3.906358 """la source
Dans la vraie vie, vous souhaitez simuler un signal avec un bruit blanc. Vous devez ajouter à votre signal des points aléatoires qui ont une distribution gaussienne normale. Si nous parlons d'un appareil dont la sensibilité est donnée en unité / SQRT (Hz), vous devez alors en déduire l'écart type de vos points. Ici, je donne la fonction "white_noise" qui fait cela pour vous, et le reste d'un code est une démonstration et vérifie s'il fait ce qu'il doit.
%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()la source
Super 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 utilisé. Le partage en cas utile aux autres également
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_awayla source
AWGN similaire à la fonction Matlab
def awgn(sinal): regsnr=54 sigpower=sum([math.pow(abs(sinal[i]),2) for i in range(len(sinal))]) sigpower=sigpower/len(sinal) noisepower=sigpower/(math.pow(10,regsnr/10)) noise=math.sqrt(noisepower)*(np.random.uniform(-1,1,size=len(sinal))) return noisela source