#!/usr/bin/env python
# coding: utf-8

# # MUSIC (MUltiple SIgnal Classification) Notebook - SIMULATION

# ***This notebook is essentially a copy of the MUSIC.ipynb, except that it shows how to insert artificial MSTID's into the data processing stream for the purpose of evaluating the processing algorithms.***
# 
# <a href="#simulate">Click here to jump to the part where the data simulator is used.</a>

# [simulate]: The MUSIC (MUltiple SIgnal Classification) algorithm is a rather generic and well-know algorithm in signal process used for the detection of multiple waves within an array of sensors [Schmidt, 1986].  MUSIC was originally used in SuperDARN by Samson et al. [1990] and Bristow et al. [1994] for the estimation of the characteristics of Medium Scale Traveling Ionospheric Disturbances (MSTIDs) detected by SuperDARN Radars.  These papers have inspired this implementation of the MUSIC algorithm.
# 
# While this code has been written with SuperDARN detected MSTIDs in mind, it may be useful for working with other wave-like perturbations moving across the field of view of some geophysical instrument.
# 
# This notebook demonstrates the use of the DaViTPy MUSIC module for MSTID parameter estimation.
# 
# The DaViTPy MUSIC module is used by creating a musicArray object, which includes built-in methods for keeping track of all steps in the processing algorithm.  As the data is processed, all preceeding variants of the data are saved within the object and a history is created.  This makes it very easy to see what has been done with the data at any step along the way.
# 
# *Written by N.A. Frissell, 4 November 2013.*

# **References**
# 
# Bristow, W. A., R. A. Greenwald, and J. C. Samson (1994), Identification of high-latitude acoustic gravity wave sources using the Goose Bay HF Radar, J. Geophys. Res., 99(A1), 319–331, doi:10.1029/93JA01470.
# 
# Samson, J. C., R. A. Greenwald, J. M. Ruohoniemi, A. Frey, and K. B. Baker (1990), Goose Bay radar observations of Earth-reflected, atmospheric gravity waves in the high-latitude ionosphere, J. Geophys. Res., 95(A6), 7693–7709, doi:10.1029/JA095iA06p07693.
# 
# Schmidt, R.O., "Multiple emitter location and signal parameter estimation," Antennas and Propagation, IEEE Transactions on, vol.34, no.3, pp.276,280, Mar 1986, doi:10.1109/TAP.1986.1143830.

# In[1]:


#Import the modules we need.
get_ipython().run_line_magic('pylab', 'inline')
import datetime

from matplotlib import pyplot as plt 
import numpy as np

from davitpy import pydarn


# # MUSIC Processing Example

# ## Loading Data and Basic Plotting

# In[2]:


#Choose the radar and time of the data that we want to look at.
#Then, connect to the data source using the standard DaViTPy
#pydarn.sdio.radDataOpen() routine.

rad='wal'
sTime = datetime.datetime(2011,5,9,8)
eTime = datetime.datetime(2011,5,9,19)

myPtr = pydarn.sdio.radDataOpen(sTime,rad,eTime=eTime)


# In[3]:


# Now create the data object.

# By creating this object, data is taken from the SuperDARN database and rearranged into a numpy.array with
# the shape (Nr Times, Nr Beams, Nr Gates).  Missing data is stored as np.nan.

# By default, this command only loads in data where the ground scatter flag is True.

# The fovModel='GS' indicates to calculate the field-of-view coordinates in using the ground-scatter mapping
# formula.  This is standard when looking at MSTIDs using ground scatter.  See Bristow et al. [1994] for details.
dataObj     = pydarn.proc.music.musicArray(myPtr,fovModel='GS')


# In[4]:


# We can list all of the data sets stored in this object.  Right now there is only one data set, but
# as we go along, each step of processing will create a new data set.

# Each data set is simply stored as an attribute of the dataObj, but the names of all of the data sets can be
# accessed through the get_data_sets() method.
dataObj.get_data_sets()


# In[5]:


# Information about the data set is saved in it's metadata dictionary.  Some of the information in the metadata
# dictionary is used to control certain plotting and processing variables in this module.

dataObj.DS000_originalFit.printMetadata()


# In[6]:


# The music module also keeps track of each step of processing in using the data set's history attribute.

dataObj.DS000_originalFit.printHistory()


# In[7]:


# The 'active' attribute of dataObj is a reference to the most recently used data set of dataObj.  It is also
# the default data set used by all DaViTPy MUSIC processing and plotting routines.  You can see this command
# produces the same result as the cell above.

dataObj.active.printHistory()


# In[8]:


# We can create an RTI plot and a fan plot from the data we loaded.
# You can see that both range gate and geographic latitude are given on the y-axis, and the solar terminator
# is also shaded in.  The terminator does not go below approximately range gate 10.  This is because the ground
# scatter mapping formula is not defined for close-in range gates.

fig = pydarn.plotting.musicRTI(dataObj)


# In[9]:


# We can also make a fan plot.
plotTime = datetime.datetime(2011,5,9,14)
fig = pydarn.plotting.musicFan(dataObj,time=plotTime)


# ## Data Selection

# In[10]:


# We want to focus on just the data that contains the MSTID we are looking for, so we can apply some limits.

# proc.music.defineLimits() does not actually modify the data, it just marks the data by putting an an entry
# into the dataObj.active.metadata dictionary.  The limits will be applied later.  Right now, you can see
# what data will be eliminated by making a new RTI plot.

pydarn.proc.music.defineLimits(dataObj,gateLimits=[30,45])
fig = pydarn.plotting.musicRTI(dataObj)


# In[11]:


# We also to restrict the amount of time we process.  Before actually running the music algorithm, we
# will be filtering the data using a FIR filter that will eat up data at the beginning and end of the filter.
# We can calculate exactly how much time that will be if we know some of the filter characteristics and
# the time resolution of the data.  This can then be used to give us new start and end times.

# For now, let's say we are going to use a filter with 101 taps and a dataset with 120 s resolution.
numtaps = 101
timeres = 120

#Let's also say that we are interested in the MSTID feature between 1400 and 1600 UT.
sTime_of_interest = datetime.datetime(2011,5,9,14)
eTime_of_interest = datetime.datetime(2011,5,9,16)

#Now calculate the new start and end times...
new_times = pydarn.proc.music.filterTimes(sTime_of_interest, eTime_of_interest, timeres, numtaps)
pydarn.proc.music.defineLimits(dataObj,timeLimits=new_times)

fig = pydarn.plotting.musicRTI(dataObj)


# In[12]:


# Now we apply the limits and replot once more.
# Note that many of the processing routines will automatically call applyLimits() before they
# run the processing algorithm.

dataObj.active.applyLimits()
fig = pydarn.plotting.musicRTI(dataObj)


# In[13]:


# Note that a new data set was created when we applied the limits.
dataObj.get_data_sets()


# In[14]:


# Also note that the history of the new object was updated.
dataObj.active.printHistory()


# ## Data Processing

# ### Interpolation and Cartesian Coordinates

# In[15]:


# Before we feed into the MUSIC Algorithm, we don't want our data to have any gaps in time or space.
# Let's start by interpolating in space along the beams.
pydarn.proc.music.beamInterpolation(dataObj)


# In[16]:


fig = plt.figure(figsize=(20,10))
ax  = fig.add_subplot(121)
pydarn.plotting.musicPlot.musicFan(dataObj,plotZeros=True,dataSet='originalFit',axis=ax,time=plotTime)
ax  = fig.add_subplot(122)
pydarn.plotting.musicPlot.musicFan(dataObj,plotZeros=True,axis=ax,time=plotTime)


# In[17]:


# We also want to interpolate in time.  timeres=120 [seconds] was set in an earlier cell.
pydarn.proc.music.timeInterpolation(dataObj,timeRes=timeres)


# In[18]:


pydarn.plotting.musicPlot.timeSeriesMultiPlot(dataObj,dataSet='timeInterpolated',dataSet2='beamInterpolated')


# In[19]:


# We also want need to calculate a local cartesian grid for each cell we are going to use.
pydarn.proc.music.determineRelativePosition(dataObj)


# In[20]:


#The black cell marks the center of the array.
fig = pydarn.plotting.plotRelativeRanges(dataObj,time=plotTime)


# ### <a name='simulate' />Data Simulation

# In[21]:


# At this point, we can inject a simulated MSTID into the processing chain.  This way, we can see exactly what
# processing is doing to ideal data.  It is really easy to inject a default simulated wave.
pydarn.proc.music.simulator(dataObj)
fig = pydarn.plotting.musicRTI(dataObj)
fig = pydarn.plotting.musicFan(dataObj,time=plotTime)


# In[22]:


#Or, you can do something more complicated.  Here we define two custom MSTIDs and add some noise.
# Each MSTID is evaluated as a cosine and then summed together.  The cosine evaluated is::
# sig  = amp * np.cos(kx*xgrid + ky*ygrid - 2.*np.pi*f*t + phi) + dc

sigs = []
#           (amp,    kx,      ky,      f, phi, dcOffset)
sigs.append((  5,  -0.01,   0.010, 0.0007,   0,       5.))
sigs.append((  10, 0.022,  -0.023, 0.0004,   0,       5.))

#And some arbitrary white gaussian noise...
noiseFactor = 4


# In[23]:


pydarn.proc.music.simulator(dataObj,sigs=sigs,noiseFactor=noiseFactor)
fig = pydarn.plotting.musicRTI(dataObj)
fig = pydarn.plotting.musicFan(dataObj,time=plotTime)


# ### Filtering

# In[24]:


# Now filter the data.
# numtaps=101 was set in an above cell.  The cutoff_frequencies are in Hz.
filt = pydarn.proc.music.filter(dataObj, numtaps=numtaps, cutoff_low=0.0003, cutoff_high=0.0012)


# In[25]:


# At this point, the data has been filtered and saved to a new data set in dataObj.
# Before we look at the data, let's look at the transfer function and impulse response of the data.

fig = filt.plotTransferFunction(xmax=0.004)
fig = filt.plotImpulseResponse()


# In[26]:


# Let's look at the filtered RTI plot.  We should set autoScale to True since the magnitudes will be much
# lower than the original data.

fig = pydarn.plotting.musicRTI(dataObj,autoScale=True)


# In[27]:


# You can see that the filter already marked off the data that you shouldn't use.
# Just applyLimits() and replot.
dataObj.active.applyLimits()
fig = pydarn.plotting.musicRTI(dataObj,autoScale=True)


# In[28]:


# Look at a fan plot.
fig = pydarn.plotting.musicFan(dataObj,time=plotTime,autoScale=True)


# ### Spectral Analysis

# In[29]:


# We have now run all of the processing needed to feed the data into the spectral analysis and MUSIC.
# Let's print the history just to recap what we have done.
dataObj.active.printHistory()


# In[30]:


# First thing to do is to calculate the FFT of every cell...
pydarn.proc.music.calculateFFT(dataObj)


# In[31]:


# We can look at the spectrum of select cells...
pydarn.plotting.spectrumMultiPlot(dataObj,xlim=(-0.0025,0.0025))
pydarn.plotting.musicPlot.spectrumMultiPlot(dataObj,plotType='magnitude',xlim=(0,0.0025))
pydarn.plotting.musicPlot.spectrumMultiPlot(dataObj,plotType='phase',xlim=(0,0.0025))


# In[32]:


# We can also plot the full spectrum.  Here, every FFT bin contains 16 slices,
# showing the data for each of the 16 radar beams from left to right.
# Range gates are shown on the y-axis.
pydarn.plotting.plotFullSpectrum(dataObj,xlim=(0,0.00175))


# ### MUSIC Processing

# In[33]:


# Calculate the cross-spectral matrix Dlm.
pydarn.proc.music.calculateDlm(dataObj)


# In[34]:


# We can look at the Dlm matrix.  This usually isn't necessary for routine processing, but
# it is good to know where the numbers are going...
pydarn.plotting.plotDlm(dataObj)


# In[35]:


# Now, we finally run detect the horizontal wave numbers.
pydarn.proc.music.calculateKarr(dataObj)


# In[36]:


# We can then plot the final result.
pydarn.plotting.musicPlot.plotKarr(dataObj)


# In[37]:


# Again, we can see all of the history gone into processing this plot.
dataObj.active.printHistory()


# ### Feature/Signal Detection

# In[38]:


# We can use the image processing library in scikit-image to automatically find the peaks representing signals.
pydarn.proc.music.detectSignals(dataObj)
pydarn.plotting.musicPlot.plotKarr(dataObj)


# In[39]:


#All of the parameters connected with these signals is located in the sigDectect attribute.
dataObj.active.sigDetect.info


# In[40]:


# To get some insight on how the signal processing algorithm works, we can plot the detected groups.
# You can adjust the threshold and neighborhood keywords on pydarn.proc.music.detectSignals() to tweak
# the autodetection.
pydarn.plotting.plotKarrDetected(dataObj)


# In[ ]: