在Python中给信号添加噪声

127

我希望在Python中对一些100个bin的信号添加一些随机噪声,使其更加逼真。

基本上,我的第一个想法是逐个bin地生成一个特定范围内的随机数,并将其加或减到信号中。

由于这是Python,我希望可能有更智能的方法通过numpy或其他模块来实现。(我认为最好每个bin添加从高斯分布中选出的数字)

非常感谢您提前的回复。


我正在计划代码,还没有什么可以展示的东西。我只是在想,也许有更复杂的方式来生成噪声。

就输出而言,如果我有10个bins,分别具有以下值:

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

我只是想知道是否有预定义的函数可以添加噪声以给我类似以下的结果:

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

如果没有,我将逐个bin地添加从高斯分布中选择的数字。

谢谢您。


实际上,我正在模拟来自射电望远镜的信号。我希望最终能够选择模拟的信噪比。


2
请展示你尝试过的代码或者你遇到的具体问题。同时提供样例输入和期望输出也会有很大帮助。 - g.d.d.c
3
它是什么类型的信号?你想引入什么类型的噪音?“真实”取决于信号的类型。例如,音频噪声与图像噪声不同。 - Diego Basch
9个回答

168

您可以生成噪声数组,并将其添加到信号中

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

24
在某些情境下,将信号乘以一个噪声数组(以1为中心)可能比添加噪声数组更有意义,但这当然取决于您试图模拟的噪声的性质。 - Edward Loper
我们如何生成非高斯随机噪声? - Dhiraj Dhakal
2
@EdwardLoper你能举个例子,说明何时乘噪声数组比加噪声数组更有意义吗? - Abe Brandsma

104
对于那些试图将SNR与由numpy生成的正常随机变量联系起来的人们:
[1] SNR ratio,这里需要牢记P是平均功率。
或者用dB表示:
[2] SNR dB2 在这种情况下,我们已经有一个信号,我们想生成噪声以给我们所需的SNR。
噪声的类型因所建模型而异,但对于这个射电望远镜的例子来说,一个很好的起点是加性白高斯噪声(AWGN)。如前面的回答所述,要建模 AWGN,需要将一个零均值高斯随机变量添加到原始信号中。该随机变量的方差将影响平均噪声功率。
对于高斯随机变量 X,平均功率Ep,也称为第二矩, 是
[3] Ex 因此,对于白噪声,Ex,平均功率等于方差Ex
当在Python中对此进行建模时,您可以选择:
1. 根据所需的SNR和一组现有测量值计算方差。如果您希望测量具有相当一致的幅度值,则此方法适用。
2. 或者,您可以将噪声功率设置为已知水平,以匹配接收器噪声之类的内容。可以通过将望远镜对准自由空间并计算平均功率来测量接收器噪声。

无论哪种方式,重要的是要确保向信号添加噪声,并在线性空间而不是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()

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


83

对于那些像我一样刚接触numpy的人,

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

13

如果想要向 Pandas DataFrame 或者 numpy 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

在现实生活中,您希望模拟具有白噪声的信号。您应该向信号中添加具有正态分布的随机点。如果我们谈论一个灵敏度以单位/根号赫兹给出的设备,则需要从中设计出您的点的标准差。这里我提供了函数“ white_noise”,它可以为您完成此操作,其余代码是演示和检查它是否做到了应该做的事情。

%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


4

AWGN类似于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

2
你可以尝试这样做:
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


2

Akavall和Noel提供了非常棒的答案(对我很有用)。同时,我看到了一些关于不同发行版的评论。我也尝试了一种解决方案,即对我的变量进行测试,并找出它更接近哪个分布。

numpy.random

numpy.random有不同的分布可以使用,这些可以在其文档中查看:

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])

我希望这篇翻译能够帮助那些寻找与原问题相关分支的人。

0

以上的回答很棒。最近我有需要生成模拟数据,这是我使用的方法。分享一下,希望对其他人有所帮助。

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

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接