(* ::Package:: *)

(* Mathematica Package *)

(* Created by Joaquin Mur Amada, July 2008    *)
(* Modified by Joaquin Mur Amada, August 2010 *)


BeginPackage["ColorHueGray2`"]
ColorHueGray::usage = "ColorHueGray[hue,saturation,brightness] gives a ColorFunction so that it is plotted nice in color and \
black&white printers. The hue and the lightness of the color are changed linearly. It can be thought as a gray level \
ColorFunction coloured to increase level identification whereas the graylevel is mantained so that it prints well in B&W. The \
conjecture is that this sequence always terminates.
ColorFunctions: ColorizeFromBlack and ColorizeFromWhite.
Typical Usage: \
ColorFunction \[Rule] RGBColor[ ColorHueGray[0.55- 2(1-z),0.25+0.75z]]"; 
ColorizeFromBlack::usage = "ColorFunctions -> ColorizeFromBlack and ColorizeFromWhite.
Typical Usage: ContourPlot[x^2 \
+y^2,{x,-2,2},{y,-2,2},ColorFunction\[Rule]ColorizeFromBlack]"; 
ColorizeFromBlackLight::usage = "ColorFunctions -> ColorizeFromBlack and ColorizeFromWhite.
Typical Usage: ContourPlot[x^2 \
+y^2,{x,-2,2},{y,-2,2},ColorFunction\[Rule]ColorizeFromBlack]"; 
ColorizeFromWhite::usage = "ColorFunctions -> ColorizeFromBlack and ColorizeFromWhite.
Typical Usage: ContourPlot[x^2 \
+y^2,{x,-2,2},{y,-2,2},ColorFunction\[Rule]ColorizeFromWhite]"; 
ColorizeFromWhiteLight::usage = "ColorFunctions -> ColorizeFromBlack and ColorizeFromWhite.
Typical Usage: ContourPlot[x^2 \
+y^2,{x,-2,2},{y,-2,2},ColorFunction\[Rule]ColorizeFromWhite]"; 
EstiloGrafica::usage = 
   "Variable containing default options for graphs. Example: 
Plot[f[x], {x,-1,1},  Evaluate@EstiloGrafica[[2]] ]"; 
EstiloLogPlot::usage = 
   "Variable containing default options for LogPlot. Example: 
LogPlot[f[x], {x,-1,1},  Evaluate@EstiloLogPlot ]"; 
DetrendSignal::usage = "Remove linear tendency of the signal. Example: 
DetrendSignal[DetrendedSignal]"; 


SmoothCyclicRepetition::usage = "SmoothCyclicRepetition makes a soft transition to smooth out the replication of the signal in \
the next period.
This is done to full use the data, specially in short samples where it is desired to full use the \
data.
SmoothCyclicRepetition weighs with a sinusoidal ramp. 
The cyclic replication does not generate bias in Fourier \
transform because the starting and ending weithing factors sum unity, but it introduces the following bias in \
PSD:
PSDbiasDueToSmoothCyclicRepetition = 1-\!\(\*FractionBox[\"HalfOverlapSamples\", 
RowBox[{\"2\", \" \", \"nSignal\"}]]\) \

Therefore, this bias is corrected in EstimatePSD.
Example: \

SmoothCyclicRepetition[DetrendedSignal,HalfOverlapSamples,HalfJoinSamples]"; 


EstimatePSD::usage = "Estimate PSD, first detrending the signal and then making a SmoothCyclicRepetition, and finally, smoothing \
the auto power spectrum. It returns two elements: the PSD estimation and the phase of the spectrum, since it is computed in the FFT and it can be necessary for some analysis.
Example: 
{autopowerspectrum,spectrumphase} = EstimatePSD[Signal,MinSamples,MaxSamples,HalfOverlapSamples,HalfJoinSamples]
Signal = time \
series

MinSamples = Spectral lines used  to decrease PSD variance at low frequency. A recommended value is about 10. \

Variance at low frequencies is about 1/MinSamples. If MinSamples is too large, the frequency resolution can be very poor at \
low frequencies. If MinSamples is too small, the PSD can experience excesive uncertainty.

MaxSamples = Spectral lines used  \
to decrease PSD variance at low frequency. A recommended value is about 200 for stochastic signals where there is not needed \
high frequency resolution at the end of the spectrum.
At higher frequencies, the spectral lines used in the running average \
increase geometicly with frequency up to MaxSamples. If MaxSamples is too high, the resolution can be poor at high \
frequencies.
Note: 2 > MinSamples and MinSamples +2 >= MaxSamples. If MinSamples == MaxSamples, the number of spectral lines \
used along all frequencies are the same.

HalfOverlapSamples = number of points used to smooth the jump between the start and \
the end of the signal. This lowers the spectral leackage due to the cyclic repetition of the signal implicit in Fourier \
Analysis.
A recommended value is HalfOverlapSamples ~ fsampling/fmain, where fsampling is the sampling frequency and fmain is \
the fundamental frequency (lowest frequency with high energy content).

HalfJoinSamples = number of points used to smooth \
aditionally the very jump between the start and the end of the signal. An order 3 polinomial is fitted to first and last \
HalfJoinSamples to avoid very sharp derivatives in junction.
A recommended value is HalfJoinSamples ~ 5 samples, wich limits \
the frequency about fmax ~ fNyquist/5 efectivelly only in the 5 samples of the cyclic repetition.

(c) Joaquin Mur Amada, \
2009
"; 


DecimatedPSD::usage = "Decimate FullPSD vector geometrically, obtaining a vector of nSpectralLines distributed geometrically. \

Returns a table of two columns {freq, PSD}.
 Sometimes, the variance between the spectral lines is desired instead of the \
average value at such frequency
 In such cases, set the Accumulated to True \
*)
Example:
DecimatedPSD[FullPSD,MinSamples,\[CapitalDelta]f,Accumulated]"; 


Spectrogram::usage = "Spectrogram[data_,opts]: Computes the spectrogram of a timeseries along its periodogram and rms \
value;

The spectrogram is computed applying EstimatePSD to chunks of the time series.

Spectrogram returns a list of vector and matrixes: {SpectogramData, PeriodogramData, WindowQuantilesData, spectrumphase} 

Recommended usage: {Spectogram, Periodogram, WindowQuantiles, Spectrumphase} = Spectrogram[data_,opts]

The spectrogram is contained in the matrix Spectogram. 
The FFT phase is contained in matrix Spectrumphase.
The peridogram is contained in the vector Periodogram.
WindowQuantiles is a matrix containing the median (50% quantile) of each data chunk in the first row and the quantiles 25%, 75%,5%, 95% in the following rows.

Data is the timeseries whose spectrogram is to be computed.

Spectrogram Options 
(default option is shown first and the elsewere option is shown last);

WindowType\[Rule] \
Automatic (rectangular) or \"hanning\". The data is linearly detrended and the data is partially mixed at the extremes of the \
chunk to make smoother the cyclic extension of the signal portion and avoid spectral leakage. Therefore, even with rectangular \
window, there is a litlle windowing at the signal end and begining.

WindowLength\[Rule] Automatic or a number (for example, 32, 64, \
128...). Is the number of samples in each chunk.

WindowOverlap\[Rule] Automatic or a number between 0 (no overlap in contiguous \
chunks of data) and 0.9 (the overlap is 90% and the spectrogram runs very slowly). The automatic value corresponds to 50 % \
overlapping (contiguous chunks share 50% of their data). In fact, overlappings above 50% does not increase time resolution and \
they slow the spectrogram calculus.

LowFreqResolution\[Rule] Automatic or a number.  Spectral lines used  to decrease PSD variance \
at low frequency. A recommended value is about 10. 
Variance at low frequencies is about 1/MinSamples. If MinSamples is too \
large, the frequency resolution can be very poor at low frequencies. If MinSamples is too small, the PSD can experience excesive \
uncertainty.

HighFreqResolution\[Rule] Automatic or a number. Spectral lines used  to decrease PSD variance at low frequency. A \
recommended value is about 200 for stochastic signals where there is not needed high frequency resolution at the end of the \
spectrum.
At higher frequencies, the spectral lines used in the running average increase geometicly with frequency up to \
MaxSamples. If MaxSamples is too high, the resolution can be poor at high frequencies.
Note: 2 > MinSamples >= MaxSamples. If \
MinSamples == MaxSamples, the number of spectral lines used along all frequencies are the same.

CiclicRepetitionOverlap\[Rule] \
Automatic or a number. Number of points used to smooth the jump between the start and the end of the signal. This lowers the \
spectral leackage due to the cyclic repetition of the signal implicit in Fourier Analysis.
A recommended value is \
HalfOverlapSamples ~ fsampling/fmain, where fsampling is the sampling frequency and fmain is the fundamental frequency (lowest \
frequency with high energy content).

CiclicRepetitionJoinSamples\[Rule] Automatic or a number. Number of points used to smooth \
aditionally the very jump between the start and the end of the signal. An order 3 polinomial is fitted to first and last \
HalfJoinSamples to avoid very sharp derivatives in junction.
A recommended value is HalfJoinSamples ~ 5 samples,which limits \
the frequency about fmax ~ fNyquist/5 only in the 5 samples of the cyclic repetition.

Example:
{Spectogram, \
Periodogram, WindowQuantiles, spectrumphase} = Spectrogram[chirp,WindowType\[Rule] \"hamming\", WindowLength\[Rule] 128, WindowOverlap\[Rule] 0.5, \
LowFreqResolution\[Rule]1, HighFreqResolution\[Rule]1, CiclicRepetitionOverlap\[Rule] 5, CiclicRepetitionJoinSamples\[Rule] 2]
"; 



PlotSpectrogram::usage = "PlotSpectrogram[data_,opts]: Plots the spectrogram of a timeseries along its periodogram and rms \
value;

The spectrogram, periodogram and rms value is computed internally calling to Spectrogram. Thus, PlotSpectrogram \
supports the options of Spectrogram. 

Data is the timeseries sampled at SampleRate frequency;

PlotSpectrogram Options \

(default option is shown first and the elsewere option is shown last);

SampleRate\[Rule] (Automatic or a number) indicates the \
sampling rate of the data. If Automatic is specified, the sampling rate is asumed to be 1 Hz;

SpectrogramType \[Rule] Automatic and \
\"Logarithmic\" shows the peridogram, the spectrogram and the signal RMS. \"Logarithmic\" uses logarithmic frequency scale and \
the periodogram shows f.PSD so that the shadowed area corresponds to the variance of the signal. For other specifications, only \
the spectrogram is shown with linear frequency axis;

TimeUnits\[Rule] Automatic or a string. Time is measured in seconds by \
default. This is only used in plot labels.;
FrequencyUnits\[Rule] Automatic or a string. Frequency is measured in Herzs by default. \
This is only used in plot labels.;

SpectrogramColors\[Rule] 
      Automatic: Colorize with lowest value white and biggest value \
black with varying hue to get an acceptable representation in grayscale devices such as monochrome printers. This is equivalent \
to ColorFunction ->ColorizeFromWhite;
      \"Lighter\": Colormap is a bit lighter than the previous one. This is equivalent to \
ColorFunction -> ColorizeFromWhiteLight.
      \"Darker\": Colorize with biggest value white and lowest value black with \
varying hue compatible with grayscale devices. This is equivalent to ColorFunction ->ColorizeFromBlack.
      The default case \
is colorize from dark red up (smallest value) to white (biggest value). This is equivalent to ColorFunction \
->ColorizeFromBlack;

DynamicRange\[Rule] Automatic or a number. If a number is specified, the density plot color is graded acording \
to the logarithm of the PSD estimate. The dynamic range by default is 10 Log10[10^-3] = 30 dB. Logarithmic option increases the \
dynamic range of the plot at low values, but it distorts the density plot if verly high dynamic range is used. Automatic grades \
linearly according to PSD estimate and it produces more realistic plots, eventhough some details of little vibrations can be \
scarcelly visible.

SignalName -> Automatic or a string. If a string is provided, it is used as vertical plot label in RMS \
plot

PeriodogramUnits -> Automatic or a string. Periodogram units are used in legend scale and periodogram. If linear \
frequency axis is used, it is the signal units squared divided by frequency units. If logatithmic frequency axis is used, it is \
the signal units squared.

DateStart -> Automatic or a date in any Date Mathematica format (see AbsoluteTime). If set to a \
date (in format of absolute time or in Mathematica Date format) and DateEnd set to a date, then the RMS plot uses date format \
instead of number of time units (seconds, days...). Otherwise, a customary time axis in TimeUnits is used.
DateEnd -> Automatic \
or a date in any Date Mathematica format. If set to a date together with DateStart, date format is used in time axis of RMS \
plot..

SpectogramData -> Automatic or a matrix with the same conventions as the Spectogram function.
PeriodogramData -> \
Automatic or a vector with the same conventions as the Spectogram function.
WindowQuantilesData -> Automatic or a is a matrix containing the median (50% quantile) of each data chunk in the first row and the quantiles 25%, 75%,5%, 95% in the following rows (the \
same conventions as the Spectogram function).

TextExplanation -> Automatic or a expression to be shown in the bottom left \
corner.

LogPeriodogramPlot -> Automatic, True or False. If set to True, the periodogram is plot in logarithmic magnitude scale

MultipliedByFrequency -> Automatic, True or False. Indicates if periodogram and spectrogram must be plotted as f*PSD(f) and f*PSD(f,t), respectivelly.

PeridogramLabel -> Automatic or a string. Allows to change the default label for the magnitude of the periodogram.

ShadowInterQuantiles->Automatic, True or False. Indicates if it is expected to receive in windowquantiles the median and quantiles 5%, 25%, 75% and 95% in that order.
A shadow between the simmetric quanttiles to indicate signal envelop. If ShadowInterQuantiles is False, no shadow is added and references to the quantiles in the legend is removed..

TimePlotOptions->Automatic or a list with the additional options of the interquantil time plot.

SpectrogramLimitFactor\[Rule] Automatic or a number. The spectrogram color table is the minimum of periodogram divided by SpectrogramLimitFactor up to the maximum of the periodogram multiplied by SpectrogramLimitFactor. By default, the automatic value corresponds to 2.

If SpectogramData is Automatic or it is not set, the periodogram, \
the spectogram and the windowquantiles is computed with the function {spect,periodogram,windowquantiles,spectrumphase} = \
Spectrogram[data];

Example:
PlotSpectrogram[chirp,SampleRate\[Rule] fsampling, SpectrogramType \[Rule] \"Logarithmic\", \
SpectrogramColors -> \"Lighter\", DynamicRange\[Rule] 20 (*Dynamic range in dB*), TimeUnits\[Rule] \"days\", FrequencyUnits\[Rule] \
\"cycles/day\",SpectrogramColors\[Rule] Lighter ,WindowType\[Rule] \"hamming\" ,WindowLength\[Rule] 128,WindowOverlap\[Rule] \
0.5,LowFreqResolution\[Rule]1,HighFreqResolution\[Rule]1,CiclicRepetitionOverlap\[Rule] 5,CiclicRepetitionJoinSamples\[Rule] 2]
"; 



PlotSpectrogram2::usage = "PlotSpectrogram2[spectogram_,periodogram_:Automatic,TimeSeries_:Automatic,opts]: Plots the \
spectrogram of a timeseries along its periodogram and rms value;

This function is provided to plots spectogram where the spectogram has been computed previously.

The options are the same than PlotSpectrogram and refefer to this function for further information.

TimeSeries can be substituted with a WindowQuantiles matrix containing the median (50% quantile) of each data chunk in the first row and the quantiles 25%, 75%,5%, 95% in the following rows.


The only difference is that spectrogram, periodogram and rms time series are passed to the function instead of the raw data. 
It is usefull when the \
spectrogram has been computed previously.
"; 
SafeDivide::usage = "SafeDivide[x1_,x2_,maxratio_: (10^12)] computes the vector (x2[[i]]/x1[[i]]) or the matrix \
(x2[[i,j]]/x1[[i,j]]) elemen-wise division avoiding zero division errors and outliers due to divisions by a number close to zero \
limiting the maximum possible ratio x2/x1 before computing it.
The maximum default ratio is 10^12, but any value is \
possible.
Example: 
* Compute vector element wise division with a maximum ratio of 100 -> \
SafeDivide[{1,2,3,4},{2,4,6,0.000001}, 100];
* Compute matrix element wise division with a maximum ratio of 100 -> \
SafeDivide[{{1,2},{3,4}},{{2,4},{6,0.000001}}, 100]
"; 


LogartithmicRunningAverage::usage=
"Produces a running average of the signal with a varying number of samples.
It is expected MinSamples <= MaxSamples. If MinSamples == MaxSamples, then the signal is filtered with a constant number of samples. 
For example, if MinSamples == MaxSamples == 3, then the filtered signal consist in the average of original points with the adjacent ones.

However, it is intended for filtering a signal with variable strength.It is usefull for filtering spectrum, where the first samples have lower uncertainty (or more information) and \
last samples have bigger uncertainty (or less information).

It is also usefull for logarithmic plots, since the samples of the beginning are less filtered and the end ones are more strongly filtered.
Since the strength of the filter is geometric, the filtered signal appears with constant variance when plotted in logaritmic plots.

Exaple: 
LogartithmicRunningAverage[Signal,MinSamples (*by default, MinSamples==3 *) ,MaxSamples (* by default, MaxSamples == 9 *)]
A running average of the signal is returned. At the start of the signal, only 4 samples are used in the samples

LogartithmicRunningAverage[Signal]
";


NumberToString::usage="Writes a number in CForm using only the number of characters especified in precision and a point. Sometimes uses the 1.3E8 format and others use the 13.31 or 0.012 formats. It returns a string.
Usage: NumberToString[-1.2345678*10^-8]  =>  \"-1.E-8\"
NumberToString[-1.2345678*10^-8,4]  =>  \"-1.24E-8\" ";



Options[Spectrogram]= {WindowType-> Automatic ,WindowLength-> Automatic,WindowOverlap-> Automatic,LowFreqResolution-> Automatic,HighFreqResolution-> Automatic,CiclicRepetitionOverlap-> Automatic,CiclicRepetitionJoinSamples-> Automatic};
Options[PlotSpectrogram]= {SampleRate-> Automatic, SpectrogramType -> Automatic, TimeUnits-> Automatic, FrequencyUnits-> Automatic, SpectrogramColors-> Automatic, DynamicRange-> Automatic,SignalName -> Automatic, PeriodogramUnits -> Automatic, DateStart -> Automatic, DateEnd -> Automatic,SpectogramData -> Automatic, PeriodogramData -> Automatic, WindowQuantilesData -> Automatic, SignalTimeDuration-> Automatic, TextExplanation -> Automatic, LogPeriodogramPlot -> Automatic, MultipliedByFrequency -> Automatic, PeridogramLabel -> Automatic,ShadowInterQuantiles->Automatic,TimePlotOptions->Automatic, SpectrogramLimitFactor-> Automatic};
Options[PlotSpectrogram2]= {SampleRate-> Automatic, SpectrogramType -> Automatic, TimeUnits-> Automatic, FrequencyUnits-> Automatic, SpectrogramColors-> Automatic, DynamicRange-> Automatic,SignalName -> Automatic, PeriodogramUnits -> Automatic, DateStart -> Automatic, DateEnd -> Automatic,SpectogramData -> Automatic, PeriodogramData -> Automatic, WindowQuantilesData -> Automatic, SignalTimeDuration-> Automatic, TextExplanation -> Automatic, LogPeriodogramPlot -> Automatic, MultipliedByFrequency -> Automatic, PeridogramLabel -> Automatic,ShadowInterQuantiles->Automatic,TimePlotOptions->Automatic, SpectrogramLimitFactor-> Automatic};


UnwrapPhase1::usage="UnwrapPhase1[x_, gap_: \[Pi], jump_: 2 \[Pi] ] 
(* gap must be greater than average difference between x elements in order to avoid detecting fake 2 pi phase jumps *)
(* divide matrix in row vectors *)

(* Guess a phase jump when the phase difference between adjacent elemts is bigger than gap. It only considers 2 points, no model of the phase vector.*)
(* Phase is in radians if gap \[TildeTilde] \[Pi] and jump = 2 \[Pi]
 Phase is in degrees if gap \[TildeTilde] 180 and jump = 360 *)

Unwrap phase\[LineSeparator]General model without asuming any pattern in phase vector, using 2 samples.
Gap must be greater than average difference between adjacent elements in order to detect phase jumps not observed only looking to the next sample.
Phase is in radians if gap \[TildeTilde] \[Pi] and jump = 2 \[Pi]
 Phase is in degrees if gap \[TildeTilde] 180 and jump = 360 

Example:
wrappedphase =UnwrapPhase1[phase];

See: UnwrapPhase1, UnwrapPhase2, UnwrapPhase3 and UnwrapLeastSquares
";



UnwrapPhase2::usage="UnwrapPhase2[x_,gap_:\[Pi], jump_: 2 \[Pi] ]
Unwrap phase\[LineSeparator]a) General model without asuming any pattern in phase vector, using 4 samples.
Gap must be greater than average difference between 4 elements in order to detect phase jumps not observed only looking to the next sample.
Phase is in radians if gap \[TildeTilde] \[Pi] and jump = 2 \[Pi]
 Phase is in degrees if gap \[TildeTilde] 180 and jump = 360 

See: UnwrapPhase1, UnwrapPhase2, UnwrapPhase3 and UnwrapLeastSquares

Example:
wrappedphase =UnwrapPhase2[phase];

";


UnwrapPhase3::usage="UnwrapPhase3[x_,gap_:\[Pi], jump_: 2 \[Pi] ]
a) General model without asuming any pattern in phase vector, using 6 samples.
Gap must be greater than average difference between 6 elements in order to detect phase jumps not observed only looking to the next sample.
Phase is in radians if gap \[TildeTilde] \[Pi] and jump = 2 \[Pi]
 Phase is in degrees if gap \[TildeTilde] 180 and jump = 360 

See: UnwrapPhase1, UnwrapPhase2, UnwrapPhase3 and UnwrapLeastSquares

Example:
wrappedphase =UnwrapPhase3[phase];
";


UnwrapPhase::usage="UnwrapPhase3[phase_]=UnwrapLeastSquares[UnwrapPhase1[phase]]

Unwrap phase combines two models that complements each other:\[LineSeparator]a) Simple and general model without asuming any pattern in phase vector using 2 samples.
b) Model using 6 samples specially devoted when the phase is expected to change at a approximatelly constant pace:
* the phase of a signal of more or less constant frequency.
* phase of senoidal signals, angular phase of motors turning at a speed that can not change instantaneously...
* the FFT of a timeseres 
=> a shift T in time domain increases the argument of the FFT in 2\[Pi] f T (thus, the argument of FFT is expected to vary more or less linearly with frequency). 
=> an signal step in time, not centered in t = 0, produces an FFT argument proportional to frequency.
* the phase of the complex coherence or the complex power spectrum 

See: UnwrapPhase1, UnwrapPhase2, UnwrapPhase3 and UnwrapLeastSquares
"


PhaseIncrement::usage=" PhaseIncrement[x_] (* Differenciate Integrate vector x[[i]]-x[[i-1]], assumming inital value x[[0]]=x[[1]] *)";


IntegratePhase::usage=" IntegratePhase[x_,x0_:0]
Integrate vector y[[i]] = Sum[x[[j]],{i,1,i}]
";


UnwrapLeastSquares::usage="UnwrapLeastSquares[x_] 

Unwrap phase\[LineSeparator]a) General model without asuming any pattern in phase vector, using 2 samples, 4 and 6.
b) Model using 6 samples specially devoted when the phase is expected to change at a approximatelly constant pace:
* the phase of a signal of more or less constant frequency.
* phase of senoidal signals, angular phase of motors turning at a speed that can not change instantaneously...
* the FFT of a timeseres 
=> a shift T in time domain increases the argument of the FFT in 2\[Pi] f T (thus, the argument of FFT is expected to vary more or less linearly with frequency). 
=> an signal step in time, not centered in t = 0, produces an FFT argument proportional to frequency.
* the phase of the complex coherence or the complex power spectrum 
Use: 
you can use first the general model using 2 and 4 samples and then, fine tune the result running the model specific for unwarp FFT of timeseries.
The same result has been obtained also only applying the specific, 6 samples FFT

(* Phase must be in radians *)

Example:
For best results, use first the general model using 2 samples 
wrappedphase =UnwrapPhase1[phase];

Or even with 4 samples:

Then, fine tune the result running the model specific for unwarp FFT.
wrappedphase = UnwrapLeastSquares[wrappedphase ];


See: UnwrapPhase1, UnwrapPhase2, UnwrapPhase3 and UnwrapLeastSquares
";


MatrixPlotWithScale::usage = 
"MatrixPlotWithScale[colorMatrix, contourMatrix ,opts]: Plots a matrix with a legend of the color scale. 
The color of each point is interpollated from colorMatrix and the contour lines are drawn from contourMatrix. If contourMatrix is not a numeric matrix, the contour lines are derived from colorMatrix.

There are several options intended for plotting probability matrix. Aditionnaly, the options addmited by ListDensityPlot and ListContourPlot can be used (for example, Frame -> False elimitane the naming of the matrix rows and columns).

StateNames -> The rows and columns of the matrix can be specified.
    The StateNames of an equally spaced state discretization can be defined as {fistStateValue, lastStateValue}. If StateNames is a vector, then the DataRange option of the Mathematica built-in function MatrixPlot is DataRange-> {option,option}.
If StateNames is not a vector, then its value is passed directly to the Mathematica built-in function MatrixPlot as DataRange option (i.e, DataRange-> {option,option})

Units\[Rule] The units of the elements of the matrix are only used in the color legend;

MinRange -> The minimum value represented in the colormap. If MinRange is not specified, then the matrix minimum is used. DynamicRange can modify this value

MaxRange -> The maximum value represented in the colormap. If MinRange is not specified, then the matrix minimum is used.

Colors\[Rule] 
      Automatic: Colorize with lowest value white and biggest value \
black with varying hue to get an acceptable representation in grayscale devices such as monochrome printers. This is equivalent \
to ColorFunction ->ColorizeFromWhite;
      \"Lighter\": Colormap is a bit lighter than the previous one. This is equivalent to \
ColorFunction -> ColorizeFromWhiteLight.
      \"Darker\": Colorize with biggest value white and lowest value black with \
varying hue compatible with grayscale devices. This is equivalent to ColorFunction ->ColorizeFromBlack.
      The default case \
is colorize from dark red up (smallest value) to white (biggest value). This is equivalent to ColorFunction \
->ColorizeFromBlack;

DynamicRange\[Rule] Automatic or a number. If a number is specified, the color scale is logarithmic. The dynamic range by default is 10 Log10[10^-3] = 30 dB. Logarithmic option increases the \
dynamic range of the plot at low values, but it distorts the plot if verly high dynamic range is used. Automatic grades \
linearly according  and it produces more realistic plots, eventhough some details of very low values can be \
scarcelly visible.

Scaling \[Rule] \"Logarithmic\" uses logarithmic color scale. For other specifications, only the spectrogram is shown with linear frequency axis. 

LimitFactor\[Rule] Automatic or a number (for example, 1. or 2.). The color table is the minimum matrix entry divided by LimitFactor up to the maximum of the periodogram multiplied by LimitFactor. By default, the automatic value corresponds to 2.

ShowColorScale is True by default. If ShowColorScale-> False, then the legend is not shown.

ShowContours \[Rule] True by default. If ShowContours \[Rule] False, then the contour plots are ommited.

==================================================================================================
Examples:
MatrixPlotWithScale[(* Matrix for drawing color levels *) 100 ProbabilityTransitionMatrix, (* Matrix for drawing contours and its text *) Automatic, StateLabels\[Rule] {3.6,38.88},DynamicRange\[Rule]15, Units\[Rule] \"%\",LimitFactor\[Rule] 1,Scaling \[Rule] \"Logarithmic\",MinRange\[Rule] 0,MaxRange\[Rule] 100,FrameLabel\[Rule]{\"previous state\",\"next state\"},PlotLabel\[Rule]\"Transition Probability Matrix\", ShowColorScale\[Rule] True, ShowContours \[Rule] True]


MatrixPlotWithScale[(* Matrix for drawing color levels *) 100 ProbabilityTransitionMatrix, (* Matrix for drawing contours and its text *) 
Automatic, StateLabels\[Rule] {3.6,38.88},DynamicRange\[Rule]Automatic, Units\[Rule] \"%\",LimitFactor\[Rule] 1,Scaling \[Rule] \"Linear\",MinRange\[Rule] 0,MaxRange\[Rule] 100,FrameLabel\[Rule]{\"previous state\",\"next state\"},PlotLabel\[Rule]\"Transition Probability Matrix\", ShowColorScale\[Rule] True, ShowContours \[Rule] False]


MatrixPlotWithScale[(* Matrix for drawing color levels *)   ProbabilityTransitionMatrix,(* Matrix for drawing contours and its text *) 100CumulativeProbabilityTransitionMatrix, StateLabels\[Rule] {3.6,38.88},DynamicRange->Automatic, Units\[Rule] \"%\",LimitFactor\[Rule] 1,Scaling \[Rule] \"Linear\",FrameLabel\[Rule]{\"previous state\",\"next state\"},PlotLabel\[Rule]\"Transition Probability Matrix\", ShowColorScale\[Rule] False,
     ShowContours \[Rule] True,Contours\[Rule]19]

MatrixPlotWithScale[
    (* Matrix for drawing color levels *)   CumulativeProbabilityTransitionMatrix,
    (* Matrix for drawing contours and its text *) Automatic, 
    StateLabels\[Rule] {3.6,38.88},DynamicRange->Automatic, Units\[Rule] \"%\",LimitFactor\[Rule] 1,Scaling \[Rule] \"Linear\",FrameLabel\[Rule]{\"previous state\",\"next state\"},PlotLabel\[Rule]\"Transition Probability Matrix\", ShowColorScale\[Rule] False,
     ShowContours \[Rule] True,Contours\[Rule]9]"; 


Options[MatrixPlotWithScale]= {StateLabels->Automatic,Units->Automatic,MinRange-> Automatic,MaxRange-> Automatic,Colors->Automatic,DynamicRange->Automatic, Scaling->Automatic,LimitFactor-> Automatic, ShowColorScale -> True, ShowContours -> True};


EstimateCumulativeProbabilityTransitionMatrixBiased::usage = "EstimateCumulativeProbabilityTransitionMatrixBiased[ProbabilityTransitionMatrix] : Estimates (with bias) the cumulative transition matrix from the probability transition matrix of a Markov chain. The cumulative probability transition matrix is computed with a bias of half sate.

In fact, aggregates the row elements, starting on the left. The last element is forced to uniti (1) since roundoff errors can make the element on the right deviate from unity.

Agregates the jump probability of a Markov ProbabilityTransitionMatrix. Plots a matrix with a legend of the color scale. Aditionnaly, the options addmited by PlotMatrix can be used (for example, Frame -> False elimitane the naming of the matrix rows and columns).


==================================================================================================
Examples:
ProbabilityTransitionMatrix=
{{0.6811,0.2916,0.0232,0.0036,0.0004},
{0.186,0.5852,0.2146,0.0131,0.0011},
{0.0116,0.2493,0.5571,0.1755,0.0065},
{0.0013,0.0119,0.2421,0.5755,0.1693},
{0.0002,0.0006,0.0159,0.2569,0.7264}};

EstimateCumulativeProbabilityTransitionMatrix[ProbabilityTransitionMatrix]

EstimateCumulativeProbabilityTransitionMatrixBiased[
{{a11,a12,a13},
{a21,a22,a23},
{a31,a32,a33}}
]
 
Result: 
{{a11,a11+a12,1},
{a21,a21+a22,1},
{a31,a31+a32,1}}

TableForm[%]
\!\(\*TagBox[GridBox[{
   {\"a11\", 
    RowBox[{\"a11\", \"+\", \"a12\"}], \"1\"},
   {\"a21\", 
    RowBox[{\"a21\", \"+\", \"a22\"}], \"1\"},
   {\"a31\", 
    RowBox[{\"a31\", \"+\", \"a32\"}], \"1\"}
  },\nGridBoxAlignment->{\"Columns\" -> {{Left}}, \"ColumnsIndexed\" -> {}, \"Rows\" -> {{Baseline}}, \"RowsIndexed\" -> {}},\nGridBoxSpacings->{\"Columns\" -> {Offset[0.27999999999999997`], {Offset[0.7]}, Offset[0.27999999999999997`]}, \"ColumnsIndexed\" -> {}, \"Rows\" -> {Offset[0.2], {Offset[0.4]}, Offset[0.2]}, \"RowsIndexed\" -> {}}],
 Function[BoxForm`e$, TableForm[BoxForm`e$]]]\)

"; 


EstimateCumulativeProbabilityTransitionMatrix::usage = "EstimateCumulativeProbabilityTransitionMatrix[ProbabilityTransitionMatrix] : Estimates the cumulative transition matrix from the probability transition matrix of a Markov chain.
 To avoid biases in the cumulative probability transition matrix, the cumulative probability is computed in the median of each bin. It is asummed that half of the probability lies on the right and on the left of the state centroid.
This is the reason for accounting half of the probability for each element.
This algorithm is recommended, at least for ploting the contours of the ProbabilityTransitionMatrix.

In fact, aggregates the row elements, starting on the left. The last element is forced to uniti (1) since roundoff errors can make the element on the right deviate from unity.

Agregates the jump probability of a Markov ProbabilityTransitionMatrix. Plots a matrix with a legend of the color scale. Aditionnaly, the options addmited by PlotMatrix can be used (for example, Frame -> False elimitane the naming of the matrix rows and columns).


==================================================================================================
Examples:
ProbabilityTransitionMatrix=
{{0.6811,0.2916,0.0232,0.0036,0.0004},
{0.186,0.5852,0.2146,0.0131,0.0011},
{0.0116,0.2493,0.5571,0.1755,0.0065},
{0.0013,0.0119,0.2421,0.5755,0.1693},
{0.0002,0.0006,0.0159,0.2569,0.7264}};

EstimateCumulativeProbabilityTransitionMatrix[ProbabilityTransitionMatrix]
TableForm[%]
Result:
0.34055      0.8269       0.9843       0.9977       0.9997
0.093        0.4786       0.8785       0.99235      0.99945
0.0058       0.13625      0.53945      0.90575      0.99675
0.00065      0.00725      0.13425      0.54305      0.91545
0.0001       0.0005       0.00875      0.14515      0.6368

EstimateCumulativeProbabilityTransitionMatrix[
{{a11,a12,a13},
{a21,a22,a23},
{a31,a32,a33}}
]

Result: 
{{0.5 a11,0.5 a11+0.5 (a11+a12),0.5 a11+0.5 (a11+a12)+0.5 (a12+a13)},
{0.5 a21,0.5 a21+0.5 (a21+a22),0.5 a21+0.5 (a21+a22)+0.5 (a22+a23)},
{0.5 a31,0.5 a31+0.5 (a31+a32),0.5 a31+0.5 (a31+a32)+0.5 (a32+a33)}}

TableForm[Simplify[%]]

Result: 
0.5 a11      a11 + 0.5 a12     a11 + a12 + 0.5 a13
0.5 a21      a21 + 0.5 a22     a21 + a22 + 0.5 a23
0.5 a31      a31 + 0.5 a32     a31 + a32 + 0.5 a33
"; 


Begin["`Private`"];

(*Clear[ColorHueGray];*)
(* La funci\[OAcute]n de color toma las coordenadas x e y del punto y la elevaci\[OAcute]n (valor de la funci\[OAcute]n representada) en z. Con ello, s\[OAcute]lo utiliza el valor de la funci\[OAcute]n representada para dibujar *)
(* with only one number, define a function that maps well in color and grayscale spaces colorHueGray[h_]:=RGBColor[ colorHueGray[h,h]]*) 
ColorizeFromBlack[z_]:=ColorHueGray[0.55- 2(1-z),Min[1, 0.25+ 1.15 Abs[Cos[\[Pi] z]-0.4]],0.25+0.75z];
ColorizeFromBlack[y_,z_]:=ColorizeFromBlack[z];
ColorizeFromBlack[x_,y_,z_]:=ColorizeFromBlack[z];
ColorizeFromWhite[z_]:=ColorizeFromBlack[1-z];
ColorizeFromWhite[y_,z_]:=ColorizeFromWhite[z];
ColorizeFromWhite[x_,y_,z_]:=ColorizeFromWhite[z];

ColorizeFromBlackLight[z_]:=ColorHueGray[(1.35-z)4/6,Min[1, 0.25+  1.15Abs[Cos[\[Pi] (z-0.15)]-0.4]],Min[1,0.45+0.60z]];
ColorizeFromBlackLight[y_,z_]:=ColorizeFromBlackLight[z];
ColorizeFromBlackLight[x_,y_,z_]:=ColorizeFromBlackLight[z];
ColorizeFromWhiteLight[z_]:=ColorizeFromBlackLight[1-z];
ColorizeFromWhiteLight[y_,z_]:=ColorizeFromWhiteLight[z];
ColorizeFromWhiteLight[x_,y_,z_]:=ColorizeFromWhiteLight[z];

ColorHueGray[hue_,saturation_,brightness_]:=Module[{hi,f,p,q,t,rgb={1,1,1},limitingcolors={1,2,3},weighting= {0.299,0.587,0.114} },
Which[brightness <=  0, rgb={0,0,0},brightness >= 1, rgb={1,1,1},0< brightness < 1,
(* Default case *)
hi = Mod[Floor[6 hue],6];
f = 6 hue - Floor[6 hue];
(* saturation -> 1 ; *)
p = (1-saturation);
q = (1 - f saturation);
t = (1-(1-f) saturation);
(* In first instance, try saturated colors saturation -> 1 ;
p =0;
q = (1 - f );
t = f; *)
(* Compute the hue with value 1 and then, adjust value from brightness *)
rgb =Switch[hi,0, {1,t,p}, 1,{q,1,p},2,{p,1,t},3,{p,q,1},4,{t,p,1},5,{1,p,q},_, {1,t,p}];
rgb = rgb (brightness /( weighting.rgb));
(* Para niveles de gris oscuros, hay que bajar la saturaci\[OAcute]n si se quiere mantener el tono *)
If[Max[rgb]>1, (* find largest color to compute the saturation*)
limitingcolors = Ordering[rgb];
(* more saturated channel *)
rgb[[limitingcolors[[3]] ]]=1;
(* try to adjust channel with intermediate content proportionally*)
rgb[[limitingcolors[[2]] ]] = Min[1,(brightness-  weighting[[limitingcolors[[3]] ]]* rgb[[limitingcolors[[3]] ]])*If[rgb[[limitingcolors[[2]] ]] <= 0, 1/(weighting[[limitingcolors[[1]] ]]  +weighting[[limitingcolors[[2]] ]] ),rgb[[limitingcolors[[2]] ]]/(weighting[[limitingcolors[[1]] ]]*rgb[[limitingcolors[[1]] ]]+weighting[[limitingcolors[[2]] ]]*rgb[[limitingcolors[[2]] ]])
]
];
(* Adjust the lightest channel to obtain the desired graylevel *)
(* despeja la saturaci\[OAcute]n necesaria para que el color tenga un tono -hue- y un tono de gris equivalente *)
rgb[[limitingcolors[[1]] ]] = Min[1,Max[0,(brightness -  weighting[[limitingcolors[[3]] ]] * rgb[[limitingcolors[[3]] ]]-  weighting[[limitingcolors[[2]] ]] rgb[[limitingcolors[[2]] ]])/weighting[[limitingcolors[[1]] ]]
]
];
]; (*endif*)
]; (*end which *)
rgb
];(*end module *)


NivelGris={0,1/2,1/3,5/6,2/3,1/6,1};
tamano = {{500},{400}};
tamano = Medium;
axes = Automatic;
joined = True;
PuntosYrayas={{},Tiny,{0.02},{0.04,0.03},
{0.08,0.02,0.01,0.02},{0.08,0.02},{0.01,0.02},{0.005,0.02,0.05,0.02},{0.005,0.02,0.02,0.02},{0.01},{0.06,0.02,0.01,0.02,0.01,0.02}};
EstiloGrafica ={
{GridLines-> Automatic,Frame->True, Axes-> False,AxesOrigin->axes,PlotStyle->Table[{Hue[(-y)/6], Thickness[0.01],Dashing[PuntosYrayas[[y+1]]]},{y,0,6}], Filling-> Axis,Background->None, ImageSize->tamano},
{GridLines-> Automatic,Frame->True, Axes-> False,AxesOrigin->axes,PlotStyle->Table[{Hue[(-y)/6], Thickness[0.006*If[y>0,1.5,2.5]],Dashing[PuntosYrayas[[y+1]]]},{y,0,6}], Filling-> Axis,Background->None, ImageSize->tamano},
{GridLines-> Automatic,Frame->True, Axes-> False,AxesOrigin->axes,PlotStyle->Table[{Hue[(-y)/6], Thickness[0.006*If[y>0,1.5,2.5]],Dashing[{}]},{y,0,6}], Filling-> Axis,Background->None, ImageSize->tamano},
{GridLines-> None,Frame->True, Axes-> False,AxesOrigin->axes,PlotStyle->Table[{Hue[-y/10], Thickness[0.01],Dashing[PuntosYrayas[[y+1]]]},{y,0,10}], Filling-> Axis,Background->None},
{GridLines-> Automatic,Frame->True, Axes-> False,AxesOrigin->axes,PlotStyle->Table[{GrayLevel[NivelGris[[y+1]]], Thickness[0.01],Dashing[PuntosYrayas[[y+1]]]},{y,0,6}], Filling-> Bottom,Background->None, ImageSize->tamano},
{GridLines-> Automatic,Frame->True, Axes-> False,AxesOrigin->axes,PlotStyle->Table[{GrayLevel[NivelGris[[y+1]]], Thickness[0.006*If[y>0,1.5,2.5]],Dashing[PuntosYrayas[[y+1]]]},{y,0,6}], Filling-> Bottom,Background->None, ImageSize->tamano}
};
EstiloLogPlot ={Frame->True, Axes-> False, AxesOrigin->Automatic, PlotStyle->Table[{GrayLevel[NivelGris[[y+1]]], Thickness[0.006*If[y>0,1.5,2.5]],Dashing[PuntosYrayas[[y+1]]]},{y,0,6}], Filling-> None, Background->None, GridLines->{Automatic,Automatic}, GridLinesStyle->Directive[Yellow], ImageSize->tamano};
EstiloListLogPlot = Join[{Joined->joined}, EstiloLogPlot];




(SetOptions[Plot, #1] &) /@ EstiloGrafica[[1]];
(SetOptions[ListPlot, #1] &) /@ EstiloGrafica[[1]];
(SetOptions[ListLinePlot, #1] &) /@ EstiloGrafica[[1]];
(SetOptions[DateListPlot, #1] &) /@ EstiloGrafica[[1]];
(SetOptions[LogPlot, #1] &) /@  EstiloLogPlot;
(SetOptions[ListLogPlot, #1] &) /@ EstiloListLogPlot;
(SetOptions[LogLogPlot, #1] &) /@ EstiloLogPlot;
(SetOptions[ListLogLogPlot, #1] &) /@ EstiloListLogPlot;
(SetOptions[ListLogLinearPlot, #1] &) /@ EstiloListLogPlot;
(SetOptions[LogLinearPlot, #1] &) /@ EstiloLogPlot;
(SetOptions[ListLogLogPlot, #1] &) /@ EstiloListLogPlot;



(*
(SetOptions[Plot3D, #1] &) /@ {ColorFunction->(RGBColor@ColorizeFromBlackLight[#3]&), MeshFunctions->{#1&,#2&,#3&}, Mesh->{5,5,9},MeshStyle->{Gray,Gray,Darker@Brown}};
(SetOptions[ContourPlot, #1] &) /@ {ColorFunction->(RGBColor@ColorizeFromWhiteLight[#3]&),ContourLabels -> Automatic };
(SetOptions[DensityPlot, #1] &) /@ {ColorFunction->(RGBColor@ColorizeFromWhiteLight[#3]&), MeshFunctions->{#1&,#2&,#3&}, Mesh-> {5,5,9},MeshStyle->{Gray,Gray,Darker@Brown}};
*)

SetOptions[Plot3D, ColorFunction->(ColorHueGray[0.9` -0.6666666666666666`#3,Min[1,0.25` +1.15` Abs[-0.4`+Cos[\[Pi] (-0.15`+#3)]]],Min[1,0.45` +0.6` #3]]&), MeshFunctions->{#1&,#2&,#3&}, Mesh->{5,5,9},MeshStyle->{LightGray,LightGray,Darker@Brown}];
SetOptions[ContourPlot, ColorFunction->(RGBColor@ColorHueGray[0.2333333333333334` +0.6666666666666666` #1,Min[1,0.25` +1.15` Abs[-0.4`+Cos[2.670353755551324` -\[Pi] #1]]],Min[1,1.05` -0.6` #1]]&),ContourLabels -> All ];
SetOptions[DensityPlot, ColorFunction->(RGBColor@ColorHueGray[0.2333333333333334` +0.6666666666666666` #1,Min[1,0.25` +1.15` Abs[-0.4`+Cos[2.670353755551324` -\[Pi] #1]]],Min[1,1.05` -0.6` #1]]&), MeshFunctions->{#1&,#2&,#3&}, Mesh->{5,5,9}, MeshStyle->{LightGray,LightGray,Darker@Brown}];
SetOptions[ListDensityPlot, ColorFunction->(RGBColor@ColorHueGray[0.2333333333333334` +0.6666666666666666` #1,Min[1,0.25` +1.15` Abs[-0.4`+Cos[2.670353755551324` -\[Pi] #1]]],Min[1,1.05` -0.6` #1]]&), MeshFunctions->{#1&,#2&,#3&}, Mesh->{5,5,9}, MeshStyle->{LightGray,LightGray,Darker@Brown}];



(**********************************************************************************************************************)
DetrendSignal[DetrendedSignal_]:=
LinearModelFit[DetrendedSignal,
(* SignalTrendModel *)  x (* {x,x^2,x^3,x^4,x^5} or x; *),
(* Variable in Signal Trend *)x
][ (* Returned vector DetrendedSignal *)
"FitResiduals"]


(**********************************************************************************************************************)
SmoothCyclicRepetition[DetrendedSignal_,HalfOverlapSamples_:30,HalfJoinSamples_:5]:=
Module[{a,x1,x2,CyclicJoint,WindowedVector, residualVector,SmoothCyclicRepetitionWeithing,nSignal,output},
(* Make a soft transition to smooth out the replication of the signal in the next period to prevent introducing spectral leackage *)
(* This is done to full use the data, specially in short samples where it is desired to full use the data *)
(* SmoothCyclicRepetition weighs with a sinusoidal ramp *)
SmoothCyclicRepetitionWeithing[x_]=(1+Sin[ \[Pi] x /(2HalfOverlapSamples) ])/2;
nSignal = Dimensions[DetrendedSignal][[1]];
(* DetrendedSignal is a vector passed by value, not by ref *)
(* The value of the vector cannot be changed *)
(* Thus, the vector cannot be modified in-line *)
(* The vector is sized here *)
output =DetrendedSignal;
(* make a weighted transition at the start/end *)
Do[a = SmoothCyclicRepetitionWeithing[i];
x1 = a*DetrendedSignal[[i+1]]+(1-a)*DetrendedSignal[[nSignal-i]];
x2 =(1- a)*DetrendedSignal[[i+1]]+a*DetrendedSignal[[nSignal-i]];
output[[i+1]] = x1;
output[[nSignal-i]] = x2;
,{i,0,HalfOverlapSamples}
];
(* To avoid very sharp derivatives in junction, smooth out some samples in the joint *)
WindowedVector = Join[output[[(nSignal-HalfJoinSamples);;nSignal]],output[[1;;HalfJoinSamples+1]] ];
residualVector=LinearModelFit[WindowedVector,{x,x^2,x^3},x]["FitResiduals"];

If[HalfJoinSamples>3, (* smooth joint *)
(* Now fit a smoother version in the HalfJoinSamples *)
Do[a =Sin[\[Pi] i /(2 HalfJoinSamples)]^2; (* Windowing weighting *)
output[[i+1]] =WindowedVector[[HalfJoinSamples+i+2]] -(1-a)* residualVector[[HalfJoinSamples+i+2]];
output[[nSignal-i]] =WindowedVector[[HalfJoinSamples-i+1]]-(1-a)* residualVector[[HalfJoinSamples-i+1]];
,{i,0,Round[HalfJoinSamples]}
];
]; (* end of smooth joint *)
(* Returned vector DetrendedSignal *)
output
]


(**********************************************************************************************************************)
(* Accelerated version with subfunctions substituted and comments removed *)
EstimatePSD[Signal_,MinSamples_:3,MaxSamples_:11,HalfOverlapSamples_:30,HalfJoinSamples_:5]:=
Module[{DetrendedSignal,nSignal, nFreqs,MinLateralSamples,MaxLateralSamples,nLateralSamples,factor, f1, f2,df,PSDfactor, RunningAverageHalfLegthScale, AutoPowerSpectrum, EstimatedPSD  ,geometricRatio, RunningAverage, iStart,iEnd, spectrumphase},

(*DetrendedSignal=DetrendSignal[Signal];*)DetrendedSignal=LinearModelFit[Signal,
(* SignalTrendModel *)  x (* {x,x^2,x^3,x^4,x^5} or x; *),
(* Variable in Signal Trend *)x
][ (* Returned vector DetrendedSignal *)
"FitResiduals"];

(*DetrendedSignal=SmoothCyclicRepetition[DetrendedSignal, 
(* Round[(* fblade in samples *)   0.85/\[Tau]] *) HalfOverlapSamples ,(* Round[(* fmax in joint, in samples *) 0.02/\[Tau]] *) HalfJoinSamples 
];*)
Module[{a,x1,x2,CyclicJoint,WindowedVector, residualVector,SmoothCyclicRepetitionWeithing},
(* Make a soft transition to smooth out the replication of the signal in the next period *)
(* This is done to full use the data, specially in short samples where it is desired to full use the data *)
(* SmoothCyclicRepetitionWeithing with a ramp *)
SmoothCyclicRepetitionWeithing[x_]=Simplify[ (1+Sin[ \[Pi] x /(2HalfOverlapSamples) ])/2];
nSignal = Dimensions[DetrendedSignal][[1]];
(* make a weighted transition at the start/end *)
(* this wieghting does not introduce a bias in the PSD
since \[CapitalSigma]w^2 = N *)
Do[
a =SmoothCyclicRepetitionWeithing[i];
x1 = a*DetrendedSignal[[i+1]]+(1-a)*DetrendedSignal[[nSignal-i]];
x2 =(1- a)*DetrendedSignal[[i+1]]+a*DetrendedSignal[[nSignal-i]];
DetrendedSignal[[i+1]] = x1;
DetrendedSignal[[nSignal-i]] = x2;
,{i,0,HalfOverlapSamples}
];
If[HalfJoinSamples>=3, (* smooth joint *)
(* To avoid very sharp derivatives in junction, smooth out some samples in the joint *)
WindowedVector = Join[DetrendedSignal[[(nSignal-HalfJoinSamples);;nSignal]],DetrendedSignal[[1;;HalfJoinSamples+1]] ];
residualVector=LinearModelFit[WindowedVector,{x,x^2,x^3},x]["FitResiduals"];
(* Now fit a smoother version in the HalfJoinSamples *)
Do[a =Sin[\[Pi] i /(2 HalfJoinSamples)]^2; (* Windowing weighting *)
DetrendedSignal[[i+1]] =WindowedVector[[HalfJoinSamples+i+2]] -(1-a)* residualVector[[HalfJoinSamples+i+2]];
DetrendedSignal[[nSignal-i]] =WindowedVector[[HalfJoinSamples-i+1]]-(1-a)* residualVector[[HalfJoinSamples-i+1]];
,{i,0,Round[HalfJoinSamples]}
];
] (* end of smooth joint *)
]; (* End of SmoothCyclicRepetition module *)

(* Compute AutoPowerSpectrum from signal *)
nSignal = Dimensions[DetrendedSignal][[1]];
nFreqs = Floor[(nSignal-1)/2];
EstimatedPSD=Fourier[DetrendedSignal][[2;;nFreqs+1]];
MaxLateralSamples = Round[nFreqs/16];
MinLateralSamples = Min[MaxLateralSamples,Max[0,Round[(MinSamples-1)/2]]];
MaxLateralSamples = Min[MaxLateralSamples,Max[MinLateralSamples,Round[(MaxSamples-1)/2]]];
(* Use only nonzero frequencies up to Nyquist Freq, neglect DC term *)
(* Measure one-sided PSD+ in RMS value, then use a 2 factor *)
(* The cyclic replication does not generate bias in Fourier transform because the starting and ending weithing factors sum unity, but it introduces the bias PSDbiasDueToSmoothCyclicRepetition = 1-HalfOverlapSamples/(2 nSignal) in PSD. Therefore, this bias is corrected *)
(* PSDbiasDueToSmoothCyclicRepetition =1-HalfOverlapSamples/(2 nSignal) *)

PSDfactor = N[1/(1-HalfOverlapSamples/(2 nSignal))];


AutoPowerSpectrum = (Re[EstimatedPSD])^2+(Im[EstimatedPSD ])^2;
spectrumphase = Arg[EstimatedPSD];


(*(* For debugging *) AutoPowerSpectrum = Table[2,{i,1,nFreqs}];*)

If[MaxLateralSamples > 0,
geometricRatio =Quiet@N[ ((nFreqs-MaxLateralSamples-1)/(MinLateralSamples+3))^(1/(MaxLateralSamples-MinLateralSamples))];
If[!NumberQ[geometricRatio],geometricRatio = nFreqs];
(*
(* First samples -> negative frequencies corresponds to the reflected spectrum*)
RunningAverage = Sum[ AutoPowerSpectrum[[i]],{i,2,MinLateralSamples+2}]+ Sum[ AutoPowerSpectrum[[i]],{i,1,MinLateralSamples-1}];
(*If[debugging,Print["iStart -> ",1,"\tiEnd -> ",MinLateralSamples+2,"\t","nLateralSamples -> ",MinLateralSamples,"\t",{i,-(MinLateralSamples-1),MinLateralSamples-1}];
];*)
*)
If[MinLateralSamples <= 0,
(* Since MinLateralSamples == 0, there is no averaging and data is passed without modifications *)
(* only auxiliary variables are modified to cope with the unmodified start since original data is copyed in filtered signal *)
iStart = 1;
iEnd = geometricRatio;
iEnd = Min[Ceiling[iEnd],nFreqs];
(* Just scale the data not averaged *)
Do[
EstimatedPSD[[i]] = AutoPowerSpectrum[[i]]*PSDfactor;
,{i,iStart,iEnd}
]; (* End of Do *)
(* The rest of the samples (if it is not reached the end of data) are averaged! *)
(* Compute the rest of the samples averaging some data *)
MinLateralSamples = 1;
, (* else MinLateralSamples > 0 and first samples must be computed from the only available data *)
nLateralSamples =MinLateralSamples;
RunningAverage = Sum[ AutoPowerSpectrum[[i]],{i,1,1+nLateralSamples-1}];
Do[
(*RunningAverage +=  AutoPowerSpectrum[[i+nLateralSamples]]-AutoPowerSpectrum[[
(* -(i-nLateralSamples-1) *) nLateralSamples-i+3]];
(*If[debugging, Print["nLateralSamples -> ",nLateralSamples,"\t i -> ",i,"\t",RunningAverage];
];*)
EstimatedPSD[[i]] =RunningAverage*PSDfactor/(2*nLateralSamples+1);*)
RunningAverage +=  AutoPowerSpectrum[[i+nLateralSamples]];
EstimatedPSD[[i]] =RunningAverage*PSDfactor/(i+nLateralSamples);
,{i,iStart=1,iEnd = nLateralSamples}
]; (* End of Do *)
]; (* End of If[MinLateralSamples <= 0,...] *)
(* Medium samples -> Spectrum at frequencies beyond Nyquist is assumed to be null*)
(* iEnd =MinLateralSamples; *)

If[iEnd < nFreqs-MinLateralSamples,

If[(MaxLateralSamples < MinLateralSamples) || (!NumberQ[geometricRatio]),
geometricRatio = (*iEnd/iStart*) (* End of samples *)(nFreqs-MaxLateralSamples) / iEnd; (*  iEnd del paso anterior *)
]; (* If[MaxLateralSamples <= MinLateralSamples, then smooth spectrum evenly with MinLateralSamples... *)

(* Clear[nLateralSamples]; *)
Do[
iStart = iEnd+1;
iEnd *= geometricRatio;
iEnd = Min[Ceiling[iEnd],nFreqs-MaxLateralSamples];
If[iEnd>=iStart,
RunningAverage = Sum[ AutoPowerSpectrum[[i]],{i,iStart-nLateralSamples1-1,iStart+nLateralSamples1-1}];
factor =PSDfactor/(2*nLateralSamples1+1) ;
Do[
RunningAverage += AutoPowerSpectrum[[i+nLateralSamples1]]- AutoPowerSpectrum[[i-nLateralSamples1-1]];
(*RunningAverage = Sum[ AutoPowerSpectrum[[i]],{i,i+1-nLateralSamples1-2,i+1+nLateralSamples1-1}];*)
(*If[debugging, Print["nLateralSamples -> ",nLateralSamples1,"\t i -> ",i,"\t",RunningAverage];
];
*)
EstimatedPSD[[i]] =RunningAverage*factor;
,{i,iStart,iEnd}
]; (* End of Do *)
(*Print[nLateralSamples1,"samples from: ",iStart," to ",iEnd];
*)]; (* End of If[iEnd>=iStart,... *)
(*If[debugging, Print["nLateralSamples -> ",nLateralSamples1,"\t iStart -> ",iStart,"\t","iEnd -> ",iEnd,"\tRunningAverage -> ",RunningAverage];
];*)
,{nLateralSamples1,MinLateralSamples,MaxLateralSamples }
]; (* End of Do *)

]; (* End of If[iEnd < nFreqs-MinLateralSamples,...*)

(* End -> at the end, only available samples are considered for the average*)
iStart =nFreqs-MaxLateralSamples+1;
iEnd = nFreqs;
(*If[debugging, Print["iStart -> ",iStart,"\tiEnd -> ",iEnd,"\t","nLateralSamples -> ",MaxLateralSamples,"\t",{i,iStart-MaxLateralSamples-1,iStart+MaxLateralSamples-1}]; 
];*)
RunningAverage = Sum[ AutoPowerSpectrum[[i]],{i,iStart-MaxLateralSamples-1,iEnd(*iStart+MaxLateralSamples*)}];
nLateralSamples =MaxLateralSamples;
(*factor =PSDfactor/(2*nLateralSamples+1) ;*)
Do[
RunningAverage -=  AutoPowerSpectrum[[i-nLateralSamples-1]]; (*beyond Nyquist frequency the spectrum content is considered null*)
EstimatedPSD[[i]] =RunningAverage*PSDfactor/(iEnd-i+nLateralSamples+1);
(*RunningAverage +=  AutoPowerSpectrum[[2*nFreqs-i-nLateralSamples-1 (*ReflectedFrequency-> nFreqs-(i+nLateralSamples+1-nFreqs)*)]]
-AutoPowerSpectrum[[i-nLateralSamples-1]];
EstimatedPSD[[i]] =RunningAverage*factor;*)
,{i,iStart,iEnd}
]; (* End of Do *)

, (* If MaxLateralSamples == 0, then return AutoPowerSpectrum after applying the PSDfactor *)
EstimatedPSD = AutoPowerSpectrum * PSDfactor;
]; (* End of If[MaxLateralSamples > 0... *)
(*Returned vector: EstimatedPSD and phase *)
{ EstimatedPSD, spectrumphase }
];


(**********************************************************************************************************************)
DecimatedPSD[FullPSD_,MinSamples_:10,\[CapitalDelta]f_:1 (*Default value *),Accumulated_:False]:=
Module[{nFreqs,nSpectralLines,geometricRatio, iStart,iEnd,i,SummarizedPSD},
(* Decimate FullPSD vector geometrically, obtaining a vector of nSpectralLines distributed geometrically  *)
nFreqs =  Dimensions[FullPSD][[1]];
nSpectralLines =  Round[nFreqs/MinSamples];
geometricRatio =N[ nFreqs^(1/(nSpectralLines-1))];
(* Sizing reduced PSD by excess due to Celing upper integer rounding. This rounding has a lot of impact in the very first samples, since geometricRatio is bellow 2 (a typical value is 1,3) 
*) 
SummarizedPSD = Table[{0,0},{Round[nSpectralLines]}];
(* Allways store first and last elements of FullPSD *)
i=iStart =iEnd = 1;
SummarizedPSD[[i]]={ \[CapitalDelta]f,FullPSD[[iStart]]};
While[iEnd<nFreqs ,
iStart = iEnd;
iEnd *= geometricRatio;
iEnd = Min[Ceiling[Max[iEnd,iStart+1]],nFreqs];
(* Store decimated values as the average between (iStart+1) and iEnd *)
SummarizedPSD[[++i]]={(iStart+1+iEnd)/2 \[CapitalDelta]f,Sum[FullPSD[[j]],{j,iStart+1,iEnd}]};
(* Sometimes, the variance between the spectral lines is desired instead of the average value at such frequency*)
(* In such cases, set the Accumulated to True *)
If[Accumulated,SummarizedPSD[[i,2]]/=(iEnd-iStart)];
(*Sum[1,{j,iStart+1,iEnd}]==(iEnd-iStart)*)
]; (* Return this table of two columns {freq, PSD} *)
SummarizedPSD[[1;;i]]
];


LogartithmicRunningAverage[Signal_,MinSamples_:3,MaxSamples_:9]:=
Module[{FilteredSignal,nFreqs,MinLateralSamples,MaxLateralSamples,factor,geometricRatio, RunningAverage, iStart,iEnd},
(* Dimensionize the filtered signal *)
FilteredSignal = Signal;
(* Set auxiliar variables *)
nFreqs = Dimensions[Signal][[1]];
MaxLateralSamples = Round[nFreqs/16];
MinLateralSamples = Min[MaxLateralSamples,Max[0,Round[(MinSamples-1)/2]]];
MaxLateralSamples = Min[MaxLateralSamples,Max[MinLateralSamples,Round[(MaxSamples-1)/2]]];

If[MaxLateralSamples > 0,
geometricRatio =Quiet@N[ ((nFreqs-MaxLateralSamples-1)/(MinLateralSamples+3))^(1/(MaxLateralSamples-MinLateralSamples))];
If[!NumberQ[geometricRatio],geometricRatio = nFreqs];
(* First samples -> negative frequencies corresponds to the reflected spectrum*)
If[MinLateralSamples <= 0,
(* Since MinLateralSamples == 0, there is no averaging and data is passed without modifications *)
(* only auxiliary variables are modified to cope with the unmodified start since original data is copyed in filtered signal *)
iStart = 1;
iEnd = geometricRatio;
iEnd = Min[Ceiling[iEnd],nFreqs];
(* The rest of the samples (if it is not reached the end of data) are averaged! *)
(* Compute the rest of the samples averaging some data *)
MinLateralSamples = 1;
, (* else MinLateralSamples > 0 and first samples must be computed from the only available data *)
RunningAverage = Sum[ Signal[[i]],{i,1,1+MinLateralSamples-1}];
Do[
RunningAverage +=  Signal[[i+MinLateralSamples]];
FilteredSignal[[i]] =RunningAverage/(i+MinLateralSamples);
,{i,iStart=1,iEnd = MinLateralSamples}
]; (* End of Do *)
]; (* End of If[MinLateralSamples <= 0,...*)

If[iEnd < nFreqs-MinLateralSamples,
(* Medium samples -> Spectrum at frequencies beyond Nyquist is assumed to be null*)
(* iEnd =MinLateralSamples; *)
If[(MaxLateralSamples < MinLateralSamples) || (!NumberQ[geometricRatio]),
geometricRatio = (*iEnd/iStart*) (* End of samples *)(nFreqs-MaxLateralSamples) / iEnd; (*  iEnd del paso anterior *)
]; (* If[MaxLateralSamples <= MinLateralSamples, then smooth spectrum evenly with MinLateralSamples... *)
Do[
iStart = iEnd+1;
iEnd *= geometricRatio;
iEnd = Min[Ceiling[iEnd],nFreqs-MaxLateralSamples];
If[iEnd>=iStart,
RunningAverage = Sum[ Signal[[i]],{i,iStart-nLateralSamples1-1,iStart+nLateralSamples1-1}];
factor =1/(2*nLateralSamples1+1) ;
Do[
RunningAverage += Signal[[i+nLateralSamples1]]- Signal[[i-nLateralSamples1-1]];
FilteredSignal[[i]] =RunningAverage*factor;
,{i,iStart,iEnd}
]; (* End of Do *)
(*Print[nLateralSamples1,"samples from: ",iStart," to ",iEnd];
*)]; (* End of If[iEnd>=iStart,... *)
,{nLateralSamples1,MinLateralSamples,MaxLateralSamples }
]; (* End of Do *)
]; (* End of If[iEnd < nFreqs-MinLateralSamples,...*)
(* End -> at the end, only available samples are considered for the average*)
iStart =nFreqs-MaxLateralSamples+1;
iEnd = nFreqs;
RunningAverage = Sum[ Signal[[i]],{i,iStart-MaxLateralSamples-1,iEnd(*iStart+MaxLateralSamples*)}];
(*factor =1/(2*nLateralSamples+1) ;*)
Do[
RunningAverage -=  Signal[[i-MaxLateralSamples-1]]; (*beyond Nyquist frequency the spectrum content is considered null*)
FilteredSignal[[i]] =RunningAverage/(iEnd-i+MaxLateralSamples+2);
,{i,iStart,iEnd}
]; (* End of Do *)
(* If MaxLateralSamples == 0, then return the input signal without modifying it *)
]; (* End of If[MaxLateralSamples > 0... *)
(* Return the filtered signal *)
FilteredSignal
]


Spectrogram[data_,opts:OptionsPattern[{WindowType-> Automatic ,WindowLength-> Automatic,WindowOverlap-> Automatic,LowFreqResolution-> Automatic,HighFreqResolution-> Automatic,CiclicRepetitionOverlap-> Automatic,CiclicRepetitionJoinSamples-> Automatic}]] := 
Module[{spect,windowquantiles,periodogram, signaltimeduration,tmp, n, n2, s, win, windowed,MinSamples,MaxSamples,HalfOverlapSamples,HalfJoinSamples,spectphase},
(* opts2 = opts;*)
tmp = WindowLength /. Flatten[{opts, Options[Spectrogram]}];
n = If[tmp === Automatic, 2^Round[Log[2, Sqrt[Length[data]]] + 1], tmp];
n2 = IntegerPart[n/2];
tmp = WindowOverlap /. Flatten[{opts, Options[Spectrogram]}];
(* Allowed overlaps are between 0% and 100% (less 8 samples to avoid infinity data chunks) with a default value of 50% *)
s =Min[n, Round@If[tmp === Automatic,n2, Max[n (1- tmp),8]]];
tmp = WindowType /. Flatten[{opts, Options[Spectrogram]}];
(* Scaled window to make the variance unmodified *)
If[tmp === Automatic,
windowed = False; win = 1;, 
windowed =True;  win = Table[(*0.54/0.3974*)1.3588324+ (*0.46/0.3974*)1.1575239*Cos[(2.*Pi)/n*i], {i, -(n - 1)/2, (n - 1)/2}];
];
tmp = LowFreqResolution /. Flatten[{opts, Options[Spectrogram]}];
tmp =If[tmp === Automatic,n/500 , tmp];
MinSamples=Round@Max[1,tmp];
tmp =HighFreqResolution /. Flatten[{opts, Options[Spectrogram]}];
tmp = If[tmp === Automatic,n/200, tmp];
MaxSamples=Round@Max[MinSamples+0,tmp];
tmp = CiclicRepetitionOverlap /. Flatten[{opts, Options[Spectrogram]}];
HalfOverlapSamples = If[tmp === Automatic,Max[5,Round[n/20]], tmp];
tmp = CiclicRepetitionJoinSamples /. Flatten[{opts, Options[Spectrogram]}];
HalfJoinSamples = If[tmp === Automatic,Max[3,Round[n/50]], tmp];
spect = EstimatePSD[If[windowed,#,win*#],MinSamples,MaxSamples,HalfOverlapSamples,HalfJoinSamples]& /@ 
             Partition[data, n, s ];
periodogram = Mean[spect[[All,1,All]]];
(* windowquantiles =((If[windowed,Total[#],Total[win*#]]&) /@ Partition[data, n, s])/ If[windowed, 1.3588324 n, n];
(*
windowquantiles =((If[windowed,Total[#],Total[win*#]]&) /@ Partition[data, n, s])/ If[windowed, 1.3588324 n, n];
*)
*)

(* WindowQuantiles is a matrix containing the median (50% quantile) of each data chunk in the first row and the quantiles 25%, 75%,5%, 95% in the following rows.*)
(* Compute quantiles using linear interpolation (hydrologist method) for better accuracy on small chunks. Deactivated due to extrapollation with wind nulls, giving negative values *)
windowquantiles = Transpose[(Quantile[#,{1/2,1/4,3/4,1/20,19/20}(*,{{1/2,0},{0,1}}*)]&) /@ Partition[data, n, s]];
(* Return values *)
{spect[[All,1,All]] (* auto power spectrum modulus *),periodogram,windowquantiles,spect[[All,2,All]] (* spectrum phase *)}
];


NumberToString[x_,precision_:5]:= Module[
{exponent ,  mantissa, effprecision },
(* Rest one space due to the minus sign of the mantissa if applicable *)
If[Abs[x]<10^-99, 
"0.0",
exponent = Floor[Log10[Abs[x]]];
effprecision =precision-If[x<0,1,0];
(* Rest one space due to E character if applicable *)
effprecision -=If[-4<exponent< effprecision, 
(* Number X, 0.X (exponent == -1) or 0.0X (exponent == -2)*)
If [-4<exponent< 0, (* number of zeros*) 1-exponent, (* number without 0.x or mantissaEexponent *)1] ,
(* number of type mantissa E exponent *)
(* E letter*) 1 + (* exponent sign *) If [exponent<0,1,0]+(* exponent digits *) Floor[1+Log10[Max[1,Abs[exponent]]]] ] ;
(* at least, write two significant digit*)
effprecision = Max[effprecision,1];
If[-4<exponent<= effprecision,
mantissa = N@Round[x ,10^Round[exponent-effprecision]];
 ToString[mantissa(* If NumberForm[mantissa, effprecision] is used, sometimes Mathematica adds 0's as in 1.000E10*)]
 , 
mantissa = N@Round[x*10^(-exponent),10^Round[-effprecision]];
ToString[NumberForm[mantissa ,effprecision]] <> "E"<>ToString[exponent]]
]
];


Needs["PlotLegends`"];
PlotSpectrogram2[spect_,periodogram_:Automatic,windowquantiles_:Automatic,opts:OptionsPattern[{SampleRate->Automatic,SpectrogramType->Automatic,TimeUnits->Automatic,FrequencyUnits->Automatic,SpectrogramColors->Automatic,DynamicRange->Automatic,SignalName->Automatic,PeriodogramUnits->Automatic,DateStart->Automatic,DateEnd->Automatic,SpectogramData->Automatic,PeriodogramData->Automatic,WindowQuantilesData->Automatic,SignalTimeDuration->Automatic,TextExplanation->Automatic,LogPeriodogramPlot->Automatic,MultipliedByFrequency->Automatic,PeridogramLabel->Automatic,ShadowInterQuantiles->Automatic,TimePlotOptions->Automatic, SpectrogramLimitFactor-> Automatic}]]:=
Module[{signaltimeduration,srate,n,tmp,freqs,timeunits,frequencyunits,linearcolormap,logcolormap,colorize,fpsd,linmapping,logmapping,nfreqs,logaritmicspectrogram,grading,logarithmicstep,dynamicrange,fperiodogram,periodogramdata,ntimes,ListDensityPlotOptions,PeriodogramPlotOptions,TimePlotOptions,ShowLegendOptions,TimeseriesLabel,periodogramunit,startdate,enddate,stringspaces,textexplanation,logperiodogramplot,multipliedbyfrequency,peridogramlabel,spectrogramtype,shadowinterquantiles,additionaltimeplotoptions,spectrogramlimitfactor},
nfreqs=Dimensions[spect][[2]];
n=2(nfreqs+1);
ntimes=Dimensions[spect][[1]];
tmp=TimeUnits/.Flatten[{opts,Options[PlotSpectrogram]}];
timeunits=If[tmp===Automatic,"s",tmp];
tmp=FrequencyUnits/.Flatten[{opts,Options[PlotSpectrogram]}];
frequencyunits=If[tmp===Automatic,"Hz",tmp];
srate=SampleRate/.Flatten[{opts,Options[PlotSpectrogram]}];

spectrogramtype=SpectrogramType/.Flatten[{opts,Options[PlotSpectrogram]}];
tmp = !((spectrogramtype === Automatic) || (spectrogramtype === "Logarithmic"));
If[srate===Automatic,
srate=1;
ListDensityPlotOptions={(*Frame->False, *)ImagePadding-> None,Frame->True};
PeriodogramPlotOptions={Frame->{False,False},GridLines->None};
TimePlotOptions={Frame->{False,False},GridLines->None};
,
ListDensityPlotOptions=If[tmp,{Frame->tmp},{ImagePadding-> None,Frame->True}];
PeriodogramPlotOptions={Frame->{True,True}};
TimePlotOptions={Frame->{True,True}};
];
signaltimeduration=SignalTimeDuration/.Flatten[{opts,Options[PlotSpectrogram2]}];
If[signaltimeduration===Automatic,signaltimeduration=Length[spect]-1];
grading=DynamicRange/.Flatten[{opts,Options[PlotSpectrogram]}];
(*dynamicrange=If[NumberQ[grading],N[10^(-Abs[grading]/10)],10^-3];*)
dynamicrange=Evaluate@If[Evaluate@NumberQ[grading],Evaluate@N[10^(-Abs[grading]/10)],10^-3];
linearcolormap=Evaluate@Switch[
Evaluate[SpectrogramColors/.Flatten[{opts,Options[PlotSpectrogram2]}]] 
,Automatic,
Evaluate[(ColorHueGray[0.55-2#1,Min[1,0.25+1.15Abs[0.4+Cos[\[Pi] #1]]],1-0.75#1]&)]
,"Lighter",
Evaluate[(ColorHueGray[0.2333333333333334+0.6666666666666666*#1,Min[1,0.25+1.15Abs[-0.4+Cos[2.670353755551324-\[Pi] #1]]],Min[1,1.05-0.6*#1]]&)]
,"Darker",
Evaluate[(ColorHueGray[0.55-2(1-#1),Min[1,0.25+1.15Abs[Cos[\[Pi] #1]-0.4]],0.25+0.75#1]&)]
,_,
Evaluate[(ColorHueGray[(1.35-#1)4/6,Min[1,0.25+1.15Abs[Cos[\[Pi] (#1-0.15)]-0.4]],Min[1,0.45+0.60#1]]&)]
];
(*Print[linearcolormap];
*)
(*
linearcolormap:=Switch[SpectrogramColors/.Flatten[{opts,Options[PlotSpectrogram]}]
,Automatic,
(ColorHueGray[0.55-2#1,Min[1,0.25+1.15Abs[0.4+Cos[\[Pi] #1]]],1-0.75#1]&)
,"Lighter",
(ColorHueGray[0.2333333333333334+0.6666666666666666*#1,Min[1,0.25+1.15Abs[-0.4+Cos[2.670353755551324-\[Pi] #1]]],Min[1,1.05-0.6*#1]]&)
,"Darker",
(ColorHueGray[0.55-2(1-#1),Min[1,0.25+1.15Abs[Cos[\[Pi] #1]-0.4]],0.25+0.75#1]&)
,_,
(ColorHueGray[(1.35-#1)4/6,Min[1,0.25+1.15Abs[Cos[\[Pi] (#1-0.15)]-0.4]],Min[1,0.45+0.60#1]]&)
];
*)
(*Print[dynamicrange];
*)(*logcolormap=((ReleaseHold[linearcolormap[1-Evaluate@Log[Evaluate@Min[1,Max[Evaluate@dynamicrange,#1]]]/Evaluate@Log[dynamicrange]]])&);
*)
(*logcolormap=Expand@(Evaluate[linearcolormap[1-Evaluate@Log[Evaluate@Min[1,Evaluate@Max[Evaluate[dynamicrange],#1]]]/Evaluate@Log[Evaluate[dynamicrange]]]]&);*)
(*logcolormap=linearcolormap[1-Evaluate@Log[Evaluate@Min[1,Evaluate@Max[Evaluate[dynamicrange],#1]]]/Evaluate@Log[Evaluate[dynamicrange]]]]&);*)
logcolormap= (linearcolormap[(1-Log[Min[1,Max[dynamicrange,#1]]]/Log[dynamicrange])])&;
(*logcolormap[x_]= linearcolormap[(1-Log[Min[1,Max[dynamicrange,x]]]/Log[dynamicrange])];*)
(*Print[logcolormap];*)

colorize=If[grading===Automatic,linearcolormap,logcolormap];
startdate=DateStart/.Flatten[{opts,Options[PlotSpectrogram]}];
enddate=DateEnd/.Flatten[{opts,Options[PlotSpectrogram]}];
freqs=Table[i N[srate/n],{i,1,nfreqs}];
logperiodogramplot=(LogPeriodogramPlot/.Flatten[{opts,Options[PlotSpectrogram]}]);
If[!(logperiodogramplot === True),logperiodogramplot = False];
textexplanation=(TextExplanation/.Flatten[{opts,Options[PlotSpectrogram]}]);
If[textexplanation===Automatic,textexplanation=Text[Style[If[spectrogramtype=="Logarithmic",
(* The first line is \ [Letterspace] or Esc_Esc characters with a point in the start and in the end *)
(* This line has been added to separate legend from periodogram and to dimension *)
(* the width of text, that is computed from the first line of a multiline string *)
(* The number of spaces has been adjusted experimentally from try and error since *)
(* Text uses non-monospace text. If conventional spaces are used, a starting and ending *)
(* point must be included and if the line is divided, Mathematica thows an error *)
(* If points are not included, Mathematica throws an error when computing text width *)
(* Finally, non-divisible spaces (that are considered as leters by Mathematica) are used *)
(* to dimension the width of the text *)
(* This is quite surprising, but this is an operative solution for Mathematica 7!!! *)
"\[LetterSpace]\[LetterSpace]\[LetterSpace]\[LetterSpace]\[LetterSpace]\[LetterSpace]\[LetterSpace]\[LetterSpace]\[LetterSpace]\[LetterSpace]\[LetterSpace]\[LetterSpace]\[LetterSpace]\[LetterSpace]\[LetterSpace]\[LetterSpace]\[LetterSpace]\[LetterSpace]\[LetterSpace]\[LetterSpace]\[LetterSpace]\[LetterSpace]\[LetterSpace]\[LetterSpace]\[LetterSpace]\[LetterSpace]\[LetterSpace]\[LetterSpace]\[LetterSpace]\[LetterSpace]\[LetterSpace]
Periodogram or f*PSD(f) at top left \[UpArrow],
Spectrogram or f*PSD(f,t) at top right \[UpperRightArrow]
and signal at bottom right ->
(quantiles 5%, 25%, 50%, 75%, 95%).

The shadowed area of the periodogram
and the darkness of the spectrogram 
represent the variance of the signal
(the PSD is multiplied by f in plots
with logarithmic frequency axis).",
"\[LetterSpace]\[LetterSpace]\[LetterSpace]\[LetterSpace]\[LetterSpace]\[LetterSpace]\[LetterSpace]\[LetterSpace]\[LetterSpace]\[LetterSpace]\[LetterSpace]\[LetterSpace]\[LetterSpace]\[LetterSpace]\[LetterSpace]\[LetterSpace]\[LetterSpace]\[LetterSpace]\[LetterSpace]\[LetterSpace]\[LetterSpace]\[LetterSpace]\[LetterSpace]\[LetterSpace]\[LetterSpace]\[LetterSpace]\[LetterSpace]\[LetterSpace]\[LetterSpace]\[LetterSpace]\[LetterSpace]
Periodogram or PSD(f) at top left \[UpArrow],
Spectrogram or PSD(f,t) at top right \[UpperRightArrow]
and signal quartiles at bottom right ->
(quantiles 5%, 25%, 50%, 75% and 95%).

The shadowed area of the periodogram
and the darkness of the spectrogram 
represent the variance of the signal."
],FormatType->Text,TextAlignment->Center]]
];

shadowinterquantiles = (ShadowInterQuantiles/.Flatten[{opts,Options[PlotSpectrogram]}]);
If[shadowinterquantiles===Automatic, shadowinterquantiles = True, If[!(shadowinterquantiles=== True),shadowinterquantiles= False]];

TimeseriesLabel=(SignalName/.Flatten[{opts,Options[PlotSpectrogram]}]);
If[tmp===Automatic,TimeseriesLabel=If[shadowinterquantiles === Automatic,"Signal(t) excursion","Signal(t)",tmp]];

tmp=(PeriodogramUnits/.Flatten[{opts,Options[PlotSpectrogram]}]);
periodogramunit=If[tmp===Automatic,""," "<>tmp];
tmp=StringLength[periodogramunit];
stringspaces=StringJoin[Table[" ",{10+2*If[tmp>27,tmp-27,tmp]}]];

multipliedbyfrequency=(MultipliedByFrequency/.Flatten[{opts,Options[PlotSpectrogram]}]);
If[multipliedbyfrequency===Automatic,multipliedbyfrequency=If[spectrogramtype==="Logarithmic",True,False]];

peridogramlabel=(PeridogramLabel/.Flatten[{opts,Options[PlotSpectrogram]}]);
If[peridogramlabel===Automatic,peridogramlabel=If[multipliedbyfrequency,"f\[CenterDot]PSD(f)","PSD(f)"]<>periodogramunit];

spectrogramlimitfactor= (SpectrogramLimitFactor/.Flatten[{opts,Options[PlotSpectrogram]}]);
If[spectrogramlimitfactor===Automatic,spectrogramlimitfactor=2, If[!NumberQ[spectrogramlimitfactor],spectrogramlimitfactor=2]];

periodogramdata=If[(periodogram===Automatic),Mean[spect],periodogram];
If[multipliedbyfrequency,periodogramdata *= freqs]; 

Module[{minrange = Min[periodogramdata], maxrange = Max[periodogramdata],minmaxrange, ShadowFilling},
minmaxrange = (maxrange-minrange);

If[ (minrange > 0) && (maxrange > 0) && (spectrogramtype === "Logarithmic"), 
(* For only positivie quantities  and logarithmic scales, use geometric scaling *)
ShadowFilling = minrange;
minrange /= spectrogramlimitfactor;
maxrange *= spectrogramlimitfactor;
,
(* Else, use linear scaling *)
 ShadowFilling = 0.0;
 minrange -= (spectrogramlimitfactor-1)*minmaxrange;
 maxrange += (spectrogramlimitfactor-1)*minmaxrange;
];

periodogramdata=Transpose@{Join[{ShadowFilling},periodogramdata,{ShadowFilling}],Join[{0.999Max[freqs[[1]],freqs[[1]]-(freqs[[nfreqs]]-freqs[[1]])/1000]},freqs,{Min[1.001freqs[[nfreqs]],freqs[[nfreqs]]+(freqs[[nfreqs]]-freqs[[1]])/1000]}]};

ListDensityPlotOptions=Join[ListDensityPlotOptions,{Background->None,ColorFunction->colorize,DataRange->{{0,signaltimeduration},{N[srate/n],N[srate/2]}},PlotRange->{minrange,maxrange},FrameLabel->{"Time ["<>timeunits<>"]","Frequency ["<>frequencyunits<>"]"},GridLines->None,Mesh->False,RotateLabel->True,MeshFunctions->{#1&,#2&,#3&},Mesh->{5,5,9},MeshStyle->{LightGray,LightGray,Darker@Brown},ClippingStyle->Automatic,ImageSize->{390,390},ImageMargins->{{12,0},{0,0}}}];

PeriodogramPlotOptions=Join[PeriodogramPlotOptions,{Background->None,PlotRange->All,FrameLabel->{peridogramlabel,"Frequency ("<>frequencyunits<>")"},Filling->Bottom,FillingStyle->Automatic,Axes->False,ImageSize->{{200},{450}},AspectRatio->450/200,RotateLabel->True,ClippingStyle->Automatic,Joined->True,ImageMargins->{{0,0},{0,16(*4*)}}}];
(* WindowQuantiles is a matrix containing the median (50% quantile) of each data chunk in the first row and the quantiles 25%, 75%,5%, 95% in the following rows.*)
(* It also works if WindowQuantiles is a vector or a matrix of only one row because ListLinePlot and DateListPlot accept those inputs*)
additionaltimeplotoptions=(TimePlotOptions/.Flatten[{opts,Options[PlotSpectrogram]}]);
If[additionaltimeplotoptions===Automatic,additionaltimeplotoptions={}];
TimePlotOptions=Join[additionaltimeplotoptions,TimePlotOptions,{DataRange->{0,signaltimeduration},FrameLabel->{"Time ["<>timeunits<>"]",TimeseriesLabel},
Filling->If[shadowinterquantiles,Switch[Length[windowquantiles],3,{2-> {3}},5,{2-> {3}, 4-> {5}},6,{2-> {{3},Directive[Red,Opacity[0.2]]}, 5-> {{6},Directive[Blue,Opacity[0.2]]}},7,{2-> {3}, 4-> {5}, 6-> {7}},_ ,Bottom],None],FillingStyle->Automatic,
PlotStyle->If[shadowinterquantiles,If[Length[windowquantiles]==6,{{Red,Thick,Dashing[None]},{Red,Thin,Dotted},{Red,Thin,Dotted},{Blue,Thick,Dashing[None]},{Blue,Thin,DotDashed},{Blue,Thin,DotDashed},{Hue[-(1/6)],Thin,DotDashed}},{{Hue[0],Thick,Dashing[None]},{Hue[-(2/6)],Thin,Dotted},{Hue[-(3/6)],Thin,Dotted},{Hue[-(1/6)],Thin,Dotted},{Hue[-(4/6)],Thin,Dotted},{Hue[-1],Thick,Dashing[None]},{Hue[-(2/6)],Thin,DotDashed},{Hue[-(3/6)],Thin,DotDashed},{Hue[-(1/6)],Thin,DotDashed},{Hue[-(4/6)],Thin,DotDashed},{Hue[-(5/6)],Thick,Dashing[None]}}],Automatic],
Axes->False,ImageSize->{{470},{200}},AspectRatio->200/460,RotateLabel->True,ClippingStyle->Automatic,Joined->True,ImageMargins->{{0,30},{0,0}}}(*{{0,20},{0,0}}}*)];
ShowLegendOptions={RGBColor[colorize[1-#1]]&,50,StringJoin[NumberToString[maxrange],periodogramunit,stringspaces],StringJoin[stringspaces,NumberToString[minrange],periodogramunit],LegendPosition->{1.04,-0.7},LegendShadow->None,LegendBorder->None,LegendSize->{0.15,1.4},LegendTextDirection->{0,1},LegendTextOffset->{0,+1},LegendBorderSpace->0,LegendLabelSpace->Automatic,ImageMargins->0,Framed->None,Background->None,ShadowBackground->Transparent,ShadowForeground->Transparent,ShadowOffset->0};

];
If[(spectrogramtype === Automatic) || (spectrogramtype === "Logarithmic"),
Grid[{{
If[(periodogram === Automatic),"",
If[(spectrogramtype === Automatic),
If[logperiodogramplot,
(* "Logarithmic" magnitude axis*) ListLogLinearPlot[periodogramdata,PeriodogramPlotOptions],
(* both axis are linear *) ListLinePlot[periodogramdata,PeriodogramPlotOptions]]
, (* "Logarithmic" frequency axis and magnitude *)
If[logperiodogramplot,ListLogLogPlot[periodogramdata,ImageMargins->{{0,0},{0,0}},PeriodogramPlotOptions],ListLogPlot[periodogramdata,PeriodogramPlotOptions]]
]]
,
ShowLegend[
If[(spectrogramtype === Automatic)
,(* Linear frequency axis *)
ListDensityPlot[Transpose[spect],FilterRules[{opts},Options[ListDensityPlot]],ListDensityPlotOptions]
, (* "Logarithmic" frequency axis*)
linmapping=Log[freqs];logarithmicstep=(linmapping[[nfreqs]]-linmapping[[1]])/Min[8nfreqs,300];
logmapping=Range[linmapping[[1]],linmapping[[nfreqs]],logarithmicstep];
logaritmicspectrogram=Table[
fpsd=ListInterpolation[If[multipliedbyfrequency, freqs*spect[[i]],spect[[i]]],{linmapping}];
fpsd/@logmapping,
{i,1,ntimes}];
ListDensityPlot[Transpose@logaritmicspectrogram,FilterRules[{opts},Options[ListDensityPlot]],ListDensityPlotOptions]
]
,ShowLegendOptions]
},
If[(windowquantiles===Automatic),{"",""},
{textexplanation,
(* Compute variance of each data chunk from spectrogram *)
If[(startdate===Automatic)||(enddate===Automatic),
Quiet@ListLinePlot[windowquantiles,TimePlotOptions],
Quiet@DateListPlot[windowquantiles,{startdate,enddate},TimePlotOptions]
]
}]
},Spacings->{0,0},Spacings->0,Alignment->{Right,Top}]
, (* Other types of spectrogramtype apart from Automatic and "Logarithmic"*)
ShowLegend[
ListDensityPlot[Transpose[spect],Frame->True,ImageMargins->0,FilterRules[{opts},Options[ListDensityPlot]],ListDensityPlotOptions],ShowLegendOptions]
] (* End of If[(spectrogramtype === Automatic) || (spectrogramtype === "Logarithmic"),... *)
];


PlotSpectrogram[data_,opts:OptionsPattern[ {SampleRate-> Automatic,SpectrogramType -> Automatic, TimeUnits-> Automatic, FrequencyUnits-> Automatic,SpectrogramColors-> Automatic,DynamicRange-> Automatic,SignalName -> Automatic, PeriodogramUnits -> Automatic, DateStart -> Automatic, DateEnd -> Automatic,SpectogramData -> Automatic, PeriodogramData -> Automatic,WindowQuantilesData -> Automatic, SignalTimeDuration-> Automatic,TextExplanation -> Automatic,LogPeriodogramPlot -> Automatic,  MultipliedByFrequency -> Automatic, PeridogramLabel -> Automatic,ShadowInterQuantiles->Automatic,TimePlotOptions->Automatic,SpectrogramLimitFactor-> Automatic}]] := 
Module[{spect,periodogram,windowquantiles,srate,spectrumphase},
spect = SpectogramData /. Flatten[{opts, Options[PlotSpectrogram]}];
periodogram = PeriodogramData /. Flatten[{opts, Options[PlotSpectrogram]}];
windowquantiles = WindowQuantilesData /. Flatten[{opts, Options[PlotSpectrogram]}];

(* Compute spectogram if it is not supplied *)
If[spect === Automatic,{spect,periodogram,windowquantiles,spectrumphase} = Spectrogram[data,FilterRules[{opts},Options[Spectrogram]]];];

srate = SampleRate /. Flatten[{opts, Options[PlotSpectrogram]}];
srate = If[srate === Automatic,1,srate];

(* Calls the function with the spectogram computed *)
PlotSpectrogram2[spect,periodogram,windowquantiles,opts,SignalTimeDuration -> N[(Length[data]-1)/srate] ]
]; (* End of PlotSpectrogram *)


(* Compute vector (x2[[i]]/x1[[i]]) and matrix (x2[[i,j]]/x1[[i,j]]) elemen-wise division avoiding zero division errors and outliers due to divisions by a number close to zero*)
SafeDivide[x1_,x2_,maxratio_: (10^12)]:=
If[ MatrixQ[x1] && MatrixQ[x2]
,
(* matrix element-wise division *)
Table[If[Abs[x2[[i,j]]]>=maxratio Abs[x1[[i,j]]],maxratio,x2[[i,j]]/x1[[i,j]]],{i,1,Min[Length[x1],Length[x2]]},{j,1,Min[Length[x1[[i]]],Length[x2[[i]]]]}]
,
(* default: vector element-wise division *)Table[If[Abs[x2[[i]]]>=maxratio Abs[x1[[i]]],maxratio,x2[[i]]/x1[[i]]],{i,1,Min[Length[x1],Length[x1]]}]
];

IntegratePhase[x_,x0_:0] :=
(* Integrate vector x, assumming inital value x0 *)
If[ MatrixQ[x,NumberQ],
(* divide matrix in row vectors *)
Table[IntegratePhase[x[[i]],x0[[i]]],{i,1,Length[x]}]
,
If[ VectorQ[x,NumberQ],
(* vector, assumed to be radians between -pi/+pi *)
Module[{former= x0},
Table[
former += x[[i]]
,{i,1,Length[x]}]
]
,
(* default: return the input. *)
x
]
];

UnwrapLeastSquares[x_ ] :=
(* gap must be greater than average difference between x elements in order to avoid detecting fake 2 pi phase jumps *)
If[ MatrixQ[x,NumberQ],
(* divide matrix in row vectors *)
Table[UnwrapLeastSquares[x[[i]]],{i,1,Length[x]}]
,
If[ VectorQ[x,NumberQ],
Module[{ now=0,previous1 =0,previous2 =0,next = x[[2]]-x[[1]], next2 = x[[3]]-x[[2]],next3 = x[[4]]-x[[3]],n=0, imax=Length[x],fiterror,fiterrorplus,fiterrorminus,jumpplus= False, jumpminus = False},
(* vector, assumed to be radians between -pi/+pi *)
Table[
(* average difference based on 2, 4 and 6 points *)
(*
{mean, slope} = LeastSquares[{{1, -2}, {1, -1}, {1, 0}, {1, 1}, {1, 2}, {1, 3}}, data = {previous2, previous1, now, next, next2, next3}];
nextfiterror = Simplify[Total[ (data - (mean + slope {-2, -1, 0, 1, 2, 3}))^2]];
{mean, slope} = LeastSquares[{{1, -2}, {1, -1}, {1, 0}, {1, 1}, {1, 2}, {1, 3}}, data = {previous2, previous1, now + 2 \[Pi], next, next2, next3}];
nextfiterrorplus = Simplify[Total[ (data - (mean + slope {-2, -1, 0, 1, 2, 3}))^2]];
Simplify@N@Simplify[nextfiterrorplus - nextfiterror];
{mean, slope} = LeastSquares[{{1, -2}, {1, -1}, {1, 0}, {1, 1}, {1, 2}, {1, 3}}, data = {previous2, previous1, now - 2 \[Pi], next, next2, next3}];
nextfiterrorminus = Simplify[Total[ (data - (mean + slope {-2, -1, 0, 1, 2, 3}))^2]];
Simplify@N@Simplify[nextfiterrorminus - nextfiterror]
*)
fiterror =0.819047619047619` next^2+0.7047619047619047` next2^2+0.47619047619047616` next3^2-0.1904761904761904` next3 now+0.819047619047619` now^2+0.09523809523809523` next3 previous1-0.4190476190476191` now previous1+0.7047619047619048` previous1^2+next (-0.419047619047619` next2-0.4761904761904762` next3-0.3047619047619048` now-0.2476190476190476` previous1-0.19047619047619047` previous2)+next2 (-0.7619047619047619` next3-0.24761904761904763` now-0.07619047619047617` previous1+0.09523809523809525` previous2)+0.38095238095238093` next3 previous2-0.4761904761904761` now previous2-0.7619047619047619` previous1 previous2+0.47619047619047616` previous2^2;
fiterrorplus = fiterror +32.334703942616564` -1.9148755221880644` next-1.5558363617778024` next2-1.1967972013675403` next3+10.292455931760847` now-2.632953843008589` previous1-2.991993003418851` previous2;
fiterrorminus = fiterror +32.334703942616564` +1.9148755221880644` next+1.5558363617778024` next2+1.1967972013675403` next3-10.292455931760847` now+2.632953843008589` previous1+2.991993003418851` previous2;

If[!((fiterror<=fiterrorplus)&&(fiterror<=fiterrorminus) ) ,
If[ (fiterrorplus<=fiterrorminus) && !jumpplus,
(* There is a maximum reduction in phase linear least square fit when guessing a jump now *)
(* Update next values for the following loop *)
now += 2 \[Pi];
n++;
(* phase jump is less than pi except occasionally *)
(* the least square error is the same if there is a constant 2 pi jump in phase *)
(* the jump in the previous loop is tracked to prevent this event *)
jumpplus= True;
jumpminus = False;
,
jumpminus = jumpplus = False;
];
If[ (fiterrorplus>fiterrorminus) && !jumpminus,
now -= 2 \[Pi];
n--;
jumpplus= False;
jumpminus =True ;
,
jumpplus= jumpminus =  False;
];
,
jumpplus= jumpminus = False;
];
(* Cascade update *)
previous2 = previous1;
previous1 = now;
now = next;
next = next2;
next2 = next3;
next3 =x[[Min[i+4,imax]]]-x[[Min[i+3,imax]]];
(* return unwarp phase without jumps computed previously *) 
x[[i]] + 2 \[Pi] n 
,{i,1,Length[x]}]
] (* End Module *)
,
(* default: return the input. *)
x
]
];

PhaseIncrement[x_] :=
(* gap must be greater than average difference between x elements in order to avoid detecting fake 2 pi phase jumps *)
If[ MatrixQ[x,NumberQ],
(* divide matrix in row vectors *)
Table[PhaseIncrement[x[[i]]],{i,1,Length[x]}]
,
If[ VectorQ[x,NumberQ],
(* vector, assumed to be radians between -pi/+pi *)
Module[{former= x[[1]] , difference},
Table[
difference =  x[[i]]-former;
former = x[[i]];
(* return difference = *)
difference
,{i,1,Length[x]}]
]
,
(* default: return the input. *)
x
]
];

UnwrapPhase1[x_, gap_: \[Pi], jump_: 2 \[Pi] ] :=
 (* gap must be greater than average difference between x elements in order to avoid detecting fake 2 pi phase jumps *)
 If[ MatrixQ[x, NumberQ],
  (* divide matrix in row vectors *)
  Table[UnwrapPhase1[x[[i]], gap], {i, 1, Length[x]}]
  ,
If[ VectorQ[x,NumberQ],
(* vector, assumed to be radians between -pi/+pi *)
    Module[{former = x[[1]] , n = 0, difference},
Table[
(* temptative value for unwarped phase *)
difference =  x[[i]] + jump n -former ;
If[difference > gap, n--, 
     If[difference < -gap, n++]
];
(* return unwarp phase without jumps *) 
former = x[[i]] + jump n 
, {i, 1, Length[x]}]
]
,
(* default: return the input. *)
x
]
];

UnwrapPhase2[x_,gap_:\[Pi], jump_: 2 \[Pi] ] :=
(* gap must be greater than average difference between x elements in order to avoid detecting fake 2 pi phase jumps *)
If[ MatrixQ[x,NumberQ],
(* divide matrix in row vectors *)
Table[UnwrapPhase3[x[[i]]],{i,1,Length[x]}]
,
If[ VectorQ[x,NumberQ],
(* vector, assumed to be radians between -pi/+pi *)
Module[{ now= x[[1]] ,previous1 = x[[1]],next = x[[2]], next2 = x[[3]],n=0, imax=Length[x],difference},
Table[
 (*average difference based on 2,4 and 6 points
Simplify[(next-previous+1/2*(next2+next-previous-previous1))/2.]
*)
(* average difference based on 2, 4 and 6 points *)
difference = 0.75` next+0.25` next2-0.75` now-0.25` previous1;
If[difference > gap,
(* There is a reduction in phase difference when guessing a jump *)
n--;
(* Update next values for the following loop *)
now -= jump;
next -= jump;
next2 -= jump;
, 
     If[difference < -gap, 
(* There is a reduction in phase difference when guessing a jump *)
n++;
(* Update next values for the following loop *)
now += jump;
next += jump;
next2 += jump;
];
];
(* Cascade update *)
previous1 = now;
now = next;
next = next2;
next2 = x[[Min[i+3,imax]]] + jump n;
(* return unwarp phase without jumps computed previously *) 
previous1
,{i,1,Length[x]}]
]
,
(* default: return the input. *)
x
]
];

UnwrapPhase3[x_,gap_:\[Pi], jump_: 2 \[Pi] ] :=
(* gap must be greater than average difference between x elements in order to avoid detecting fake 2 pi phase jumps *)
If[ MatrixQ[x,NumberQ],
(* divide matrix in row vectors *)
Table[UnwrapPhase3[x[[i]]],{i,1,Length[x]}]
,
If[ VectorQ[x,NumberQ],
(* vector, assumed to be radians between -pi/+pi *)
Module[{ now= x[[1]] ,previous1 = x[[1]],previous2 = x[[1]],next = x[[2]], next2 = x[[3]],next3 = x[[4]],n=0, imax=Length[x],difference},
Table[
 (*average difference based on 2,4 and 6 points
Simplify[(next-previous+1/2*(next2+next-previous-previous1)+1/3*(next3+next2+next-previous-previous1-previous2))/3.]
*)
(* average difference based on 2, 4 and 6 points *)
difference = 0.611111111111111` next+0.2777777777777778` next2+0.1111111111111111` next3-0.611111111111111` now-0.2777777777777778` previous1-0.1111111111111111` previous2;
If[difference > gap,
(* There is a reduction in phase difference when guessing a jump *)
n--;
(* Update next values for the following loop *)
(*now -= jump;*)
next -= jump;
next2 -= jump;
next3-=jump;
, 
     If[difference < -gap, 
(* There is a reduction in phase difference when guessing a jump *)
n++;
(* Update next values for the following loop *)
(*now += jump;*)
next += jump;
next2 += jump;
next3+=jump;
];
];
(* Cascade update *)
previous2 = previous1;
previous1 = now;
now = next;
next = next2;
next2 = next3;
next3 = x[[Min[i+4,imax]]] + jump n;
(* return unwarp phase without jumps computed previously *) 
previous1
,{i,1,Length[x]}]
]
,
(* default: return the input. *)
x
]
];


UnwrapPhase[phase_]:=UnwrapLeastSquares[UnwrapPhase1[phase]];



MatrixPlotWithScale[matrix_,LevelMatrix_:Automatic,opts:OptionsPattern[{StateLabels->Automatic,Units->Automatic,MinRange-> Automatic,MaxRange-> Automatic,Colors->Automatic,DynamicRange->Automatic, Scaling->Automatic,LimitFactor-> Automatic, ShowColorScale -> True, ShowContours -> True}]]:=
Module[{datarange,linearcolormap,logcolormap,colorize,colormap,logarithmiccolorscale,limitfactor,tmp,grading,dynamicrange,ListDensityPlotOptions,ListContourPlotOptions, ShowLegendOptions,units,stringspaces,minrange , maxrange ,minmaxrange, showlegend },

datarange=StateLabels/.Flatten[{opts,Options[MatrixPlotWithScale]}];
ListDensityPlotOptions=If[datarange===Automatic,
{}(* do not show a value in matrix indexes *),
If[VectorQ[datarange],
(* Assume a vector of centroids, or the maximum and minimum centroid of an equally spaced state discretization *)
{DataRange-> {datarange,datarange}}
,
(* If datarange is not a vector, pass the argument directly to DataRange *)
{DataRange->datarange}
]
];

tmp=(Units/.Flatten[{opts,Options[MatrixPlotWithScale]}]);
units=If[tmp===Automatic,""," "<>tmp];
tmp=StringLength[units];
stringspaces=StringJoin[Table[" ",{10+2*If[tmp>27,tmp-27,tmp]}]];

(* Color Function *)
(* Compute color scale maximums *)
tmp=(MinRange/.Flatten[{opts,Options[MatrixPlotWithScale]}]);
minrange=If[tmp===Automatic, Min[matrix],tmp];
tmp=(MaxRange/.Flatten[{opts,Options[MatrixPlotWithScale]}]);
maxrange=If[tmp===Automatic, Max[matrix],tmp];
minmaxrange = (maxrange-minrange);
grading=DynamicRange/.Flatten[{opts,Options[MatrixPlotWithScale]}];
dynamicrange=If[Evaluate@NumberQ[grading],Evaluate@N[10^(-Abs[grading]/10)],10^-3];

colormap=Evaluate@Switch[
Evaluate[Colors/.Flatten[{opts,Options[MatrixPlotWithScale]}]] 
,Automatic,
Evaluate[(ColorHueGray[0.55-2#1,Min[1,0.25+1.15Abs[0.4+Cos[\[Pi] #1]]],1-0.75#1]&)]
,"Lighter",
Evaluate[(ColorHueGray[0.2333333333333334+0.6666666666666666*#1,Min[1,0.25+1.15Abs[-0.4+Cos[2.670353755551324-\[Pi] #1]]],Min[1,1.05-0.6*#1]]&)]
,"Darker",
Evaluate[(ColorHueGray[0.55-2(1-#1),Min[1,0.25+1.15Abs[Cos[\[Pi] #1]-0.4]],0.25+0.75#1]&)]
,_,
Evaluate[(ColorHueGray[(1.35-#1)4/6,Min[1,0.25+1.15Abs[Cos[\[Pi] (#1-0.15)]-0.4]],Min[1,0.45+0.60#1]]&)]
];

limitfactor= (LimitFactor/.Flatten[{opts,Options[MatrixPlotWithScale]}]);
If[limitfactor===Automatic,limitfactor=1,
If[!NumberQ[limitfactor],limitfactor=1]];

tmp=(Scaling/.Flatten[{opts,Options[MatrixPlotWithScale]}]);
logarithmiccolorscale=If[tmp===Automatic,False,If[tmp==="Logarithmic",True,False]];

If[ (minrange > 0) && (maxrange > 0) && logarithmiccolorscale, 
(* For only positivie quantities  and logarithmic scales, use geometric scaling *)
minrange /= limitfactor;
maxrange *= limitfactor;
,
(* Else, use linear scaling *)
minrange -= (limitfactor-1)*minmaxrange;
maxrange += (limitfactor-1)*minmaxrange;
];

showlegend =  TrueQ[(ShowColorScale/.Flatten[{opts,Options[MatrixPlotWithScale]}])];
(* Bug Issue:  If ColorFunctionScaling-> False, DensityPlotOptions does not use the adequate color scale *)
(* Solution: use allways ColorFunctionScaling-> True even thought the color scale does not much with high accuracy *)
showlegend = False;
linearcolormap =If[ showlegend,
(* ColorFunctionScaling->False *)
(colormap[(*Min[1,Max[0,*)(#1-minrange)/minmaxrange(*]]*)])&,
(* ColorFunctionScaling->True *)
 (colormap[#1])&
];
logcolormap= If[ showlegend,
(* ColorFunctionScaling->False *)
(colormap[(1-Log[Min[maxrange,Max[maxrange*dynamicrange,#1]]]/Log[maxrange])])&,
(* ColorFunctionScaling->True *)
(colormap[(1-Log[Min[1,Max[dynamicrange,#1]]]/Log[dynamicrange])])&
];
colorize=If[logarithmiccolorscale,logcolormap,linearcolormap];

showlegend =  TrueQ[(ShowColorScale/.Flatten[{opts,Options[MatrixPlotWithScale]}])];


(* options for ListContourPlots, which works best with a colormap legend *)
ListDensityPlotOptions= Join[ListDensityPlotOptions,
{
(* Frame -> False, then image looks cooler but FrameLabels and FrameTicks are not shown *)Frame->{{True,True(*False*)},{True,True(*False*)}} ,FrameStyle->Opacity[1] ,LabelStyle->Opacity[1],FrameTicksStyle->Opacity[1],FrameTicks->Automatic(*All*),ImagePadding->All (* if None, no ticks and numbers are shown *),Background->None,RotateLabel->True,ClippingStyle->Automatic,ImageSize->{390,390},ImageMargins->None,(*ImageMargins->{{12,12},{12,12}},*)(*{{12,0},{0,0}}*) FrameMargins-> Tiny
}
];
(* Additional ListContourPlotOptions options if legend is shown *)
ListContourPlotOptions=Join[ {PlotRange->Full,ContourLabels->All,Contours->9,ContourShading-> None,Mesh->False} ,ListDensityPlotOptions];

(* Additional ListDensityPlotOptions options if legend is shown *)
ListDensityPlotOptions= Join[ListDensityPlotOptions,
(* options for ListDensityPlots, which works best with a colormap legend *)
{ ColorFunction->(RGBColor[colorize[#1]]&),ColorFunctionScaling->If[ showlegend,True,True],(*PlotRange->All,*)PlotRange->{Full,Full,{minrange,maxrange}}(*{minrange,maxrange}*)(*PlotRange-> {All,All,{minrange,maxrange}} gives error*),Mesh->False(*MeshFunctions->{#3&} (* Plot contours based on matrix values *),Mesh->{10},MeshStyle->{Darker@Brown}*) (* Contours color *)}
];

ShowLegendOptions={RGBColor[colorize[1-#1]]&,50,StringJoin[NumberToString[maxrange],units,stringspaces],StringJoin[stringspaces,NumberToString[minrange],units],(*ColorFunctionScaling->False,*)LegendPosition->{1.04,-0.7},LegendShadow->None,LegendBorder->None,LegendSize->{0.15,1.4},LegendTextDirection->{0,1},LegendTextOffset->{0,+1},LegendBorderSpace->0,LegendLabelSpace->Automatic,ImageMargins->0,Framed->True(* of Frame -> False, then image looks cooler but FrameLabels are not shown ,FrameStyle->Opacity[0]*),LabelStyle->Opacity[1],FrameTicksStyle->Opacity[1],FrameTicks->All,Background->None,ShadowBackground->Transparent,ShadowForeground->Transparent,ShadowOffset->0,ImageMargins->None};
(*********************************************)
If[ showlegend, 
ShowLegend[ 
(* ArrayPlot has been used instead of MatrixPlot to match colors with ColorFunctionScaling-> True *)
If[TrueQ[(ShowContours/.Flatten[{opts,Options[MatrixPlotWithScale]}])],
Show[ (* show Density + Contours*)
ListDensityPlot[matrix,FilterRules[{opts},Options[ListDensityPlot]],FilterRules[{ListDensityPlotOptions},Options[ListDensityPlot]]]
,
ListContourPlot[If[MatrixQ[LevelMatrix,NumberQ],LevelMatrix,matrix],FilterRules[{opts},Options[ListContourPlot]],FilterRules[{ListContourPlotOptions},Options[ListContourPlot]] ]
]
, (* else ShowContours -> show only Density*)
ListDensityPlot[matrix,FilterRules[{opts},Options[ListDensityPlot]],FilterRules[{ListDensityPlotOptions},Options[ListDensityPlot]]]
] (* end ShowContours *)
,ShowLegendOptions
],
If[TrueQ[(ShowContours/.Flatten[{opts,Options[MatrixPlotWithScale]}])],
Show[ (* show Density + Contours*)
ListDensityPlot[matrix,FilterRules[{opts},Options[ListDensityPlot]],FilterRules[{ListDensityPlotOptions},Options[ListDensityPlot]]]
,
ListContourPlot[If[MatrixQ[LevelMatrix,NumberQ],LevelMatrix,matrix],FilterRules[{opts},Options[ListContourPlot]],FilterRules[{ListContourPlotOptions},Options[ListContourPlot]] ]
]
, (* else ShowContours -> show only Density*)
ListDensityPlot[matrix,FilterRules[{opts},Options[ListDensityPlot]],FilterRules[{ListDensityPlotOptions},Options[ListDensityPlot]]]
] (* end ShowContours *)
] (* End showlegend *)
] (* End MatrixPlotWithScale *)


(* cumulative probability transition matrix computed with a bias of half sate *)
EstimateCumulativeProbabilityTransitionMatrixBiased[ProbabilityTransitionMatrix_] := Module[{cumsum =0, nrows = Length[ProbabilityTransitionMatrix], ncols = Dimensions[ProbabilityTransitionMatrix][[2]]}, 
Table[ 
Switch[j,
(* first column *)
1,cumsum =ProbabilityTransitionMatrix[[i,1]], 
(* last column *)
ncols,1  ,(* Avoid roundoff error *)
_,cumsum +=ProbabilityTransitionMatrix[[i,j]]
] ,
{i,nrows},
{j,ncols}
] (* Table *) 
] (* Module *)


(* To avoid biases in the cumulative probability transition matrix, the cumulative probability is computed in the median of each bin. It is asummed that half of the probability lies on the right and on the left of the state centroid.
This is the reason for accounting half of the probability for each element.
This algorithm is recommended, at least, for ploting the contours of the ProbabilityTransitionMatrix *)
EstimateCumulativeProbabilityTransitionMatrix[ProbabilityTransitionMatrix_] := Module[{cumsum =0,  nrows = Length[ProbabilityTransitionMatrix], ncols = Dimensions[ProbabilityTransitionMatrix][[2]]}, 
Table[ 
Switch[j,
1,cumsum =0.5*ProbabilityTransitionMatrix[[i,1]],
(*,ncols,1  Avoid roundoff error, *)
_,cumsum +=0.5*(ProbabilityTransitionMatrix[[i,j-1]]+ProbabilityTransitionMatrix[[i,j]])
] ,
{i,nrows},
{j,ncols}
] (* Table *) 
] (* Module *)


End[];
EndPackage[];















