r - tag - meta title länge 2018



Wie führe ich einen Hochpass- oder Tiefpassfilter für Datenpunkte in R aus? (6)

Überprüfen Sie diesen Link, wo es R-Code zum Filtern gibt (medizinische Signale). Es ist von Matt Shotwell und die Seite ist voller interessanter R / Statistik-Informationen mit einem medizinischen Hintergrund:

biostattmat.com

Das Paket fftfilt enthält viele Filteralgorithmen, die ebenfalls helfen sollten.

Ich bin ein Anfänger in R und habe versucht, Informationen über Folgendes zu finden, ohne etwas zu finden.

Das grüne Diagramm im Bild setzt sich aus den roten und gelben Diagrammen zusammen. Nehmen wir aber an, ich habe nur die Datenpunkte von so etwas wie dem grünen Graphen. Wie extrahiere ich die tiefen / hohen Frequenzen (dh ungefähr die roten / gelben Graphen) mit einem Tiefpass- / Hochpassfilter ?

Update: Die Grafik wurde mit erzeugt

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)

Update 2: Mit dem im signal enthaltenen Butterworth-Filter erhalte ich den folgenden Hinweis:

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)

Die Berechnungen waren ein bisschen arbeitsintensiv, signal.pdf gab kaum Hinweise darauf, welche Werte W sollte, aber die Original- Oktavdokumentation erwähnte zumindest das radians das mich zum Laufen brachte. Die Werte in meinem ursprünglichen Diagramm wurden nicht für eine bestimmte Frequenz ausgewählt, daher f_low = 1/500 * 2 = 1/250 Folgendes: f_low = 1/500 * 2 = 1/250 , f_high = 1/500 * 2*10 = 1/25 und die Abtastfrequenz f_s = 500/500 = 1 . Dann habe ich irgendwo zwischen den tiefen und hohen Frequenzen einen f_c für die Tief- / Hochpassfilter gewählt (1/100 bzw. 1/50).


Auf CRAN gibt es ein Paket namens FastICA , das die Approximation der unabhängigen Quellsignale berechnet. Um jedoch beide Signale zu berechnen, benötigen Sie eine Matrix von mindestens 2xn gemischten Beobachtungen (in diesem Beispiel). Dieser Algorithmus kann die beiden nicht bestimmen unabhängige Signale mit nur 1xn Vektor. Siehe das folgende Beispiel. hoffe das kann dir helfen.

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 = "")

Eine Methode ist die Verwendung der fast fourier transform die in R als fft . Hier ist ein Beispiel eines Hochpassfilters. Ausgehend von den obigen Darstellungen besteht die Idee in diesem Beispiel darin, die Serie in Gelb zu erhalten, beginnend mit der Serie in Grün (Ihre realen Daten).

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

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


Ich bin kürzlich auf ein ähnliches Problem gestoßen und fand die Antworten hier nicht besonders hilfreich. Hier ist ein alternativer Ansatz.

Beginnen wir mit der Definition der Beispieldaten aus der Frage:

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")

Die grüne Linie ist also der Datensatz, den wir für Tiefpass- und Hochpassfilter verwenden möchten.

Randnotiz: Die Linie könnte in diesem Fall als Funktion ausgedrückt werden, indem ein kubischer Spline ( spline(x,y, n = length(x)) ) verwendet wird. Bei Daten aus der realen Welt ist dies jedoch selten der Fall. Nehmen wir also an, dass dies eine Funktion ist Es ist nicht möglich, den Datensatz als Funktion auszudrücken.

Der einfachste Weg, solche Daten, auf die ich smooth.spline bin, zu glätten, ist die Verwendung von loess oder smooth.spline mit entsprechendem span / spar . Laut Statistikern ist loess / smooth.spline hier wahrscheinlich nicht der richtige Ansatz , da es in diesem Sinne kein wirklich definiertes Modell der Daten darstellt. Eine Alternative ist die Verwendung der Funktion Generalized Additive Models ( gam() aus dem Paket mgcv ). Mein Argument für die Verwendung von Löß oder geglättetem Spline ist, dass es einfacher ist und keinen Unterschied macht, da wir an dem sichtbaren resultierenden Muster interessiert sind. Reale Datasets sind komplizierter als in diesem Beispiel und es kann schwierig sein, eine definierte Funktion zum Filtern mehrerer ähnlicher Datasets zu finden. Wenn die sichtbare Anpassung gut ist, warum sollte sie dann mit R2- und p-Werten komplizierter werden? Für mich ist die Anwendung visuell, für die Löß / geglättete Splines geeignete Methoden sind. Beide Methoden gehen von Polynombeziehungen aus, mit dem Unterschied, dass Löß auch bei Polynomen höheren Grades flexibler ist, während der kubische Spline immer kubisch ist (x ^ 2). Welche verwendet werden soll, hängt von den Trends in einem Datensatz ab. Der nächste Schritt ist das smooth.spline() eines Tiefpassfilters auf den Datensatz mit smooth.spline() oder 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)

Die rote Linie ist der geglättete Spline-Filter und die blaue Linie der Lössfilter. Wie Sie sehen, unterscheiden sich die Ergebnisse geringfügig. Ich denke, ein Argument für die Verwendung von GAM wäre, die beste Anpassung zu finden, wenn die Trends in den Datensätzen wirklich so klar und konsistent wären, aber für diese Anwendung sind beide Anpassungen gut genug für mich.

Nach dem Finden eines passenden Tiefpassfilters ist die Hochpassfilterung so einfach wie das Subtrahieren der tiefpassgefilterten Werte von y :

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

Diese Antwort kommt zu spät, aber ich hoffe, sie hilft jemandem, der mit ähnlichen Problemen zu kämpfen hat.


Ich hatte auch Mühe herauszufinden, wie der W-Parameter in der Butter-Funktion auf den Filter-Cut-Off abgebildet wird, zum Teil, weil die Dokumentation für filter und filtfilt zum Zeitpunkt der Veröffentlichung nicht korrekt ist (dies deutet darauf hin, dass W = .1 eine 10 ergeben würde Hz-LP-Filter in Kombination mit Filtfilt bei einer Signalabtastrate von Fs = 100. Tatsächlich handelt es sich jedoch nur um ein 5-Hz-LP-Filter. Bei Verwendung von Filtfilt beträgt die Halbamplitudenabschaltung 5 Hz, die Halbleistungsabschaltung jedoch 5 Hz, wenn Sie den Filter nur einmal mit der Filterfunktion anwenden). Ich veröffentliche einen Demo-Code, der mir dabei geholfen hat, zu bestätigen, wie das alles funktioniert, und mit dem Sie überprüfen können, ob ein Filter das tut, was Sie wollen.

#Example usage of butter, filter, and filtfilt functions
#adapted from https://rdrr.io/cran/signal/man/filtfilt.html

library(signal)

Fs <- 100; #sampling rate

bf <- butter(3, 0.1);       
#when apply twice with filtfilt, 
#results in a 0 phase shift 
#5 Hz half-amplitude cut-off LP filter
#
#W * (Fs/2) == half-amplitude cut-off when combined with filtfilt
#
#when apply only one time, using the filter function (non-zero phase shift),
#W * (Fs/2) == half-power cut-off


t <- seq(0, .99, len = 100)   # 1 second sample

#generate a 5 Hz sine wave
x <- sin(2*pi*t*5)

#filter it with filtfilt
y <- filtfilt(bf, x)

#filter it with filter
z <- filter(bf, x)

#plot original and filtered signals
plot(t, x, type='l')
lines(t, y, col="red")
lines(t,z,col="blue")

#estimate signal attenuation (proportional reduction in signal amplitude)
1 - mean(abs(range(y[t > .2 & t < .8]))) #~50% attenuation at 5 Hz using filtfilt

1 - mean(abs(range(z[t > .2 & t < .8]))) #~30% attenuation at 5 Hz using filter

#demonstration that half-amplitude cut-off is 6 Hz when apply filter only once
x6hz <- sin(2*pi*t*6)

z6hz <- filter(bf, x6hz)

1 - mean(abs(range(z6hz[t > .2 & t < .8]))) #~50% attenuation at 6 Hz using filter


#plot the filter attenuation profile (for when apply one time, as with "filter" function):

hf <- freqz(bf, Fs = Fs);

plot(c(0, 20, 20, 0, 0), c(0, 0, 1, 1, 0), type = "l", 
 xlab = "Frequency (Hz)", ylab = "Attenuation (abs)")

lines(hf$f[hf$f<=20], abs(hf$h)[hf$f<=20])

plot(c(0, 20, 20, 0, 0), c(0, 0, -50, -50, 0),
 type = "l", xlab = "Frequency (Hz)", ylab = "Attenuation (dB)")

lines(hf$f[hf$f<=20], 20*log10(abs(hf$h))[hf$f<=20])

hf$f[which(abs(hf$h) - .5 < .001)[1]] #half-amplitude cutoff, around 6 Hz

hf$f[which(20*log10(abs(hf$h))+6 < .2)[1]] #half-amplitude cutoff, around 6 Hz

hf$f[which(20*log10(abs(hf$h))+3 < .2)[1]] #half-power cutoff, around 5 Hz

Verwenden Sie die Funktion filtfilt anstelle des Filters (Paketsignal), um die Signalverschiebung zu beseitigen.

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





frequency-analysis