-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathdemo_stationaritytest.m
118 lines (96 loc) · 3.04 KB
/
demo_stationaritytest.m
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
111
112
113
114
115
116
% Demo script for stationarity test (Hauber et al 2021) using normalized data
%from e15c064 (Sigloch et al. 2020). Each individual section should not take
%much longer than a few seconds to compute on a modern machine.
tsInit
%% Simulate data (AR[2] process) with frequency modulation
nTimePoints = 1000;
nRealizations = 1;
decay = 100;
discardLength = 5*decay;
baseline = 20;
switchTime = 300;
switchLevel = 40;
decays = decay * ones(nTimePoints + discardLength,1);
period = baseline * ones(nTimePoints + discardLength,1);
period((discardLength+switchTime):end) = baseline+switchLevel;
arPars = tsCalcAr2Pars(period, decays);
[data, ~] = tsSimulateAR(arPars, nTimePoints, nRealizations, discardLength);
figure
plot(data)
hold on
xline(switchTime, '--')
hold off
xlabel('Time')
title('Simulated data (AR[2] process) with FM')
%% Perform stationarity test
statConfig.signLevel = 0.05;
statConfig.surrN = 1000;
statConfig.surrMethod = 'AFFT';
idata = 1;
[instPhase, instPhaseSurr, lb, ub, rejection] = ...
tsStationarityTest2(data(:,idata), statConfig);
figure
plot(instPhaseSurr(:,1:5:end), 'Color', 0.7*[1 1 1], 'HandleVisibility', 'off')
hold on
plot(instPhase, 'Color', 'k', 'LineWidth',2)
plot(lb,'k-.', 'LineWidth',2)
plot(ub,'k-.','HandleVisibility','off', 'LineWidth',2)
xline(switchTime,'--')
hold off
xlabel('Time')
title('Instantaneous phase - simulated data (AR[2] process) with FM')
%% Simulate data (AR[2] process) without frequency modulation
nTimePoints = 1000;
nRealizations = 1;
decay = 100;
discardLength = 5*decay;
baseline = 20;
decays = decay;
period = baseline;
arPars = tsCalcAr2Pars(period, decays);
[data, ~] = tsSimulateAR(arPars, nTimePoints, nRealizations, discardLength);
figure
plot(data)
xlabel('Time')
title('Simulated data (AR[2] process) without FM')
%% Perform stationarity test
statConfig.signLevel = 0.05;
statConfig.surrN = 1000;
statConfig.surrMethod = 'AFFT';
idata = 1;
[instPhase, instPhaseSurr, lb, ub, rejection] = ...
tsStationarityTest2(data(:,idata), statConfig);
figure
plot(instPhaseSurr(:,1:5:end), 'Color', 0.7*[1 1 1], 'HandleVisibility', 'off')
hold on
plot(instPhase, 'Color', 'k', 'LineWidth',2)
plot(lb,'k-.', 'LineWidth',2)
plot(ub,'k-.','HandleVisibility','off', 'LineWidth',2)
xline(switchTime,'--')
hold off
xlabel('Time')
title('Instantaneous phase - simulated data (AR[2] process) without FM')
%% Experimental data
import = csvread('demodata_stationaritytest.csv');
times = import(:,1); data = import(:,2);
figure
plot(times, data)
xlabel('Time [min]')
title('Data - e15c064')
%% Perform stationarity test
rng(2)
statConfig.signLevel = 0.05;
statConfig.surrN = 1000;
statConfig.surrMethod = 'AFFT';
idata = 1;
[instPhase, instPhaseSurr, lb, ub, rejection] = ...
tsStationarityTest2(data, statConfig);
figure
plot(times',instPhaseSurr(:,1:4:end), 'Color', 0.7*[1 1 1], 'HandleVisibility', 'off')
hold on
plot(times,instPhase, 'Color', 'k', 'LineWidth',2)
plot(times,lb,'k-.', 'LineWidth',2)
plot(times,ub,'k-.','HandleVisibility','off', 'LineWidth',2)
hold off
xlabel('Time')
title('Instantaneous phase - e15c064')