Skip to content

Instantly share code, notes, and snippets.

@mike-lawrence

mike-lawrence/CEEMD_example.R

Created Dec 29, 2011
Star
Embed
What would you like to do?
Learn more about clone URLs
Download ZIP
CEEMD demo
Raw
CEEMD_example.R
#load doMC for parallelism
library(doMC)
options(cores=4) #adjust to match your system
registerDoMC()
#load plyr to use doMC implicitly
library(plyr)
#load the toy data set
source('data.R')
str(data)
#show a bit of the data
plot(data,type='l')
#set some parameters
num_iterations = 1e2 #more iterations = more accurate IMFs but longer compute time
max_imfs = 1e1 #this value should increase with larger amounts of data
#pre-compute some values to be used below
signal_length = nrow(data)
noise_sd = sd(data$value)*.1
#set the random seed for replicability
set.seed(1)
#start a timer
start = proc.time()[3]
#run the CEEMD (this is the section to be converted to CUDA)
source('EMD.R') #the code in EMD.R will need to be converted to CUDA
from_ceemd = llply(
.parallel = TRUE
, .progress = 'text' #the progress bar isn't accurate when .parallel=TRUE; it'll be done at 100/ncores %
, .data = 1:num_iterations
, .fun = function(iteration){
noise = rnorm(signal_length,0,noise_sd)
imfs1 = emd(
xt = data$value+noise
, tt = data$time
, boundary = 'symmetric'
, max.imf = max_imfs
, plot.imf = F
, print_progress = F
)$imf
imfs2 = emd(
xt = data$value-noise
, tt = data$time
, boundary = 'symmetric'
, max.imf = max_imfs
, plot.imf = F
, print_progress = F
)$imf
min_cols = min(c(ncol(imfs1),ncol(imfs2)))
return((imfs1[,1:min_cols]+imfs2[,1:min_cols])/2) #return the average of the two IMF matricies
}
)
#print the time taken
proc.time()[3]-start
#collapse the results to a mean
imfs = matrix(0,nrow = signal_length, ncol = max_imfs)
for(i in 1:length(from_ceemd)){
imfs[,1:ncol(from_ceemd[[i]])] = (imfs[,1:ncol(from_ceemd[[i]])]+from_ceemd[[i]])/(i*2-1)
}
#drop any all-zero IMFs
imfs = imfs[,colSums(imfs)!=0]
#show a "filtered" version of the data by summing all but the few fastest-varying and slowest-varying IMFs (like a bandpass)
plot(x=data$time,y=rowSums(imfs[,3:(ncol(imfs)-2)]),type='l')
Raw
data.R
data = structure(list(time = c(0, 0.0039062, 0.0078125, 0.0117188, 0.015625,
0.0195312, 0.0234375, 0.0273438, 0.03125, 0.0351562, 0.0390625,
0.0429688, 0.046875, 0.0507812, 0.0546875, 0.0585938, 0.0625,
0.0664062, 0.0703125, 0.0742188, 0.078125, 0.0820312, 0.0859375,
0.0898438, 0.09375, 0.0976562, 0.1015625, 0.1054688, 0.109375,
0.1132812, 0.1171875, 0.1210938, 0.125, 0.1289062, 0.1328125,
0.1367188, 0.140625, 0.1445312, 0.1484375, 0.1523438, 0.15625,
0.1601562, 0.1640625, 0.1679688, 0.171875, 0.1757812, 0.1796875,
0.1835938, 0.1875, 0.1914062, 0.1953125, 0.1992188, 0.203125,
0.2070312, 0.2109375, 0.2148438, 0.21875, 0.2226562, 0.2265625,
0.2304688, 0.234375, 0.2382812, 0.2421875, 0.2460938, 0.25, 0.2539062,
0.2578125, 0.2617188, 0.265625, 0.2695312, 0.2734375, 0.2773438,
0.28125, 0.2851562, 0.2890625, 0.2929688, 0.296875, 0.3007812,
0.3046875, 0.3085938, 0.3125, 0.3164062, 0.3203125, 0.3242188,
0.328125, 0.3320312, 0.3359375, 0.3398438, 0.34375, 0.3476562,
0.3515625, 0.3554688, 0.359375, 0.3632812, 0.3671875, 0.3710938,
0.375, 0.3789062, 0.3828125, 0.3867188, 0.390625, 0.3945312,
0.3984375, 0.4023438, 0.40625, 0.4101562, 0.4140625, 0.4179688,
0.421875, 0.4257812, 0.4296875, 0.4335938, 0.4375, 0.4414062,
0.4453125, 0.4492188, 0.453125, 0.4570312, 0.4609375, 0.4648438,
0.46875, 0.4726562, 0.4765625, 0.4804688, 0.484375, 0.4882812,
0.4921875, 0.4960938, 0.5, 0.5039062, 0.5078125, 0.5117188, 0.515625,
0.5195312, 0.5234375, 0.5273438, 0.53125, 0.5351562, 0.5390625,
0.5429688, 0.546875, 0.5507812, 0.5546875, 0.5585938, 0.5625,
0.5664062, 0.5703125, 0.5742188, 0.578125, 0.5820312, 0.5859375,
0.5898438, 0.59375, 0.5976562, 0.6015625, 0.6054688, 0.609375,
0.6132812, 0.6171875, 0.6210938, 0.625, 0.6289062, 0.6328125,
0.6367188, 0.640625, 0.6445312, 0.6484375, 0.6523438, 0.65625,
0.6601562, 0.6640625, 0.6679688, 0.671875, 0.6757812, 0.6796875,
0.6835938, 0.6875, 0.6914062, 0.6953125, 0.6992188, 0.703125,
0.7070312, 0.7109375, 0.7148438, 0.71875, 0.7226562, 0.7265625,
0.7304688, 0.734375, 0.7382812, 0.7421875, 0.7460938, 0.75, 0.7539062,
0.7578125, 0.7617188, 0.765625, 0.7695312, 0.7734375, 0.7773438,
0.78125, 0.7851562, 0.7890625, 0.7929688, 0.796875, 0.8007812,
0.8046875, 0.8085938, 0.8125, 0.8164062, 0.8203125, 0.8242188,
0.828125, 0.8320312, 0.8359375, 0.8398438, 0.84375, 0.8476562,
0.8515625, 0.8554688, 0.859375, 0.8632812, 0.8671875, 0.8710938,
0.875, 0.8789062, 0.8828125, 0.8867188, 0.890625, 0.8945312,
0.8984375, 0.9023438, 0.90625, 0.9101562, 0.9140625, 0.9179688,
0.921875, 0.9257812, 0.9296875, 0.9335938, 0.9375, 0.9414062,
0.9453125, 0.9492188, 0.953125, 0.9570312, 0.9609375, 0.9648438,
0.96875, 0.9726562, 0.9765625, 0.9804688, 0.984375, 0.9882812,
0.9921875, 0.9960938, 1, 1.0039062, 1.0078125, 1.0117188, 1.015625,
1.0195312, 1.0234375, 1.0273438, 1.03125, 1.0351562, 1.0390625,
1.0429688, 1.046875, 1.0507812, 1.0546875, 1.0585938, 1.0625,
1.0664062, 1.0703125, 1.0742188, 1.078125, 1.0820312, 1.0859375,
1.0898438, 1.09375, 1.0976562, 1.1015625, 1.1054688, 1.109375,
1.1132812, 1.1171875, 1.1210938, 1.125, 1.1289062, 1.1328125,
1.1367188, 1.140625, 1.1445312, 1.1484375, 1.1523438, 1.15625,
1.1601562, 1.1640625, 1.1679688, 1.171875, 1.1757812, 1.1796875,
1.1835938, 1.1875, 1.1914062, 1.1953125, 1.1992188, 1.203125,
1.2070312, 1.2109375, 1.2148438, 1.21875, 1.2226562, 1.2265625,
1.2304688, 1.234375, 1.2382812, 1.2421875, 1.2460938, 1.25, 1.2539062,
1.2578125, 1.2617188, 1.265625, 1.2695312, 1.2734375, 1.2773438,
1.28125, 1.2851562, 1.2890625, 1.2929688, 1.296875, 1.3007812,
1.3046875, 1.3085938, 1.3125, 1.3164062, 1.3203125, 1.3242188,
1.328125, 1.3320312, 1.3359375, 1.3398438, 1.34375, 1.3476562,
1.3515625, 1.3554688, 1.359375, 1.3632812, 1.3671875, 1.3710938,
1.375, 1.3789062, 1.3828125, 1.3867188, 1.390625, 1.3945312,
1.3984375, 1.4023438, 1.40625, 1.4101562, 1.4140625, 1.4179688,
1.421875, 1.4257812, 1.4296875, 1.4335938, 1.4375, 1.4414062,
1.4453125, 1.4492188, 1.453125, 1.4570312, 1.4609375, 1.4648438,
1.46875, 1.4726562, 1.4765625, 1.4804688, 1.484375, 1.4882812,
1.4921875, 1.4960938, 1.5, 1.5039062, 1.5078125, 1.5117188, 1.515625,
1.5195312, 1.5234375, 1.5273438, 1.53125, 1.5351562, 1.5390625,
1.5429688, 1.546875, 1.5507812, 1.5546875, 1.5585938, 1.5625,
1.5664062, 1.5703125, 1.5742188, 1.578125, 1.5820312, 1.5859375,
1.5898438, 1.59375, 1.5976562, 1.6015625, 1.6054688, 1.609375,
1.6132812, 1.6171875, 1.6210938, 1.625, 1.6289062, 1.6328125,
1.6367188, 1.640625, 1.6445312, 1.6484375, 1.6523438, 1.65625,
1.6601562, 1.6640625, 1.6679688, 1.671875, 1.6757812, 1.6796875,
1.6835938, 1.6875, 1.6914062, 1.6953125, 1.6992188, 1.703125,
1.7070312, 1.7109375, 1.7148438, 1.71875, 1.7226562, 1.7265625,
1.7304688, 1.734375, 1.7382812, 1.7421875, 1.7460938, 1.75, 1.7539062,
1.7578125, 1.7617188, 1.765625, 1.7695312, 1.7734375, 1.7773438,
1.78125, 1.7851562, 1.7890625, 1.7929688, 1.796875, 1.8007812,
1.8046875, 1.8085938, 1.8125, 1.8164062, 1.8203125, 1.8242188,
1.828125, 1.8320312, 1.8359375, 1.8398438, 1.84375, 1.8476562,
1.8515625, 1.8554688, 1.859375, 1.8632812, 1.8671875, 1.8710938,
1.875, 1.8789062, 1.8828125, 1.8867188, 1.890625, 1.8945312,
1.8984375, 1.9023438, 1.90625, 1.9101562, 1.9140625, 1.9179688,
1.921875, 1.9257812, 1.9296875, 1.9335938, 1.9375, 1.9414062,
1.9453125, 1.9492188, 1.953125, 1.9570312, 1.9609375, 1.9648438,
1.96875, 1.9726562, 1.9765625, 1.9804688, 1.984375, 1.9882812,
1.9921875, 1.9960938, 2, 2.0039062, 2.0078125, 2.0117188, 2.015625,
2.0195312, 2.0234375, 2.0273438, 2.03125, 2.0351562, 2.0390625,
2.0429688, 2.046875, 2.0507812, 2.0546875, 2.0585938, 2.0625,
2.0664062, 2.0703125, 2.0742188, 2.078125, 2.0820312, 2.0859375,
2.0898438, 2.09375, 2.0976562, 2.1015625, 2.1054688, 2.109375,
2.1132812, 2.1171875, 2.1210938, 2.125, 2.1289062, 2.1328125,
2.1367188, 2.140625, 2.1445312, 2.1484375, 2.1523438, 2.15625,
2.1601562, 2.1640625, 2.1679688, 2.171875, 2.1757812, 2.1796875,
2.1835938, 2.1875, 2.1914062, 2.1953125, 2.1992188, 2.203125,
2.2070312, 2.2109375, 2.2148438, 2.21875, 2.2226562, 2.2265625,
2.2304688, 2.234375, 2.2382812, 2.2421875, 2.2460938, 2.25, 2.2539062,
2.2578125, 2.2617188, 2.265625, 2.2695312, 2.2734375, 2.2773438,
2.28125, 2.2851562, 2.2890625, 2.2929688, 2.296875, 2.3007812,
2.3046875, 2.3085938, 2.3125, 2.3164062, 2.3203125, 2.3242188,
2.328125, 2.3320312, 2.3359375, 2.3398438, 2.34375, 2.3476562,
2.3515625, 2.3554688, 2.359375, 2.3632812, 2.3671875, 2.3710938,
2.375, 2.3789062, 2.3828125, 2.3867188, 2.390625, 2.3945312,
2.3984375, 2.4023438, 2.40625, 2.4101562, 2.4140625, 2.4179688,
2.421875, 2.4257812, 2.4296875, 2.4335938, 2.4375, 2.4414062,
2.4453125, 2.4492188, 2.453125, 2.4570312, 2.4609375, 2.4648438,
2.46875, 2.4726562, 2.4765625, 2.4804688, 2.484375, 2.4882812,
2.4921875, 2.4960938, 2.5, 2.5039062, 2.5078125, 2.5117188, 2.515625,
2.5195312, 2.5234375, 2.5273438, 2.53125, 2.5351562, 2.5390625,
2.5429688, 2.546875, 2.5507812, 2.5546875, 2.5585938, 2.5625,
2.5664062, 2.5703125, 2.5742188, 2.578125, 2.5820312, 2.5859375,
2.5898438, 2.59375, 2.5976562, 2.6015625, 2.6054688, 2.609375,
2.6132812, 2.6171875, 2.6210938, 2.625, 2.6289062, 2.6328125,
2.6367188, 2.640625, 2.6445312, 2.6484375, 2.6523438, 2.65625,
2.6601562, 2.6640625, 2.6679688, 2.671875, 2.6757812, 2.6796875,
2.6835938, 2.6875, 2.6914062, 2.6953125, 2.6992188, 2.703125,
2.7070312, 2.7109375, 2.7148438, 2.71875, 2.7226562, 2.7265625,
2.7304688, 2.734375, 2.7382812, 2.7421875, 2.7460938, 2.75, 2.7539062,
2.7578125, 2.7617188, 2.765625, 2.7695312, 2.7734375, 2.7773438,
2.78125, 2.7851562, 2.7890625, 2.7929688, 2.796875, 2.8007812,
2.8046875, 2.8085938, 2.8125, 2.8164062, 2.8203125, 2.8242188,
2.828125, 2.8320312, 2.8359375, 2.8398438, 2.84375, 2.8476562,
2.8515625, 2.8554688, 2.859375, 2.8632812, 2.8671875, 2.8710938,
2.875, 2.8789062, 2.8828125, 2.8867188, 2.890625, 2.8945312,
2.8984375, 2.9023438, 2.90625, 2.9101562, 2.9140625, 2.9179688,
2.921875, 2.9257812, 2.9296875, 2.9335938, 2.9375, 2.9414062,
2.9453125, 2.9492188, 2.953125, 2.9570312, 2.9609375, 2.9648438,
2.96875, 2.9726562, 2.9765625, 2.9804688, 2.984375, 2.9882812,
2.9921875, 2.9960938, 3, 3.0039062, 3.0078125, 3.0117188, 3.015625,
3.0195312, 3.0234375, 3.0273438, 3.03125, 3.0351562, 3.0390625,
3.0429688, 3.046875, 3.0507812, 3.0546875, 3.0585938, 3.0625,
3.0664062, 3.0703125, 3.0742188, 3.078125, 3.0820312, 3.0859375,
3.0898438, 3.09375, 3.0976562, 3.1015625, 3.1054688, 3.109375,
3.1132812, 3.1171875, 3.1210938, 3.125, 3.1289062, 3.1328125,
3.1367188, 3.140625, 3.1445312, 3.1484375, 3.1523438, 3.15625,
3.1601562, 3.1640625, 3.1679688, 3.171875, 3.1757812, 3.1796875,
3.1835938, 3.1875, 3.1914062, 3.1953125, 3.1992188, 3.203125,
3.2070312, 3.2109375, 3.2148438, 3.21875, 3.2226562, 3.2265625,
3.2304688, 3.234375, 3.2382812, 3.2421875, 3.2460938, 3.25, 3.2539062,
3.2578125, 3.2617188, 3.265625, 3.2695312, 3.2734375, 3.2773438,
3.28125, 3.2851562, 3.2890625, 3.2929688, 3.296875, 3.3007812,
3.3046875, 3.3085938, 3.3125, 3.3164062, 3.3203125, 3.3242188,
3.328125, 3.3320312, 3.3359375, 3.3398438, 3.34375, 3.3476562,
3.3515625, 3.3554688, 3.359375, 3.3632812, 3.3671875, 3.3710938,
3.375, 3.3789062, 3.3828125, 3.3867188, 3.390625, 3.3945312,
3.3984375, 3.4023438, 3.40625, 3.4101562, 3.4140625, 3.4179688,
3.421875, 3.4257812, 3.4296875, 3.4335938, 3.4375, 3.4414062,
3.4453125, 3.4492188, 3.453125, 3.4570312, 3.4609375, 3.4648438,
3.46875, 3.4726562, 3.4765625, 3.4804688, 3.484375, 3.4882812,
3.4921875, 3.4960938, 3.5, 3.5039062, 3.5078125, 3.5117188, 3.515625,
3.5195312, 3.5234375, 3.5273438, 3.53125, 3.5351562, 3.5390625,
3.5429688, 3.546875, 3.5507812, 3.5546875, 3.5585938, 3.5625,
3.5664062, 3.5703125, 3.5742188, 3.578125, 3.5820312, 3.5859375,
3.5898438, 3.59375, 3.5976562, 3.6015625, 3.6054688, 3.609375,
3.6132812, 3.6171875, 3.6210938, 3.625, 3.6289062, 3.6328125,
3.6367188, 3.640625, 3.6445312, 3.6484375, 3.6523438, 3.65625,
3.6601562, 3.6640625, 3.6679688, 3.671875, 3.6757812, 3.6796875,
3.6835938, 3.6875, 3.6914062, 3.6953125, 3.6992188, 3.703125,
3.7070312, 3.7109375, 3.7148438, 3.71875, 3.7226562, 3.7265625,
3.7304688, 3.734375, 3.7382812, 3.7421875, 3.7460938, 3.75, 3.7539062,
3.7578125, 3.7617188, 3.765625, 3.7695312, 3.7734375, 3.7773438,
3.78125, 3.7851562, 3.7890625, 3.7929688, 3.796875, 3.8007812,
3.8046875, 3.8085938, 3.8125, 3.8164062, 3.8203125, 3.8242188,
3.828125, 3.8320312, 3.8359375, 3.8398438, 3.84375, 3.8476562,
3.8515625, 3.8554688, 3.859375, 3.8632812, 3.8671875, 3.8710938,
3.875, 3.8789062, 3.8828125, 3.8867188, 3.890625, 3.8945312,
3.8984375, 3.9023438), value = c(-335.332553252417, -325.579494919085,
-328.529622696863, -335.545636585752, -331.999206030195, -323.181504641307,
-323.011393530195, -333.434951863529, -334.045643530195, -327.790883807973,
-323.349007419085, -333.710994919084, -337.082990752418, -330.941960196863,
-321.579497696862, -327.126799085751, -331.654594919085, -330.339489363529,
-323.455811585751, -325.512722696862, -330.559968530196, -326.35425047464,
-320.995337974641, -322.897621307974, -332.359865752418, -334.224474085751,
-330.744037974641, -326.290417141307, -332.752226863529, -332.753126863529,
-327.127282419085, -319.043900474639, -325.434953252418, -333.372468530195,
-332.215378252418, -325.17241297464, -329.32427547464, -339.435854641307,
-335.748814363529, -328.042222696863, -325.43886297464, -333.889821307973,
-333.892022696862, -328.416782419085, -321.215764363529, -329.921919919085,
-335.179329641307, -331.963644919085, -325.831231030197, -329.494412974641,
-336.949724085752, -335.052633807974, -326.407667141307, -326.669296307974,
-334.642417141307, -332.586454641308, -325.244890752418, -322.125990752418,
-332.469236585751, -336.083835196863, -332.428499085752, -328.023075474639,
-336.788682419085, -343.171972696863, -342.609490752418, -335.843854641307,
-339.158469919085, -344.622922696863, -339.954547696862, -331.644222696863,
-330.570807419084, -338.225764363529, -337.138574085751, -330.120342141307,
-325.023501863529, -333.010912974641, -333.538739363529, -326.497949085751,
-317.789147696862, -326.500901863529, -334.027447696862, -334.554806030196,
-328.659399085752, -332.592911585751, -341.838193530196, -341.283515752418,
-333.220168530196, -330.043024085752, -337.736175474641, -337.196274085751,
-332.116003252418, -327.442308807974, -338.168871307973, -341.512240752418,
-337.127721307974, -327.266132419085, -333.117636585752, -338.286093530196,
-335.125961585751, -327.596443530195, -331.173719919084, -341.312529641307,
-339.943233807973, -333.736206030196, -332.089414363529, -341.898035196863,
-341.840371307974, -335.138140752417, -328.132449085751, -337.483103252418,
-340.518328252418, -336.068264363528, -327.158939363529, -330.366781030196,
-335.523525474639, -332.762261585751, -324.595597696862, -327.470087974639,
-337.185801863529, -336.231885196862, -330.103854641306, -326.953156030195,
-337.921472696863, -341.125954641307, -336.421557419085, -329.548231030196,
-336.176665752418, -338.798690752418, -334.864249085752, -327.130714363529,
-330.242207419084, -336.85596297464, -335.480551863529, -331.116842141307,
-333.204106030196, -344.183626863529, -346.635036585751, -339.518718530196,
-334.60820047464, -343.45835047464, -346.979211585752, -343.77783797464,
-336.167997696863, -343.068593530196, -347.506990752418, -343.820393530196,
-335.568657419085, -339.22831297464, -346.487056030195, -344.816490752418,
-338.173265752418, -336.753925474641, -344.800375474641, -344.26353797464,
-338.417619919085, -331.613749085751, -340.444467141307, -343.729247696862,
-339.565626863529, -329.344668530196, -334.792554641307, -340.697133807973,
-336.876385196863, -328.094694919085, -331.424157419085, -341.145875474641,
-338.93063797464, -331.889428252418, -330.009576863528, -340.104201863529,
-342.305201863529, -335.681101863529, -328.078585196863, -337.092471307974,
-343.080397696863, -340.488833807973, -332.81781297464, -338.30123797464,
-343.387683807973, -337.936317141307, -327.916785196862, -328.454883807973,
-338.294354641307, -338.206704641307, -332.995347696862, -330.072069919085,
-339.76392547464, -338.753139363529, -332.812183807974, -325.062540752418,
-332.977019919084, -337.306956030196, -334.893785196862, -327.582568530195,
-332.530837974641, -339.816928252418, -338.202356030195, -328.079967141307,
-328.10071297464, -339.78435186353, -338.439776863529, -331.816094919084,
-329.835108807973, -342.151515752418, -345.247056030196, -339.272671307973,
-330.621157419085, -336.293403252418, -340.720553252418, -336.512268530196,
-326.219731030196, -329.36591297464, -339.381578252418, -339.007482419085,
-331.211058807973, -330.474836585751, -340.460131030195, -339.237511585752,
-331.514871307974, -324.002649085752, -332.197085196862, -334.023542141307,
-328.97147547464, -319.651978252418, -326.203136585751, -332.016567141307,
-328.500554641306, -319.691517141307, -323.030058807973, -333.807342141307,
-333.908064363529, -328.720169919085, -325.653653252418, -336.487032419084,
-338.091657419085, -333.650728252418, -326.085131030197, -335.172325474641,
-339.692781030196, -334.460192141307, -324.675851863529, -327.765651863529,
-332.538265752418, -330.836494919084, -328.461012974641, -332.744333807973,
-340.815140752418, -339.15458797464, -333.915018530197, -328.715283807974,
-335.221803252419, -336.401090752418, -332.041331030195, -325.352903252418,
-333.848103252419, -339.132468530196, -335.944543530196, -329.677172696862,
-335.740454641307, -344.029594919085, -340.953658807974, -331.657675474641,
-329.535169919084, -335.471393530196, -332.983167141307, -325.268335196862,
-320.822069919085, -333.430572696863, -337.216206030196, -332.750965752419,
-326.476585196863, -334.610683807974, -340.521747696863, -336.269624085752,
-326.032210196862, -328.951404641307, -335.167164363529, -331.886821307973,
-324.062199085751, -320.895856030196, -330.720533807974, -330.327764363529,
-323.917232419085, -319.331631030195, -329.189690752418, -333.477496307974,
-330.006618530195, -324.218785196863, -331.086819919085, -336.729660196863,
-334.155047696862, -324.813494919085, -327.216160196863, -336.344254641307,
-334.581299085751, -326.912028252418, -324.711415752417, -334.915846307974,
-335.197121307973, -329.198435196862, -322.516953252418, -331.57120047464,
-335.223156030195, -331.271781030195, -321.402407419084, -325.378924085751,
-330.938869919085, -329.553894919084, -322.737496307973, -324.003061585751,
-334.158472696862, -334.968397696862, -329.107721307974, -327.335099085752,
-338.25133797464, -339.516565752419, -334.450167141307, -326.697957419085,
-332.408443530196, -332.284790752418, -327.685846307973, -320.028282419084,
-325.584208807973, -333.892847696863, -333.454068530196, -326.711442141308,
-326.384131030196, -336.802986585752, -334.88462547464, -323.01483797464,
-315.379821307973, -327.309514363529, -332.918918530196, -331.099040752418,
-325.954918530196, -336.572486585751, -341.964910196862, -335.56171297464,
-326.668904641307, -327.428399085751, -335.059092141307, -333.257458807974,
-325.64721297464, -326.579868530195, -339.935818530195, -342.798668530196,
-335.806499085751, -328.058179641307, -335.685339363529, -337.123769919084,
-333.054331030196, -327.047819919085, -333.673631030196, -338.990536585751,
-334.833854641307, -325.447126863529, -327.674915752418, -336.87548797464,
-335.513093530196, -326.49010186353, -324.153235196862, -335.747449085751,
-340.214932419084, -335.26397547464, -327.556053252418, -336.001347696862,
-337.643742141307, -331.718407419085, -320.760057419084, -323.097251863529,
-330.061717141306, -329.50008797464, -325.28611297464, -328.289103252418,
-339.169781030196, -336.907200474641, -329.628128252418, -326.35898797464,
-336.183208807974, -337.980111585752, -332.890281030196, -326.666308807974,
-335.575132419085, -339.519174085752, -334.813443530197, -324.752239363529,
-330.147143530195, -337.446244919084, -335.172385196862, -328.031332419085,
-326.884160196863, -337.392429641307, -337.494010196863, -330.728390752419,
-324.575574085751, -333.631554641307, -335.900247696863, -332.076476863529,
-325.016990752418, -331.997007419084, -337.178467141307, -334.719262974641,
-328.230114363529, -333.420171307973, -342.084697696863, -339.583412974641,
-329.634624085751, -326.342922696862, -336.784333807974, -337.569519919085,
-331.982731030195, -327.612899085751, -339.300381030196, -342.711442141307,
-337.19452547464, -327.465342141308, -331.112865752418, -337.467947696862,
-335.264839363529, -328.359451863529, -330.97920047464, -340.216660196862,
-339.597740752418, -331.787849085751, -327.493985196863, -339.985321307973,
-340.816046307974, -332.316061585751, -326.019156030196, -334.852886585752,
-339.630711585751, -338.543051863529, -331.131129641307, -336.552936585751,
-343.001792141307, -340.726199085751, -333.14417547464, -334.834624085752,
-344.594219919085, -345.03996436353, -338.839890752418, -335.817275474641,
-347.498268530196, -349.144569919084, -344.318194919084, -337.042996307974,
-344.012585196863, -350.506115752418, -345.839889363529, -336.617222696862,
-341.731758807974, -349.143260196862, -347.56731297464, -338.472719919085,
-337.443568530196, -346.178836585752, -344.265247696863, -337.60685047464,
-331.532994919084, -340.334653252418, -345.444497696863, -341.596419919085,
-333.37895047464, -340.312939363529, -345.897197696863, -342.291735196862,
-335.23443936353, -336.936646307974, -345.905442141308, -346.21491297464,
-339.175429641307, -337.858081030197, -347.914558807974, -348.278715752418,
-343.278737974641, -336.79257547464, -345.574683807973, -347.591214363529,
-343.174983807974, -337.642853252418, -344.625007419084, -351.345108807973,
-347.083396307974, -337.052147696862, -338.651921307973, -348.306040752418,
-347.260461585751, -337.432786585752, -333.405818530197, -342.889789363529,
-344.527817141307, -339.51395186353, -331.045587974641, -339.511290752418,
-345.393272696862, -341.757008807974, -333.263932419085, -337.40365186353,
-342.746997696862, -341.27348797464, -335.19580186353, -338.378904641307,
-350.011318530195, -349.84250047464, -342.12290047464, -337.521279641306,
-347.441871307974, -347.405443530196, -341.230969919084, -334.000469919084,
-341.930557419085, -346.131979641307, -341.916731030195, -332.728361585751,
-337.408431030196, -347.159782419085, -344.614233807974, -337.578635196862,
-337.252185196862, -346.616371307974, -347.776096307974, -341.754831030196,
-334.287333807974, -341.664071307974, -341.958806030195, -336.346853252418,
-328.031293530197, -336.654951863529, -344.000056030196, -340.536083807974,
-332.565157419085, -336.695306030196, -346.016967141308, -344.482693530196,
-337.793028252418, -335.778210196862, -343.847242141307, -343.280874085751,
-337.110746307973, -331.388922696863, -341.453999085751, -345.301265752418,
-343.666294919085, -335.986561585751, -341.44398936353, -346.384136585751,
-342.773053252418, -333.285643530197, -335.549461585752, -344.908910196862,
-343.588168530196, -335.688447696862, -332.309049085752, -341.941876863529,
-344.908049085752, -341.870725474641, -337.400204641307, -347.130197696862,
-350.306021307973, -343.831661585751, -333.518261585751, -338.387133807973,
-344.027833807974, -338.960562974641, -332.61116436353, -336.12803936353,
-347.672342141307, -347.081656030196, -339.505710196862, -334.884117141306,
-343.750015752418, -343.74048797464, -338.660219919084, -332.635015752419,
-342.572487974641, -348.257417141307, -344.587297696863, -335.369847696863,
-340.068132419085, -348.024346307974, -345.639824085751, -338.516557419084,
-337.86285047464, -348.309049085751, -346.865506030195, -339.395456030195,
-334.819889363529, -342.787767141307, -345.098567141307, -341.038694919084,
-334.26566297464, -343.578118530196, -349.539117141306, -346.888528252418,
-340.056485196862, -341.928826863529, -351.764364363529, -350.21663797464,
-342.424974085751, -340.574222696862, -353.246974085751, -356.319483807974,
-348.424103252418, -339.160176863528, -348.671450474641, -352.309503252418,
-347.964031030196, -339.131972696862, -344.289928252418, -350.156269919085,
-348.378514363528, -340.680592141307, -343.675333807974, -350.920590752418,
-347.817747696863, -340.256129641307, -337.699211585751, -348.602447696863,
-349.339436585751, -343.016981030196, -335.294721307974, -342.882375474641,
-345.522149085752, -341.641968530196, -333.439261585751, -339.065526863529,
-347.723553252418, -348.653253252418, -340.697954641307, -340.395394919084,
-349.738736585752, -348.173207419085, -340.576849085751, -338.616279641307,
-348.241732419085, -349.23742686353, -343.613308807973, -336.74736297464,
-346.703169919085, -352.245636585751, -349.63193936353, -342.096767141307,
-348.529562974641, -358.600671307974, -358.012151863529, -350.619354641307,
-348.894926863529, -359.595465752418, -358.936629641307, -353.129778252418,
-346.597619919086, -354.727393530196, -357.138022696863, -352.933167141307,
-344.715258807973, -351.414456030196, -357.019531030197, -352.963544919085,
-343.695318530196, -348.231733807973, -357.530803252418, -356.12805047464,
-347.408007419085, -342.249971307973, -349.89625047464, -349.118944919084,
-340.687114363529, -335.524731030196, -347.250860196863, -351.056881030195,
-346.51307686353, -340.368943530196, -348.528178252418, -354.906271307974,
-350.868521307974, -340.766947696862, -342.813771307974, -349.915378252418,
-347.083794919085, -340.71750047464, -336.924042141307, -346.898029641307,
-349.818618530196, -343.577328252418, -336.543853252418, -346.250003252418,
-351.17105186353, -347.603364363529, -339.97052547464, -345.400599085752,
-351.713172696862, -347.088172696862, -338.269592141307, -337.720922696862,
-347.911065752418, -347.406315752418, -338.761365752418, -333.458768530196,
-343.656687974641, -346.520465752418, -341.454940752418, -333.796478252418,
-343.134549085751, -347.859442141307, -344.489215752418, -335.56472547464,
-339.229589363529, -349.018274085751, -348.730540752418, -340.539571307973,
-339.963968530196, -351.258285196863, -350.717518530196, -344.090779641308,
-340.165367141307, -349.858499085752, -352.756115752418, -348.200592141307,
-341.625469919085, -346.779511585751, -348.21057547464, -344.773942141307,
-337.735307419084, -341.70615047464, -351.258708807974, -349.514378252418,
-343.266985196863, -341.217871307974, -350.41058936353, -348.373729641307,
-342.023953252418, -335.723097696863, -345.980042141307, -350.46967547464,
-347.415447696862, -339.156322696863, -344.483507419085, -349.060397696863,
-345.454076863529, -338.661536585751, -341.343743530195, -350.273061585751,
-349.875072696863, -342.31343936353, -338.421433807973, -348.058597696862,
-348.486994919085, -341.870732419084, -334.038643530196, -343.97438797464,
-349.874193530196, -348.947987974641, -343.140211585752, -349.259103252418,
-354.266065752418, -350.453199085751, -342.06256436353, -340.206151863529,
-347.878926863529, -348.04130047464, -342.87507686353, -340.427515752418,
-349.749129641307, -353.089861585752, -347.362871307973, -337.71223797464,
-345.865861585752, -350.031261585751, -345.667143530197, -340.303833807974,
-346.798590752417, -354.226574085751, -351.079033807973, -342.436243530196,
-340.421418530196, -349.590715752418, -351.150642141307, -347.722707419085,
-344.563364363529, -353.486104641307, -354.076835196863, -347.290854641307,
-339.39193797464, -344.84330186353, -349.029968530196, -346.786529641307,
-338.334257419084, -343.071174085752, -353.385424085752, -351.060385196863,
-343.55301297464, -341.037756030196, -351.883244919085, -353.672772696862,
-348.797383807974, -342.502601863529, -349.543831030195, -352.901068530196,
-348.769179641307, -339.626321307974, -346.813358807974, -352.428408807974,
-348.129839363529, -337.360751863529, -338.735668530196, -347.892390752418,
-347.002654641308, -341.606417141307, -338.76171297464, -347.044276863528,
-347.572079641308, -342.585147696862, -337.097225474641, -347.147562974641,
-351.078168530196, -347.874651863529, -339.260032419085, -343.281240752417,
-349.430162974641, -346.528715752418, -338.70233797464, -339.721350474641,
-348.452267141307, -346.336394919084, -340.283914363528, -337.881931030196,
-350.978719919084, -353.095936585751, -347.894592141307, -340.911035196863,
-347.683154641307, -352.646746307974, -349.810836585751, -342.655418530196,
-347.179236585751, -353.918408807974, -350.397582419084, -342.224415752418,
-339.899715752417, -347.378469919085, -345.940557419085, -339.650231030195,
-335.12542547464, -346.179253252418, -351.144549085751, -347.593374085751,
-340.021283807973, -347.848519919085, -352.771726863529, -349.477492141307,
-340.435797696863, -342.101993530196, -351.236554641307, -351.497010196863,
-346.003531030197, -344.619353252418, -355.006956030195, -354.974442141307,
-347.762649085752, -341.225697696862, -349.405815752418, -354.056435196862,
-348.64240047464, -339.852874085751, -344.26473797464, -349.595072696862,
-345.584675474641, -337.634165752418, -338.131935196863, -345.804704641307,
-343.728783807974, -337.057793530196, -335.34027547464, -344.65147686353,
-344.70795186353, -340.285222696862, -334.380231030196, -344.623269919085,
-349.793442141307, -345.903710196863, -338.121565752418, -344.33462686353,
-351.233521307974, -346.898042141307, -336.129824085751, -337.605890752418,
-347.728740752418, -346.168454641307, -336.409819919085, -332.455761585751,
-347.595093530195, -351.861167141307, -344.755279641307, -334.020847696863,
-339.97741436353, -345.913212974641, -341.659322696862, -334.890215752418,
-340.77211297464, -348.590303252418, -345.319506030195, -338.332535196862,
-339.796882419084, -348.789096307974, -348.948389363529, -342.304743530196,
-338.366321307973, -351.011722696863, -356.44100047464, -351.801693530196,
-343.219208807974, -350.663174085751, -357.263917141307, -355.382874085752,
-347.338603252418, -348.990885196862, -356.976176863529, -354.694065752419,
-345.437569919084, -343.577694919085, -352.261751863529, -353.781744919085,
-348.798674085752, -342.758274085751, -351.635421307974, -354.292590752418,
-348.102501863529, -339.395894919085, -344.356346307974, -350.316451863529,
-346.624194919085, -338.811697696862, -342.269958807973, -351.454024085751,
-345.515681030196, -337.267867141307, -337.965290752418, -351.610676863529,
-353.177547696862, -345.882003252418, -339.304692141307, -347.802081030196,
-350.895426863529, -345.97835047464, -336.510861585751, -343.079846307974,
-349.347681030196, -348.059937974641, -340.549531030195, -341.804257419084,
-350.693160196862, -350.602917141307, -343.298667141307, -338.964418530195,
-349.139769919085, -352.462711585751, -349.750517141307, -344.741354641307,
-352.661028252419, -358.013082419085, -353.774382419085, -343.604215752418,
-346.927953252417, -352.552529641307, -350.802136585752, -345.936253252418,
-346.751729641307, -356.805565752418, -356.628961585751, -349.132007419085,
-342.729618530196, -350.600707419084, -352.019142141307, -346.915435196862,
-340.261306030195, -349.328561585751, -356.08598797464, -354.216194919085,
-346.036522696863, -350.386724085751, -357.033010196863, -354.384594919085,
-347.137222696862, -344.783422696862, -354.063804641307, -354.521757419084,
-349.107696307973, -342.029085196863, -350.976568530196, -353.200135196863,
-347.400669919085, -339.007418530196, -346.684894919084, -353.322094919084,
-349.685383807974, -342.477474085752, -344.592014363529, -353.41798797464,
-350.216621307973, -343.580340752418)), .Names = c("time", "value"
), row.names = c(NA, 1000L), class = "data.frame")
Raw
EMD.R
########
# This code comes from the EMD package (http://cran.r-project.org/web/packages/EMD/index.html) but incorporates optimizations that resulted from this thread: https://stat.ethz.ch/pipermail/r-help/2011-February/268550.html
########
diff1 <- function (x){
return(x[-1]-x[-length(x)])
}
rev1 <- function (x){
return(x[length(x):1L])
}
Sequence <- function(nvec) {
seq_len(sum(nvec)) - rep(cumsum(c(0L,nvec[-length(nvec)])), nvec)
}
isEndOfStrictlyIncreasingRun <- function(x) {
c(FALSE, diff1(diff1(x) > 0) < 0, FALSE)
}
isStartOfStrictlyIncreasingRun <- function(x) {
c(FALSE, diff1(diff1(x) <= 0) < 0, FALSE)
}
isEndOfStrictlyDecreasingRun <- function(x) {
isEndOfStrictlyIncreasingRun(-x)
}
isStartOfStrictlyDecreasingRun <- function(x) {
isStartOfStrictlyIncreasingRun(-x)
}
isStartOfZeroRun <- function(x) {
x==0 & c(TRUE, x[-length(x)]!=0)
}
nToEndOfCurrentFlatRun <- function(x) {
# for each element of x, how far to end of current
# run of equal values.
rev1( Sequence(rle(rev1(x))$lengths) - 1L)
}
indexOfEndOfCurrentFlatRun <- function(x) {
# should be a more direct way of doing this, but this is pretty quick
seq_len(length(x)) + nToEndOfCurrentFlatRun(x)
}
extrema <- function(x) {
f <- indexOfEndOfCurrentFlatRun(x)
isMaxStart <- isEndOfStrictlyIncreasingRun(x) & isStartOfStrictlyDecreasingRun(x)[f]
maxindex <- cbind(which(isMaxStart), f[isMaxStart], deparse.level=0)
isMinStart <- isEndOfStrictlyDecreasingRun(x) & isStartOfStrictlyIncreasingRun(x)[f]
minindex <- cbind(which(isMinStart), f[isMinStart], deparse.level=0)
# zero-crossings are bit weird: Report either the before-zero/after-zero
# pair or the start and stop of a a run of zero's (even if the run is
# not part of an actual zero-crossing). Do them separately then sort.
# Also, if the entire sequence never actually crosses 0, do
# not report the zero-touchings.
# Also, if length(x)==2, the original doesn't report a zero-crossing or
# a zero run. We do not copy that.
if (max(x) > 0 && min(x) < 0) {
indexOfZeroCrossingStart <- which(c(abs(diff1(sign(x)))==2,FALSE))
indexOfZeroCrossingEnd <- indexOfZeroCrossingStart + 1L
indexOfZeroRunStart <- which(isStartOfZeroRun(x))
indexOfZeroRunEnd <- f[indexOfZeroRunStart]
crossStart <- c(indexOfZeroCrossingStart, indexOfZeroRunStart)
cross <- cbind(crossStart, c(indexOfZeroCrossingEnd,indexOfZeroRunEnd), deparse.level=0)[order(crossStart),, drop=FALSE]
} else {
cross <- cbind(integer(), integer())
}
# extrema likes to return NULL instead of a zero-row matrix,
# so we follow it. Zero-row matrices are easier to deal with.
list(
minindex=if (nrow(minindex)) minindex else NULL,
maxindex=if (nrow(maxindex)) maxindex else NULL,
nextreme=nrow(minindex) + nrow(maxindex),
cross=if(nrow(cross)) cross else NULL,
ncross=nrow(cross) # extrema() returns numeric, not integer, here
)
}
emd <- function (xt, tt = NULL, tol = sd(xt) * 0.1^2, max.sift = 20,
stoprule = "type1", boundary = "periodic", smlevels = c(1),
sm = "none", spar = NA, weight = 20, check = FALSE, max.imf = 10,
plot.imf = TRUE, interm = NULL, print_progress = TRUE)
{
if (is.ts(xt))
xt <- as.numeric(xt)
if (is.null(tt))
tt <- 1:length(xt)
if (is.null(interm) || all(interm <= 0))
intermtest <- FALSE
else intermtest <- TRUE
ndata <- length(xt)
ndatam1 <- ndata - 1
residue <- xt
rangext <- range(residue)
imf <- NULL
j <- 1
repeat {
if (j > max.imf)
break
if ((any(j == smlevels) || smlevels == "all") & sm ==
"spline")
tmp <- extractimf(residue, tt, tol, max.sift, stoprule = stoprule,
boundary = boundary, sm = sm, spar = spar, check = check)
else tmp <- extractimf(residue, tt, tol, max.sift, stoprule = stoprule,
boundary = boundary, check = check)
if (is.null(tmp$imf)) {
break
}
if (plot.imf) {
plot(tt, residue, type = "l", xlab = "", ylab = "",
main = paste(j - 1, "-th residue=", j, "-th imf+",
j, "-th residue", sep = ""))
abline(h = 0)
}
if (intermtest && length(interm) >= j && interm[j] >
0) {
tmpimf <- tmp$imf
tmpresidue <- tmp$residue
tmpinterm <- extrema(tmpimf)
tmpncross <- tmpinterm$ncross
zerocross <- as.numeric(round(apply(tmpinterm$cross,
1, mean)))
if (abs(tt[zerocross[3]] - tt[zerocross[1]]) > interm[j]) {
tmpresidue[1:zerocross[3]] <- tmpresidue[1:zerocross[3]] +
tmpimf[1:zerocross[3]]
tmpimf[1:zerocross[3]] <- 0
}
for (k in seq(3, tmpncross - 3, by = 2)) if (abs(tt[zerocross[k +
2]] - tt[zerocross[k]]) > interm[j]) {
tmpresidue[zerocross[k]:zerocross[k + 2]] <- tmpresidue[zerocross[k]:zerocross[k +
2]] + tmpimf[zerocross[k]:zerocross[k + 2]]
tmpimf[zerocross[k]:zerocross[k + 2]] <- 0
}
if (!(tmpncross%%2)) {
if (abs(tt[zerocross[tmpncross]] - tt[zerocross[tmpncross -
1]]) > interm[j]/2) {
tmpresidue[zerocross[tmpncross - 1]:ndata] <- tmpresidue[zerocross[tmpncross -
1]:ndata] + tmpimf[zerocross[tmpncross -
1]:ndata]
tmpimf[zerocross[tmpncross - 1]:ndata] <- 0
}
}
else {
if (abs(tt[zerocross[tmpncross]] - tt[zerocross[tmpncross -
2]]) > interm[j])
tmpresidue[zerocross[tmpncross - 2]:ndata] <- tmpresidue[zerocross[tmpncross -
2]:ndata] + tmpimf[zerocross[tmpncross -
2]:ndata]
tmpimf[zerocross[tmpncross - 2]:ndata] <- 0
}
tmp$imf <- tmpimf
tmp$residue <- tmpresidue
}
imf <- cbind(imf, tmp$imf)
residue <- tmp$residue
if (plot.imf) {
plot(tt, imf[, j], type = "l", xlab = "", ylab = "",
main = paste(j, "-th imf", sep = ""))
abline(h = 0)
plot(tt, residue, type = "l", xlab = "", ylab = "",
main = paste(j, "-th residue", sep = ""))
abline(h = 0)
locator(1)
}
if(print_progress){
print(paste('IMF #',j,' obtained.',sep=''))
}
j <- j + 1
}
list(imf = imf, residue = residue, nimf = j - 1)
}
extractimf <- function (residue, tt = NULL, tol = sd(residue) * 0.1^2, max.sift = 20,
stoprule = "type1", boundary = "periodic", sm = "none", spar = NA,
weight = 20, check = FALSE)
{
if (boundary == "none")
minextrema <- 4
else minextrema <- 2
if (sm == "spline" & is.null(spar))
stop("Provide the smoothing parameter of spline smoothing.\n")
ndata <- length(residue)
ndatam1 <- ndata - 1
if (is.ts(residue))
residue <- as.numeric(residue)
if (is.null(tt))
tt <- 1:length(residue)
emin <- emax <- em <- h <- imf <- NULL
n2data <- 2 * ndata
tt2 <- 1:n2data
n3data <- 3 * ndata
n3datam1 <- n3data - 1
tt3 <- 1:n3data
input <- residue
rangext <- range(residue)
j <- 1
repeat {
tmp <- extrema(input)
if (tmp$nextreme <= minextrema)
break
if (j == 1 || boundary == "wave") {
minindex <- unique(c(t(tmp$minindex)))
minn <- length(minindex)
maxindex <- unique(c(t(tmp$maxindex)))
maxn <- length(maxindex)
extminindex <- minindex
extmaxindex <- maxindex
tmpwavefreq1 <- diff1(sort(c(tt[1], tt[minindex[1]],
tt[maxindex[1]])))
tmpwavefreq2 <- diff1(sort(c(tt[minindex[minn]], tt[maxindex[maxn]],
tt[ndata])))
if (input[1] <= input[minindex[1]] && input[1] <=
input[maxindex[1]]) {
extminindex <- c(1, extminindex)
wavefreq1 <- 2 * tmpwavefreq1[1]
}
else if (input[1] >= input[minindex[1]] && input[1] >=
input[maxindex[1]]) {
extmaxindex <- c(1, extmaxindex)
wavefreq1 <- 2 * tmpwavefreq1[1]
}
else if (input[1] >= (input[minindex[1]] + input[maxindex[1]])/2) {
wavefreq1 <- tmpwavefreq1[2] + max(tmpwavefreq1[2],
2 * tmpwavefreq1[1])
}
else {
wavefreq1 <- tmpwavefreq1[2] + max(tmpwavefreq1[2],
round(1.5 * tmpwavefreq1[1]))
}
if (input[ndata] <= input[minindex[minn]] && input[ndata] <=
input[maxindex[maxn]]) {
extminindex <- c(extminindex, ndata)
wavefreq2 <- 2 * tmpwavefreq2[2]
}
else if (input[ndata] >= input[minindex[minn]] &&
input[ndata] >= input[maxindex[maxn]]) {
extmaxindex <- c(extmaxindex, ndata)
wavefreq2 <- 2 * tmpwavefreq2[2]
}
else if (input[ndata] >= (input[minindex[minn]] +
input[maxindex[maxn]])/2) {
wavefreq2 <- tmpwavefreq2[1] + max(tmpwavefreq2[1],
2 * tmpwavefreq2[2])
}
else {
wavefreq2 <- tmpwavefreq2[1] + max(tmpwavefreq2[1],
round(1.5 * tmpwavefreq2[2]))
}
extminn <- length(extminindex)
extmaxn <- length(extmaxindex)
extttminindex <- c(tt[extminindex[1]] - 4:1 * wavefreq1,
tt[extminindex], tt[extminindex[extminn]] + 1:4 *
wavefreq2)
extttmaxindex <- c(tt[extmaxindex[1]] - 4:1 * wavefreq1,
tt[extmaxindex], tt[extmaxindex[extmaxn]] + 1:4 *
wavefreq2)
if (sm == "none") {
f <- splinefun(extttminindex, c(rep(input[extminindex[1]],
4), input[extminindex], rep(input[extminindex[extminn]],
4)))
emin <- cbind(emin, f(tt))
f <- splinefun(extttmaxindex, c(rep(input[extmaxindex[1]],
4), input[extmaxindex], rep(input[extmaxindex[extmaxn]],
4)))
emax <- cbind(emax, f(tt))
}
else if (sm == "spline") {
f <- sreg(extttminindex, c(rep(input[extminindex[1]],
4), input[extminindex], rep(input[extminindex[extminn]],
4)), lambda = spar)
llambda <- f$lambda * weight
f <- sreg(extttminindex, c(rep(input[extminindex[1]],
4), input[extminindex], rep(input[extminindex[extminn]],
4)), lambda = llambda)
llambda <- f$lambda
emin <- cbind(emin, predict(f, tt))
f <- sreg(extttmaxindex, c(rep(input[extmaxindex[1]],
4), input[extmaxindex], rep(input[extmaxindex[extmaxn]],
4)), lambda = spar)
ulambda <- f$lambda * weight
f <- sreg(extttmaxindex, c(rep(input[extmaxindex[1]],
4), input[extmaxindex], rep(input[extmaxindex[extmaxn]],
4)), lambda = ulambda)
ulambda <- f$lambda
emax <- cbind(emax, predict(f, tt))
}
em <- cbind(em, (emin[, j] + emax[, j])/2)
if (check) {
plot(tt, input, type = "l", col = 3, xlab = "",
ylab = "", main = paste("Boundary = ", boundary,
sep = ""))
points(tt[unique(c(t(tmp$minindex)))], input[unique(c(t(tmp$minindex)))],
col = 4)
points(tt[unique(c(t(tmp$maxindex)))], input[unique(c(t(tmp$maxindex)))],
col = 2)
lines(tt, emin[, j], col = 4)
lines(tt, emax[, j], col = 2)
lines(tt, em[, j])
locator(1)
}
}
else {
if (boundary == "none") {
minindex <- c(1, unique(c(t(tmp$minindex))),
ndata)
maxindex <- c(1, unique(c(t(tmp$maxindex))),
ndata)
fmin <- splinefun(tt[minindex], input[minindex])
emin <- cbind(emin, fmin(tt))
fmax <- splinefun(tt[maxindex], input[maxindex])
emax <- cbind(emax, fmax(tt))
em <- cbind(em, (emin[, j] + emax[, j])/2)
if (check) {
plot(tt, input, type = "l", col = 3, xlab = "",
ylab = "", main = paste("Boundary = ", boundary,
sep = ""))
points(tt[unique(c(t(tmp$minindex)))], input[unique(c(t(tmp$minindex)))],
col = 4)
points(tt[unique(c(t(tmp$maxindex)))], input[unique(c(t(tmp$maxindex)))],
col = 2)
lines(tt, emin[, j], col = 4)
lines(tt, emax[, j], col = 2)
lines(tt, em[, j])
locator(1)
}
}
if (boundary == "symmetric" || boundary == "periodic") {
if (boundary == "symmetric") {
inputext <- c(rev1(input[-1]), input, rev1(input)[-1])
ttext <- c(tt[1] - rev1(cumsum(diff1(tt))), tt,
tt[ndata] + cumsum(rev1(diff1(tt))))
tmp <- extrema(inputext)
}
else if (boundary == "periodic") {
inputext <- c(input[-ndata], input, input[-1])
ttext <- c(rev1(tt[1] - cumsum(rev1(diff1(tt)))),
tt, tt[ndata] + cumsum(diff1(tt)))
tmp <- extrema(inputext)
}
minindex <- unique(c(t(tmp$minindex)))
maxindex <- unique(c(t(tmp$maxindex)))
if (sm == "none") {
fmin <- splinefun(ttext[minindex], inputext[minindex])
tmpmin <- fmin(ttext[ndata:(n2data - 1)])
fmax <- splinefun(ttext[maxindex], inputext[maxindex])
tmpmax <- fmax(ttext[ndata:(n2data - 1)])
}
else if (sm == "spline") {
fmin <- sreg(ttext[minindex], inputext[minindex],
lambda = spar)
llambda <- fmin$lambda * weight
fmin <- sreg(ttext[minindex], inputext[minindex],
lambda = llambda)
tmpmin <- predict(fmin, ttext[ndata:(n2data -
1)])
fmax <- sreg(ttext[maxindex], inputext[maxindex],
lambda = spar)
ulambda <- fmax$lambda * weight
fmax <- sreg(ttext[maxindex], inputext[maxindex],
lambda = ulambda)
tmpmax <- predict(fmax, ttext[ndata:(n2data -
1)])
}
emin <- cbind(emin, tmpmin[1:ndata])
emax <- cbind(emax, tmpmax[1:ndata])
tmpm <- (tmpmin + tmpmax)/2
em <- cbind(em, tmpm[1:ndata])
if (check) {
plot(ttext[ndata:(n2data - 1)], inputext[ndata:(n2data -
1)], type = "l", col = 3, xlab = "", ylab = "",
main = paste("Boundary = ", boundary, sep = ""))
points(ttext[minindex[minindex > ndata & minindex <
n2data]], inputext[minindex[minindex > ndata &
minindex < n2data]], col = 4)
points(ttext[maxindex[maxindex > ndata & maxindex <
n2data]], inputext[maxindex[maxindex > ndata &
maxindex < n2data]], col = 2)
lines(ttext[ndata:(n2data - 1)], tmpmin, col = 4)
lines(ttext[ndata:(n2data - 1)], tmpmax, col = 2)
lines(ttext[ndata:(n2data - 1)], tmpm)
locator(1)
}
}
else if (boundary == "evenodd") {
inputeven <- c(input, rev1(input), input)
ttext <- c(tt, tt[ndata] + tt[ndata] - tt[ndatam1],
tt[ndata] + tt[ndata] - tt[ndatam1] + cumsum(rev1(diff1(tt))))
tmp <- extrema(inputeven)
minindex <- unique(c(t(tmp$minindex)))
minindex <- minindex[minindex <= n2data + minindex[1]]
maxindex <- unique(c(t(tmp$maxindex)))
maxindex <- maxindex[maxindex <= n2data + maxindex[1]]
if (sm == "none") {
f <- splinefun(ttext[minindex], inputeven[minindex])
emineven <- f(ttext[1:ndata])
f <- splinefun(ttext[maxindex], inputeven[maxindex])
emaxeven <- f(ttext[1:ndata])
}
else if (sm == "spline") {
fmin <- sreg(ttext[minindex], inputeven[minindex],
lambda = spar)
llambda <- fmin$lambda * weight
fmin <- sreg(ttext[minindex], inputeven[minindex],
lambda = llambda)
emineven <- predict(fmin, ttext[1:ndata])
fmax <- sreg(ttext[maxindex], inputext[maxindex],
lambda = spar)
ulambda <- fmax$lambda * weight
fmax <- sreg(ttext[maxindex], inputext[maxindex],
lambda = ulambda)
emaxeven <- predict(fmax, ttext[1:ndata])
}
inputodd <- c(input, -rev1(input), input)
tmp <- extrema(inputodd)
minindex <- unique(c(t(tmp$minindex)))
minindex <- minindex[minindex <= n2data + minindex[1]]
maxindex <- unique(c(t(tmp$maxindex)))
maxindex <- maxindex[maxindex <= n2data + maxindex[1]]
if (sm == "none") {
f <- splinefun(ttext[minindex], inputodd[minindex])
eminodd <- f(ttext[1:ndata])
f <- splinefun(ttext[maxindex], inputodd[maxindex])
emaxodd <- f(ttext[1:ndata])
}
else if (sm == "spline") {
fmin <- sreg(ttext[minindex], inputodd[minindex],
lambda = spar)
llambda <- fmin$lambda * weight
fmin <- sreg(ttext[minindex], inputodd[minindex],
lambda = llambda)
eminodd <- predict(fmin, ttext[1:ndata])
fmax <- sreg(ttext[maxindex], inputext[maxindex],
lambda = spar)
ulambda <- fmax$lambda * weight
fmax <- sreg(ttext[maxindex], inputext[maxindex],
lambda = ulambda)
emaxodd <- predict(fmax, ttext[1:ndata])
}
emin <- cbind(emin, (emineven + eminodd)/2)
emax <- cbind(emax, (emaxeven + emaxodd)/2)
em <- cbind(em, (emin[, j] + emax[, j])/2)
if (check) {
plot(tt, (inputeven + inputodd)/2[1:ndata],
type = "l", col = 3, ylab = "")
lines(tt, emin, col = 4)
lines(tt, emax, col = 2)
lines(tt, (emin + emax)/2)
locator(1)
}
}
}
h <- cbind(h, input - em[, j])
if (stoprule == "type1" && (all(abs(em[, j]) < tol) ||
j >= max.sift)) {
imf <- h[, j]
residue <- residue - imf
break
}
if (stoprule == "type2" && j >= 2) {
if (sum((h[2:ndatam1, j - 1] - h[2:ndatam1, j])^2/h[2:ndatam1,
j - 1]^2) < tol || j >= max.sift) {
imf <- h[, j]
residue <- residue - imf
break
}
}
input <- h[, j]
j <- j + 1
}
if (check)
list(emin = emin, emax = emax, em = em, h = h, imf = imf,
residue = residue, niter = j)
else list(imf = imf, residue = residue, niter = j)
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment