diff --git a/Libs/cbrewer2/cbrewer2.m b/Libs/cbrewer2/cbrewer2.m new file mode 100644 index 0000000..dcdc2ec --- /dev/null +++ b/Libs/cbrewer2/cbrewer2.m @@ -0,0 +1,260 @@ +%CBREWER2 Interpolated versions of Cynthia Brewer's ColorBrewer colormaps +% CBREWER2(CNAME, NCOL) returns the colour scheme CNAME with the number +% of colours equal to NCOL. If there is a ColorBrewer scheme with exactly +% this number of colours, the color scheme is returned as-is. If NCOL +% larger (or smaller) than the designed colormaps for this scheme, the +% largest (smallest) one is interpolated to provide enough colours, +% unless the requested colour scheme CNAME is a qualitative palette. For +% a qualitative scheme, the colours are repeated, cycling from the +% beginning again, to output the requested NCOL colours. +% +% CBREWER2(CNAME) without an NCOL input will use the same number of +% colours as the current colormap. +% +% CBREWER2(CNAME, NCOL, INTERP_METHOD) allows you to change the method +% used for the interpolation. The default is 'cubic'. +% +% CBREWER2(CNAME, NCOL, INTERP_METHOD, INTERP_SPACE) allows you to +% change the colorspace used for the interpolation. By default, this is +% in the CIELAB colorspace, which is approximately perceptually uniform. +% Options for INTERP_SPACE are +% 'rgb' : interpolation in sRGB (as used in original CBREWER) +% 'lab' : interpolation in CIELAB (default) +% 'lch' : interpolation in CIELCH_ab (not recommended due to the +% discontinuities at C=0 and H=0) +% Anything else supported by COLORSPACE will also function. +% +% The input format CBREWER2(TYPE, ...) can also be used, where TYPE is +% one of 'seq', 'div', 'qual'. This input is redandant and will be +% ignored. This input format is provided for backwards compatibility with +% the original CBREWER. +% +% Example 1 (sequential heatmap): +% C = [0 2 4 6; 8 10 12 14; 16 18 20 22]; +% imagesc(C); +% colorbar; +% colormap(cbrewer('YlOrRd', 256); +% +% Example 2 (line plot): +% x = 0:0.01:2; +% sc = [0.5; 1; 2]; +% t0 = [0; 0.2; 0.4]; +% t = bsxfun(@rdivide, bsxfun(@plus, x, t0), sc); +% y = sin(t * 2 * pi); +% cmap = cbrewer2('Set1', numel(sc)); +% axes('ColorOrder', cmap, 'NextPlot', 'ReplaceChildren'); +% plot(x, y); +% +% Example 3 (divergent heatmap): +% [X,Y,Z] = peaks(30); +% surfc(X,Y,Z); +% colormap(cbrewer2('RdBu')); +% +% This product includes color specifications and designs developed by +% Cynthia Brewer (http://colorbrewer.org/). For more information on +% ColorBrewer, please visit http://colorbrewer.org/. +% +% CBREWER2 uses a cached copy of the Cynthia Brewer color schemes which +% was converted to .mat format by Charles Robert for use with CBREWER. +% CBREWER is available from the MATLAB FileExchange under the MIT license. +% +% See also CBREWER, BREWERMAP, COLORSPACE, INTERP1. + + +% Copyright (c) 2016 Scott Lowe +% +% Licensed under the Apache License, Version 2.0 (the "License"); +% you may not use this file except in compliance with the License. +% You may obtain a copy of the License at +% +% http://www.apache.org/licenses/LICENSE-2.0 +% +% Unless required by applicable law or agreed to in writing, software +% distributed under the License is distributed on an "AS IS" BASIS, +% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +% See the License for the specific language governing permissions and +% limitations under the License. + + +function colormap = cbrewer2(... + cname, ncol, interp_method, interp_space, varargin) + +% Definitions ------------------------------------------------------------- + +% List of all of Cynthia Brewer's colormaps and their types +% seq: sequential +% div: divergent +% qual: qualitative +cbdict = {... + 'Blues', 'seq'; ... + 'BuGn', 'seq'; ... + 'BuPu', 'seq'; ... + 'GnBu', 'seq'; ... + 'Greens', 'seq'; ... + 'Greys', 'seq'; ... + 'Oranges', 'seq'; ... + 'OrRd', 'seq'; ... + 'PuBu', 'seq'; ... + 'PuBuGn', 'seq'; ... + 'PuRd', 'seq'; ... + 'Purples', 'seq'; ... + 'RdPu', 'seq'; ... + 'Reds', 'seq'; ... + 'YlGn', 'seq'; ... + 'YlGnBu', 'seq'; ... + 'YlOrBr', 'seq'; ... + 'YlOrRd', 'seq'; ... + 'BrBG', 'div'; ... + 'PiYG', 'div'; ... + 'PRGn', 'div'; ... + 'PuOr', 'div'; ... + 'RdBu', 'div'; ... + 'RdGy', 'div'; ... + 'RdYlBu', 'div'; ... + 'RdYlGn', 'div'; ... + 'Spectral', 'div'; ... + 'Accent', 'qual'; ... + 'Dark2', 'qual'; ... + 'Paired', 'qual'; ... + 'Pastel1', 'qual'; ... + 'Pastel2', 'qual'; ... + 'Set1', 'qual'; ... + 'Set2', 'qual'; ... + 'Set3', 'qual'; ... + }; + + +% Input handling ---------------------------------------------------------- + +narginchk(1, 5); + +% Initialise variables if not supplied +if nargin<2 + ncol = []; +end +if nargin<3 + interp_method = []; +end +if nargin<4 + interp_space = []; +end +if nargin<5 + varargin = {[]}; +end + +% Check if the colormap type was unnecessarily input +types = unique(cbdict(:, 2)); +if nargin > 1 && ischar(cname) && ischar(ncol) + LI = ismember({cname ncol}, types); + if ~any(LI); error('Number of colors cannot be a string'); end; + if all(LI); error('Incorrect colormap name'); end; + if LI(1) + vgn = {cname; ncol; interp_method; interp_space}; + cname = vgn{2}; + ncol = vgn{3}; + interp_method = vgn{4}; + interp_space = varargin{1}; + ctype_input = vgn{1}; + elseif LI(2) + vgn = {cname; ncol; interp_method; interp_space}; + cname = vgn{1}; + ncol = vgn{3}; + interp_method = vgn{4}; + interp_space = varargin{1}; + ctype_input = vgn{2}; + end +else + ctype_input = ''; +end + +% Default values +if isempty(ncol) + % Number of colours in the colormap + ncol = size(get(gcf,'colormap'), 1); +end +if isempty(interp_method) + interp_method = 'pchip'; +end +if isempty(interp_space) + interp_space = 'lab'; +end + + +% Load colorbrewer data --------------------------------------------------- +Tmp = load('colorbrewer.mat'); +colorbrewer = Tmp.colorbrewer; + +[TF, idict] = ismember(lower(cname), lower(cbdict(:, 1))); + +if ~TF + error('%s is not a recognised Brewer colormap',cname); +end + +cname = cbdict{idict, 1}; +ctype = cbdict{idict, 2}; + +if (~isfield(colorbrewer.(ctype), cname)) + error('Colormap %s is not present in loaded data',cname); +end + + +% Main script ------------------------------------------------------------- + +if ncol > length(colorbrewer.(ctype).(cname)) + % If we specified too many colours, we take the maximum and interpolate + colormap = colorbrewer.(ctype).(cname){length(colorbrewer.(ctype).(cname))}; + colormap = colormap ./ 255; +elseif isempty(colorbrewer.(ctype).(cname){ncol}) + % If we specified too few colours, we take the minimum and interpolate + nmin = find(~cellfun(@isempty, colorbrewer.(ctype).(cname)), 1); + colormap = colorbrewer.(ctype).(cname){nmin}; + colormap = colormap./255; +else + % If we specified a number of colours in the pre-designed range, no + % need to interpolate + colormap = (colorbrewer.(ctype).(cname){ncol}) ./ 255; + return; +end + +% Don't interpolate if qualitative type +if strcmp(ctype,'qual') + if size(colormap, 1) >= ncol + colormap = colormap(1:ncol, :); + return; + end + warning('CBREWER2:QualTooManyColors', ... + ['Too many colors requested: cannot interpolate a qualitative' ... + ' colorscheme']); + % Cycle the colours from the beginning again, so we have enough to + % return + colormap = repmat(colormap, ceil(ncol / size(colormap, 1)), 1); + colormap = colormap(1:ncol, :); + return; +end + +% Make sure we have colorspace downloaded from the FEX +if ~strcmpi(interp_space, 'rgb') && ~exist('colorspace.m', 'file') + P = requireFEXpackage(28790); + if isempty(P); + error(... + ['You need to download COLORSPACE from the MATLAB FEX and' ... + ' add it to the MATLAB path.']); + end; +end + +% Move to perceptually uniform space +if ~strcmpi(interp_space,'rgb') + colormap = colorspace(['rgb->' interp_space], colormap); +end + +% Linearly interpolate +X = linspace(0, 1, size(colormap, 1)); +XI = linspace(0, 1, ncol); +colormap = interp1(X, colormap, XI, interp_method); + +% Move from perceptually uniform space back to sRGB +if ~strcmpi(interp_space,'rgb') + colormap = colorspace(['rgb<-' interp_space], colormap); +end + +end \ No newline at end of file diff --git a/Libs/cbrewer2/colorbrewer.mat b/Libs/cbrewer2/colorbrewer.mat new file mode 100644 index 0000000..ec59ef4 Binary files /dev/null and b/Libs/cbrewer2/colorbrewer.mat differ diff --git a/Libs/cbrewer2/requireFEXpackage.m b/Libs/cbrewer2/requireFEXpackage.m new file mode 100644 index 0000000..86838cc --- /dev/null +++ b/Libs/cbrewer2/requireFEXpackage.m @@ -0,0 +1,194 @@ +function AddedPath = requireFEXpackage(FEXSubmissionID) +%Function requireFEXpackage - +%installs Matlab Central File Exchange (FEX) submission +%with given ID into the directory chosen by the user. +%A new FEX submissions may use previous FEX submissions as its part. +%The function 'requireFEXpackage' helps in adding those previous +%submissions to the user's MATLAB installation. +% +%This function is a part of File Exchange submission 31069. +%Download the entire submission: +%http://www.mathworks.com/matlabcentral/fileexchange/31069 +% +% SYNTAX: +% AddedPath = requireFEXpackage(FEXSubmissionID) +% +% INPUT: +% ID of the required submission to File Exchange +% +% OUTPUT: +% the path to that submission added to the user's MATLAB path. +% +% HOW TO CALL: +% The command +% P = requireFEXpackage(8277) +% will download and install the package with ID 8277 +% (namely, nice 'fminsearchbnd' by John D'Errico) +% +% EXAMPLES -- HOW TO USE: +% +% EXAMPLE 1 (using 'exist' command): +% +% % first, somewhere in the very beginning of your code, +% % check if the function 'fminsearchbnd' from the FEX package 8277 +% % is on your MATLAB path, and if it is not there, +% % require the FEX package 8277: +% if ~(exist('fminsearchbnd', 'file') == 2) +% P = requireFEXpackage(8277); % fminsearchbnd is part of 8277 +% end +% +% % Then just use 'fminsearchbnd' where you need it: +% syms x +% RosenbrockBananaFunction = @(x) (1-x(1)).^2 + 100*(x(2)-x(1).^2).^2; +% x = fminsearchbnd(RosenbrockBananaFunction,[3 3]) + +% EXAMPLE 2 (using 'try-catch' command): +% +% syms x +% RosenbrockBananaFunction = @(x) (1-x(1)).^2+100*(x(2)-x(1).^2).^2; +% try +% % if function 'fminsearchbnd' already exists in your MATLAB +% % installation, just use it: +% x = fminsearchbnd(RosenbrockBananaFunction,[3 3]) +% catch +% % if function 'fminsearchbnd' is not present in your MATLAB +% % installation, first get the package 8277 (to which it belongs) +% % from the MATLAB Central File Exchange (FEX) +% P = requireFEXpackage(8277); % fminsearchbnd is part of 8277 +% % and then use that function: +% x = fminsearchbnd(RosenbrockBananaFunction,[3 3]) +% end +% +% +% NOTE: on Mac platform, the title of the dialog box for +% choosing the directory for installing the required FEX package +% is not shown; this is not a bug, this is how UIGETDIR works on Macs -- +% see the documentation for UIGETDIR +% http://www.mathworks.com/help/techdoc/ref/uigetdir.html +% +% (C) Igor Podlubny, 2011 + +% Copyright (c) 2011, Igor Podlubny +% All rights reserved. +% +% Redistribution and use in source and binary forms, with or without +% modification, are permitted provided that the following conditions are +% met: +% +% * Redistributions of source code must retain the above copyright +% notice, this list of conditions and the following disclaimer. +% * Redistributions in binary form must reproduce the above copyright +% notice, this list of conditions and the following disclaimer in +% the documentation and/or other materials provided with the distribution +% +% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +% POSSIBILITY OF SUCH DAMAGE. + + +ID = num2str(FEXSubmissionID); + +% Ask user for the confirmation of the installation +% of the required FEX package +yes = ['YES, Install package ' ID]; +no = 'NO, do not install'; +userchoice = questdlg(['The Matlab function/toolbox, which you are running, ' ... + 'requires the presence of the package ' ID ... + ' from Matlab Central File Exchange.' ... + sprintf('\n\n') ... + 'Would you like to install the FEX package ' ID ' now?'] , ... + ['Required package ' ID], ... + yes, no, yes); + +% Handle response +switch userchoice + case yes, + install = 1; + case no, + install = 0; + otherwise, + install = 0; +end + + +if install == 1 + baseURL = 'http://www.mathworks.com/matlabcentral/fileexchange/'; + query = '?download=true'; + + location = uigetdir(pwd, ['Select the directory for installing the required FEX package' ID ]); + if location ~= 0 + % download package 'ID' from Matlab Central File Exchange + filetosave = [location filesep ID '.zip']; + FEXpackage = [baseURL ID query]; + [f,status] = urlwrite(FEXpackage,filetosave); + if status==0 + warndlg(['No connection to Matlab Central File Exchange,' ' or package ' ID ' does not exist.' ... + ' Package ' ID ' has not been installed. ' ... + ' Check you internet settings and the ID of the required package, and try again. '] , ... + ['No connection to Matlab Central File Exchange' ' or package ' ID ' does not exist'], ... + 'modal'); + AddedPath = ''; + return + end + % unzip the downloaded file to the subdirectory 'ID' + todir = [location filesep ID]; + % if the directory 'ID' doesn't exist at given location, create it + if ~(exist([location filesep ID], 'dir') == 7) + mkdir(location, ID); + end + try + unzip(filetosave, todir); + % after unzipping, delete the downloaded ZIP file + delete(filetosave); + % prepend the paths to the downloaded package to the MATLAB path + P = genpath([location filesep ID]); + path(P,path); + catch + % if the FEX package is not ZIP, then it is a single m-file + % just move the file to the ID directory + [pathstr, name, ext] = fileparts(filetosave); + movefile(filetosave, [todir filesep name '.m']); + P = genpath([location filesep ID]); + path(P,path); + end + else + P = ''; + end + + AddedPath = P; +else + AddedPath = ''; +end + + +if install == 1, + % Ask user about reviewing and saving the modified MATLAB path, + % and take him to PATHTOOL, if the user wants to save the modified path + yes = 'YES, I want to review and save the MATLAB path'; + no = 'NO, I don''t want to save the path permanently'; + userchoice = questdlg(['After adding the package ' ID ... + ' from Matlab Central File Exchange to your MATLAB installation,' ... + ' the MATLAB path has been modified accordingly. ', ... + 'Would you like to review and save the modified MATLAB path?'] , ... + 'Review and save the modified MATLAB path for future use?', ... + yes, no, yes); + + % Handle response + switch userchoice + case yes, + pathtool; + case no, + otherwise, + end +end + + + diff --git a/Libs/linspecer.m b/Libs/linspecer.m new file mode 100644 index 0000000..434a398 --- /dev/null +++ b/Libs/linspecer.m @@ -0,0 +1,261 @@ +% function lineStyles = linspecer(N) +% This function creates an Nx3 array of N [R B G] colors +% These can be used to plot lots of lines with distinguishable and nice +% looking colors. +% +% lineStyles = linspecer(N); makes N colors for you to use: lineStyles(ii,:) +% +% colormap(linspecer); set your colormap to have easily distinguishable +% colors and a pleasing aesthetic +% +% lineStyles = linspecer(N,'qualitative'); forces the colors to all be distinguishable (up to 12) +% lineStyles = linspecer(N,'sequential'); forces the colors to vary along a spectrum +% +% % Examples demonstrating the colors. +% +% LINE COLORS +% N=6; +% X = linspace(0,pi*3,1000); +% Y = bsxfun(@(x,n)sin(x+2*n*pi/N), X.', 1:N); +% C = linspecer(N); +% axes('NextPlot','replacechildren', 'ColorOrder',C); +% plot(X,Y,'linewidth',5) +% ylim([-1.1 1.1]); +% +% SIMPLER LINE COLOR EXAMPLE +% N = 6; X = linspace(0,pi*3,1000); +% C = linspecer(N) +% hold off; +% for ii=1:N +% Y = sin(X+2*ii*pi/N); +% plot(X,Y,'color',C(ii,:),'linewidth',3); +% hold on; +% end +% +% COLORMAP EXAMPLE +% A = rand(15); +% figure; imagesc(A); % default colormap +% figure; imagesc(A); colormap(linspecer); % linspecer colormap +% +% See also NDHIST, NHIST, PLOT, COLORMAP, 43700-cubehelix-colormaps +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% by Jonathan Lansey, March 2009-2013 – Lansey at gmail.com % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +%% credits and where the function came from +% The colors are largely taken from: +% http://colorbrewer2.org and Cynthia Brewer, Mark Harrower and The Pennsylvania State University +% +% +% She studied this from a phsychometric perspective and crafted the colors +% beautifully. +% +% I made choices from the many there to decide the nicest once for plotting +% lines in Matlab. I also made a small change to one of the colors I +% thought was a bit too bright. In addition some interpolation is going on +% for the sequential line styles. +% +% +%% + +function lineStyles=linspecer(N,varargin) + +if nargin==0 % return a colormap + lineStyles = linspecer(128); + return; +end + +if ischar(N) + lineStyles = linspecer(128,N); + return; +end + +if N<=0 % its empty, nothing else to do here + lineStyles=[]; + return; +end + +% interperet varagin +qualFlag = 0; +colorblindFlag = 0; + +if ~isempty(varargin)>0 % you set a parameter? + switch lower(varargin{1}) + case {'qualitative','qua'} + if N>12 % go home, you just can't get this. + warning('qualitiative is not possible for greater than 12 items, please reconsider'); + else + if N>9 + warning(['Default may be nicer for ' num2str(N) ' for clearer colors use: whitebg(''black''); ']); + end + end + qualFlag = 1; + case {'sequential','seq'} + lineStyles = colorm(N); + return; + case {'white','whitefade'} + lineStyles = whiteFade(N);return; + case 'red' + lineStyles = whiteFade(N,'red');return; + case 'blue' + lineStyles = whiteFade(N,'blue');return; + case 'green' + lineStyles = whiteFade(N,'green');return; + case {'gray','grey'} + lineStyles = whiteFade(N,'gray');return; + case {'colorblind'} + colorblindFlag = 1; + otherwise + warning(['parameter ''' varargin{1} ''' not recognized']); + end +end +% *.95 +% predefine some colormaps + set3 = colorBrew2mat({[141, 211, 199];[ 255, 237, 111];[ 190, 186, 218];[ 251, 128, 114];[ 128, 177, 211];[ 253, 180, 98];[ 179, 222, 105];[ 188, 128, 189];[ 217, 217, 217];[ 204, 235, 197];[ 252, 205, 229];[ 255, 255, 179]}'); +set1JL = brighten(colorBrew2mat({[228, 26, 28];[ 55, 126, 184]; [ 77, 175, 74];[ 255, 127, 0];[ 255, 237, 111]*.85;[ 166, 86, 40];[ 247, 129, 191];[ 153, 153, 153];[ 152, 78, 163]}')); +set1 = brighten(colorBrew2mat({[ 55, 126, 184]*.85;[228, 26, 28];[ 77, 175, 74];[ 255, 127, 0];[ 152, 78, 163]}),.8); + +% colorblindSet = {[215,25,28];[253,174,97];[171,217,233];[44,123,182]}; +colorblindSet = {[215,25,28];[253,174,97];[171,217,233]*.8;[44,123,182]*.8}; + +set3 = dim(set3,.93); + +if colorblindFlag + switch N + % sorry about this line folks. kind of legacy here because I used to + % use individual 1x3 cells instead of nx3 arrays + case 4 + lineStyles = colorBrew2mat(colorblindSet); + otherwise + colorblindFlag = false; + warning('sorry unsupported colorblind set for this number, using regular types'); + end +end +if ~colorblindFlag + switch N + case 1 + lineStyles = { [ 55, 126, 184]/255}; + case {2, 3, 4, 5 } + lineStyles = set1(1:N); + case {6 , 7, 8, 9} + lineStyles = set1JL(1:N)'; + case {10, 11, 12} + if qualFlag % force qualitative graphs + lineStyles = set3(1:N)'; + else % 10 is a good number to start with the sequential ones. + lineStyles = cmap2linspecer(colorm(N)); + end + otherwise % any old case where I need a quick job done. + lineStyles = cmap2linspecer(colorm(N)); + end +end +lineStyles = cell2mat(lineStyles); + +end + +% extra functions +function varIn = colorBrew2mat(varIn) +for ii=1:length(varIn) % just divide by 255 + varIn{ii}=varIn{ii}/255; +end +end + +function varIn = brighten(varIn,varargin) % increase the brightness + +if isempty(varargin), + frac = .9; +else + frac = varargin{1}; +end + +for ii=1:length(varIn) + varIn{ii}=varIn{ii}*frac+(1-frac); +end +end + +function varIn = dim(varIn,f) + for ii=1:length(varIn) + varIn{ii} = f*varIn{ii}; + end +end + +function vOut = cmap2linspecer(vIn) % changes the format from a double array to a cell array with the right format +vOut = cell(size(vIn,1),1); +for ii=1:size(vIn,1) + vOut{ii} = vIn(ii,:); +end +end +%% +% colorm returns a colormap which is really good for creating informative +% heatmap style figures. +% No particular color stands out and it doesn't do too badly for colorblind people either. +% It works by interpolating the data from the +% 'spectral' setting on http://colorbrewer2.org/ set to 11 colors +% It is modified a little to make the brightest yellow a little less bright. +function cmap = colorm(varargin) +n = 100; +if ~isempty(varargin) + n = varargin{1}; +end + +if n==1 + cmap = [0.2005 0.5593 0.7380]; + return; +end +if n==2 + cmap = [0.2005 0.5593 0.7380; + 0.9684 0.4799 0.2723]; + return; +end + +frac=.95; % Slight modification from colorbrewer here to make the yellows in the center just a bit darker +cmapp = [158, 1, 66; 213, 62, 79; 244, 109, 67; 253, 174, 97; 254, 224, 139; 255*frac, 255*frac, 191*frac; 230, 245, 152; 171, 221, 164; 102, 194, 165; 50, 136, 189; 94, 79, 162]; +x = linspace(1,n,size(cmapp,1)); +xi = 1:n; +cmap = zeros(n,3); +for ii=1:3 + cmap(:,ii) = pchip(x,cmapp(:,ii),xi); +end +cmap = flipud(cmap/255); +end + +function cmap = whiteFade(varargin) +n = 100; +if nargin>0 + n = varargin{1}; +end + +thisColor = 'blue'; + +if nargin>1 + thisColor = varargin{2}; +end +switch thisColor + case {'gray','grey'} + cmapp = [255,255,255;240,240,240;217,217,217;189,189,189;150,150,150;115,115,115;82,82,82;37,37,37;0,0,0]; + case 'green' + cmapp = [247,252,245;229,245,224;199,233,192;161,217,155;116,196,118;65,171,93;35,139,69;0,109,44;0,68,27]; + case 'blue' + cmapp = [247,251,255;222,235,247;198,219,239;158,202,225;107,174,214;66,146,198;33,113,181;8,81,156;8,48,107]; + case 'red' + cmapp = [255,245,240;254,224,210;252,187,161;252,146,114;251,106,74;239,59,44;203,24,29;165,15,21;103,0,13]; + otherwise + warning(['sorry your color argument ' thisColor ' was not recognized']); +end + +cmap = interpomap(n,cmapp); +end + +% Eat a approximate colormap, then interpolate the rest of it up. +function cmap = interpomap(n,cmapp) + x = linspace(1,n,size(cmapp,1)); + xi = 1:n; + cmap = zeros(n,3); + for ii=1:3 + cmap(:,ii) = pchip(x,cmapp(:,ii),xi); + end + cmap = (cmap/255); % flipud?? +end + + + diff --git a/construct_classes.m b/construct_classes.m deleted file mode 100644 index 99a8091..0000000 --- a/construct_classes.m +++ /dev/null @@ -1,2 +0,0 @@ - - diff --git a/loop_simulation.m b/loop_simulation.m deleted file mode 100644 index 99a8091..0000000 --- a/loop_simulation.m +++ /dev/null @@ -1,2 +0,0 @@ - - diff --git a/process_simulation.m b/process_simulation.m deleted file mode 100644 index 5fa3831..0000000 --- a/process_simulation.m +++ /dev/null @@ -1,37 +0,0 @@ - - - - -figure(22) -clf -subplot(2,1,1) -hold on -plot(reference.signal(10000:end-20)); -plot(rx_series(10000:end-20)); -hold off -subplot(2,1,2) -hold on -plot(reference.signal(4150:4175)); -plot(rx_series(4150:4175)); -hold off - -figure() -hold on -stem(X.signal(1,:)) -stem(bitpattern(:,1)') -hold off - -a = fftshift(xcorr(X.signal(1,:),circshift(bitpattern(:,1)',0))); - -% BER -if (length(X.signal) ~= length(bitpattern')) - warning("TX and RX bitstreams have different length...") -end - -[bits,errors,BER] = calc_ber(X.signal(:,10000:end-20),bitpattern(10000:end-19,:)',0); - -disp(X.logbook); -disp(['BER: ', sprintf('%2E',BER)]); - - - diff --git a/python_version/amp_mode.py b/python_version/amp_mode.py new file mode 100644 index 0000000..65da429 --- /dev/null +++ b/python_version/amp_mode.py @@ -0,0 +1,7 @@ +from enum import IntEnum + +class amp_mode(IntEnum): + ideal_no_noise = 1 + edfa_increase_nase = 2 + edfa_replace_nase = 3 + # edfa_set_osnr = 4 TODO: implement mode to achieve desired OSNR \ No newline at end of file diff --git a/python_version/ampclass.py b/python_version/ampclass.py new file mode 100644 index 0000000..6d1789b --- /dev/null +++ b/python_version/ampclass.py @@ -0,0 +1,86 @@ +import numpy as np +import scipy.constants as Constant +import amp_mode +import gain_mode +import nase_mode + +class Amplifier: + def __init__(self, options=None): + self.amp_mode = options.amp_mode if options and hasattr(options, 'amp_mode') else amp_mode.ideal_no_noise + self.gain_mode = options.gain_mode if options and hasattr(options, 'gain_mode') else gain_mode.output_power + self.nase_mode = options.nase_mode if options and hasattr(options, 'nase_mode') else nase_mode.pass_ase + self.amplification_db = options.amplification_db if options and hasattr(options, 'amplification_db') else 0 + self.noifig = options.noifig if options and hasattr(options, 'noifig') else 0 + self.fsimu = options.fsimu if options and hasattr(options, 'fsimu') else None + + def process(self, signalclass_in): + signalclass_in = self.process_(signalclass_in) + lbdesc = 'Amp ' + signalclass_in = signalclass_in.logbookentry(lbdesc) + signalclass_out = signalclass_in + return signalclass_out + + def process_(self, X_in): + a_lin = self.calculateGain(X_in.signal) + X_in.signal = a_lin * X_in.signal + + if isinstance(X_in, Opticalsignal): + if self.amp_mode == amp_mode.ideal_no_noise: + X_in.nase = self.onlyAmplifyAse(X_in.nase, a_lin) + elif self.amp_mode == amp_mode.edfa_increase_nase: + X_in.nase = self.increaseAse(X_in.nase, a_lin, self.noifig, X_in.lambda_) + elif self.amp_mode == amp_mode.edfa_replace_nase: + X_in.nase = self.replaceAse(X_in.nase, a_lin, self.noifig, X_in.lambda_) + + if self.nase_mode == nase_mode.generate_ase: + nase = X_in.nase + fs = X_in.fs + dimension = X_in.signal.shape + nase_numeric = self.generateAseNoise(nase, fs, dimension) + X_in.signal, osnr = self.applyAseToSignal(X_in.signal, nase_numeric) + X_in.nase = 0 + elif self.nase_mode == nase_mode.pass_ase: + X_in.nase = X_in.nase + + X_out = X_in + return X_out + + def calculateGain(self, xin): + if self.gain_mode == gain_mode.output_power: + pow_in = np.mean(np.abs(xin) ** 2) + pow_out = 10 ** (self.amplification_db / 10 - 3) + a_lin = np.sqrt(pow_out / pow_in) + elif self.gain_mode == gain_mode.gain: + a_lin = 10 ** (self.amplification_db / 20) + return a_lin + + def onlyAmplifyAse(self, nase_old, a_lin): + nase_amped = a_lin ** 2 * nase_old + return nase_amped + + def calculateAseFromThisAmp(self, a_lin, noisefig, lambda_): + h = Constant.Planck + c = Constant.speed_of_light + nase_new = 0.5 * 10 ** (noisefig / 10) * h * c / lambda_ * (a_lin ** 2 - 1) + return nase_new + + def increaseAse(self, nase_old, a_lin, noisefig, lambda_): + nase_fromThisAmp = self.calculateAseFromThisAmp(a_lin, noisefig, lambda_) + nase_old = a_lin ** 2 * nase_old + ase_increased = nase_old + nase_fromThisAmp + return ase_increased + + def replaceAse(self, nase_old, a_lin, noisefig, lambda_): + nase_fromThisAmp = self.calculateAseFromThisAmp(a_lin, noisefig, lambda_) + nase_new = nase_fromThisAmp + return nase_new + + def generateAseNoise(self, nase, fs, dimension): + nase_numeric = (np.random.randn(*dimension) + 1j * np.random.randn(*dimension)) * np.sqrt(nase / 2 * fs) + return nase_numeric + + def applyAseToSignal(self, opticalsignal, nase_numeric): + noisy_sig = opticalsignal + nase_numeric + osnr = pow(10 * np.log10(np.mean(np.abs(opticalsignal) ** 2) / np.mean(np.abs(nase_numeric) ** 2))) + print(osnr) + return noisy_sig, osnr diff --git a/python_version/gain_mode.py b/python_version/gain_mode.py new file mode 100644 index 0000000..6fa928d --- /dev/null +++ b/python_version/gain_mode.py @@ -0,0 +1,5 @@ +from enum import IntEnum + +class gain_mode(IntEnum): + gain = 1 + output_power = 2 \ No newline at end of file diff --git a/python_version/nase_mode.py b/python_version/nase_mode.py new file mode 100644 index 0000000..f88e2f7 --- /dev/null +++ b/python_version/nase_mode.py @@ -0,0 +1,5 @@ +from enum import IntEnum + +class nase_mode(IntEnum): + generate_ase = 1 + pass_ase = 2 \ No newline at end of file diff --git a/python_version/python_trial.py b/python_version/python_trial.py new file mode 100644 index 0000000..26a5866 --- /dev/null +++ b/python_version/python_trial.py @@ -0,0 +1,3 @@ +from ampclass import Amplifier + +a = Amplifier()