Update README.md

Signed-off-by: David Rotermund <54365609+davrot@users.noreply.github.com>

pywavelet/README.md

@@ -302,7 +302,121 @@ plt.show()
 
 This looks already better...
 
-
+### Fixing the problems -- Cone of influence
+
+If the look at the edges of the 2d plot, we see that the power tapers of. There regions are invalid results because part of the wavelet hangs outside of the time series. The larger the frequency, the larger the region. 
+
+```python
+import numpy as np
+import matplotlib.pyplot as plt
+import pywt
+
+
+# Calculate the wavelet scales we requested
+def calculate_wavelet_scale(
+    number_of_frequences: int,
+    frequency_range_min: float,
+    frequency_range_max: float,
+    dt: float,
+) -> np.ndarray:
+    s_spacing: np.ndarray = (1.0 / (number_of_frequences - 1)) * np.log2(
+        frequency_range_max / frequency_range_min
+    )
+    scale: np.ndarray = np.power(2, np.arange(0, number_of_frequences) * s_spacing)
+    frequency_axis_request: np.ndarray = frequency_range_min * np.flip(scale)
+
+    return 1.0 / (frequency_axis_request * dt)
+
+
+def calculate_cone_of_influence(dt: float, frequency_axis: np.ndarray):
+    wave_scales = 1.0 / (frequency_axis * dt)
+    cone_of_influence: np.ndarray = np.ceil(np.sqrt(2) * wave_scales).astype(np.int64)
+    return cone_of_influence
+
+
+def get_y_ticks(
+    reduction_to_ticks: int, frequency_axis: np.ndarray, round: int
+) -> tuple[np.ndarray, np.ndarray]:
+    output_ticks = np.arange(
+        0,
+        frequency_axis.shape[0],
+        int(np.floor(frequency_axis.shape[0] / reduction_to_ticks)),
+    )
+    if round < 0:
+        output_freq = frequency_axis[output_ticks]
+    else:
+        output_freq = np.round(frequency_axis[output_ticks], round)
+    return output_ticks, output_freq
+
+
+def get_x_ticks(
+    reduction_to_ticks: int, dt: float, number_of_timesteps: int, round: int
+) -> tuple[np.ndarray, np.ndarray]:
+    time_axis = dt * np.arange(0, number_of_timesteps)
+    output_ticks = np.arange(
+        0, time_axis.shape[0], int(np.floor(time_axis.shape[0] / reduction_to_ticks))
+    )
+    if round < 0:
+        output_time_axis = time_axis[output_ticks]
+    else:
+        output_time_axis = np.round(time_axis[output_ticks], round)
+    return output_ticks, output_time_axis
+
+
+f_test: float = 50  # Hz
+number_of_test_samples: int = 1000
+
+# The wavelet we want to use
+mother = pywt.ContinuousWavelet("cmor1.5-1.0")
+
+# Parameters for the wavelet transform
+number_of_frequences: int = 25  # frequency bands
+frequency_range_min: float = 15  # Hz
+frequency_range_max: float = 200  # Hz
+dt: float = 1.0 / 1000  # sec
+
+t_test: np.ndarray = np.arange(0, number_of_test_samples) * dt
+test_data: np.ndarray = np.sin(2 * np.pi * f_test * t_test)
+
+wave_scales = calculate_wavelet_scale(
+    number_of_frequences=number_of_frequences,
+    frequency_range_min=frequency_range_min,
+    frequency_range_max=frequency_range_max,
+    dt=dt,
+)
+
+complex_spectrum, frequency_axis = pywt.cwt(
+    data=test_data, scales=wave_scales, wavelet=mother, sampling_period=dt
+)
+
+cone_of_influence = calculate_cone_of_influence(dt, frequency_axis)
+
+plt.imshow(abs(complex_spectrum) ** 2, cmap="hot", aspect="auto")
+plt.plot(cone_of_influence, np.arange(0, cone_of_influence.shape[0]), "g")
+plt.plot(
+    complex_spectrum.shape[1] - 1 - cone_of_influence,
+    np.arange(0, cone_of_influence.shape[0]),
+    "g",
+)
+plt.colorbar()
+
+y_ticks, y_labels = get_y_ticks(
+    reduction_to_ticks=10, frequency_axis=frequency_axis, round=1
+)
+
+x_ticks, x_labels = get_x_ticks(
+    reduction_to_ticks=10, dt=dt, number_of_timesteps=complex_spectrum.shape[1], round=2
+)
+
+plt.yticks(y_ticks, y_labels)
+plt.xticks(x_ticks, x_labels)
+
+plt.xlabel("Time [sec]")
+plt.ylabel("Frequency [Hz]")
+plt.show()
+```
+
+![figure 7](image7.png)