使用R语言计算功率谱密度函数R(w)

以下是使用R语言编写的程序，用于计算 R(w) 的值：

calculate_R <- function(y, w, n) {
  # 执行FFT变换
  fft_result <- fft(y)

  # 计算 ∑y[t]e^(-itw)
  sum_y <- sum(fft_result * exp(-1i * seq_along(y) * w))

  # 计算模的平方
  modulus_squared <- Mod(sum_y) ^ 2

  # 计算 R(w)
  R <- (1 / (2 * pi * n)) * modulus_squared

  return(R)
}

# 示例输入
y <- c(1+2i, 2-1i, 3+0i)  # y[t]序列
n <- length(y)  # y[t]序列的长度

# 计算R(w)
w_values <- seq(0, pi, length.out = 100)  # 在[0,π]上均匀取100个w值
R_values <- numeric(length(w_values))  # 存储对应的R(w)值

for (i in seq_along(w_values)) {
  R_values[i] <- calculate_R(y, w_values[i], n)
}

# 打印R(w)值
print(R_values)

此程序首先定义了一个计算 R(w) 的函数 calculate_R，函数参数包括 y[t] 序列、w 值和 y[t] 序列的长度 n。函数内部使用FFT变换进行频域计算，计算 ∑y[t]e^(-itw) 的值，并得到模的平方。然后将模的平方除以 2πn，得到 R(w) 的值。接下来，在给定的范围 [0,π] 上均匀取100个 w 值，通过循环遍历计算每个 w 值对应的 R(w) 值，并将结果存储在一个向量中。最后打印出 R(w) 值。您可以根据需要修改和扩展输入的 y[t] 序列。