43

私はRの初心者で、何も見つからずに次の情報を見つけようとしました。

写真の緑のグラフは、赤と黄色のグラフで構成されています。しかし、緑のグラフのようなデータ ポイントしか持っていないとしましょう。ローパス/ハイパスフィルターを使用して低/高周波数 (つまり、ほぼ赤/黄色のグラフ) を抽出するにはどうすればよいですか?

変調された高周波正弦曲線を持つ低周波正弦曲線

更新：グラフはで生成されました

number_of_cycles = 2
max_y = 40

x = 1:500
a = number_of_cycles * 2*pi/length(x)

y = max_y * sin(x*a)
noise1 = max_y * 1/10 * sin(x*a*10)

plot(x, y, type="l", col="red", ylim=range(-1.5*max_y,1.5*max_y,5))
points(x, y + noise1, col="green", pch=20)
points(x, noise1, col="yellow", pch=20)

更新 2: パッケージでバターワース フィルターを使用するsignalと、次の結果が得られることが示唆されました。

フィルタリングされたグラフが追加された元の画像

library(signal)

bf <- butter(2, 1/50, type="low")
b <- filter(bf, y+noise1)
points(x, b, col="black", pch=20)

bf <- butter(2, 1/25, type="high")
b <- filter(bf, y+noise1)
points(x, b, col="black", pch=20)

計算はちょっとした作業でした。signal.pdf には、どのような値Wが必要かについてのヒントはほとんどありませんでしたが、元のオクターブのドキュメントには、少なくともラジアンが記載されていたため、うまくいきました。元のグラフの値は、特定の周波数を念頭に置いて選択されていないため、単純ではない次の周波数になりました:f_low = 1/500 * 2 = 1/250とf_high = 1/500 * 2*10 = 1/25サンプリング周波数f_s = 500/500 = 1. 次に、低域/高域フィルターの低域と高域の間の f_c を選択しました (それぞれ 1/100 と 1/50)。

4

8 に答える 8

8

1 つの方法はfast fourier transform、R に実装されている asを使用することfftです。ハイパスフィルタの例を次に示します。上記のプロットから、この例で実装されているアイデアは、緑色のセリエ (実際のデータ) から黄色のセリエを取得することです。

# I've changed the data a bit so it's easier to see in the plots
par(mfrow = c(1, 1))
number_of_cycles = 2
max_y = 40
N <- 256

x = 0:(N-1)
a = number_of_cycles * 2 * pi/length(x)

y = max_y * sin(x*a)
noise1 = max_y * 1/10 * sin(x*a*10)
plot(x, y, type="l", col="red", ylim=range(-1.5*max_y,1.5*max_y,5))
points(x, y + noise1, col="green", pch=20)
points(x, noise1, col="yellow", pch=20)

### Apply the fft to the noisy data
y_noise = y + noise1
fft.y_noise = fft(y_noise)


# Plot the series and spectrum
par(mfrow = c(1, 2))
plot(x, y_noise, type='l', main='original serie', col='green4')
plot(Mod(fft.y_noise), type='l', main='Raw serie - fft spectrum')

y ノイズと fft スペクトル

### The following code removes the first spike in the spectrum
### This would be the high pass filter
inx_filter = 15
FDfilter = rep(1, N)
FDfilter[1:inx_filter] = 0
FDfilter[(N-inx_filter):N] = 0
fft.y_noise_filtered = FDfilter * fft.y_noise

ここに画像の説明を入力

par(mfrow = c(2, 1))
plot(x, noise1, type='l', main='original noise')
plot(x, y=Re( fft( fft.y_noise_filtered, inverse=TRUE) / N ) , type='l', 
     main = 'filtered noise')

ここに画像の説明を入力

于 2016-06-24T04:06:30.027 に答える
7

OPのリクエストごと：

信号パッケージには、信号処理用のあらゆる種類のフィルターが含まれています。そのほとんどは、Matlab/Octaveの信号処理機能に匹敵する/互換性があります。

于 2011-08-19T08:26:11.387 に答える
3

フィルタリング用のRコード（医療信号）があるこのリンクをチェックしてください。これは Matt Shotwell によるもので、このサイトには興味深い R/stats 情報が満載で、医学的な傾向があります。

biostattmat.com

パッケージ fftfilt には、役立つフィルタリング アルゴリズムが多数含まれています。

于 2011-08-18T10:35:33.270 に答える
2

CRAN には という名前のパッケージがありFastICA、これは独立したソース信号の近似を計算しますが、両方の信号を計算するには、少なくとも 2xn の混合観測の行列が必要です (この例では)。このアルゴリズムは 2 つの独立したものを決定できません。 1xn ベクトルだけの信号。以下の例を参照してください。これがあなたを助けることを願っています。 

number_of_cycles = 2
max_y = 40

x = 1:500
a = number_of_cycles * 2*pi/length(x)

y = max_y * sin(x*a)
noise1 = max_y * 1/10 * sin(x*a*10)

plot(x, y, type="l", col="red", ylim=range(-1.5*max_y,1.5*max_y,5))
points(x, y + noise1, col="green", pch=20)
points(x, noise1, col="yellow", pch=20)
######################################################
library(fastICA)
S <- cbind(y,noise1)#Assuming that "y" source1 and "noise1" is source2
A <- matrix(c(0.291, 0.6557, -0.5439, 0.5572), 2, 2) #This is a mixing matrix
X <- S %*% A 

a <- fastICA(X, 2, alg.typ = "parallel", fun = "logcosh", alpha = 1,
method = "R", row.norm = FALSE, maxit = 200,
tol = 0.0001, verbose = TRUE)

par(mfcol = c(2, 3))
plot(S[,1 ], type = "l", main = "Original Signals",
xlab = "", ylab = "")
plot(S[,2 ], type = "l", xlab = "", ylab = "")
plot(X[,1 ], type = "l", main = "Mixed Signals",
xlab = "", ylab = "")
plot(X[,2 ], type = "l", xlab = "", ylab = "")
plot(a$S[,1 ], type = "l", main = "ICA source estimates",
xlab = "", ylab = "")
plot(a$S[, 2], type = "l", xlab = "", ylab = "")
于 2014-11-08T21:55:25.473 に答える
1

どのフィルターがあなたにとって最良の方法であるかはわかりません. その目的のためのより便利な手段は、高速フーリエ変換です。

于 2013-03-01T12:43:30.063 に答える