silas/imdd_silas
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

WIP büro

Browse Source
This commit is contained in:
Silas Oettinghaus
2023-05-24 15:49:14 +02:00
parent 81ccb80240
commit 9e0f298962
12 changed files with 821 additions and 41 deletions
Download Patch File Download Diff File Expand all files Collapse all files

260
Libs/cbrewer2/cbrewer2.m Normal file
View File

@@ -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

BIN
Libs/cbrewer2/colorbrewer.mat Normal file
View File

Binary file not shown.

194
Libs/cbrewer2/requireFEXpackage.m Normal file
View File

@@ -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

261
Libs/linspecer.m Normal file
View File

@@ -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

2
construct_classes.m
View File

@@ -1,2 +0,0 @@

2
loop_simulation.m
View File

@@ -1,2 +0,0 @@

37
process_simulation.m
View File

@@ -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)]);

7
python_version/amp_mode.py Normal file
View File

@@ -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

86
python_version/ampclass.py Normal file
View File

@@ -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

5
python_version/gain_mode.py Normal file
View File

@@ -0,0 +1,5 @@
from enum import IntEnum
class gain_mode(IntEnum):
gain = 1
output_power = 2

5
python_version/nase_mode.py Normal file
View File

@@ -0,0 +1,5 @@
from enum import IntEnum
class nase_mode(IntEnum):
generate_ase = 1
pass_ase = 2

3
python_version/python_trial.py Normal file
View File

@@ -0,0 +1,3 @@
from ampclass import Amplifier
a = Amplifier()