R中数字信号处理的卷积

Posted 2023-04-12

【中文标题】R中数字信号处理的卷积

【英文标题】：Convolution for Digital Signal Processing in R

【发布时间】：2013-01-17 23:16:28

【问题描述】：

我有一个简单的数字系统，它有一个输入 x(n) = u(n) - u(n-4)。

我正在尝试使用“signal”包中的 conv() 函数或“stats”包中的 convolve() 函数查找输出 y(n)，并绘制 y(n) 与 n 的 -10 ≤ n ≤ 10。

到目前为止，我有以下代码：

library(signal)

n <- c(-10:10)                           # Time index
x <- c(rep(0, 10), rep(1, 4), rep(0, 7)) # Input Signal
h1 <- c(rep(0, 11), 0.5, rep(0, 9))      # Filter 1
h2 <- 0.8^n                              # Filter 2
h2[0:11] <- 0                            #

system <- data.frame(n, x, h1, h2)


y <- conv(x + conv(x, h1), h2)           # Output Signal

system <- transform(system, y=y[1:21]) 

plot(system$n, system$y)  

我检查了这个情节，这是非常错误的。我认为在进行卷积时会回收一些向量，并且 conv() 函数的输出似乎与原始时间索引不一致。我似乎无法弄清楚如何在这里解决我的逻辑。我意识到 conv(n, m) 函数返回一个长度为 (m+n)-1 的向量，有没有一种好方法可以轻松地将这个向量与时间索引向量匹配？

这将需要一些数字信号处理知识以及 R 中的编码知识，如果有人有为此目的使用 R 的经验并可以提供一些指导，那就太好了。提前致谢。

【参考方案1】：

我想通了。conv() 函数的输出中心与时间索引向量的中心对齐。因此：

library(signal)

n <- c(-10:10)                           # Time index
x <- c(rep(0, 10), rep(1, 4), rep(0, 7)) # Input Signal, square pulse
h1 <- c(rep(0, 11), 0.5, rep(0, 9))      # Filter 1
h2 <- 0.8^n                              # Filter 2
h2[1:10] <- 0                            #

system <- data.frame(n, x, h1, h2)

y <- conv(x + conv(x, h1)[11:31], h2)    # Output Signal

system <- transform(system, y=y[11:31]) 

plot(system$n, system$y)

我将使用通用表单来完成此操作，因为我会定期执行此操作，并且不想每次都手动执行此操作。如果有人打败我，请分享。 :)

更新

创建了一个通用形式的 conv() 函数来自动排列输入和输出向量的索引。这是以无法获得完整卷积为代价的，因此您必须首先将输入设置为描绘整个感兴趣区域。

library(signal) # Should this be inside the func. with attach(), detach()?

conv2 <- function(x, y)
    conv(x, y)[ceiling(length(x)/2):(length(x)+floor(length(x)/2))]


# so 
y <- conv2(x + conv2(x, h1), h2)

更新 2

我想要一个函数来与 FFT 进行比较。我对这个版本不太满意，我想使用 sapply()，但它确实有效。现在，它会做......我会努力改进。

conv3 <- function(x, h)
m <- length(x)
n <- length(h)
X <- c(x, rep(floor(n/2), 0, floor(n/2)))   
H <- c(h, rep(floor(m/2), 0, floor(m/2)))
   Y <- vector()

for(i in 1:n+m-1)
    Y[i] <- 0 
    for(j in 1:m)
        Y[i] <- ifelse(i-j+1>0, Y[i] + X[j]*H[i-j+1], 0)
    

Y[is.na(Y)] <- 0
Y[ceiling(m/2):(m+floor(m/2))]

接下来，我认为我需要努力使其具有多维性。