Как запустить фильтр высоких частот или низких частот на точках данных в R?

Я новичок в R, и я попытался найти информацию о следующем, ничего не найдя.

Зеленый график на картинке состоит из красных и желтых графов. Но позвольте сказать, что у меня есть только точки данных чего-то вроде зеленого графа. Как извлечь низкие/высокие частоты (т.е. Примерно красные/желтые графики) с помощью low pass/фильтр высоких частот?

Обновление: график был сгенерирован с помощью

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, но оригинальная октавная документация по крайней мере упомянутый radians, который заставил меня идти. Значения в моем исходном графе не были выбраны с учетом какой-либо конкретной частоты, поэтому я получил следующие не очень простые частоты: 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 соответственно).

Ответ 1

По запросу OP:

пакет сигналов содержит все виды фильтров для обработки сигналов. Большинство из них сопоставимо/совместимо с функциями обработки сигналов в Matlab/Octave.

Ответ 2

Недавно я столкнулся с подобной проблемой и не нашел здесь здесь особенно полезными ответы. Вот альтернативный подход.

Давайте начнем с определения данных примера из вопроса:

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)
y <- y + noise1

plot(x, y, type="l", ylim=range(-1.5*max_y,1.5*max_y,5), lwd = 5, col = "green")

Итак, зеленая линия - это набор данных, который мы хотим использовать для низкочастотного и высокочастотного фильтров.

Боковое примечание: строка в этом случае может быть выражена как функция с помощью кубического сплайна (spline(x,y, n = length(x))), но с данными реального мира это редко бывает так, поэтому предположим, что невозможно выразить набор данных как функция.

Самый простой способ сгладить такие данные, с которыми я столкнулся, - это использовать loess или smooth.spline с соответствующим span/spar. По мнению статистиков leess/smooth.spline, вероятно, здесь не правильный подход, поскольку на самом деле он не представляет определенную модель данных в этом смысле. Альтернативой является использование обобщенных аддитивных моделей (gam()функция из пакета mgcv). Мой аргумент в пользу использования лесса или сглаженного сплайна здесь заключается в том, что это проще и не имеет никакого значения, поскольку мы заинтересованы в видимом результирующем шаблоне. Наборы данных реального мира сложнее, чем в этом примере, и найти определенную функцию для фильтрации нескольких похожих наборов данных может быть затруднительным. Если видимая подгонка хороша, почему это усложняет значения R2 и p? Для меня приложение является визуальным, для которого лесс/сглаженные сплайны являются подходящими методами. Оба метода предполагают полиномиальные отношения с той разницей, что лесс более гибкий, также используя полиномы более высокой степени, тогда как кубический сплайн всегда кубический (x ^ 2). Какой из них использовать, зависит от тенденций в наборе данных. Тем не менее, следующим шагом является применение фильтра нижних частот в наборе данных с помощью loess() или smooth.spline():

lowpass.spline <- smooth.spline(x,y, spar = 0.6) ## Control spar for amount of smoothing
lowpass.loess <- loess(y ~ x, data = data.frame(x = x, y = y), span = 0.3) ## control span to define the amount of smoothing

lines(predict(lowpass.spline, x), col = "red", lwd = 2)
lines(predict(lowpass.loess, x), col = "blue", lwd = 2)

Красная линия - это сглаженный сплайн-фильтр и синий лессовый фильтр. Как видите, результаты немного отличаются. Я полагаю, что одним из аргументов использования GAM было бы найти наилучшее соответствие, если бы на самом деле тренды были такими ясными и последовательными среди наборов данных, но для этого приложения оба этих припадка достаточно хороши для меня.

После поиска подходящего фильтра нижних частот фильтрация верхних частот будет такой же простой, как вычесть значения фильтра низких частот от y:

highpass <- y - predict(lowpass.loess, x)
lines(x, highpass, lwd =  2)

Этот ответ приходит поздно, но я надеюсь, что это поможет кому-то, кто борется с подобной проблемой.

Ответ 3

Используйте функцию filterfilt вместо фильтра (пакетный сигнал), чтобы избавиться от сдвига сигнала.

library(signal)
bf <- butter(2, 1/50, type="low")
b1 <- filtfilt(bf, y+noise1)
points(x, b1, col="red", pch=20)

Ответ 4

Один метод использует fast fourier transform, реализованный в R как fft. Вот пример фильтра высоких частот. Из приведенных выше сюжетов идея, реализованная в этом примере, - получить серию желтого цвета, начиная с серии зеленым цветом (ваши реальные данные).

# I've changed the data a bit so it 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')

### 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')

Ответ 5

Проверьте эту ссылку, где есть R-код для фильтрации (медицинские сигналы). Это Мэтт Шотвелл, и сайт полон интересной информации о R/stats с медицинским изгибом:

biostattmat.com

Пакет fftfilt содержит множество алгоритмов фильтрации, которые также должны помочь.

Ответ 6

есть пакет на CRAN с именем FastICA, это вычисляет приближение независимых сигналов источника, однако для вычисления обоих сигналов вам нужна матрица не менее 2xn смешанных наблюдений (для этого примера), этот алгоритм может 'Определите два независимых сигнала с одним вектором 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 = "")

Ответ 7

Я не уверен, что какой-либо фильтр - лучший способ для вас. Более полезным инструментом для этой цели является быстрое преобразование Фурье.