Commit | Line | Data |
---|---|---|
1c6f223e BA |
1 | library("Rwave") |
2 | ||
d7d55bc1 | 3 | #Entrée : courbes synchrones, soit après étape 1 itérée, soit après chaqure étape 1 |
d03c0621 BA |
4 | step2 = function(conso) |
5 | { | |
db6fc17d BA |
6 | n <- nrow(conso) |
7 | delta <- ncol(conso) | |
8 | #TODO: automatic tune of all these parameters ? (for other users) | |
d03c0621 | 9 | nvoice <- 4 |
d7d55bc1 BA |
10 | # noctave = 2^13 = 8192 half hours ~ 180 days ; ~log2(ncol(conso)) |
11 | noctave = 13 | |
12 | # 4 here represent 2^5 = 32 half-hours ~ 1 day | |
db6fc17d BA |
13 | #NOTE: default scalevector == 2^(0:(noctave * nvoice) / nvoice) * s0 (?) |
14 | scalevector <- 2^(4:(noctave * nvoice) / nvoice) * 2 | |
15 | #condition: ( log2(s0*w0/(2*pi)) - 1 ) * nvoice + 1.5 >= 1 | |
16 | s0=2 | |
17 | w0=2*pi | |
18 | scaled=FALSE | |
19 | s0log = as.integer( (log2( s0*w0/(2*pi) ) - 1) * nvoice + 1.5 ) | |
20 | totnoct = noctave + as.integer(s0log/nvoice) + 1 | |
21 | ||
22 | # (normalized) observations node with CWT | |
23 | Xcwt4 <- lapply(seq_len(n), function(i) { | |
24 | ts <- scale(ts(conso[i,]), center=TRUE, scale=scaled) | |
25 | totts.cwt = Rwave::cwt(ts,totnoct,nvoice,w0,plot=0) | |
26 | ts.cwt = totts.cwt[,s0log:(s0log+noctave*nvoice)] | |
27 | #Normalization | |
28 | sqs <- sqrt(2^(0:(noctave*nvoice)/nvoice)*s0) | |
29 | sqres <- sweep(ts.cwt,MARGIN=2,sqs,'*') | |
30 | sqres / max(Mod(sqres)) | |
31 | }) | |
3ccd1e39 | 32 | |
db6fc17d BA |
33 | Xwer_dist <- matrix(0., n, n) |
34 | fcoefs = rep(1/3, 3) #moving average on 3 values (TODO: very slow! correct?!) | |
35 | for (i in 1:(n-1)) | |
1c6f223e | 36 | { |
db6fc17d | 37 | for (j in (i+1):n) |
d03c0621 | 38 | { |
db6fc17d BA |
39 | #TODO: later, compute CWT here (because not enough storage space for 32M series) |
40 | # 'circular=TRUE' is wrong, should just take values on the sides; to rewrite in C | |
41 | num <- filter(Mod(Xcwt4[[i]] * Conj(Xcwt4[[j]])), fcoefs, circular=TRUE) | |
42 | WX <- filter(Mod(Xcwt4[[i]] * Conj(Xcwt4[[i]])), fcoefs, circular=TRUE) | |
43 | WY <- filter(Mod(Xcwt4[[j]] * Conj(Xcwt4[[j]])), fcoefs, circular=TRUE) | |
44 | wer2 <- sum(colSums(num)^2) / sum( sum(colSums(WX) * colSums(WY)) ) | |
45 | Xwer_dist[i,j] <- sqrt(delta * ncol(Xcwt4[[1]]) * (1 - wer2)) | |
46 | Xwer_dist[j,i] <- Xwer_dist[i,j] | |
d03c0621 | 47 | } |
1c6f223e | 48 | } |
d03c0621 | 49 | diag(Xwer_dist) <- numeric(n) |
c6556868 | 50 | Xwer_dist |
1c6f223e | 51 | } |