Добавление шума к сигналу в python
Я хочу добавить некоторый случайный шум к некоторому 100 бин-сигналу, который я имитирую в Python, чтобы сделать его более реалистичным.
На базовом уровне моя первая мысль заключалась в том, чтобы запустить bin bin и просто создать случайное число между определенным диапазоном и добавить или вычесть это из сигнала.
Я надеялся (поскольку это python), что может быть более разумный способ сделать это через numpy или что-то еще. (Я полагаю, что в идеале число, полученное из гауссовского распределения и добавленное к каждому бинду, было бы лучше.)
Благодарим вас за любые ответы.
Я просто нахожусь в стадии планирования моего кода, поэтому мне нечего показать. Я просто думал, что может быть более сложный способ генерации шума.
В терминах вывода, если у меня было 10 бункеров со следующими значениями:
Бин 1:1
Бин 2: 4
Бин 3: 9
Бин 4: 16
Бин 5: 25
Бин 6: 25
Бин 7: 16
Бин 8: 9
Бин 9: 4
Бин 10: 1
Я просто задавался вопросом, была ли предустановленная функция, которая могла бы добавить шум, чтобы дать мне что-то вроде:
Бин 1:1,13
Бин 2: 4,21
Бин 3: 8,79
Бин 4: 16.08
Бин 5: 24,97
Бин 6: 25,14
Бин 7: 16.22
Бин 8: 8,90
Бин 9: 4.02
Бин 10: 0.91
Если нет, я просто перейду по отдельности и добавлю число, выбранное из gaussian-дистрибутива для каждого из них.
Спасибо.
На самом деле это сигнал от радиотелескопа, который я имитирую. Я хочу, чтобы в конечном итоге выбрать отношение сигнал/шум моего моделирования.
Ответы
Ответ 1
Вы можете создать массив шума и добавить его к вашему сигналу
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
Ответ 2
... И для тех, кто, как и я, - очень рано в своей кривой обучения numpy,
import numpy as np
pure = np.linspace(-1, 1, 100)
noise = np.random.normal(0, 1, pure.shape)
signal = pure + noise
Ответ 3
Для тех, кто пытается установить связь между SNR и обычной случайной величиной, сгенерированной numpy:
[1]
, где важно помнить, что P - это средняя мощность.
Или в дБ:
[2] ![SNR dB2]()
В этом случае у нас уже есть сигнал, и мы хотим генерировать шум, чтобы получить желаемое SNR.
В то время как шум может проявляться в разных вариантах в зависимости от того, что вы моделируете, хорошее начало (особенно для этого примера радиотелескопа) - это аддитивный белый гауссов шум (AWGN). Как указывалось в предыдущих ответах, для моделирования AWGN необходимо добавить гауссовскую случайную переменную с нулевым средним значением к исходному сигналу. Дисперсия этой случайной величины будет влиять на среднюю мощность шума.
Для гауссовой случайной величины X средняя мощность
, также известная как второй момент, равна
[3] ![Ex]()
Таким образом, для белого шума
и средняя мощность тогда равны дисперсии
.
При моделировании этого в Python вы можете либо
1. Рассчитайте дисперсию на основе требуемого SNR и набора существующих измерений, которые сработают, если вы ожидаете, что ваши измерения будут иметь достаточно постоянные значения амплитуды.
2. В качестве альтернативы, вы можете установить мощность шума на известном уровне, чтобы соответствовать чему-то вроде шума приемника. Шум приемника можно измерить, направив телескоп в свободное пространство и рассчитав среднюю мощность.
В любом случае важно убедиться, что вы добавляете шум к своему сигналу и берете средние значения в линейном пространстве, а не в единицах дБ.
Вот некоторый код для генерации сигнала и графика напряжения, мощности в ваттах и мощности в дБ:
# 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()
![Generated Signal]()
Вот пример для добавления AWGN на основе желаемого 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()
![Signal with target SNR]()
А вот пример добавления AWGN на основе известной мощности шума:
# 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()
![Signal with target noise level]()
Ответ 4
Для тех, кто хочет добавить шум к многомерному набору данных, загруженному в фрейм данных pandas или даже крошечный ndarray, вот пример:
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
"""
Ответ 5
Потрясающие ответы выше. Недавно у меня возникла необходимость генерировать смоделированные данные, и именно этим я и воспользовался. Поделиться делом полезно для других,
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