如何用numpy产生频率范围内的噪声？

Question

我有一个主信号，例如周期为 200 个样本的正弦信号。

我想向该信号添加噪声。 “噪声信号部分”的周期应在例如 5-30 个样本的范围内。

我认为这足以在这个范围内生成具有不同随机选择幅度的多个正弦：

noise = np.sin(np.array(range(N))/0.7)*np.random.random(1) + np.sin(np.array(range(N))/1.1)*np.random.random(1) + np.sin(np.array(range(N))/1.5)*np.random.random(1) 

但对于我的目的来说，这个解决方案仍然过于“确定性”。

如何生成幅度和周期随机变化的噪声？

Answer 1

在 MathWorks 的文件交换中：

fftnoise - generate noise with a specified power spectrum

，您可以找到 Aslak Grinsted 的 matlab 代码，该代码创建具有指定功率谱的噪声。它可以很容易地移植到Python：

def fftnoise(f):
    f = np.array(f, dtype='complex')
    Np = (len(f) - 1) // 2
    phases = np.random.rand(Np) * 2 * np.pi
    phases = np.cos(phases) + 1j * np.sin(phases)
    f[1:Np+1] *= phases
    f[-1:-1-Np:-1] = np.conj(f[1:Np+1])
    return np.fft.ifft(f).real

您可以将其用于您的案例，如下所示：

def band_limited_noise(min_freq, max_freq, samples=1024, samplerate=1):
    freqs = np.abs(np.fft.fftfreq(samples, 1/samplerate))
    f = np.zeros(samples)
    idx = np.where(np.logical_and(freqs>=min_freq, freqs<=max_freq))[0]
    f[idx] = 1
    return fftnoise(f)

据我所知似乎有效。聆听您新创建的噪音：

from scipy.io import wavfile

x = band_limited_noise(200, 2000, 44100, 44100)
x = np.int16(x * (2**15 - 1))
wavfile.write("test.wav", 44100, x)

Answer 2

您应该使用随机相位，而不是使用具有不同幅度的多个正弦：

import numpy as np
from functools import reduce

def band_limited_noise(min_freq, max_freq, samples=44100, samplerate=44100):
    t = np.linspace(0, samples/samplerate, samples)
    freqs = np.arange(min_freq, max_freq+1, samples/samplerate)
    phases = np.random.rand(len(freqs))*2*np.pi
    signals = [np.sin(2*np.pi*freq*t + phase) for freq,phase in zip(freqs,phases)]
    signal = reduce(lambda a,b: a+b,signals)
    signal /= np.max(signal)
    return signal

背景：白噪声意味着功率谱包含每个频率，因此如果您想要带限噪声，您可以将频带内的每个频率加在一起。噪声部分来自随机相位。由于 DFT 是离散的，因此您只需考虑给定采样率时实际出现的离散频率。

Answer 3

这是对这个“标题”的答案的改进：

“您应该使用随机相位，而不是使用具有不同幅度的多个正弦：”

这种编写方式的问题在于它会消耗大量的 RAM。在这里我已经修复了它。它仍然使用大量的CPU：

from scipy.io import wavfile
import numpy as np

# It takes my CPU about 50 sec to generate 1 sec of "noise".
def band_limited_noise(min_freq, max_freq, freq_step, samples=44100, samplerate=44100):
    t = np.linspace(0, samples/samplerate, samples)
    freqs = np.arange(min_freq, max_freq+1, freq_step)
    no_of_freqs = len(freqs)
    pi_2 = 2*np.pi
    phases = np.random.rand(no_of_freqs)*pi_2 # I am not sure why it is necessary to add randomness?
    signals = np.zeros(samples)
    for i in range(no_of_freqs):
        signals = signals + np.sin(pi_2*freqs[i]*t + phases[i])
    max_sig = np.max(signals)
    signals /= max_sig
    return signals


if __name__ == '__main__':

    # CONFIGURE HERE:
    seconds = 3
    samplerate = 44100
    freq_step = 0.25

    x = band_limited_noise(10, 20000, freq_step, seconds * samplerate, samplerate)
    wavfile.write("add all frequencies 3 sec freq_step 0,25.wav", 44100, x)

我分析了生成的“噪声”，并将其与使用以下更短、更快的代码创建的“普通”白噪声进行了比较，结果并没有更好：

from scipy.io import wavfile
from scipy import stats

sample_rate = 44100
length_in_seconds = 3
noise = stats.truncnorm(-3, 3, scale = 0.2).rvs(size = sample_rate * length_in_seconds)
wavfile.write('noise float64 -3 +3 0.2.wav', sample_rate, noise)

我所说的结果不是更好是指，当您使用 scipy.FFT 生成 2 个 WAV 文件（均为 3 秒）的频谱时，不同频率的信号“量”之间的变化或多或少相同。并且使 freq_step 变得更小也没有帮助。我或多或少地使用以下方法分析了变化：

distfit.distfit.fit_transform(numpy.abs(scipy.FFT.rfft(data_from_wav_file)))

这是详细信息：

import soundfile as sf
from scipy.fft import rfft, rfftfreq
import numpy as np
from matplotlib import pyplot as plot
import distfit

if __name__ == '__main__':

    # CONFIGURE HERE:
    soundfile = "add all frequencies 3 sec freq_step 0,125.wav"
    print("Input file: ", soundfile)
    with sf.SoundFile(soundfile, 'r') as f:
        new_samplerate = f.samplerate
        data = f.read()
    cutofflen = len(data)
    xf = rfftfreq(cutofflen, 1 / new_samplerate)
    yf = rfft(data)
    yf = np.abs(yf) # Length of a vector...
    print("Fit variations in the FFT of WAV-file to statistical distribution:")
    FFT_dist = distfit.distfit()  # Create the distfit object
    FFT_dist.fit_transform(yf)  # Call the fit transform function
    FFT_dist.plot() # dfit.plot_summary()  # Plot the summary

您可以使用以下方法看到频率响应的巨大变化：

plot.plot(xf, yf)
plot.show()

根据提到的 FFT 分析，这两种类型的噪声实际上都没有接近我想要的：特定频带中所有频率的数量完全相同，例如从 10 Hz 到 20 KHz。

Answer 4

产生全频谱白噪声然后过滤它就像你想把你房子的墙漆成白色，所以你决定把整个房子漆成白色，然后把除了墙壁之外的所有房子都漆回来。是白痴。 （但在电子学中有意义）。

我制作了一个小型 C 程序，可以在任何频率和任何带宽下生成白噪声（假设中心频率为 16kHz，“宽”为 2kHz）。 不涉及过滤。

我所做的很简单： 在主（无限）循环内，我生成一个中心频率+/--半带宽和+半带宽之间的随机数的正弦曲线，然后我为任意数量的样本（粒度）保留该频率，这就是结果：

16kHz 中心频率下 2kHz 宽的白噪声

伪代码：

while (true)
{

    f = center frequency
    r = random number between -half of bandwidth and + half of bandwidth

<secondary loop (for managing "granularity")>

    for x = 0 to 8 (or 16 or 32....)

    {

        [generate sine Nth value at frequency f+r]

        output = generated Nth value
    }


}