-
Notifications
You must be signed in to change notification settings - Fork 2
/
README.Rmd
110 lines (93 loc) · 3.76 KB
/
README.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
---
output:
md_document:
variant: gfm
---
```{r Setup, include=FALSE}
library(rwavelet)
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
fig.path = "inst/doc/readme_img/"
)
```
# rwavelet
[![CRAN_Status_Badge](http://www.r-pkg.org/badges/version/rwavelet)](http://cran.r-project.org/package=rwavelet)
![CRAN Downloads](http://cranlogs.r-pkg.org/badges/rwavelet)
![CRAN Downloads](http://cranlogs.r-pkg.org/badges/grand-total/rwavelet)
Wavelet Analysis
## Download and Install
Install the devtools package if you haven't already.
```{r, eval=FALSE}
install.packages("devtools")
```
To install the development package, type the following at the R command line:
```{r, eval=FALSE}
devtools::install_github("fabnavarro/rwavelet")
library(rwavelet)
```
To install the CRAN version of the package, type the following:
```{r, eval=FALSE}
install.packages("rwavelet")
```
## General features
* 1d, 2d and 3d MRA Forward/Inverse Wavelet Transforms Periodized Orthogonal (FWT_PO, IFWT_PO)
* 1d, 2d tensor Forward/Inverse Wavelet Transform Periodized Orthogonal (FTWT_PO, ITWT_PO)
* 1d, 2d Translation Invariant Forward/Inverse Wavelet Transform (FWT_TI, IWT_TI, also know as redundant wavelet transform, maximal overlap wavelet transform, undecimated wavelet transform or stationary wavelet transform)
* Linear wavelet estimation/approximation (using 2 fold cross validation, CVLinear)
* Non linear wavelet estimation/approximation (hard and soft thresholding, can be easily extended to other thresholding rules)
* 1d Wavelet Block Thresholding estimation/approximation
* ...
To obtain the complete list of package functions, simply type
```{r, eval=FALSE}
help(package = "rwavelet")
```
## Getting Started
Here is an example of denoising of an experimental nuclear magnetic resonance (NMR) spectrum. We start by loading the data:
```{r}
data("RaphNMR")
Y <- RaphNMR
n <- length(Y)
t <- seq(0, 1, length = n)
```
Then we specify the coarse decomposition scale $j_0$, the wavelets we want to use (here, Symmlet with 6 vanishing moments) and we perform a fast wavelet transform to get the noisy wavelet coefficients (Ywd):
```{r}
j0 <- 0
J <- log2(n)
qmf <- MakeONFilter('Symmlet', 6)
Ywd <- FWT_PO(Y, j0, qmf)
Ywnoise <- Ywd
```
We estimate $\sigma$ the standard deviation of the noise using the maximum absolute deviation (with only the finest scale coefficients). We apply a hard thresholding rule (with a universal threshold) to the coefficient estimators and obtain the estimator by applying an inverse transform:
```{r}
hatsigma <- MAD(Ywd[(2^(J-1)+1):2^J])
lambda <- sqrt(2*log(n))*hatsigma
Ywd[(2^(j0)+1):n] <- HardThresh(Ywd[(2^(j0)+1):n], lambda)
fhat <- IWT_PO(Ywd, j0, qmf)
```
Finally, we plot the resulting estimator:
```{r, NMR, fig.width=8, fig.height=5}
par(mfrow=c(2,2), mgp = c(1.2, 0.5, 0), tcl = -0.2,
mar = .1 + c(2.5,2.5,1,1), oma = c(0,0,0,0))
plot(t,Y,xlab="", ylab="", main="Observations")
plot(t,Y,xlab="", ylab="", main="Observations and Estimator")
matlines(t, fhat, lwd=2, col="blue", lty=1)
plot(Ywnoise, ylim=c(-20, 20), xlab="", ylab="", main = "Noisy Coefficients")
matlines(rep(lambda, n), lwd=2,col="red",lty=1)
matlines(-rep(lambda, n), lwd=2,col="red",lty=1)
plot(Ywd, ylim=c(-20,20), xlab="", ylab="", main = "Estimated Coefficients")
```
See the [package vignette](http://fnavarro.perso.math.cnrs.fr/rpackage/rwaveletvignette.html) or [documentation](https://fnavarro.perso.math.cnrs.fr/rpackage/rwavelet.pdf) for more details.
You could also build and see the vignette associated with the package using the following lines of code
```{r, eval=FALSE}
devtools::install_github("fabnavarro/rwavelet", build_vignettes = TRUE)
library(rwavelet)
```
Then, to view the vignette
```{r, eval=FALSE}
vignette("rwaveletvignette")
```
## How to cite
```{r Cite me}
citation("rwavelet")
```