### Setup Python Path and Imports

Source: https://github.com/kpetras/wavespace/blob/main/docs/source/tutorials/CircLinCorr.md

Ensures the project root is on the Python path and imports necessary modules for WaveSpace analysis and plotting.

```python
import sys
import os
path = os.path.abspath(os.path.join(os.path.dirname(__file__), '..'))
sys.path.insert(0, path )
print(path)

from WaveSpace.PlottingHelpers import Plotting
from WaveSpace.Utils import HelperFuns as hf
from WaveSpace.Utils import ImportHelpers
from WaveSpace.WaveAnalysis import DistanceCorrelation
from WaveSpace.SpatialArrangement import SensorLayout as sensors

import time
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
import matplotlib.gridspec as gridspec
from matplotlib import colormaps
```

--------------------------------

### Setup WaveSpace Project Path

Source: https://github.com/kpetras/wavespace/blob/main/docs/source/tutorials/Sensor_layout.md

Add the project root directory to the Python path for working with source code. This is not necessary when the package is installed.

```python
import sys
import os
path = os.path.abspath(os.path.join(os.path.dirname(__file__), '..'))
sys.path.insert(0, path )
print(path)

from WaveSpace.Utils import ImportHelpers
from WaveSpace.SpatialArrangement import SensorLayout as sensors
from WaveSpace.PlottingHelpers import Plotting
from WaveSpace.Utils import WaveData as wd

import vtk
import numpy as np
import matplotlib.pyplot as plt
import pickle
```

--------------------------------

### Setup for WaveSpace Project

Source: https://github.com/kpetras/wavespace/blob/main/docs/source/tutorials/Create_WaveData.md

Ensures the project root is on the Python path for WaveSpace modules. Imports necessary libraries like sys, os, numpy, and WaveData.

```python
import sys
import os
path = os.path.abspath(os.path.join(os.path.dirname(__file__), '..'))
sys.path.insert(0, path )
print(path)

from WaveSpace.Utils import WaveData as wd
import numpy as np
```

--------------------------------

### Install WaveSpace Package

Source: https://github.com/kpetras/wavespace/blob/main/docs/index.md

Install the WaveSpace package using pip. Replace 'x.x.x' with the specific version number.

```console
pip install WaveSpace-x.x.x-py3-none-any.whl
```

--------------------------------

### Install WaveSpace using pip

Source: https://github.com/kpetras/wavespace/blob/main/README.md

Install the latest version of WaveSpace by downloading the wheel file and using pip in your terminal.

```bash
pip install WaveSpace-1.1.8-py3-none-any.whl
```

--------------------------------

### Setup Python Environment for WaveSpace

Source: https://github.com/kpetras/wavespace/blob/main/docs/source/tutorials/optical_flow.md

Ensures the project root is on the Python path and imports necessary libraries for WaveSpace analysis.

```python
import sys
import os
import time
path = os.path.abspath(os.path.join(os.path.dirname(__file__), '..'))
sys.path.insert(0, path )
print(path)

from WaveSpace.PlottingHelpers import Plotting
from WaveSpace.Utils import HelperFuns as hf
from WaveSpace.Utils import ImportHelpers
from WaveSpace.WaveAnalysis import OpticalFlow

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
import matplotlib.gridspec as gridspec
from matplotlib import colormaps
```

--------------------------------

### Setup Grid and Distance Matrix

Source: https://github.com/kpetras/wavespace/blob/main/docs/source/tutorials/CircLinCorr.md

Sets up a regular grid and distance matrix using channel positions from the loaded WaveData object. This assumes the data is already on a regular grid.

```python
sensors.regularGrid(waveData)
```

--------------------------------

### Setup WaveSpace Environment

Source: https://github.com/kpetras/wavespace/blob/main/docs/source/tutorials/2DFFT.md

Adds the project root to the Python path and imports necessary WaveSpace modules and libraries. Ensure the project root is on the Python path when working with source code.

```python
# Add the project root directory to the Python path when working with source code,
# not necessary when package is installed
import sys
import os
path = os.path.abspath(os.path.join(os.path.dirname(__file__), '..'))
sys.path.insert(0, path )
print(path)

from WaveSpace.Utils import ImportHelpers
from WaveSpace.PlottingHelpers import Plotting
from WaveSpace.WaveAnalysis import WaveAnalysis as wa

import numpy as np
import matplotlib.pyplot as plt
```

--------------------------------

### Empirical Mode Decomposition (EMD) Example

Source: https://github.com/kpetras/wavespace/blob/main/docs/source/tutorials/Frequency_Decomposition.md

Applies EMD to a signal when FFT-based methods are unsuitable. Note that EMD is significantly slower than Filter + Hilbert. This example truncates data for faster execution.

```python
# If we cannot expect the signal to be well behaved for FFT based approaches, we can use EMD
# note that this is A LOT slower than Filter + Hilbert
#We cut down the data to a small region to speed up the example
tempWaveData = copy.deepcopy(waveData)
tempWaveData.DataBuckets["SimulatedData"].set_data(waveData.get_data("SimulatedData")[0:2,10:14,10:14,:], "trl_posx_posy_time")

emd.EMD(tempWaveData,
        siftType = 'masked_sift',
        nIMFs=7,
        dataBucketName="SimulatedData",
        noiseVar = 0.05,
        n_noiseChans = 10,
        ndir=None,
        stp_crit ='stop',
        sd=0.075,
        sd2=0.75,
        tol=0.075,
        stp_cnt=2)

#plot imfs
TrialOfInterest = 0
SelectedChannel = (1,1)
IMFOfInterest = 4
dataInds = (slice(None), TrialOfInterest, SelectedChannel[0], SelectedChannel[1])
Plotting.plot_imfs(tempWaveData, dataInds, IMFOfInterest)
```

--------------------------------

### Add Project Root to Python Path

Source: https://github.com/kpetras/wavespace/blob/main/docs/source/tutorials/Wave_Activity.md

Add the project root directory to the Python path. This is necessary when working with source code and WaveSpace is not yet installed as a package.

```python
import sys
import os
path = os.path.abspath(os.path.join(os.path.dirname(__file__), '..'))
sys.path.insert(0, path )
print(path)
```

--------------------------------

### Wavelet Convolution (Time-Domain)

Source: https://github.com/kpetras/wavespace/blob/main/docs/source/tutorials/Frequency_Decomposition.md

Performs wavelet convolution using the Morlet wavelet. Frequencies should typically be logarithmically spaced within the range of interest. This example uses specific frequencies for comparison.

```python
# 4.1 Time-Domain
# frequencies are the centre frequencies of the wavelets and would normally be logarithmically spaced within your frequency range of interest.
# for comparison with other methods we use the same two frequencies as above
frequencies = [10,17]
Morlet.wavelet_convolution(waveData, dataBucketName="SimulatedData", n_cycles=3, frequencies=frequencies)

#plot.
complexData = waveData.DataBuckets["complexData"].get_data()[0,0,18,19,:] #dimord is freq_trl_posx_posy_time
fig, axs = plt.subplots(2, 1, figsize=(10, 6), sharex=True)
fig.suptitle(f"Wavelet Convolution")
# real part and envelope
axs[0].plot(waveData.get_time(), np.real(complexData), label='Real part')
axs[0].plot(waveData.get_time(), np.abs(complexData), label='Envelope', linestyle='--')
axs[0].set_ylabel('Amplitude')
axs[0].set_title('Real part and Envelope of Analytic Signal')
axs[0].legend()
axs[0].grid()
# phase
axs[1].plot(waveData.get_time(), np.angle(complexData), color='tab:orange')
axs[1].set_ylabel('Phase (radians)')
axs[1].set_xlabel('Time (s)')
axs[1].set_title('Phase of Analytic Signal')
axs[1].grid()
plt.tight_layout()
plt.show()
```

--------------------------------

### Import WaveSpace and Utility Packages

Source: https://github.com/kpetras/wavespace/blob/main/docs/source/tutorials/Wave_Activity.md

Import necessary modules from WaveSpace for wave analysis and utility functions, along with numpy and matplotlib for data handling and visualization.

```python
from WaveSpace.WaveAnalysis import WaveActivity as wa
from WaveSpace.Utils import ImportHelpers
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
```

--------------------------------

### WaveData Object Creation and Population

Source: https://github.com/kpetras/wavespace/blob/main/docs/source/tutorials/Create_WaveData.md

Creates an empty WaveData object, sets its sample rate, and adds a DataBucket containing the generated data and metadata. Populates optional fields like channel positions, trial info, and distance matrix.

```python
# create empty waveData object
waveData = wd.WaveData()
waveData.set_sample_rate(sampleRate)
# Create Databucket:
dataBucketName = "fakeEEG"
fakeEEG = wd.DataBucket(data, dataBucketName, dimord=dimord, chanNames=channelNames, time=time, unit=unit)
waveData.add_data_bucket(fakeEEG)

waveData.set_channel_positions(chanpos)
waveData.set_trialInfo(trialInfo)
waveData.set_distMat(distMat)
```

--------------------------------

### Run Unit Tests Locally

Source: https://github.com/kpetras/wavespace/blob/main/CONTRIBUTING.md

Execute the project's unit tests using Python's built-in unittest framework. Ensure all tests pass before submitting changes.

```bash
python -m unittest discover UnitTest
```

--------------------------------

### Load Simulated Wave Data

Source: https://github.com/kpetras/wavespace/blob/main/docs/source/tutorials/optical_flow.md

Loads simulated wave data using the ImportHelpers utility from a specified output path.

```python
# Load some simulated data
dataPath  = os.path.join(path, "Examples/ExampleData/Output")
waveData = ImportHelpers.load_wavedata_object(dataPath + "/complexData")
```

--------------------------------

### Map and Plot Wave Motifs

Source: https://github.com/kpetras/wavespace/blob/main/docs/source/tutorials/optical_flow.md

Creates a motif map using pcolormesh and plots individual motifs as quiver plots with colored spines.

```python
conds = waveData.get_trialInfo()
uniques = np.unique(conds)
trial_dict = {}
for trial_idx, condition in enumerate(conds):
    if condition not in trial_dict:
        trial_dict[condition] = []
    trial_dict[condition].append(trial_idx)

motifMap = np.full((len(conds),len(waveData.get_time())), -1)
for ind, motif  in enumerate(motifs):
    trial_frames_list = motif['trial_frames']
    for trial_frame in trial_frames_list:
        (trial, (start_timepoint, end_timepoint)) = trial_frame
        motifMap[trial, start_timepoint:end_timepoint] = ind

cmap = mcolors.ListedColormap(['grey', "#8F43D1", "#c50069",'#d67258', '#416ae4', '#378b8c', "#0f3200" ,'#a05195', "#4e2f13", "#3900AB","#b3ff00", "#ff0015", "#0d15ff"])
bounds = [-1, 0, 1, 2, 3, 4, 5, 7, 8, 9, 10, 11, 12]
norm = mcolors.BoundaryNorm(bounds, cmap.N)
im =  plt.pcolormesh(motifMap, cmap=cmap, norm=norm)
plt.colorbar()

# One column per condition
fig, axs = plt.subplots(1, len(motifs[0:6]), figsize=(12, 6), gridspec_kw={'wspace': 0.3})
for motifInd, motif  in enumerate(motifs[0:6]):
    # Quiver plot
    axs[motifInd].quiver(-np.real(motif['average']), -np.imag(motif['average']), color='black')
    axs[motifInd].set_facecolor('white')
    axs[motifInd].set_aspect('equal')
    for spine in axs[motifInd].spines.values():
        spine.set_edgecolor(cmap(motifInd+1))
        spine.set_linewidth(2)
```

--------------------------------

### Load Simulated Data

Source: https://github.com/kpetras/wavespace/blob/main/docs/source/tutorials/CircLinCorr.md

Loads simulated data from a specified path. Ensure the 'os' and 'ImportHelpers' modules are available.

```python
dataPath  = os.path.join(path, "Examples/ExampleData/Output")
waveData = ImportHelpers.load_wavedata_object(dataPath + "/complexData")
```

--------------------------------

### Create Surface and Interpolate Sensor Data

Source: https://github.com/kpetras/wavespace/blob/main/docs/source/tutorials/Sensor_layout.md

This snippet demonstrates creating a surface from sensor points, calculating distances along the surface, and interpolating sensor positions to a regular grid. It's useful for spatial analysis and visualization of sensor data.

```python
#interpolate the whole thing back into a regular grid
chanInds=True
surface, polySurface = sensors.create_surface_from_points(waveData,
                                                            type='channels',
                                                            num_points=1000)

sensors.distance_along_surface(waveData, surface, tolerance=0.1, get_extent = chanInds, plotting= True)
sensors.distmat_to_2d_coordinates_Isomap(waveData) #can also use MDS here

grid_x, grid_y, mask =sensors.interpolate_pos_to_grid(
    waveData,
    dataBucketName = "EEGLayout",
    numGridBins=18,
    return_mask = True,
    mask_stretching = True)

distMat = sensors.regularGrid(waveData)

original_data_bucket = "EEGLayout"
interpolated_data_bucket = "EEGLayoutInterpolated"
OrigInd  = (0,slice(None),202)
InterpInd =(0,slice(None),slice(None),202)
fig = Plotting.plot_interpolated_data(waveData, original_data_bucket, interpolated_data_bucket,
                                        grid_x, grid_y, OrigInd, InterpInd,  type='')

```

--------------------------------

### Placeholder Data Generation

Source: https://github.com/kpetras/wavespace/blob/main/docs/source/tutorials/Create_WaveData.md

Generates placeholder data for trials, channels, and timepoints, along with essential metadata like time vector, sample rate, dimension order, and channel names. Includes optional fields for unit, channel positions, trial info, and distance matrix.

```python
# Placeholder data
nTrials = 3
nChannels = 64
nTimepoints = 500

data = np.random.rand(nTrials, nChannels, nTimepoints)  # 3 trials, 64 channels, 500 time points

# Obligatory: Create time-vector
duration = 1 # in seconds
start = -0.5
time = np.linspace(start, start+duration, nTimepoints)

# Obligatory: Sample rate
sampleRate = nTimepoints / duration

# Obligatory: dimension order
dimord = "trl_chan_time"  # trials x channels x time

# Obligatory: channel Names
channelNames = [f"Channel {i+1}" for i in range(nChannels)]

# optional: unit
unit = "uV"

# optional (but highly recommended): channel Positions
chanpos = np.random.rand(nChannels, 3)  # Random 3D positions for each channel

# optional: trial info
trialInfo = ['condition_1', 'condition_2', 'condition_3']

# optional: distance matrix (can also be calculated from positions)
distMat   = np.random.rand(nChannels, nChannels)
```

--------------------------------

### Run 2D FFT Analysis

Source: https://github.com/kpetras/wavespace/blob/main/docs/source/tutorials/2DFFT.md

Performs a 2D Fast Fourier Transform on the provided WaveSpace data. It restricts the analysis to a specified frequency range and stores the results. Ensure the 'SimulatedData' bucket exists in the WaveSpace object.

```python
# restrict to (temp) frequencies between lower and upper bound:
lowerBound = 2
upperBound = 40

wa.FFT_2D(waveData, sourcePointsDiagonal, lowerBound, upperBound, DataBucketName="SimulatedData")

result = waveData.get_data("Result")
```

--------------------------------

### Load Wave Data Object

Source: https://github.com/kpetras/wavespace/blob/main/docs/source/tutorials/Wave_Activity.md

Load a pre-existing wave data object from a specified file path. This function is part of the ImportHelpers utility.

```python
waveData = ImportHelpers.load_wavedata_object("ExampleData/Output/complexData")
```

--------------------------------

### Load Simulated WaveData

Source: https://github.com/kpetras/wavespace/blob/main/docs/source/tutorials/Frequency_Decomposition.md

Loads simulated neural data from a specified path using the ImportHelpers.load_wavedata_object function. Ensure the 'path' variable is defined and points to the correct directory.

```python
# Load some simulated data
dataPath  = os.path.join(path, "Examples/ExampleData/Output")
waveData = ImportHelpers.load_wavedata_object(dataPath + "/SimulatedData")
```

--------------------------------

### Simulate Rotating Wave Data

Source: https://github.com/kpetras/wavespace/blob/main/docs/source/tutorials/Simulate_WaveData.md

Generates rotating wave data with specified parameters and mixes it with spatial pink noise. Ensure `SimulationFuns` is imported and necessary variables like `nTrials`, `MatrixSize`, `SampleRate`, `SimDuration`, `TemporalFrequency`, `SpatialFrequency`, `WaveOnset`, `WaveDuration`, `CenterX`, and `CenterY` are defined.

```python
Type = "RotatingWave"
WaveDirection = [1, 1, -1, -1]  # 1 = CW, -1 = CCW

spiralWave = SimulationFuns.simulate_signal(
    Type, nTrials, MatrixSize, SampleRate, SimDuration,
    SimLayout="grid",
    TemporalFrequency=TemporalFrequency,
    SpatialFrequency=SpatialFrequency,
    WaveDirection=WaveDirection,
    WaveOnset=WaveOnset,
    WaveDuration=WaveDuration,
    CenterX=CenterX,
    CenterY=CenterY
)

spiralNoise = SimulationFuns.simulate_signal("SpatialPinkNoise", nTrials, MatrixSize, SampleRate, SimDuration, SimLayout="grid")
spiralWaveData = SimulationFuns.SNRMix(spiralWave, spiralNoise, 0.8, SimLayout="grid")
```

--------------------------------

### Simulate Local Oscillation Data

Source: https://github.com/kpetras/wavespace/blob/main/docs/source/tutorials/Simulate_WaveData.md

Generates local oscillation data with random oscillator phases and mixes it with spatial pink noise. Ensure `SimulationFuns` is imported and necessary variables like `nTrials`, `MatrixSize`, `SampleRate`, `SimDuration`, `WaveOnset`, `WaveDuration`, `TemporalFrequency` are defined.

```python
Type = "LocalOscillation"

localOscillators = SimulationFuns.simulate_signal(
    Type=Type, ntrials=nTrials,
    MatrixSize=MatrixSize, SampleRate=SampleRate, SimDuration=SimDuration,
    SimLayout="grid",
    WaveOnset=WaveOnset,
    WaveDuration=WaveDuration,
    OscillatoryPhase="Random",
    TemporalFrequency=TemporalFrequency,
    OscillatorProportion=0.4
)

oscillatorNoise = SimulationFuns.simulate_signal("SpatialPinkNoise", nTrials, MatrixSize, SampleRate, SimDuration, SimLayout="grid")
oscillatorWaveData = SimulationFuns.SNRMix(localOscillators, oscillatorNoise, 0.8, SimLayout="grid")
```

--------------------------------

### Animate Simulated Wave Data

Source: https://github.com/kpetras/wavespace/blob/main/docs/source/tutorials/Simulate_WaveData.md

Animates simulated wave data across a grid and saves each animation as a GIF. Requires the waveData object, a data bucket name, and a data index.

```python
for trl in range(waveData.get_data("SimulatedData").shape[0]):
    ani = Plotting.animate_grid_data(
        waveData,
        DataBucketName="SimulatedData",
        dataInd=trl,
        probepositions=[(0,15), (5,15), (10,15), (15,15), (19,15), (19,15)]
    )
    plot_file = output_path + f"/SimulationAnimation_{trl}.gif"
    ani.save(plot_file)
```

--------------------------------

### Load WaveData Object

Source: https://github.com/kpetras/wavespace/blob/main/docs/source/tutorials/Sensor_layout.md

Load a WaveData object from a specified file path. Ensure the 'saveFolder' variable points to the correct directory containing the data.

```python
#load data from file:
saveFolder = "ExampleData/Output/"
waveData = ImportHelpers.load_wavedata_object(saveFolder + "SimulatedData")
```

--------------------------------

### Combine Simulated Data

Source: https://github.com/kpetras/wavespace/blob/main/docs/source/tutorials/Simulate_WaveData.md

Combines multiple simulated datasets (plane wave, target wave, spiral wave, oscillator wave) into a single dataset with specified conditions. Requires `os` module for directory creation and saving. Ensure `SimulationFuns` is imported and `Conditions` is defined.

```python
simCondList = []
for item in Conditions:
    simCondList.append(item)
    simCondList.append(item)

waveData = SimulationFuns.combine_SimData(
    [planeWaveData, targetWaveData, spiralWaveData, oscillatorWaveData],
    dimension='trl',
    SimCondList=simCondList
)

output_path = "ExampleData/Output"
if not os.path.exists(output_path):
    os.makedirs(output_path)
waveData.save_to_file(os.path.join(output_path, "SimulatedData"))
```

--------------------------------

### Load Simulated Data

Source: https://github.com/kpetras/wavespace/blob/main/docs/source/tutorials/2DFFT.md

Loads a WaveSpace wavedata object from a specified folder. This function is used to retrieve pre-existing simulated data for analysis.

```python
saveFolder = "ExampleData/Output/"
waveData = ImportHelpers.load_wavedata_object(saveFolder + "SimulatedData")
```

--------------------------------

### Detect Wave Motifs

Source: https://github.com/kpetras/wavespace/blob/main/docs/source/tutorials/optical_flow.md

Identifies wave motifs based on specified thresholds for magnitude, pixel density, and duration. Requires frequency and data bucket information.

```python
foi = 10
cycleLength = waveData.get_sample_rate()/ foi
freqInd = 0
motifs = hf.find_wave_motifs(waveData,
                        dataBucketName="UV",
                        threshold = 0.8,
                        nTimepointsEdge=cycleLength,
                        mergeThreshold = 0.8,
                        minFrames=cycleLength,
                        pixelThreshold = 0.6,
                        magnitudeThreshold=.1,
                        dataInds = (freqInd, slice(None), slice(None), slice(None), slice(None)),
                        Mask = False)
```

--------------------------------

### Define Simulation Conditions

Source: https://github.com/kpetras/wavespace/blob/main/docs/source/tutorials/Simulate_WaveData.md

Define a list of strings representing different types of wave simulations to be performed. These include plane waves, target waves, rotating waves, and local oscillations.

```python
Conditions = [
  "PlaneWave_45", "PlaneWave_135",
  "TargetWave_in", "TargetWave_out",
  "RotatingWave_CW", "RotatingWave_CCW",
  "LocalOscillationRandom", "LocalOscillationSynched"
]
```

--------------------------------

### Extract Basis Functions from All Data

Source: https://github.com/kpetras/wavespace/blob/main/docs/source/tutorials/Wave_Activity.md

Calculate spatial basis functions using all available data at once for a more global view of wave activity. The results can be interpreted as linear combinations of all included waves.

```python
nBases = 5
dataInd = None
wa.find_wave_activity(waveData, dataBucketName="complexData", dataInd=dataInd, nBases=nBases)

bases = waveData.get_data('Bases')
```

```python
fig, axs = plt.subplots(1, nBases, figsize=(nBases*6, 6))
if nBases == 1:
    axs = [axs]
for b in range(nBases):
    im = axs[b].imshow(
        np.angle(bases[:, :, b]),
        cmap='hsv',
        vmin=-np.pi,
        vmax=np.pi,
        origin='lower',
        aspect='auto'
    )
    axs[b].set_title(f'wave map {b+1}')
    axs[b].set_xlabel('posy')
    axs[b].set_ylabel('posx')
    fig.colorbar(im, ax=axs[b], fraction=0.046, pad=0.04, label='Phase (rad)')

plt.tight_layout()
plt.show()
```

--------------------------------

### Simulate Target Wave Data

Source: https://github.com/kpetras/wavespace/blob/main/docs/source/tutorials/Simulate_WaveData.md

Generates target wave data with specified center coordinates and wave properties. Noise is added using `SNRMix`. Wave direction can be inward (-1) or outward (1).

```python
Type = "TargetWave"
CenterX, CenterY = 2, 2

SpatialFrequency = [0.6, 0.6, 0.6, 0.6]
TemporalFrequency = [10, 10, 10, 10]
WaveDirection = [-1, -1, 1, 1]  # -1 = inward, 1 = outward

WaveOnset = 300
WaveDuration = 1000

targetWave = SimulationFuns.simulate_signal(
    Type, nTrials, MatrixSize, SampleRate, SimDuration,
    SimLayout="grid",
    TemporalFrequency=TemporalFrequency,
    SpatialFrequency=SpatialFrequency,
    WaveDirection=WaveDirection,
    WaveOnset=WaveOnset,
    WaveDuration=WaveDuration,
    CenterX=CenterX,
    CenterY=CenterY
)

targetNoise = SimulationFuns.simulate_signal("SpatialPinkNoise", nTrials, MatrixSize, SampleRate, SimDuration, SimLayout="grid")
targetWaveData = SimulationFuns.SNRMix(targetWave, targetNoise, 0.8, SimLayout="grid")
```

--------------------------------

### Create a New Git Branch

Source: https://github.com/kpetras/wavespace/blob/main/CONTRIBUTING.md

Use this command to create a new branch for your feature or bugfix before making changes.

```bash
git checkout -b feature/your-feature-name
```

--------------------------------

### Apply Filter-Hilbert Approach for Analytic Signal

Source: https://github.com/kpetras/wavespace/blob/main/docs/source/tutorials/Frequency_Decomposition.md

Applies a narrow-band filter around frequencies of interest (10Hz and 17Hz) and then computes the analytic signal using the Hilbert transform. This method is suitable when a narrowband oscillation at a known frequency is expected. The code then plots the real part, envelope, and phase of the analytic signal for one of the extracted complex time series.

```python
# we filter the data narrowly around our frequency of interest (10Hz) and then apply the Hilbert transform to get the analytic signal.
# Note that this only makes sense if we **already know** that there is a narrowband oscillation at the frequency of interest.
# To demonstrate this, we will filter the data at 17Hz as well.
for freqInd, freq in enumerate([10, 17]):
    filt.filter_narrowband(waveData, dataBucketName = "SimulatedData", LowCutOff=freq-1, HighCutOff=freq+1, type = "FIR", order=100, causal=False)
    waveData.DataBuckets[str(freq)] =  waveData.DataBuckets.pop("NBFiltered") #Rename

temp = np.stack((waveData.DataBuckets["10"].get_data(), waveData.DataBuckets["17"].get_data()),axis=0) #Stack into single filtered databucket
waveData.add_data_bucket(wd.DataBucket(temp, "NBFiltered", "freq_trl_posx_posy_time", sampleRate=waveData.get_sample_rate(), chanNames=waveData.get_channel_names()))
#remove the individual filtered data
waveData.delete_data_bucket("10")
waveData.delete_data_bucket("17")
# get complex timeseries
hilb.apply_hilbert(waveData, dataBucketName = "NBFiltered")

#plot. Try both frequencies and see for which one the phase makes sense
complexData = waveData.DataBuckets["complexData"].get_data()[0,0,18,19,:] #dimord is freq_trl_posx_posy_time
fig, axs = plt.subplots(2, 1, figsize=(10, 6), sharex=True)
# real part and envelope
axs[0].plot(waveData.get_time(), np.real(complexData), label='Real part')
axs[0].plot(waveData.get_time(), np.abs(complexData), label='Envelope', linestyle='--')
axs[0].set_ylabel('Amplitude')
axs[0].set_title('Real part and Envelope of Analytic Signal')
axs[0].legend()
axs[0].grid()
# phase
axs[1].plot(waveData.get_time(), np.angle(complexData), color='tab:orange')
axs[1].set_ylabel('Phase (radians)')
axs[1].set_xlabel('Time (s)')
axs[1].set_title('Phase of Analytic Signal')
axs[1].grid()
plt.tight_layout()
plt.show()

waveData.save_to_file(os.path.join(dataPath, "complexData"))
```

--------------------------------

### Assign Spatial Layout to Channels and Interpolate Data

Source: https://github.com/kpetras/wavespace/blob/main/docs/source/tutorials/Sensor_layout.md

Assigns a spatial layout to channels by normalizing and mapping 3D positions to a regular grid. It then interpolates the original data onto this grid and subsamples it at the channel positions. This is useful for demonstrating interpolation from 3D positions to a regular grid.

```python
# pick some channels and assign a spatial layout to them (this makes no sense at all for real data and
#is only to demonstrate the interpolation to a reagular grid from 3d positions)
gridshape = waveData.get_data('SimulatedData').shape[1:3]
chanpos = np.load( saveFolder + 'exampleChanpos.npy')
waveData.set_channel_positions(chanpos)

x_normalized = (chanpos[:, 0] - np.min(chanpos[:, 0])) / (np.max(chanpos[:, 0]) - np.min(chanpos[:, 0]))
y_normalized = (chanpos[:, 1] - np.min(chanpos[:, 1])) / (np.max(chanpos[:, 1]) - np.min(chanpos[:, 1]))
z_normalized = (chanpos[:, 2] - np.min(chanpos[:, 2])) / (np.max(chanpos[:, 2]) - np.min(chanpos[:, 2]))

x_indices = np.round(x_normalized * (gridshape[0] - 1)).astype(int)
y_indices = np.round(y_normalized * (gridshape[1] - 1)).astype(int)
z_indices = np.round(z_normalized * (gridshape[1] - 1)).astype(int)

original_data = waveData.get_data('SimulatedData')
grid_data = np.repeat(original_data[:, :, :, np.newaxis, :], gridshape[1], axis=3)

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(x_indices, y_indices, z_indices)
X, Y, Z = np.meshgrid(np.arange(gridshape[0]), np.arange(gridshape[1]), np.arange(gridshape[1]))
ax.scatter(X, Y, Z, alpha=0.1)
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
plt.show()
#now we simply subsample the original grid at the channel positions
newdata = grid_data[:, y_indices, x_indices, z_indices, :]
dataBucket = wd.DataBucket(newdata, "EEGLayout",'trl_chan_time', [] )
waveData.add_data_bucket(dataBucket)
```

--------------------------------

### Visualize 2D FFT Results

Source: https://github.com/kpetras/wavespace/blob/main/docs/source/tutorials/2DFFT.md

Visualizes the results of the 2D FFT analysis, including channel activity over time, FFT magnitude plots, and comparisons of 'Max Along Power' and 'Max Reverse Power'. This snippet iterates through different experimental conditions to display relevant plots.

```python
# Get the number of trials
n_trials = waveData.get_data("SimulatedData").shape[0]
trialInfo = waveData.get_trialInfo()
conditions = np.unique(trialInfo)

for condition in conditions:
    indices = [i for i, cond in enumerate(trialInfo) if cond == condition]
    logRatios=np.mean(np.log(result["Max Along Power"][indices]/result["Max Reverse Power"][indices]))
    newlineseries = np.zeros((len(sourcePointsDiagonal),waveData.get_data("SimulatedData").shape[3]))

    for ind, position in enumerate(sourcePointsDiagonal):
        newlineseries[ind] = np.mean(waveData.get_data("SimulatedData")[indices, position[0], position[1], :], axis=0)
    plt.figure()
    plt.imshow(newlineseries, aspect=4)
    plt.title(f"Channels over time (Condition {condition})")
    plt.show()

    plot = Plotting.plotfft_zoomed(np.mean(waveData.get_data("FFT_ABS")[indices,:,:],axis=0), waveData.get_sample_rate(), -20, 20, "fft abs", scale='log')
    plot.show()

    x_labels = np.arange(1)
    plt.figure()
    plt.bar(x_labels, [np.mean(result["Max Along Power"][indices])], color='b', width = 0.25 )
    plt.bar(x_labels + 0.25, [np.mean(result["Max Reverse Power"][indices])] , color='r', width = 0.25 )
    plt.legend(labels=["Along", "Reverse"])
    plt.xticks(x_labels + 0.125, ["0 degree"])
    plt.title(f"Max Power (Condition {condition})")
    plt.show()
```

--------------------------------

### Broadband Filtering and Generalized Phase Analysis

Source: https://github.com/kpetras/wavespace/blob/main/docs/source/tutorials/Frequency_Decomposition.md

Applies broadband filtering to wave data and then computes the generalized phase. Includes plotting of the real part, envelope, and phase of the analytic signal.

```python
lowerCutOff = 1
higherCutOff = 40
filt.filter_broadband(waveData, "SimulatedData", lowerCutOff, higherCutOff, 5)
GenPhase.generalized_phase(waveData, "BBFiltered")
#plot
complexSignal = waveData.DataBuckets["complexData"].get_data()[0,0,0,:] #dimord is freq_trl_posx_posy_time
origSignal = waveData.DataBuckets["SimulatedData"].get_data()[0,0,0,:]

fig, axs = plt.subplots(2, 1, figsize=(10, 6), sharex=True)
fig.suptitle(f"Generalized Phase")
# real part and envelope
axs[0].plot(waveData.get_time(), np.real(complexSignal), label='Real part')
axs[0].plot(waveData.get_time(), np.abs(complexSignal), label='Envelope', linestyle='--')
axs[0].set_ylabel('Amplitude')
axs[0].set_title('Real part and Envelope of Analytic Signal')
axs[0].legend()
axs[0].grid()
# phase
axs[1].plot(waveData.get_time(), np.angle(complexSignal), color='tab:orange')
axs[1].set_ylabel('Phase (radians)')
axs[1].set_xlabel('Time (s)')
axs[1].set_title('Phase of Analytic Signal')
axs[1].grid()
plt.tight_layout()
plt.show()

#alternative plot closer to the figure shown in # https://github.com/mullerlab/generalized-phase
time = waveData.get_time()
xw = np.real(origSignal)
xgp = complexSignal
phase = np.angle(xgp)
fig = plt.figure(figsize=(12.5, 4.2))
fig.suptitle(f"Generalized phase")
ax1 = fig.add_axes([0.08, 0.15, 0.7, 0.75])
ax1.plot(time, xw, linewidth=4, color='k', label='wideband signal')
# Colored phase line
points = np.array([time, np.real(xgp)]).T.reshape(-1, 1, 2)
segments = np.concatenate([points[:-1], points[1:]], axis=1)
norm = plt.Normalize(-np.pi, np.pi)
lc = LineCollection(segments, cmap='hsv', norm=norm)
lc.set_array(phase)
lc.set_linewidth(5)
ax1.add_collection(lc)
# Normal axes
ax1.set_xlim([time[0], time[-1]])
ax1.set_xlabel('Time (s)')
ax1.set_ylabel('Amplitude (a.u.)')
ax1.spines['top'].set_visible(False)
ax1.spines['right'].set_visible(False)
ax2 = fig.add_axes([0.1116, 0.6976, 0.0884, 0.2000], polar=True)
theta = np.linspace(-np.pi, np.pi, 100)
for i in range(len(theta)-1):
    ax2.plot(theta[i:i+2], [1, 1], color=plt.cm.hsv(norm(theta[i])), linewidth=6)
ax2.set_yticklabels([])
ax2.set_xticklabels([])
ax2.set_axis_off()
plt.show()
```

--------------------------------

### Explore Basis Functions for Data Subset

Source: https://github.com/kpetras/wavespace/blob/main/docs/source/tutorials/Wave_Activity.md

Extract and visualize spatial basis functions for a specific subset of trials and a particular frequency. This involves specifying data indices and the number of bases to find.

```python
nBases = 3
dataInd = (slice(0,1), slice(10,12), slice(None), slice(None), slice(None))
wa.find_wave_activity(waveData, dataBucketName="complexData", dataInd=dataInd, nBases=nBases)

bases = waveData.get_data('Bases')
```

```python
fig, axs = plt.subplots(1, nBases, figsize=(nBases*6, 6))
if nBases == 1:
    axs = [axs]
for b in range(nBases):
    im = axs[b].imshow(
        np.angle(bases[:, :, b]),
        cmap='hsv',
        vmin=-np.pi,
        vmax=np.pi,
        origin='lower',
        aspect='auto'
    )
    axs[b].set_title(f'wave map {b+1}')
    axs[b].set_xlabel('posy')
    axs[b].set_ylabel('posx')
    fig.colorbar(im, ax=axs[b], fraction=0.046, pad=0.04, label='Phase (rad)')

plt.tight_layout()
plt.show()
```

--------------------------------

### Calculate and Plot Power Spectral Density (PSD)

Source: https://github.com/kpetras/wavespace/blob/main/docs/source/tutorials/Frequency_Decomposition.md

Calculates and plots the Power Spectral Density (PSD) for each trial condition in the waveData. This helps confirm the dominant frequencies in the simulated data. The plot displays power on a logarithmic scale and highlights a specific frequency (10Hz) with a vertical line.

```python
#waves were simulated at 10Hz. We can confirm that by plotting the PSD (for each trial-type)
trialInfo = waveData.get_trialInfo() #this contains the condition name for each trial
unique_conds = np.unique(trialInfo)
for cond in unique_conds:
    trial_indices = [i for i, info in enumerate(waveData.get_trialInfo()) if info == cond]
    f, psd = welch(waveData.get_data("SimulatedData")[trial_indices], fs=waveData.get_sample_rate(), nperseg=256)
    #average over trials and grid positions (channels)
    psd = np.mean(psd, axis=(0,1,2))
    plt.semilogy(f, psd)
    plt.title(f"Power Spectral Density - {cond}")
    plt.xlabel("Frequency (Hz)")
    plt.ylabel("Power/Frequency (dB/Hz)")
    plt.axvline(10, color='red', linestyle=':', linewidth=2)
    plt.xlim(0, 50)
    plt.grid()
    plt.show()
```

--------------------------------

### Visualize Optical Flow

Source: https://github.com/kpetras/wavespace/blob/main/docs/source/tutorials/optical_flow.md

Generates and saves an animated GIF of the optical flow, plotting angle and normalizing vector length.

```python
trialToPlot = 4
waveData.set_active_dataBucket('UV')
ani = Plotting.plot_optical_flow(waveData,
                                UVBucketName = 'UV',
                                PlottingDataBucketName = 'complexData',
                                dataInds = (0, trialToPlot, slice(None), slice(None), slice(None)),
                                plotangle=True,
                                normVectorLength = True)
output_path = os.path.join(path, "Examples/ExampleData/Output/")
ani.save( output_path + 'OpticalFlowAfterFilter_Hilbert.gif')
```

--------------------------------

### Define Source Points for 2D FFT

Source: https://github.com/kpetras/wavespace/blob/main/docs/source/tutorials/2DFFT.md

Creates sample source points for 2D FFT analysis. Points are selected along the diagonal of the channel array. This function is used to specify regions of interest for the FFT.

```python
# Create 10 sample points along the diagonal of the channel array
gridSize = waveData.get_data("SimulatedData").shape[1]
nPoints = range(0,gridSize,2)
sourcePointsDiagonal = []
for i in nPoints:
    sourcePointsDiagonal.append([i,i])
```

--------------------------------

### Project 3D Coordinates to 2D Space

Source: https://github.com/kpetras/wavespace/blob/main/docs/source/tutorials/Sensor_layout.md

Projects 3D sensor coordinates to a 2D space while preserving distances using MDS. It also visualizes both the original 3D and the embedded 2D contact positions.

```python
#project 3D coordinates to 2D space, preserving distanes between them as good as possible
sensors.distmat_to_2d_coordinates_MDS(waveData)
#plot 3D and 2D contact positions:
fig = plt.figure()
ax = fig.add_subplot(projection='3d')
ax.scatter(waveData.get_channel_positions()[:,0], waveData.get_channel_positions()[:,1],waveData.get_channel_positions()[:,2])
plt.title('Contact position 3D ')
plt.figure()
plt.scatter(waveData.get_2d_coordinates()[:,0], waveData.get_2d_coordinates()[:,1])
plt.title('Contact position 2D embedding preserving inter-contact distances. Arbitrary units')
```

--------------------------------

### Simulate Plane Wave Data

Source: https://github.com/kpetras/wavespace/blob/main/docs/source/tutorials/Simulate_WaveData.md

Generates unidirectional plane waves mixed with spatial noise. Configure wave properties like frequency, direction, and onset. Noise is added using `SNRMix`.

```python
Type = "PlaneWave"
nTrials = 4
MatrixSize = 20
SampleRate = 500
SimDuration = 1.6

SpatialFrequency = [0.6, 0.6, 0.6, 0.6]
TemporalFrequency = [10, 10, 10, 10]
WaveDirection = [45, 45, 135, 135]
SimLayout = "grid"

WaveOnset = [300, 300, 300, 300]
WaveDuration = 1000

planeWave = SimulationFuns.simulate_signal(
    Type, nTrials, MatrixSize, SampleRate, SimDuration,
    SimLayout=SimLayout,
    TemporalFrequency=TemporalFrequency,
    SpatialFrequency=SpatialFrequency,
    WaveDirection=WaveDirection,
    WaveOnset=WaveOnset,
    WaveDuration=WaveDuration
)

planeWaveNoise = SimulationFuns.simulate_signal(
    Type="SpatialPinkNoise",
    ntrials=nTrials,
    MatrixSize=MatrixSize,
    SampleRate=SampleRate,
    SimLayout=SimLayout,
    SimDuration=SimDuration
)

SNR = 0.8
planeWaveData = SimulationFuns.SNRMix(planeWave, planeWaveNoise, SNR, SimLayout="grid")
```

--------------------------------

### Load and Plot Saved Correlation Data

Source: https://github.com/kpetras/wavespace/blob/main/docs/source/tutorials/CircLinCorr.md

Loads previously saved correlation data and visualizes the phase-distance correlation over time for different conditions. Assumes data is saved in a specific format.

```python
waveData = ImportHelpers.load_wavedata_object("ExampleData/Output/DistanceCorrelation")
pointRange = (20,20)
sourcePoints = []
for x in range(pointRange[0]):
    for y in range(pointRange[1]):
        sourcePoints.append((x,y))

phaseDistCorr= waveData.get_data("PhaseDistanceCorrelation")
conditions = waveData.get_trialInfo()[::2]

shape = waveData.get_data("complexData").shape
selectedTrial = 4

rho = np.zeros((8,20,20))
for condInd, condition in enumerate(conditions):
    for i, (x,y) in enumerate(sourcePoints):
        phaseDistCorrOverTime = phaseDistCorr.loc[(phaseDistCorr["trialInd"] == condInd*2) &
                                                (phaseDistCorr["sourcePointX"] == x) &
                                                (phaseDistCorr["sourcePointY"] == y)]
        rho[condInd,x,y] = np.mean(phaseDistCorrOverTime["rho"][300:500])
    fig, ax = plt.subplots(figsize=(8,6))
    im = ax.imshow(rho[condInd], origin="lower", )
    ax.set_title(condition)
    fig.colorbar(im, ax=ax)
    plt.show()
```

--------------------------------

### Calculate and Plot Full Grid Correlation

Source: https://github.com/kpetras/wavespace/blob/main/docs/source/tutorials/CircLinCorr.md

Calculates and plots the average phase-distance correlation for a specified time range and all source points. Requires waveData object and a data bucket name.

```python
# Only do if you have too much time on your hands:
# Calculate and plot average phase-distance correlation for 600 to 1000 ms for all points
pointRange = (20,20)
sourcePoints = []
for x in range(pointRange[0]):
    for y in range(pointRange[1]):
        sourcePoints.append((x,y))

DistanceCorrelation.calculate_distance_correlation(waveData, dataBucketName = "complexData", sourcePoints=sourcePoints, pixelSpacing=1)

output_path = os.path.join(path, "Examples/ExampleData/Output")
waveData.save_to_file(os.path.join(output_path, "DistanceCorrelation"))
```

--------------------------------

### Calculate and Plot Sensor Distance Matrix

Source: https://github.com/kpetras/wavespace/blob/main/docs/source/tutorials/Sensor_layout.md

Calculates the sensor-to-sensor distance matrix for a regular grid of contacts and visualizes it using a heatmap. The 'waveData' object must contain channel position information.

```python
#regular grid
#calculate sensor to sensor distance, where chanpos has x and y coordinates
sensors.regularGrid(waveData) #adds a distance matrix to the data object

#plot
plt.imshow(waveData.get_distMat(), origin= 'lower')
plt.colorbar()
plt.title('Contact-to-Contact distance')
plt.xlabel('Contact')
plt.ylabel('Contact')
```

--------------------------------

### Compute Optical Flow

Source: https://github.com/kpetras/wavespace/blob/main/docs/source/tutorials/optical_flow.md

Computes optical flow vectors (uv) for wave data. Parameters control Gaussian smoothing, flow type, and iteration.

```python
tStart = time.time()
print("OpticalFlow started")
OpticalFlow.create_uv(waveData,
        applyGaussianBlur=False,
        type = "angle",
        Sigma=1,
        alpha = 0.1,
        maxIter = 200,
        dataBucketName="complexData",
        is_phase = False)
print('optical flow took: ', time.time()-tStart)
```

--------------------------------

### Calculate Generalized Phase Distance Correlation

Source: https://github.com/kpetras/wavespace/blob/main/docs/source/tutorials/CircLinCorr.md

Calculates the generalized phase distance correlation for the WaveData object. Requires specifying a data bucket name, evaluation angle, and tolerance.

```python
DistanceCorrelation.calculate_distance_correlation_GP(waveData, dataBucketName = "complexData", evaluationAngle=np.pi, tolerance=0.2)
dataFrame = waveData.get_data("PhaseDistanceCorrelation")
```

--------------------------------

### Calculate Distance Correlation for Selected Source Points

Source: https://github.com/kpetras/wavespace/blob/main/docs/source/tutorials/CircLinCorr.md

Calculates distance correlation for a specified range of source points. Requires defining the point range, source points list, and pixel spacing.

```python
pointRange = range(0,20,2)
sourcePoints = []
for i in pointRange:
    sourcePoints.append((i,i))
DistanceCorrelation.calculate_distance_correlation(waveData, dataBucketName = "complexData", sourcePoints=sourcePoints, pixelSpacing=1)
```