|
| 1 | +# Find Periodicity Using Autocorrelation |
| 2 | + |
| 3 | +([Example](https://www.mathworks.com/help/signal/ug/find-periodicity-using-autocorrelation.html) from Matlab documentation) |
| 4 | + |
| 5 | +The autocorrelation sequence of a periodic signal has the same cyclic characteristics as the signal itself. Thus, autocorrelation can help verify the presence of cycles and determine their durations. |
| 6 | + |
| 7 | +Consider a set of temperature data collected by a thermometer inside an office building. |
| 8 | +The device takes a reading every half hour for four months. Load the data and plot it. |
| 9 | +Subtract the mean to concentrate on temperature fluctuations. Convert the temperature to degrees Celsius. |
| 10 | +Measure time in days. The sample rate is thus 2 measurements/hour × 24 hours/day = 48 measurements/day. |
| 11 | + |
| 12 | + |
| 13 | +\begin{examplefig}{} |
| 14 | +```julia |
| 15 | +using GMT |
| 16 | + |
| 17 | +D = gmtread(TESTSDIR * "assets/temps.dat"); |
| 18 | +tempC = (D-32)*5/9; |
| 19 | +tempnorm = tempC - mean(tempC); |
| 20 | +fs = 2*24; |
| 21 | +t = (0:length(tempnorm) - 1)/fs; |
| 22 | +plot(t,tempnorm, xlabel="Time (days)", ylabel="Temperature (°C)", title="Temperature time series", figsize=(15, 6), show=true) |
| 23 | +``` |
| 24 | +\end{examplefig} |
| 25 | + |
| 26 | +The temperature does seem to oscillate, but the lengths of the cycles cannot be read out easily. |
| 27 | + |
| 28 | +Compute the autocorrelation of the temperature such that it is unity at zero lag. Restrict the |
| 29 | +positive and negative lags to three weeks. Note the double periodicity of the signal. |
| 30 | + |
| 31 | + |
| 32 | +\begin{examplefig}{} |
| 33 | +```julia |
| 34 | +using GMT # Hide |
| 35 | +D = gmtread(TESTSDIR * "assets/temps.dat"); # Hide |
| 36 | +tempC = (D-32)*5/9; # Hide |
| 37 | +tempnorm = tempC - mean(tempC);# Hide |
| 38 | +fs = 2*24;# Hide |
| 39 | + |
| 40 | +lags = 0:3*7*fs; |
| 41 | +ac = GMT.autocor(tempnorm, lags); |
| 42 | +locs = findpeaks(ac); |
| 43 | +plot(lags/fs, ac, figsize=(15, 7)) |
| 44 | +scatter!(lags[locs]/fs,ac[locs], show=1) |
| 45 | +``` |
| 46 | +\end{examplefig} |
| 47 | + |
| 48 | +Determine the short and long periods by finding the peak locations and determining the average time differences between them. |
| 49 | + |
| 50 | +To find the long period, restrict findpeaks to look for peaks separated by more than the short period and with a minimum height of 0.3. |
| 51 | + |
| 52 | +\begin{examplefig}{} |
| 53 | +```julia |
| 54 | +using GMT # Hide |
| 55 | + |
| 56 | +D = gmtread(TESTSDIR * "assets/temps.dat"); # Hide |
| 57 | +tempC = (D-32)*5/9; # Hide |
| 58 | +tempnorm = tempC-mean(tempC);# Hide |
| 59 | +fs = 2*24;# Hide |
| 60 | +lags = 0:3*7*fs; # Hide |
| 61 | +ac = GMT.autocor(tempnorm, lags);# Hide |
| 62 | +locs = findpeaks(ac);# Hide |
| 63 | + |
| 64 | +shortT = mean(diff(sort(locs))) / fs |
| 65 | +locs_long = findpeaks(ac, min_dist=ceil(shortT)*fs, min_height=0.3); |
| 66 | +longT = mean(diff(locs_long))/fs |
| 67 | +plot(lags/fs, ac, figsize=(15, 7)) |
| 68 | +scatter!(lags[locs]/fs,ac[locs], fill=:red) |
| 69 | +scatter!(lags[locs_long]/fs, ac[locs_long], fill=:blue, show=1) |
| 70 | +``` |
| 71 | +\end{examplefig} |
| 72 | + |
| 73 | +To a very good approximation, the autocorrelation oscillates both daily and weekly. |
| 74 | +This is to be expected, since the temperature in the building is higher when people |
| 75 | +are at work and lower at nights and on weekends. |
| 76 | + |
0 commit comments