如何在 R 中模拟粉红噪声

Posted 2023-02-14

【中文标题】如何在 R 中模拟粉红噪声

【英文标题】：How to simulate pink noise in R

【发布时间】：2012-01-31 15:00:00

【问题描述】：

我知道可以通过将rnorm() 的输出视为时间序列来实现白噪声。关于如何模拟粉红噪声有什么建议吗？

【参考方案1】：

包tuneR 具有noise 函数，可以生成白噪声或粉红噪声的波对象：

require(tuneR)
w <- noise(kind = c("white"))
p <- noise(kind = c("pink"))
par(mfrow=c(2,1))
plot(w,main="white noise")
plot(p,main="pink noise")

编辑：我意识到上面的方法不会生成向量 (doh)。将其转换为向量的残酷方法是添加以下代码：

writeWave(p,"p.wav")#writes pink noise on your hard drive
require(audio)#loads `audio` package to use `load.wave` function
p.vec <- load.wave("path/to/p.wav")#this will load pink noise as a vector

【讨论】：

p@left 还不够制作矢量吗？ （由于 CRAN 故障，我无法检查。）

只是出于兴趣，如何编写一个广义的“颜色”噪声函数，即抑制带宽的任意区域？对于一些 R-nerd 来说，这可能是一个令人愉快的新年项目 :-)

@Carl：您生成高斯白噪声，然后通过滤波器运行样本以生成所需的功率谱。粉红噪声被定义为具有“1/f”功率谱的噪声，因此您需要设计一个具有“1/sqrt(f)”频率响应的滤波器。通常，您会设计一个 FIR（有限脉冲响应）滤波器，以在某些感兴趣的频带中近似期望响应。

【参考方案2】：

正如@mbq 所说，您可以只使用 p@left 来获取向量，而不是保存和读取 wav 文件。另一方面，您可以直接使用 tuneR 中生成时间序列的函数：

TK95 <- function(N, alpha = 1) 
    f <- seq(from=0, to=pi, length.out=(N/2+1))[-c(1,(N/2+1))] # Fourier frequencies
    f_ <- 1 / f^alpha # Power law
    RW <- sqrt(0.5*f_) * rnorm(N/2-1) # for the real part
    IW <- sqrt(0.5*f_) * rnorm(N/2-1) # for the imaginary part
    fR <- complex(real = c(rnorm(1), RW, rnorm(1), RW[(N/2-1):1]), 
                  imaginary = c(0, IW, 0, -IW[(N/2-1):1]), length.out=N)
     # Those complex numbers that are to be back transformed for Fourier Frequencies 0, 2pi/N, 2*2pi/N, ..., pi, ..., 2pi-1/N 
     # Choose in a way that frequencies are complex-conjugated and symmetric around pi 
     # 0 and pi do not need an imaginary part
    reihe <- fft(fR, inverse=TRUE) # go back into time domain
    return(Re(reihe)) # imaginary part is 0

这很完美：

par(mfrow=c(3,1))
replicate(3,plot(TK95(1000,1),type="l",ylab="",xlab="time"))