您可以使用所需的任何噪声频谱生成噪声序列（包括$1/f$，也称为粉红噪声）通过在频谱空间中生成噪声系数。应选择系数的大小以给出所需的频谱，并且应随机选择相位。然后，您只需执行傅里叶逆变换即可给出序列值。以下算法将产生 1/f 噪声：

1. 确定点数，$n$, 和长度,$T$, 你的序列。这也确定了具有波数的光谱空间$-k_{\max}$到$k_{\max}$对应于频率值$f_k=kT/(2\pi)$.
2. 设置频谱系数的大小：$C_k = 1/\left|f_k\right|$. 放$C_0 = 0$给噪声序列零均值。
3. 将频谱系数的相位设置为随机值，如果需要实值噪声，请确保对称。
   为了$k=0..k_{\max}$
   $~~~\phi_k=\left(~\operatorname{rand\in[0,2\pi)}~\right)$,
   $~~~C_k = C_ke^{i\varphi_k}$,
   $~~~C_{-k} = C_{-k}e^{-i\varphi_k}$
4. 对频谱系数进行傅里叶逆变换，得到你的噪声序列，$\{y_i\} = \operatorname{ifft}(\{C_k\})$

如果您想要更高或更低的噪声功率，您可以简单地将结果序列乘以一个常数（或将所有$C_k$值由一个常数，这是相同的效果）。

根据定义，以这种方式生成的噪声具有完全指定的频谱，因为您是通过指定其频谱来创建噪声的。

您还可以生成噪声谱作为随机微分方程的解。这提供了一种逐个时间步生成噪声的迭代方法，而不是首先求解整个频谱（即，您只需再求解 SDE 的一个步骤，这是您的下一个噪声项）。本文针对您的噪声类型对一些 SDE 进行了很好的处理。

这是在 C# 中创建粉红噪声的一种非常简单的方法，它只是将许多波（以对数间隔）加在一起！这些函数输出一个值从 -1 到 1 的双精度数组。这分别表示波形中的最低点和最高点。

“质量”参数表示为发出声音而产生的波数。我发现 5000 波差不多是我无法检测到更高值的任何明显改进的阈值，但为了安全起见，您可以（可选）将其增加到大约 10,000 或更高。此外，根据 Wikipedia 的说法，就我们能听到的声音而言，20 赫兹大约是人类感知的下限，但如果你愿意，你也可以改变它。

请注意，由于技术原因，声音变得更安静，质量值更高，因此您可能（可选）希望通过 volumeAdjust 参数调整音量。

public double[] createPinkNoise(double seconds, int quality=5000, double lowestFrequency=20, double highestFrequency = 44100, double volumeAdjust=1.0)
{
    long samples = (long)(44100 * seconds);
    double[] d = new double[samples];
    double[] offsets = new double[samples];
    double lowestWavelength = highestFrequency / lowestFrequency;
    Random r = new Random();
    for (int j = 0; j < quality; j++)
    {
        double wavelength = Math.Pow(lowestWavelength, (j * 1.0) / quality)  * 44100 / highestFrequency;
        double offset = r.NextDouble() * Math.PI*2;     // Important offset is needed, as otherwise all the waves will be almost in phase, and this will ruin the effect!
        for (int i = 0; i < samples; i++)
        {
            d[i] += Math.Cos(i * Math.PI * 2 / wavelength + offset) / quality * volumeAdjust;
        }
    }
    return d;
}

我想我找到了一个很好的解决方案：

import numpy as np
import math
import cmath
from scipy.io import wavfile


def pink_noise(f_ref, f_min, f_max, length, f_sample):
    aliasfil_len = 10000
    fil_Time = aliasfil_len * 1/f_sample
    L = length + 2 * fil_Time
    f_low = 1 / L
    f_high = f_sample
    T = f_low * 2 * np.pi
    k_max = int(f_high / f_low / 2) + 1
    print(k_max)

    # Create frequencies
    f = np.array([(k * T)/(2 * np.pi) for k in range(0, k_max)])

    # Create 1/f noise amplitude in band
    C = np.array([(1 / f[k] if (f[k] >= f_min and f[k] <= f_max) else 0)
                  for k in range(0, k_max)], dtype='complex')
    C[0] = 0
    # Create random phase in band
    Phase = np.array([(np.random.uniform(0, np.pi)
                       if (f[k] >= f_min and f[k] <= f_max)
                       else 0)
                      for k in range(0, k_max)])

    Clist_neg = list()
    Clist_pos = list()
    for k in range(-k_max + 1, -1):
        Clist_neg.append(C[-k] * cmath.exp(-1j * Phase[-k]))
    for k in range(0, k_max):
        Clist_pos.append(C[k]  * cmath.exp( 1j * Phase[k] ))

    CC = np.array(Clist_pos + Clist_neg, dtype='complex')

    # Scale to max amplitude
    maxampl = max(abs(CC))
    CC /= maxampl

    tsig = np.fft.ifft(CC)
    sig = np.real(np.sign(tsig)) * np.abs(tsig)

    # Filter aliassing
    sig = sig[aliasfil_len:-aliasfil_len]

    # clip to maximum signal and
    # correct for amplitude at reference frequency
    if f_ref > ((f_max + f_min) / 2):
        print("WARNING: f_ref ({} Hz) should be smaller or equal to the mid "
              "between {} Hz and {} Hz "
              "to prevent clipping.\n"
              "f_ref changed to {} Hz"
              .format(f_ref,
                      f_min,
                      f_max,
                      ((f_max + f_min) / 2)))
        f_ref = ((f_max + f_min) / 2)
    maxampl = max(np.abs(sig))
    sig = sig / maxampl * f_ref / ((f_max + f_min) / 2)

    halfway = int(len(sig) / 2)
    # sign invert second part for a good connection,
    # it is the mirror of the first half
    sig2nd = -1 * sig[halfway:]
    sigc = np.concatenate((sig[0:halfway], sig2nd))
    # average middle point, but second point sign inverted
    sigc[halfway] = (sig[halfway-1] - sig[halfway+1])/2

    return(sigc)


length = 10.0  # seconds
f_sample = 22050  # Hz

f_ref = 1000  # Hz, The frequency for max amplitude

f_min = 1000  # Hz
f_max = 2000  # hz

sig = pink_noise(f_ref, f_min, f_max, length, f_sample)



p = 0
for x in sig:
    p += abs(x)
rms = math.sqrt(p/len(sig))
print("rms = {:f}, {:5.1f} dB".format(rms, 20 * math.log10(rms)))

print("Length of time signal: {} samples".format(len(sig)))
print("Time signal: ", sig)
x = sig * (2**15 - 1)

wavfile.write("test.wav", f_sample, np.array(x, dtype=np.int16))