How to generate noise in frequency range with numpy?

Question

I have a main signal, for example sinus with period of 200 samples.

I would like to add a noise to this signal. The periods of "noise signal parts" should be in range for example 5-30 samples.

I thought that will be enough to generate multiple sinuses in this range with different randomly chosen amplitudes:

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) 

But this solution is still too much "deterministic" for my purpose.

How could I generate noise with randomly changing amplitude and period?

Overdone

A typical approach would be to generate some white noise (e.g. using np.random.randn), then bandpass filter it in order to give it the desired frequency characteristics before adding it to your signal. – Biennial

@Biennial Yes, that is typical and completely correct approach. You are right. But I would like to avoid filtering if possible. So solution I want should be something simple like I suggest, but with better result (less deterministic). – Overdone

Why do you want to "avoid filtering"? – Biennial

@Biennial I want to use this noise to test a filter. According to my experience, filters do not remove all noise out of bandpass, or it delay the data, or it also suppress the frequencies around the bandpass border. Maybe I am wrong, but I believe that for relatively short data I will get cleaner result with some "cheating solution" than with proper filtering. – Overdone

I'm only talking about bandpass filtering the noise before you add it to your signal, so I don't see how phase shift could possibly be an issue. Your main concern seems to be that the noise will leak out into other spectral bands, but that really just depends on selecting an appropriate bandpass filter. If you want to generate something resembling band-limited white noise using individual random sinusoids then in you would need a lot of sinusoids (in principle, an infinite number of them). It would help if you could explain your exact needs in your question. – Biennial

@Biennial It seems that you are right, your solution seems to be the less complicated at the end. Do you want to create an answer, so I can mark it as solution? – Overdone

Sure, no worries. I will come back to it later today. – Biennial

Answer 1

In MathWorks' File Exchange: fftnoise - generate noise with a specified power spectrum you find matlab code from Aslak Grinsted, creating noise with a specified power spectrum. It can easily be ported to 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

You can use it for your case like this:

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)

Seems to work as far as I see. For listening to your freshly created noise:

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)

Suprematism

This behaves very strangely! For 'x = band_limited_noise(10, 20000, seconds*samplerate, samplerate)' and seconds=1 the min. and max. values in the x array are -0.0219 and 0.0189. (It should be around -1 and 1.) But with seconds=10 they are -0.0067 and 0.0067. A larger amount of samples should get bigger variations, not significantly smaller variation...!? – Bebop

Hi @henrik-r I think this arises from Numpy's definition of fft and ifft. Its definition of ifft includes the term 1/n. In contrast, in some mathematical contexts, both the fft and ifft have a term 1/sqrt(n). You'll probably want to renormalize your output by multiplying by sqrt(n) (n is the number of samples) – Bregenz

This approach also introduces a large spike at the zeroth sample of its output. I suggest this could be fixed by using (np.random.rand(Np) - 0.5) * 2 * np.pi on line 4. – Bregenz

Answer 2

Instead of using multiple sinuses with different amplitudes, you should use them with random phases:

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

Background: White noise means that the power spectrum contains every frequency, so if you want band limited noise you can add together every frequency within the band. The noisy part comes from the random phase. Because DFT is discrete, you only need to consider the discrete frequencies that actually occur given a sampling rate.

Wolfgang

HI! Can you explain in what way random phases is better? You also say: "The noisy part comes from the random phase." In reality I am not interested in any randomness... My ideal is a signal with the same 'amount' of every frequency within a certain frequency-band. Thank you. :-) – Bebop

# The following 2 lines are extremely heavy on both RAM and CPU! # In debug-mode: About 7 GB RAM for 441000 samples, but 5 GB for 44100 samples...!? signals = [np.sin(2*np.pi*freq*t + phase) for freq,phase in zip(freqs,phases)] signal = reduce(lambda a,b: a+b,signals) – Bebop

Answer 3

(Edited March 17th)

Here comes an improvement of the answer with this 'headline':

"Instead of using multiple sinuses with different amplitudes, you should use them with random phases:"

The problem with the way that was written is that it consumes enormous amounts of RAM. Here I have fixed it. It still uses a large amount of 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])
    peak_value = np.max(np.abs(signals))
    signals /= peak_value
    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)

EDITED March 17th:

I have now done a thorough analysis of the above generated "noise", compared to more 'ordinary' white noise, created with the following, much shorter and much faster code. The above has less variation of the 'amount' of signal at each frequency:

from scipy.io import wavfile
from scipy import stats
import numpy

sample_rate = 44100
length_in_seconds = 3
no_of_repetitions = 10
noise = numpy.array( stats.truncnorm(-1, 1, loc=0, scale=1).rvs(size = sample_rate * length_in_seconds * no_of_repetitions) )
if no_of_repetitions > 1:
    slices = []
    for i in range(no_of_repetitions):
        slices.append(
            numpy.array(noise[i * sample_rate * length_in_seconds:(i + 1) * sample_rate * length_in_seconds - 1]))
    added_up = numpy.zeros(len(slices[0]))
    for i in range(no_of_repetitions):
        added_up = added_up + slices[i]
else:
    added_up = noise
peak_value = numpy.max(numpy.abs(added_up))
numpy.abs(peak_value2))
added_up /= peak_value

wavfile.write('Quick noise float64 -1 +1 scale 1 rep 10.wav', sample_rate, added_up)

In order to make the comparison of the variation of the 'amount' of signal at each frequency, I had to do a 10 Hz high pass and a 20000 Hz low pass filtering of the "Quick noise" in Audacity, before comparing it's frequency spectrum (FFT) with that of the 'noise' from the more complex Python script that adds the sinuses of a large amount of frequencies. To generate the frequency spectrum of the 2 WAV-files (both 3 seconds) I used scipy.FFT.

I analysed the variation using more or less this:

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

Here comes the details of that:

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(distr=['norm', 't', 'dweibull', 'uniform'])  # Create the distfit object and choose interesting distributions.
    FFT_dist.fit_transform(yf)  # Call the fit transform function
    FFT_dist.plot() # dfit.plot_summary()  # Plot the summary

You can see the large variations of the frequency response using:

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

According to the mentioned FFT and distfit analysis, none of these 2 types of noise actually come near to what I want: An exact equal amount of ALL frequencies in a certain frequency band, for instance from 10 Hz to 20 KHz. The variation (=the standard deviation divided by the mean value, in the Normal distribution, computed by distfit) is = 0,538. That's not impressing. I would have liked 0,2. :-) Making the freq_step even smaller does not make it better. Actually worse.

Bebop

Answer 4

Producing full spectrum white noise and then filtering it is like you want to paint a wall of your house white, so you decide to paint the whole house white and then paint back all house except the wall. Is idiotic. (But has sense in electronics).

I made a small C program that can generate white noise at any frequency and any bandwidth (let's say at 16kHz central frequency and 2 kHz "wide"). No filtering involved.

What I did is simple: inside the main (infinite) loop I generate a sinusoid at center frequency +/- a random number between -half bandwidth and +halfbandwidth, then I keep that frequency for an arbitrary number of samples (granularity) and this is the result:

White noise 2kHz wide at 16kHz center frequency

Pseudo code:

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
    }


}

Inconveniency

This answer could have made sense, if it was written in Python... ;-) – Bebop