6

入力x（n）= u（n）-u（n-4）を持つ単純なデジタルシステムがあります。システム1

'signal'パッケージのconv（）関数または'stats'パッケージのconvolve（）関数を使用して出力y（n）を見つけ、-10≤n≤の場合のy（n）対nをプロットしようとしています。 10.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を使用した経験があり、いくつかの指針を与えることができれば素晴らしいと思います。前もって感謝します。

4

1 に答える 1

4

私はそれを理解しました..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))]
}

次に、多次元化に取り組む必要があると思います。

于 2013-02-02T22:28:41.437 に答える