|
@@ -0,0 +1,70 @@
|
|
|
+import numpy as np
|
|
|
+
|
|
|
+
|
|
|
+def get_c_p(angles, p):
|
|
|
+ return np.sum(np.cos(p * angles))
|
|
|
+
|
|
|
+
|
|
|
+def get_mean_c_p(angles, p):
|
|
|
+ return get_c_p(angles, p) / np.cumprod(np.shape(angles))
|
|
|
+
|
|
|
+
|
|
|
+def get_s_p(angles, p):
|
|
|
+ return np.sum(np.sin(p * angles))
|
|
|
+
|
|
|
+
|
|
|
+def get_mean_s_p(angles, p):
|
|
|
+ return get_s_p(angles, p) / np.cumprod(np.shape(angles))
|
|
|
+
|
|
|
+
|
|
|
+def get_r_p(angles, p):
|
|
|
+ return np.sqrt(get_c_p(angles, p) ** 2 + get_s_p(angles, p) ** 2)
|
|
|
+
|
|
|
+
|
|
|
+def get_t_p(angles, p):
|
|
|
+ mean_c_p = get_mean_c_p(angles, p)
|
|
|
+ mean_s_p = get_mean_s_p(angles, p)
|
|
|
+ t_p = np.arctan(mean_s_p / mean_c_p)[0]
|
|
|
+
|
|
|
+ if mean_c_p < 0:
|
|
|
+ t_p = t_p + np.pi
|
|
|
+ elif mean_s_p < 0 and mean_c_p > 0:
|
|
|
+ t_p = t_p + 2 * np.pi
|
|
|
+ return t_p
|
|
|
+
|
|
|
+
|
|
|
+def get_circular_correlation_measure(angles_1, angles_2):
|
|
|
+ t_1 = get_t_p(angles_1, 1)
|
|
|
+ t_2 = get_t_p(angles_2, 1)
|
|
|
+ r_c = np.sum(np.sin(angles_1 - t_1) * np.sin(angles_2 - t_2)) / np.sqrt(
|
|
|
+ np.sum(np.sin(angles_1 - t_1) ** 2) * np.sum(np.sin(
|
|
|
+ angles_2 - t_2) ** 2))
|
|
|
+ return r_c
|
|
|
+
|
|
|
+
|
|
|
+N = 1000
|
|
|
+
|
|
|
+random_angles_1 = np.random.uniform(-np.pi, np.pi, N)
|
|
|
+random_angles_2 = np.random.uniform(-np.pi, np.pi, N)
|
|
|
+angles_close_to_pi_half = np.random.vonmises(np.pi / 2.0, 20, N)
|
|
|
+angles_close_to_minus_pi_half = angles_close_to_pi_half-np.pi*np.random.uniform(0.9, 1.1, N)
|
|
|
+
|
|
|
+angles_description = ["random", "von_mises_pi_half", "von_mises_minus_pi_half"]
|
|
|
+angles_collection = [random_angles_1, angles_close_to_pi_half, angles_close_to_minus_pi_half]
|
|
|
+
|
|
|
+for desc, angles in zip(angles_description, angles_collection):
|
|
|
+ print(desc)
|
|
|
+ print("C_1 = {:.2f}".format(get_c_p(angles, 1)))
|
|
|
+ print("mean R_1 = {:.2f}".format(get_r_p(angles, 1) / N))
|
|
|
+ print("T_1 = {:.2f}".format(get_t_p(angles, 1)))
|
|
|
+ print()
|
|
|
+
|
|
|
+
|
|
|
+print(get_circular_correlation_measure(random_angles_1, angles_close_to_pi_half))
|
|
|
+print(get_circular_correlation_measure(angles_close_to_minus_pi_half, angles_close_to_pi_half))
|
|
|
+import matplotlib.pyplot as plt
|
|
|
+plt.hist(random_angles_1)
|
|
|
+plt.hist(angles_close_to_pi_half)
|
|
|
+plt.hist(angles_close_to_minus_pi_half)
|
|
|
+
|
|
|
+plt.show()
|