Chcę dodać losowy szum do około 100 bin sygnału, który symuluję w Pythonie - żeby był bardziej realistyczny.

Na podstawowym poziomie, moją pierwszą myślą było przejście do segmentu segment po segmencie i po prostu wygenerowanie liczby losowej z określonego zakresu i dodanie lub odjęcie tej liczby od sygnału.

Miałem nadzieję (ponieważ to jest Python), że może być bardziej inteligentny sposób na zrobienie tego za pomocą numpy lub czegoś podobnego. (Przypuszczam, że najlepiej byłoby, gdyby liczba pobrana z rozkładu Gaussa i dodana do każdego pojemnika byłaby lepsza).

Z góry dziękuję za wszelkie odpowiedzi.

Jestem dopiero na etapie planowania kodu, więc nie mam nic do pokazania. Właśnie pomyślałem, że może istnieć bardziej wyrafinowany sposób generowania hałasu.

Pod względem wydajności, gdybym miał 10 pojemników o następujących wartościach:

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

Zastanawiałem się tylko, czy istnieje wstępnie zdefiniowana funkcja, która mogłaby dodać szum, aby dać mi coś takiego:

Pojemnik 1: 1,13 Pojemnik 2: 4,21 Pojemnik 3: 8,79 Pojemnik 4: 16,08 Pojemnik 5: 24,97 Pojemnik 6: 25,14 Pojemnik 7: 16,22 Pojemnik 8: 8,90 Pojemnik 9: 4,02 Pojemnik 10: 0,91

Jeśli nie, po prostu przejdę bin-by-bin i dodam do każdego liczbę wybraną z rozkładu Gaussa.

Dziękuję Ci.

Właściwie to sygnał z radioteleskopu, który symuluję. Chcę móc ostatecznie wybrać stosunek sygnału do szumu mojej symulacji.

Możesz wygenerować tablicę szumów i dodać ją do swojego sygnału

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 noise

... a dla tych, którzy - jak ja - są na bardzo wczesnym etapie uczenia się odrętwienia,

import numpy as np
pure = np.linspace(-1, 1, 100)
noise = np.random.normal(0, 1, 100)
signal = pure + noise

Dla tych, którzy próbują nawiązać połączenie między SNR z normalną zmienną losową generowaną przez numpy:

[1] , gdzie należy pamiętać, że P to moc średnia .

Lub w dB:
[2]

W tym przypadku mamy już sygnał i chcemy wygenerować szum, aby uzyskać pożądany SNR.

Chociaż szum może mieć różne smaki w zależności od tego, co modelujesz, dobrym początkiem (szczególnie w przypadku tego przykładu radioteleskopu) jest Addytywny Biały Szum Gaussowski (AWGN) . Jak stwierdzono w poprzednich odpowiedziach, aby modelować AWGN, należy dodać zerową średnią gaussowską zmienną losową do oryginalnego sygnału. Wariancja tej zmiennej losowej wpłynie na średnią moc szumów.

Dla losowej zmiennej Gaussa X, średnia moc , znana również jako drugi moment , wynosi
[3]

Więc dla białego szumu, a średnia moc jest równa wariancji .

Podczas modelowania tego w Pythonie możesz:
1. Obliczyć wariancję na podstawie pożądanego SNR i zestawu istniejących pomiarów, co zadziała, jeśli spodziewasz się, że pomiary będą miały dość spójne wartości amplitudy.
2. Alternatywnie, możesz ustawić moc szumu na znanym poziomie, aby dopasować coś takiego jak szum odbiornika. Szum odbiornika można zmierzyć, kierując teleskop w wolną przestrzeń i obliczając średnią moc.

Tak czy inaczej, ważne jest, aby upewnić się, że dodajesz szum do sygnału i bierzesz średnie w przestrzeni liniowej, a nie w jednostkach dB.

Oto kod generujący sygnał i wykreślający napięcie, moc w watach i moc w 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()

Oto przykład dodawania AWGN na podstawie pożądanego SNR:

# 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()

A oto przykład dodawania AWGN na podstawie znanej mocy szumów:

# 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()

Dla tych, którzy chcą dodać szum do wielowymiarowego zestawu danych załadowanego w ramce danych pandy, a nawet numpy ndarray, oto przykład:

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
""" 

W prawdziwym życiu chcesz zasymulować sygnał z białym szumem. Powinieneś dodać do swojego sygnału losowe punkty, które mają normalny rozkład Gaussa. Jeśli mówimy o urządzeniu, które ma czułość podaną w jednostkach / SQRT (Hz), to musisz obliczyć z niego odchylenie standardowe punktów. Tutaj podaję funkcję "white_noise", która robi to za Ciebie, reszta kodu to demonstracja i sprawdzam, czy robi to, co powinna.

%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()

Świetne odpowiedzi powyżej. Niedawno miałem potrzebę generowania symulowanych danych i właśnie z tego skorzystałem. Dzielenie się w razie potrzeby również innym,

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

AWGN Podobny do funkcji 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 noise