adding noise to a signal in python

I want to add some random noise to some 100 bin signal that I am simulating in Python - to make it more realistic.

On a basic level, my first thought was to go bin by bin and just generate a random number between a certain range and add or subtract this from the signal.

I was hoping (as this is python) that there might a more intelligent way to do this via numpy or something. (I suppose that ideally a number drawn from a gaussian distribution and added to each bin would be better also.)

Thank you in advance of any replies.


I'm just at the stage of planning my code, so I don't have anything to show. I was just thinking that there might be a more sophisticated way of generating the noise.

In terms out output, if I had 10 bins with the following values:

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

I just wondered if there was a pre-defined function that could add noise to give me something like:

Bin 1: 1.13 Bin 2: 4.21 Bin 3: 8.79 Bin 4: 16.08 Bin 5: 24.97 Bin 6: 25.14 Bin 7: 16.22 Bin 8: 8.90 Bin 9: 4.02 Bin 10: 0.91

If not, I will just go bin-by-bin and add a number selected from a gaussian distribution to each one.

Thank you.


It's actually a signal from a radio telescope that I am simulating. I want to be able to eventually choose the signal to noise ratio of my simulation.

307786 次浏览

您可以生成一个噪声阵列,并将其添加到您的信号

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

对于那些像我一样还处于麻木学习阶段的人来说,

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

对于那些想要在熊猫数据框中或者甚至是数字图像中为多维数据集添加噪音的人,这里有一个例子:

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

对于那些试图在 SNR和由 numpy 产生的普通随机变量之间建立联系的人:

[1] < img src = “ https://chart.googleapis.com/chart? cht = tx & amp; chl =% 7BSNR% 7D% 20 =% 20% 5Cfrac% 7BP _% 5Cmathrm% 7B 信号% 7D% 7D% 7BP _% 5Cmathrm% 7B 噪音% 7D% 7D”alt = “ SNR ratio”> ,其中重要的是要记住 P 是 一般功率。

或分贝:
[2] < img src = “ https://chart.googleapis.com/chart? cht = tx & amp; chl =% 5Cmathrm% 7BSNR _% 7BdB% 7D% 20 =% 20% 7Bmathrm% 7B 信号% 2CdB% 7D% 20 -% 20P _% 5Cmathrm% 7B 噪音% 2CdB% 7D% 7D”alt = “ SNR dB2”>

在这种情况下,我们已经有了一个信号,我们希望产生噪声,给我们一个理想的信噪比。

虽然噪音可以来在不同的 口味取决于你建模,一个良好的开端(特别是对于这个射电望远镜的例子)是 Additive White Gaussian Noise (AWGN)。如前面的答案所述,为了建立 AWGN 模型,需要在原始信号中添加一个零均值高斯随机变量。该随机变量的方差将影响 一般噪声功率。

对于一个高斯随机变量 X,平均功率 < img src = “ https://chart.googleapis.com/chart? cht = tx & amp; chl = E% 5BX% 5E2% 5D”alt = “ Ep”> ,也称为第二个 片刻
[3] < img src = “ https://chart.googleapis.com/chart? cht = tx & amp; chl = E% 5BX% 5E2% 5D% 3D% 5Cmu% 5E2% 2B% 5Csigma% 5E2”alt = “ Ex”>

因此对于白噪声,平均功率等于方差。

当在 python 中对此进行建模时,您可以
1.计算方差基于一个理想的信噪比和一组现有的测量,这将工作,如果你期望你的测量有相当一致的振幅值。
2. Alternatively, you could set noise power to a known level to match something like receiver noise. Receiver noise could be measured by pointing the telescope into free space and calculating average power.

无论哪种方式,重要的是要确保您添加噪声到您的信号,并采取平均线性空间,而不是分贝单位。

下面是一些产生信号并绘制电压、功率(瓦特)和功率(分贝)的代码:

# 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 的例子:

# 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

上面的答案很棒。我最近需要生成模拟数据,这就是我最终使用的。分享对他人也有帮助,

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

类似 Matlab 函数的 AWGN

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

在现实生活中,你希望用白噪声来模拟信号。你应该在你的信号中加入具有正常正态分布的随机点。如果我们谈论的设备具有灵敏度单位/SQRT (赫兹) ,那么你需要设计出你的观点的标准差。在这里,我给出了一个函数“ white _ sound”,它可以为您完成这个任务,代码的其余部分是演示,并检查它是否完成了应该完成的任务。

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

Signal with noise

PSD

Akavall 和诺埃尔给出了很棒的回答(这对我很有效)。此外,我还看到了一些关于不同发行版的注释。我也尝试过的一个解决方案是对我的变量进行测试,找出它更接近的分布。

numpy.random

可以使用不同的发行版,可以在其文档中看到:

documentation numpy.random

作为一个来自不同发行版的例子(来自 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])

我希望这可以帮助有人从原来的问题中寻找这个特定的分支。

你可以试试这个:

import numpy as np


x = np.arange(-5.0, 5.0, 0.1)


y = np.power(x,2)


noise = 2 * np.random.normal(size=x.size)


ydata = y + noise


plt.plot(x, ydata,  'bo')
plt.plot(x,y, 'r')
plt.ylabel('y data')
plt.xlabel('x data')
plt.show()

enter image description here