Introducere
Acest post descrie interpolarea în domeniul de frecvență (IFD) și o îmbunătățire a IFD. IFD este o metodă de interpolare ușor de utilizat și produce rezultate bune pe cele mai multe date din seria de timp (eșantioane de date la fel de distanțate în timp). Este ușor de utilizat, deoarece nu necesită nicio cunoaștere prealabilă a datelor. De asemenea, deoarece IFD folosește transformarea rapidă Fourier (FFT), este rapid (chiar și pe seturi de date mari). Dezavantajul către IFD este că uneori produce o interpolare care nu este bună la începutul și sfârșitul secvenței (denumită „efecte finale”). IFD îmbunătățit reduce efectele finale. Deși efectele finale nu sunt întotdeauna o problemă, este întotdeauna mai bine să nu le ai.
Prima secțiune introduce pe scurt interpolarea. Următoarea secțiune definește algoritmul IFD standard și arată exemple de utilizare a acestuia. Secțiunea finală arată un IFD îmbunătățit care reduce efectele finale.
Exemplu de interpolare
Să zicem că verificăm un ecartament de ploaie la 7:00 și apoi o dată pe oră timp de cinci ore, rezultând următoarele măsurători: 0,01 cm, 0,19 cm, 0,31 cm, 0,40 cm, 0,52 cm și 0,59 cm. Nu am măsurat ploaia la 7:30, dar am putea estima precipitațiile la 7:30 folosind un spline liniar pentru a interpola (conectarea punctelor cu linii drepte).
Figura de mai jos arată cele cinci valori de precipitații de sus și o ilustrare a interpolării unei valori de precipitații la 730. Figura arată cele cinci măsurători ca cercuri negre, linia albastră interpolantă care leagă punctele 700 și 800 și punctul interpolat la 0730 prezentat ca un cerc roșu. Cercul roșu arată că precipitațiile la 730 sunt de 0,1 cm.
Interpolare în domeniul frecvenței
Interpolarea liniară spline este eficientă dacă datele de bază pot fi modelate de o linie dreaptă. IFD funcționează bine indiferent de tipul de date.
Pașii din algoritm pentru IFD sunt enumerate mai jos.
- Găsiți transformarea discretă Fourier (DFT) a datelor folosind funcția încorporată
fftcare este implementarea R a FFT. Datele au început în domeniul timpului (date indexate în timp). După transformare, datele sunt în domeniul frecvenței (date indexate de frecvență). - Adăugați zero la date. Aceasta se numește zero Padding și forțează DFT invers (IDFT) în următorul pas pentru a evalua datele în mai multe puncte. Când adăugați zerouri la datele domeniului de frecvență, zerourile trebuie să meargă la mijloc. Acest pas este ușor diferit pentru secvențe de lungime uniformă și ciudată.
- Efectuați IDFT pe secvența DFT cu căptușeală zero.
Acești pași sunt explicați mai jos.
Găsiți DFT -ul secvenței
Pentru a ilustra acest pas, folosim valorile de cădere a ploii de sus.
În primul rând, folosim fft Funcție pentru a găsi DFT -ul secvenței. Mai jos, voi înceta să spun: „Folosiți fft
Funcție pentru a găsi DFT. ” Voi spune doar FFT.
N <- length(x) X <- fft(x) round(X,digits = 2) ## (1) 1.43+0.00i -0.34+0.37i -0.34+0.11i -0.34-0.11i -0.34-0.37i
Zero plasează datele
În acest caz, dorim să dublăm numărul de puncte, așa că adăugăm cinci zerouri la X. Zerourile trebuie adăugate în „mijlocul” secvenței. Dacă sunteți familiarizați cu teoria Transformatului Fourier, veți observa că aceasta este între frecvențele pozitive și negative.
Xzp <- c(X(1),X(2:3),rep(0,5),X(4:5)) round(Xzp,digits = 2) ## (1) 1.43+0.00i -0.34+0.37i -0.34+0.11i 0.00+0.00i 0.00+0.00i 0.00+0.00i ## (7) 0.00+0.00i 0.00+0.00i -0.34-0.11i -0.34-0.37i
Notă dacă eșantioanele DFT X ar avea un număr egal de eșantioane, de exemplu 6, atunci zero DFT XZP ar fi C (x (1), x (2: 3), x (4)/2, rep (0, n), x (4)/2, x (5: 6)). Eșantionul X (N/2+1) (în acest caz, X (4)) este utilizat de două ori, dar se taie la jumătate pentru a menține energia la fel. Arăt un exemplu în acest sens în secțiunea IFD.
DFT invers
Acum, găsim FFT invers. FFT invers este scalat cu unul peste numărul de probe. Observați că împărțim numărul de puncte de date reale (cinci) și nu numărul total de puncte din secvența captusită (zece). Dacă împărțim la zece, prejudiciem rezultatul scăzut, deoarece zerourile nu adaugă energie la secvență.
xinterp <- fft(Xzp,inverse = TRUE)/N round(xinterp,digits=2) ## (1) 0.01+0i 0.00+0i 0.19+0i 0.33+0i 0.31+0i 0.29+0i 0.40+0i 0.55+0i 0.52+0i ## (10) 0.26+0i
De când am început cu o secvență reală, ar trebui să ne încheiem cu o secvență reală. Cu toate acestea, fft ieșirea este complexă, așa că convertim secvența în reală folosind Re funcţie.
xinterp <- Re(xinterp) round(xinterp,digits = 2) ## (1) 0.01 0.00 0.19 0.33 0.31 0.29 0.40 0.55 0.52 0.26
Figura de mai jos conține două parcele. Secvența interpolată în Blue X’s și secvența originală în cercurile negre (cercurile sunt greu de văzut, deoarece au X pe ele). Rețineți că interpolarea (în special punctele de la 7:30 și 12:30) nu este grozavă. Această metodă suferă de efecte finale și funcționează slab cu secvențe scurte. Arătăm o metodă care atenuează efectele finale în secțiunea funcției IFD îmbunătățită.
Exemplu IFD cu o secvență mai lungă
Acum, să încercăm IFD cu o secvență mai lungă. Mai jos este o rampă de 49 de probe.
N <- 49 x <-1:N
Întrucât acum știm algoritmul, voi parcurge pașii cu mai puțin explicații. În primul rând, găsiți FFT.
X <-fft(x)
În al doilea rând, zero pad cu 2.
Xzp <- c(X(1:ceiling(N/2)),rep(0,N),X((ceiling(N/2)+1):N))
În cele din urmă, IDFT.
xInterp <- Re(fft(Xzp,inverse = TRUE)/N)
Figura de mai jos conține două parcele. Secvența interpolată este prezentată ca o linie albastră, iar secvența rampei este prezentată ca cercuri negre. Puteți vedea efectele finale la începutul și sfârșitul complotului. După efectele finale, punctele interpolate sunt apropiate de rampă.
Figura de mai jos conține o vedere zoom-în cea de mai sus. O linie albastră arată secvența interpolată, iar cercurile negre arată probele originale de rampă. Puteți vedea că după efectele de la început, interpolarea este foarte bună.
Funcția IFD
Aici, punem pașii IFD în funcție findIfdcare introduce o secvență de date dataSet și factorul de interpolare M. Pentru a menține funcția simplă, M este limitat la o valoare întreagă. Funcția verifică mai întâi pentru a vedea dacă factorul de interpolare este un număr întreg, nu mai puțin de 1. În continuare, găsește FFT de dataSet și seturi N la lungimea
dataSet. Funcția determină dacă dataSet
Are un număr uniform sau ciudat de eșantioane și apoi plasează zerourile la mijloc. După zero căptușeală, găsește IFFT-ul secvenței zero. În cele din urmă, returnează secvența interpolată, convertind secvența în reală dacă intrarea dataSet a fost real. Rețineți funcția R is.numeric Returnează adevărat dacă secvența este reală.
findIfd <- function(dataSet,M=2) {
if (round(M)!=M || M<1) {
stop("M must be an integer larger than zero.")
}
X <- fft(dataSet)
N <- length(X)
if(floor(N/2)==N/2){
Xzp <- c(X(1:(N/2)),X(N/2+1)/2,rep(0,(M-1)*N-1),X(N/2+1)/2,X((N/2+2):N))
} else {
Xzp <- c(X(1:ceiling(N/2)),rep(0,((M-1)*N)),X((ceiling(N/2)+1):N))
}
xint <- fft(Xzp,inverse = TRUE)/N
if (is.numeric(dataSet))
Re(xint(1:(M*N)))
else
xint(1:(M*N))
}
Exemplu: Funcția IFD și rampa
Mai jos, folosim findIfd pentru a interpola rampa. Puteți vedea că am început cu o rampă de 50 de probe și am folosit un factor de interpolare de 8, astfel încât lungimea secvenței interpolate este de 400.
N <- 50 x <-1:N xInterp <- findIfd(x,M=8)
Graficul de mai jos arată secvența interpolată în albastru (800 de puncte fac ca graficul să fie greu de citit), iar secvența rampei este prezentată în puncte negre. După cum s -a arătat mai sus, putem vedea că interpolarea este bună departe de mijloc.
Exemplu: interpolarea unei valuri sinusoidale
IFD funcționează mai bine atunci când funcția este periodică și lină. În acest exemplu, interpolăm o secvență sinusoidală și, deoarece este periodică și netedă, ar trebui să obținem un rezultat mai bun.
N <- 50 n <- 0:(N-1) f <- 0.02 x <- sin(2*pi*f*n) M <- 4 xInterp <- findIfd(x,M) Ndx <- seq(0,N-1/M,1/M)
Figura de mai jos conține două parcele. Cercurile negre sunt secvența sinusoidală originală, iar linia albastră arată interpolarea. Puteți vedea că aceasta este o interpolare foarte bună, fără efecte finale.
Funcția IFD s -a îmbunătățit
Versiunea îmbunătățită a funcției IFD „oglindește” secvența înainte de a efectua interpolarea. Prin oglindă, mă refer la apariția unei versiuni inversate a secvenței la sfârșitul secvenței originale. Vezi mai jos pentru un exemplu de oglindire.
x <- 1:5 c(x,rev(x)) ## (1) 1 2 3 4 5 5 4 3 2 1
Funcția findIfdImproved intră o secvență
x și factorul de interpolare M și apoi iese o secvență interpolată. Diferența dintre această funcție și findIfd Această funcție reflectă secvența înainte de a implementa algoritmul IFD. Probele suplimentare la capătul secvenței interpolate (cauzate de oglindire) sunt îndepărtate înainte de a returna secvența interpolată. Notă Oglindirea dublează secvența, făcând -o uniformă. Deci, trebuie să folosim doar metoda de căptușeală zero pentru secvențe de lungime uniformă. Rețineți că findIfdImproved
dublează secvența, dar nu durează de două ori mai mult din cauza proprietăților FFT.
findIfdImproved <- function(x,M) {
if (round(M)!=M || M<1 ) {
stop("M must be an integer larger than zero.")
}
X <- fft(c(x,rev(x)))
N <- length(X)
Xzp <- c(X(1:(N/2)),X(N/2+1)/2,rep(0,(M-1)*N-1),X(N/2+1)/2,X((N/2+2):N))
xint <- fft(Xzp,inverse = TRUE)/N
if (is.numeric(x)) {
xint <- Re(xint(1:(M*length(x)))) }
else{
xint <- xint(1:(M*length(x))) }
}
Mai jos, folosim findIfdImproved pentru a interpola rampa de 50 de puncte, cu un factor de interpolare de 8.
N <- 50 x <-1:N xInterp <- findIfdImproved(x,M=8)
Graficul de mai jos arată secvența interpolată ca linie albastră și secvența rampei ca puncte negre. Vedem findIfdImproved
produce o interpolare îmbunătățită, fără efecte finale aproape fără efecte (puteți vedea efecte finale ușoare la sfârșitul secvenței).
