text:
% start : starting time of audio signal in seconds(e.g. 0 second)
% finish : e.g. ending time of audio signal in seconds (e.g. 30 second)
% filename: the input file (wavfile)
%
%
% segments: the signal segmentation with 20 msec accuracy based on RMS,
% if segments(i) = 1 then there is a change, otherwize segments(i) = 0
% The command find(segments == 1)/50 yields the signal segmentation times in sec
%
% classif: final classification of signal. classif(i) is the class
% of the i^{th} (20 msec) interval of the given signal (see also Fig 2 (b))
%
% Fig 1 shows the propability of each class per segment
%
% Fig 2 (a) shows the  segmentation of the audio signal
% Fig 2 (b) shows the  classification of the audio signal
% Fig 2 (c) shows the  RMS of the audio signal
%
% example to run: [classif, segments]=trmszc(0,60,'test3.wav')

function [classif, segments]=trmszc(start,finish,filename)

segments = [];
classif = [];

[Y1,FS,NBITS,OPTS]=wavread(filename,[1 5]);
duration = finish-start;
kanalia = size(Y1,2);
arxh = start*FS;
m = 100;
timeread = 1;

Vt = 0*[1:50];
Vta = [];
Vtz = 0*[1:50];
Vtz1 = 0*[1:50];

Vtaz = [];
Vtaz1 = [];

shma = [];
temp = [];
y = 0;
py = [1:1:50];

scp = 0;
sc = 0;
Vtp = [1:1:50];
i=1;
[PYs,FS,NBITS,OPTS]=wavread(filename,[arxh+1+(i-1)*timeread*FS  arxh+i*timeread*FS]);
sf11 = 0;
aVt = [];

%READING - Computation of RMS-ZC-RMS*ZC per 20 msec

for i=1:floor(duration/timeread),
   [Ys,FS,NBITS,OPTS]=wavread(filename,[arxh+1+(i-1)*timeread*FS  arxh+i*timeread*FS]);
   if kanalia==2,
      Ys = Ys(:,1)+Ys(:,2);
   end
   for j=1:timeread*50,
      y = Ys(round((j-1)*0.02*FS+1):floor(j*0.02*FS));
      rms = sqrt(sum(y.*y));
      Vtp = Vt;
      Vt = [Vt(2:50) rms];
      l = length(y);
      y1 = y(1:l-1);
      y2 = y(2:l);
      mm1 = y1.*y2;
      zc = 0.5*sum((abs(sign(mm1))-sign(mm1)));
      Vtz1 = [Vtz1(2:50) zc];
      zc = zc*rms;
      Vtz = [Vtz(2:50) zc];
   end

   m = mean(Vt);
   v = var(Vt);
   Vta = [Vta Vt];
   Vtaz = [Vtaz Vtz];
   Vtaz1 = [Vtaz1 Vtz1];

   if (v ~= 0 & m ~= 0)
      b(i) = v/m;
      a(i) = (m/b(i))-1;
   else
      b(i) = 1000000;
      a(i) = -1;
   end
end

%Computation of  KVRMS

l = length(Vta);
ftm(1) = mean(Vta(1:50));
ftv(1) = var(Vta(1:50));
for i=2:l-50,
  ftm(i) = ftm(i-1)+(Vta(i+49)-Vta(i-1))/50;
  ftv(i) = ftv(i-1) + ftm(i-1)*ftm(i-1) + ((Vta(i+49)*Vta(i+49)-Vta(i-1)*Vta(i-1))/50) - ftm(i)*ftm(i);

  if  ftm(i)~= 0,
     kvrms(i) =   ftv(i)/(ftm(i)*ftm(i));
  else
     kvrms(i) = 0;
  end
end

l = floor(duration/timeread);
m = mean(kvrms);
v = var(kvrms);
b1111 = v/m;
a111 = (m/b1111)-1;

%computation of similarity (omiot) between windows  based on RMS distribution
for i=2:l,
   k = (a(i)+a(i-1)+2)/2;
   omiot1(i-1) = ((2/(b(i-1)+b(i)))^k)*(b(i-1)^((a(i)+1)/2))*(b(i)^((a(i-1)+1)/2))*gamma(k)/(sqrt(gamma(a(i-1)+1)*gamma(a(i)+1)));
   if omiot1(i-1) ~= omiot1(i-1),
      omiot1(i-1) = 0;
   end
end

for i=3:l,
   k = (a(i)+a(i-2)+2)/2;
   omiot2(i-2) = ((2/(b(i-2)+b(i)))^k)*(b(i-2)^((a(i)+1)/2))*(b(i)^((a(i-2)+1)/2))*gamma(k)/(sqrt(gamma(a(i-2)+1)*gamma(a(i)+1)));
   if omiot2(i-2) ~= omiot2(i-2),
      omiot2(i-2) = 0;
   end
end

P = [];

%Segmentation
%Computation of propability  P of change  (at window i)

for i=2:l-1,
   P(i-1) =(1-omiot2(i-1))*(1-omiot1(i-1)+1-omiot1(i))*(1-omiot2(i-1));
end

for i=2:l-1,
   P(i-1) =(1-omiot2(i-1));
end

P = [0 P 0];

M = mean(P);
V1 = mean(abs(P(1:l-3)-P(2:l-2)));
V2 = mean(abs(P(1:l-4)-P(3:l-2)));

V = (V1+V2)*0.25;

V = 0.25*median(abs((P(1:l-4)-P(3:l-2))));

Pd(1) = P(1);
Pd(2) = P(2);

for i=3:l-2,
   m = mean(P([i-2:i-1 i+1:i+2]));
   m1 = max(P([i-2:i-1 i+1:i+2]));

   if (m < 0.1*M)
      m1 = 0.1*M;
   end
   d = (P(i)-m)/m1;

   if d < 0,
      d = 0;
   end
   Pd(i) = P(i)*d/V;
end
Pd(l) = P(l);
Pd(l-1) = P(l-1);

%window size = 1 sec

%Selection of candicate windows where there is a change

j = 0;
for i=3:l-3,
   if (Pd(i) > 5 & P(i) > P(i+1) & P(i) > P(i-1)),
      j = j+1;
      pos(j) = i;
   end
end

if j == 0,
   return;
end

%Computation of change with 20 msec accurancy

k2 = 1;
c = 25;
for k1=1:length(pos),
   i = 50*(pos(k1)-2);  % msec
   for j=-c:50+c,
      if (ftv(i+j) ~= 0 & ftm(i+j) ~= 0),
         b1 = ftv(i+j)/ftm(i+j);
         a1 = ftm(i+j)/b1-1;
      else
         b1 = 10000;
         a1 = -1;
      end

      if (ftv(i+j+50) ~= 0 & ftm(i+j+50) ~= 0)
         b2 = ftv(i+j+50)/ftm(i+j+50);
         a2 = ftm(i+j+50)/b2-1;
      else
         b2 = 10000;
         a2 = -1;
      end
      k = (a2+a1+2)/2;
      Om(j+c+1) = ((2/(b1+b2))^k)*(b1^((a2+1)/2))*(b2^((a1+1)/2))*gamma(k)/(sqrt(gamma(a1+1)*gamma(a2+1)));
      if Om(j+c+1)~= Om(j+c+1),
         Om(j+c+1) = 0;
      end
   end
   [x posms(k1)] = min(Om);
   h = 1-Om;
   if mean(h(max(posms(k1)-25,1):min(2*c+51,posms(k1)+25))) >  0.1,
      msec(k2) = posms(k1)-1-c;
      mpos(k2) = pos(k1);
      k2 = k2+1;
   end
end

posms = [];
posms = msec;
pos = [];
pos = mpos;

fpos = start+pos-1+0.02*posms;

pfpos = floor(50*pos+posms-50);

bvoise = 0.14339738781836;
bmusic = 0.04399754659151;
amusic = 1.66349817725076;
avoise = 2.32677887950291;

pfpos = [1 pfpos length(kvrms)];
l = length(pfpos)-1;
V = [];
Pzero = [];
Fsp1 = [];

%Classification for each segment

for i=1:l,

   d = 0;

   x1 = mean(kvrms([pfpos(i):pfpos(i+1)]));
   x = x1;
   y = Vtaz([pfpos(i)+d:pfpos(i+1)-d]);
   y = y/(2*max(Vta(pfpos(i)+d:pfpos(i+1)-d)) - min(Vta(pfpos(i)+d:pfpos(i+1)-d))-median(Vta(pfpos(i)+d:pfpos(i+1)-d)));

   thor(i) = 50*mean(exp(-y));

   Pmusic(i) = (x^amusic)*exp(-x/bmusic)/(gamma(amusic+1)*bmusic^(amusic+1));
   Pvoise(i) = 0.5*(x^avoise)*exp(-x/bvoise)/(gamma(avoise+1)*bvoise^(avoise+1));
   sm = 0.7*median(ftm(pfpos(i):pfpos(i+1)))+0.3*mean(ftm(pfpos(i):pfpos(i+1)));
   Pspace(i) =  6*exp((-sm*sm)/(2*0.6*0.6));

   zpos(i) = sum(exp(-10*Vtaz1([pfpos(i)+d:pfpos(i+1)-d])))/length(Vtaz1([pfpos(i)+d:pfpos(i+1)-d]));

   if zpos(i) > 0.08 | x > 3
     Pvoise(i)  = 10;
  else
      Pmusic(i) = (x^amusic)*exp(-x/bmusic)/(gamma(amusic+1)*bmusic^(amusic+1));
      Pvoise(i)  = 0.5*(x^avoise)*exp(-x/bvoise)/(gamma(avoise+1)*bvoise^(avoise+1));
   end

   V(i,:) = [Pmusic(i)  Pvoise(i)];
   [Pmax(i) type(i)] = max(V(i,:));

  thor2(i) = 50*mean(exp(-4*y));
  Pmusicg(i) = 0;

   if thor2(i) > 3.5,
      Pvoiseg(i) = 10;
   else
      Pmusicg(i) = (x^amusic)*exp(-x/bmusic)/(gamma(amusic+1)*bmusic^(amusic+1));
      Pvoiseg(i)  = 0.5*(x^avoise)*exp(-x/bvoise)/(gamma(avoise+1)*bvoise^(avoise+1));
   end

   Vg(i,:) = [Pmusicg(i)  Pvoiseg(i)];
   [Pmaxg(i) typeg(i)] = max(Vg(i,:));

   Vk = (sign(Vtaz1([pfpos(i):pfpos(i+1)])));
   tempV = [];
   tempV = Vta([pfpos(i):pfpos(i+1)]);
   rr = max(tempV)+0.001;
   tempV1 = tempV/rr;
   Vk0 = Vk;

   for j = 1:length(Vk)
      if (tempV1(j) < 0.1 & tempV(j) < 1.5) | tempV(j) < 1
         Vk(j) = 0;
      end
   end
   Vk1 = Vk;

   for j = 1:length(Vk)
      if (tempV1(j) < 0.4 | tempV(j) < 2 | tempV(j) < mean(tempV(j)))
         Vk1(j) = 0;
      end
   end
   Zca = [];
   Zca = Vk1.*(Vtaz1([pfpos(i):pfpos(i+1)]));

   Pol = (ones(1,length(Vk)) -  sign(Vk)).*Vk0;

   VZC = Pol.*Vtaz1([pfpos(i):pfpos(i+1)]);

   dVk = abs(Vk(2:length(Vk))-Vk(1:length(Vk)-1));

   Freq(i,1) = 50*sum(dVk)/(2*length(Vk)); %FSP%
   Freq(i,2) = length(Vk)/50;%duration in sec of the segment
   Freq(i,3) = sum(Freq(:,2));
   Freq(i,4) = zpos(i);
   Freq(i,5) = thor2(i);
   Freq(i,6) =   max(Zca);

   Pmusicg(i) = 0;
   Pvoiseg(i) = 0;
   if Freq(i,1) < 0.59 & Freq(i,2) > 2.5,
      Pmusicg(i) = 10;
   else
      Pmusicg(i) = (x^amusic)*exp(-x/bmusic)/(gamma(amusic+1)*bmusic^(amusic+1));
      Pvoiseg(i)  = 0.5*(x^avoise)*exp(-x/bvoise)/(gamma(avoise+1)*bvoise^(avoise+1));
      if x > 3,
         Pvoiseg(i) = Pvoiseg(i)+0.1;

      end

      if x > 300 ,
         Pmusicg(i) = 9;
     end

      if thor2 > 3.5
         Pvoiseg(i) = 11;
      end
      if max(Zca) > 280,
         Pmusicg(i) = 11;
         max(Zca)

      end
      if Freq(i,1) > 4.62 & Freq(i,2) > 2.5 & (thor2(i) > 0.1 | zpos(i) > 0.05),
         Pvoiseg(i) = 12;
      end

      if zpos(i) > 0.15
         Pvoiseg(i)  = 13;
      end
  end
   Vg(i,:) = [Pmusicg(i)  Pvoiseg(i)];
   [Pmaxg(i) typeg(i)] = max(Vg(i,:));

	Fsp1 = [Fsp1 Vk];

   if (Pspace(i) > 4 & thor2(i) < 1)  | (Pspace(i) > 4.5),
      type(i) = 3; %silence
      typeg(i) = 3;
   end

   if (pfpos(i+1) - pfpos(i))/50 < 0.5 & i >= 1, % short segments
      type(i) = type(i-1);
   end
end

   g = [];
   g1 = [];
   tt = [];
   pt = [];
   pm = [];
   ps = [];
   pv = [];
   pmf = [];

for i=1:l,
   tt = [tt type(i)*ones(1,pfpos(i+1)-pfpos(i))];
   pt = [pt Pmax(i)*ones(1,pfpos(i+1)-pfpos(i))];
   pm = [pm Pmusic(i)*ones(1,pfpos(i+1)-pfpos(i))];
   pv = [pv Pvoise(i)*ones(1,pfpos(i+1)-pfpos(i))];
   ps = [ps Pspace(i)*ones(1,pfpos(i+1)-pfpos(i))];
end

PosV = 100*sum(exp(-10*abs(tt-2*ones(1,length(tt)))))/length(tt);
PosM = 100*sum(exp(-10*abs(tt-ones(1,length(tt)))))/length(tt);

l1 = l;
l = length(tt);

l = l1;
 g = [];
   g1 = [];
   tt = [];
   pt = [];
   pm = [];
   ps = [];
   pv = [];
   pmf = [];

for i=1:l,
   tt = [tt typeg(i)*ones(1,pfpos(i+1)-pfpos(i))];
   pt = [pt Pmaxg(i)*ones(1,pfpos(i+1)-pfpos(i))];
   pm = [pm Pmusicg(i)*ones(1,pfpos(i+1)-pfpos(i))];
   pv = [pv Pvoiseg(i)*ones(1,pfpos(i+1)-pfpos(i))];
   ps = [ps Pspace(i)*ones(1,pfpos(i+1)-pfpos(i))];
end
l = length(tt);

PosVg = 100*sum(exp(-10*abs(tt-2*ones(1,length(tt)))))/length(tt);
PosMg = 100*sum(exp(-10*abs(tt-ones(1,length(tt)))))/length(tt);

figure(1);
subplot(5,1,1);
plot([1:l]*(finish-start-1)/l+start,tt);
title('classification type (1 : Music, 2: Voise, 3: Silence)');
subplot(5,1,2);
plot([1:l]*(finish-start-1)/l+start,pt);
title('~probability posibility of selected type');
subplot(5,1,3);
plot([1:l]*(finish-start-1)/l+start,pm);
title('~probability MUsic');
subplot(5,1,4);
plot([1:l]*(finish-start-1)/l+start,pv);
title('~probability Voise');
subplot(5,1,5);
plot([1:l]*(finish-start-1)/l+start,ps);
title('~probability silence');

for k = 1:length(Vta),
   if Vta(k) < 2 | Vta(k) < 0.1*max(Vta),
      Vtaz1(k) = 0;
   end
end
m = 100;
temp = [];
h1 = [];
temp = Vtaz1;
h1 = hist(temp,m);
h1(1) = 0;
h1 = h1/sum(h1);
a1 = [min(temp):(max(temp)-min(temp))/m:max(temp)-(max(temp)-min(temp))/m];

segbar = [];
j = 1;
i = 1;
while i <= l,
   if j <= length(fpos) & round(50*fpos(j)) == i,
      segbar(i) = 1;
      j = j+1;
   else
      segbar(i) = 0;
   end
   i = i+1;
end

segments = segbar;

figure(2);
subplot(3,1,1);
stem([1:l]*(finish-start-1)/l+start,segbar);
axis([0 duration 0 1]);
title('Segmentation');
xlabel('sec');
subplot(3,1,2);
plot([1:l]*(finish-start-1)/l+start,tt);
axis([0.5 duration 0 3.5]);
title('classification type (1 : Music, 2: Voise, 3: Silence)');
xlabel('sec');
subplot(3,1,3);
plot(start+([1:length(Vta)]*duration/length(Vta)),Vta);
axis([0 duration 0 (floor(max(Vta))+1)]);
title('RMS');
xlabel('sec');

classif = tt;

author: Sable

repo: mcbench-benchmarks


text:
% ========================================================================
% USAGE: [Coeff]=LLC_coding_appr(B,X,knn,lambda)
% Approximated Locality-constraint Linear Coding
%
% Inputs
%       B       -M x d codebook, M entries in a d-dim space
%       X       -N x d matrix, N data points in a d-dim space
%       knn     -number of nearest neighboring
%       lambda  -regulerization to improve condition
%
% Outputs
%       Coeff   -N x M matrix, each row is a code for corresponding X
%
% Jinjun Wang, march 19, 2010
% URL: http://www.ifp.illinois.edu/~jyang29/codes/CVPR10-LLC.rar
%
% Minor edits by Aditya Khosla
% ========================================================================

function [Coeff] = LLC_coding_appr(B, X, knn, beta)

if ~exist('knn', 'var') || isempty(knn),
    knn = 5;
end

if ~exist('beta', 'var') || isempty(beta),
    beta = 3e-2;
end

nframe=size(X,1);
nbase=size(B,1);

% find k nearest neighbors
D = sp_dist2(X, B);
[~, sort_idx] = sort(D, 2, 'ascend');
IDX = sort_idx(:, 1:knn);

% llc approximation coding
II = eye(knn, knn);
Coeff = zeros(nframe, nbase);
for i=1:nframe
   idx = IDX(i,:);
   z = bsxfun(@minus, B(idx,:), X(i,:));           % shift ith pt to origin
   C = z*z';                                        % local covariance
   C = C + II*beta*trace(C);                        % regularlization (K>D)
   w = C\ones(knn,1);
   w = w/sum(w);                                    % enforce sum(w)=1
   Coeff(i,idx) = w';
end

author: adikhosla

repo: feature-extraction


text:
function y=add_CP(x,Ncp)
    y = [x(:,end-Ncp+1:end) x];
end

author: LyricYang


text:
function out = DN_cv(x,k)
% DN_cv     Coefficient of variation
%
% Coefficient of variation of order k is sigma^k / mu^k (for sigma, standard
% deviation and mu, mean) of a data vector, x
%
%---INPUTS:
%
% x, the input data vector
% k, the order of coefficient of variation (k = 1 is default)

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This function is free software: you can redistribute it and/or modify it under
% the terms of the GNU General Public License as published by the Free Software
% Foundation, either version 3 of the License, or (at your option) any later
% version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
% details.
%
% You should have received a copy of the GNU General Public License along with
% this program. If not, see <http://www.gnu.org/licenses/>.
% ------------------------------------------------------------------------------

%-------------------------------------------------------------------------------
% Check inputs
%-------------------------------------------------------------------------------
if nargin < 2 || isempty(k)
    k = 1; % Do standard CV by default
end

if (rem(k,1) ~= 0) || (k < 0)
    warning('k should probably be a positive integer');
    % Carry on with just this warning, though
end

% Compute the coefficient of variation (of order k) of the data

out = (std(x))^k / (mean(x))^k;

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [varargout] = ft_plot_box(position, varargin)

% FT_PLOT_BOX plots the outline of a box that is specified by its lower
% left and upper right corner
%
% Use as
%   ft_plot_box(position, ...)
% where the position of the box is specified as is [x1, x2, y1, y2].
%
% Optional arguments should come in key-value pairs and can include
%   'facealpha'       = transparency value between 0 and 1
%   'facecolor'       = color specification as [r g b] values or a string, for example 'brain', 'cortex', 'skin', 'red', 'r'
%   'edgecolor'       = color specification as [r g b] values or a string, for example 'brain', 'cortex', 'skin', 'red', 'r'
%   'tag'             = string, the name assigned to the object. All tags with the same name can be deleted in a figure, without deleting other parts of the figure.
%
% It is possible to plot the object in a local pseudo-axis (c.f. subplot), which is specfied as follows
%   'hpos'            = horizontal position of the center of the local axes
%   'vpos'            = vertical position of the center of the local axes
%   'width'           = width of the local axes
%   'height'          = height of the local axes
%   'hlim'            = horizontal scaling limits within the local axes
%   'vlim'            = vertical scaling limits within the local axes
%   'parent'          = handle which is set as the parent for all plots
%
% Example
%   ft_plot_box([-1 1 2 3], 'facecolor', 'b')
%   axis([-4 4 -4 4])
%
% See also FT_PLOT_LINE, FT_PLOT_CROSSHAIR

% Copyrights (C) 2009-2011, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

ws = warning('on', 'MATLAB:divideByZero');

% get the optional input arguments
hpos        = ft_getopt(varargin, 'hpos');
vpos        = ft_getopt(varargin, 'vpos');
width       = ft_getopt(varargin, 'width');
height      = ft_getopt(varargin, 'height');
hlim        = ft_getopt(varargin, 'hlim');
vlim        = ft_getopt(varargin, 'vlim');
facealpha   = ft_getopt(varargin, 'facealpha', 1);
facecolor   = ft_getopt(varargin, 'facecolor', 'none');
edgecolor   = ft_getopt(varargin, 'edgecolor', 'k');
tag         = ft_getopt(varargin, 'tag',       '');
parent      = ft_getopt(varargin, 'parent', []);

% color management
if ischar(facecolor) && exist([facecolor '.m'], 'file')
	facecolor = eval(facecolor);
end
if ischar(edgecolor) && exist([edgecolor '.m'], 'file')
	edgecolor = eval(edgecolor);
end

% convert the two cornerpoints into something that the patch function understands
% the box position is represented just like the argument to the AXIS function
x1 = position(1);
x2 = position(2);
y1 = position(3);
y2 = position(4);
X = [x1 x2 x2 x1 x1];
Y = [y1 y1 y2 y2 y1];

if isempty(hlim) && isempty(vlim) && isempty(hpos) && isempty(vpos) && isempty(height) && isempty(width)
  % no scaling is needed, the input X and Y are already fine
  % use a shortcut to speed up the plotting

else
  % use the full implementation
  if isempty(hlim)
    hlim = get(gca, 'XLim');
  end

  if isempty(vlim)
    vlim = get(gca, 'YLim');
  end

  if isempty(hpos)
    hpos = (hlim(1)+hlim(2))/2;
  end

  if isempty(vpos)
    vpos = (vlim(1)+vlim(2))/2;
  end

  if isempty(width)
    width = hlim(2)-hlim(1);
  end

  if isempty(height)
    height = vlim(2)-vlim(1);
  end

  % first shift the horizontal axis to zero
  X = X - (hlim(1)+hlim(2))/2;
  % then scale to length 1
  X = X ./ (hlim(2)-hlim(1));
  % then scale to the new width
  X = X .* width;
  % then shift to the new horizontal position
  X = X + hpos;

  % first shift the vertical axis to zero
  Y = Y - (vlim(1)+vlim(2))/2;
  % then scale to length 1
  Y = Y ./ (vlim(2)-vlim(1));
  % then scale to the new width
  Y = Y .* height;
  % then shift to the new vertical position
  Y = Y + vpos;

end % shortcut

% use an arbitrary color, which will be replaced by the correct color a few lines down
C = 0;

h = patch(X, Y, C);
set(h, 'FaceAlpha', facealpha)
set(h, 'FaceColor', facecolor)
set(h, 'EdgeColor', edgecolor)
set(h, 'tag', tag);

if ~isempty(parent)
  set(h, 'Parent', parent);
end

% the (optional) output is the handle
if nargout == 1
  varargout{1} = h;
end

warning(ws); % revert to original state

author: spm


text:
function q_out=qmult(q1,q2)
% QMULT(Q1,Q2) calculates the product of two quaternions Q1 and Q2.
%    Inputs can be vectors of quaternions, but they must either have the
%    same number of component quaternions, or one input must be a single
%    quaternion.  QMULT will determine whether the component quaternions of
%    the inputs are row or column vectors according to ISQ.
%
%    The output will have the same shape as Q1.  If the component
%    quaternions of either Q1 or Q2 (but not both) are of indeterminate
%    shape (see ISQ), then the shapes will be assumed to be the same for
%    both inputs.  If both Q1 and Q2 are of indeterminate shape, then both
%    are assumed to be composed of row vector quaternions.
%
% See also ISQ.

% Release: $Name: quaternions-1_2_2 $
% $Revision: 1.12 $
% $Date: 2001/05/01 20:20:31 $

% Copyright (C) 2001, Jay A. St. Pierre.  All rights reserved.

if nargin~=2
  error('qmult() requires two input arguments');
else
  q1type = isq(q1);
  if ( q1type == 0 )
    error(['Invalid input: q1 must be a quaternion or a vector of' ...
          ' quaternions'])
  end
  q2type = isq(q2);
  if ( q2type == 0 )
    error(['Invalid input: q2 must be a quaternion or a vector of' ...
          ' quaternions'])
  end
end

% Make sure q1 is a column of quaternions (components are rows)
if ( q1type==1 | (q1type==3 & q2type==1) )
  q1=q1.';
end

% Make sure q2 is a column of quaternions (components are rows)
if ( q2type==1 | (q2type==3 & q1type==1) )
  q2=q2.';
end

num_q1=size(q1,1);
num_q2=size(q2,1);

if (  num_q1~=num_q2 & num_q1~=1 & num_q2~=1 )
  error(['Inputs do not have the same number of elements:', 10, ...
         '   number of quaternions in q1 = ', num2str(num_q1), 10,...
         '   number of quaternions in q2 = ', num2str(num_q2), 10,...
         'Inputs must have the same number of elements, or', 10, ...
         'one of the inputs must be a single quaternion (not a', 10, ...
         'vector of quaternions).'])
end

% Build up full quaternion vector if one input is a single quaternion
if ( num_q1 ~= num_q2 )
  ones_length = ones(max(num_q1,num_q2),1);
  if ( num_q1 == 1 )
    q1 = [q1(1)*ones_length ...
          q1(2)*ones_length ...
          q1(3)*ones_length ...
          q1(4)*ones_length ];
  else % num_q2 == 1
    q2 = [q2(1)*ones_length ...
          q2(2)*ones_length ...
          q2(3)*ones_length ...
          q2(4)*ones_length ];
  end
end

% Products

% If q1 and q2 are not vectors of quaternions, then:
%
%   q1*q2 = q1*[ q2(4) -q2(3)  q2(2) -q2(1)
%                q2(3)  q2(4) -q2(1) -q2(2)
%               -q2(2)  q2(1)  q2(4) -q2(3)
%                q2(1)  q2(2)  q2(3)  q2(4) ]
%
% But to deal with vectorized quaternions, we have to use the ugly
% commands below.

prod1 = ...
    [ q1(:,1).*q2(:,4) -q1(:,1).*q2(:,3)  q1(:,1).*q2(:,2) -q1(:,1).*q2(:,1)];
prod2 = ...
    [ q1(:,2).*q2(:,3)  q1(:,2).*q2(:,4) -q1(:,2).*q2(:,1) -q1(:,2).*q2(:,2)];
prod3 = ...
    [-q1(:,3).*q2(:,2)  q1(:,3).*q2(:,1)  q1(:,3).*q2(:,4) -q1(:,3).*q2(:,3)];
prod4 = ...
    [ q1(:,4).*q2(:,1)  q1(:,4).*q2(:,2)  q1(:,4).*q2(:,3)  q1(:,4).*q2(:,4)];

q_out = prod1 + prod2 + prod3 + prod4;

% Make sure output is same format as q1
if ( q1type==1 | (q1type==3 & q2type==1) )
  q_out=q_out.';
end

% NOTE that the following algorithm proved to be slower than the one used
% above:
%
% q_out = zeros(size(q1));
%
% q_out(:,1:3) = ...
%     [q1(:,4) q1(:,4) q1(:,4)].*q2(:,1:3) + ...
%     [q2(:,4) q2(:,4) q2(:,4)].*q1(:,1:3) + ...
%     cross(q1(:,1:3), q2(:,1:3));
%
% q_out(:,4) = q1(:,4).*q2(:,4) - dot(q1(:,1:3), q2(:,1:3), 2);

author: christianwengert


text:
function [params, names] = rbfKernExtractParam(kern)

% RBFKERNEXTRACTPARAM Extract parameters from the RBF kernel structure.
% FORMAT
% DESC Extract parameters from the radial basis function kernel
% structure into a vector of parameters for optimisation.
% ARG kern : the kernel structure containing the parameters to be
% extracted.
% RETURN param : vector of parameters extracted from the kernel. If
% the field 'transforms' is not empty in the kernel matrix, the
% parameters will be transformed before optimisation (for example
% positive only parameters could be logged before being returned).
%
% FORMAT
% DESC Extract parameters and parameter names from the radial basis
% function kernel structure.
% ARG kern : the kernel structure containing the parameters to be
% extracted.
% RETURN param : vector of parameters extracted from the kernel. If
% the field 'transforms' is not empty in the kernel matrix, the
% parameters will be transformed before optimisation (for example
% positive only parameters could be logged before being returned).
% RETURN names : cell array of strings giving names to the parameters.
%
% SEEALSO rbfKernParamInit, rbfKernExpandParam, kernExtractParam, scg, conjgrad
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005, 2006

% KERN

params = [kern.inverseWidth kern.variance];
if nargout > 1
  names={'inverse width', 'variance'};
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [tp_score,fp_score,totalPos,tp_points,fp_points,tp_idxs,fp_idxs,fp_dist,missing_recall] = computeTFP(dist,score,points,rectidxs_gt,rect_ignore,thresh_gt)

if (nargin < 5)
    rect_ignore = cell(length(dist),1);
end

if (nargin < 6)
    thresh_gt = 0.5;
end

assert(length(rectidxs_gt) == length(dist));

N = sum(cellfun('length',score));
fp_score = nan(N,1);
tp_score = nan(N,1);

fp_points = nan(N,2);
tp_points = nan(N,2);

fp_idxs = nan(N,1);
tp_idxs = nan(N,1);

fp_dist = nan(N,1);
missing_recall = cell(length(dist),1);

last_fpidx = 0;
last_tpidx = 0;
totalPos = 0;
nskipped = 0;

for imgidx = 1:length(dist)

%     if (imgidx == 212)
%         fprintf('blabla\n');
%     end
    assert(size(dist{imgidx},1) == size(score{imgidx},1));

    [minDist,minIdx] = min(dist{imgidx},[],2);
    match_gt = minDist <= thresh_gt;

    tp_img = nan(1,size(dist{imgidx},2));
    tp_score_img = -inf(1,size(dist{imgidx},2));

    fp_img = find(match_gt == 0);

    % exclude fp on not annotated people
    idxs_ignore = false(length(fp_img),1);
    if (~isempty(rect_ignore{imgidx}))
        for i=1:size(rect_ignore{imgidx},1)
            r = rect_ignore{imgidx}(i,:);
            condX1 = points{imgidx}(fp_img,1) >= repmat(r(1),length(fp_img),1);
            condX2 = points{imgidx}(fp_img,1) <= repmat(r(3),length(fp_img),1);
            condY1 = points{imgidx}(fp_img,2) >= repmat(r(2),length(fp_img),1);
            condY2 = points{imgidx}(fp_img,2) <= repmat(r(4),length(fp_img),1);
            s = condX1 + condX2 + condY1 + condY2;
            idxs_ignore(s == 4) = true;
        end
    end
    fp_img(idxs_ignore) = [];

    fp_score_img = score{imgidx}(fp_img);

    matchidx = find(match_gt == 1);

    nskipped_img = 0;
    for ridx = matchidx'

        % ignore tp w.r.t other rectangles
        if (ismember(minIdx(ridx),rectidxs_gt{imgidx}))

            if (isinf(tp_score_img(minIdx(ridx))))
                tp_img(minIdx(ridx)) = ridx;
                tp_score_img(minIdx(ridx)) = score{imgidx}(ridx);
            elseif (tp_score_img(minIdx(ridx)) < score{imgidx}(ridx))
                fp_img = [fp_img; tp_img(minIdx(ridx))];
                fp_score_img = [fp_score_img; tp_score_img(minIdx(ridx))];
                tp_img(minIdx(ridx)) = ridx;
                tp_score_img(minIdx(ridx)) = score{imgidx}(ridx);
            else
                fp_img = [fp_img; ridx];
                fp_score_img = [fp_score_img; score{imgidx}(ridx)];
            end
        else
            nskipped_img = nskipped_img + 1;
%             fprintf('not a member!\n');
        end
    end % rectangles

    % rectangles with the part which is not present in the image
    idxs_absent = find(sum(isnan(dist{imgidx}),1) == size(dist{imgidx},1));

    % relevant rectangles with present body part
    idxs_present = setdiff(rectidxs_gt{imgidx},idxs_absent);
%     idxs_present = rectidxs_gt{imgidx};
    totalPos = totalPos + length(idxs_present);

    idxs = find(isinf(tp_score_img));
    tp_score_img(idxs) = [];
    tp_img(idxs) = [];

    % relevant not found rectangles
    missing_recall{imgidx} = intersect(idxs,idxs_present);

    assert(length(tp_img) + length(fp_img) + nskipped_img + sum(idxs_ignore) == length(match_gt));

    idxsNaN = find(isnan(fp_score_img));
    if (~isempty(idxsNaN))
        fp_img(idxsNaN) = [];
        fp_score_img(idxsNaN) = [];
    end
    tp_score(last_tpidx+1:last_tpidx+length(tp_img)) = tp_score_img;
    fp_score(last_fpidx+1:last_fpidx+length(fp_img)) = fp_score_img;
    tp_points(last_tpidx+1:last_tpidx+length(tp_img),:) = points{imgidx}(tp_img,:);
    fp_points(last_fpidx+1:last_fpidx+length(fp_img),:) = points{imgidx}(fp_img,:);
    tp_idxs(last_tpidx+1:last_tpidx+length(tp_img)) = imgidx;
    fp_idxs(last_fpidx+1:last_fpidx+length(fp_img)) = imgidx;
    fp_dist(last_fpidx+1:last_fpidx+length(fp_img)) = minDist(fp_img);

    last_tpidx = last_tpidx+length(tp_img);
    last_fpidx = last_fpidx+length(fp_img);

    nskipped = nskipped + nskipped_img;
end

idxs = find(isnan(tp_score));
tp_score(idxs) = [];
tp_points(idxs,:) = [];
tp_idxs(idxs) = [];

idxs = find(isnan(fp_score));
fp_score(idxs) = [];
fp_points(idxs,:) = [];
fp_idxs(idxs) = [];
fp_dist(idxs) = [];

% assert(length(tp_score) + length(fp_score) + nskipped == N);

end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
function jac = p00_jac ( problem, option, nvar, x )

%*****************************************************************************80
%
%% P00_JAC evaluates the jacobian for any problem.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 October 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer PROBLEM, the problem index.
%
%    Input, integer OPTION, the option index.
%
%    Input, integer NVAR, the number of variables.
%
%    Input, real X(NVAR), the argument of the jacobian.
%
%    Output, real JAC(NVAR-1,NVAR), the jacobian matrix evaluated
%    at X.  The NVAR-th row is not set by this routine.
%
  if ( problem == 1 )
    jac = p01_jac ( option, nvar, x );
  elseif ( problem == 2 )
    jac = p02_jac ( option, nvar, x );
  elseif ( problem == 3 )
    jac = p03_jac ( option, nvar, x );
  elseif ( problem == 4 )
    jac = p04_jac ( option, nvar, x );
  elseif ( problem == 5 )
    jac = p05_jac ( option, nvar, x );
  elseif ( problem == 6 )
    jac = p06_jac ( option, nvar, x );
  elseif ( problem == 7 )
    jac = p07_jac ( option, nvar, x );
  elseif ( problem == 8 )
    jac = p08_jac ( option, nvar, x );
  elseif ( problem == 9 )
    jac = p09_jac ( option, nvar, x );
  elseif ( problem == 10 )
    jac = p10_jac ( option, nvar, x );
  elseif ( problem == 11 )
    jac = p11_jac ( option, nvar, x );
  elseif ( problem == 12 )
    jac = p12_jac ( option, nvar, x );
  elseif ( problem == 13 )
    jac = p13_jac ( option, nvar, x );
  elseif ( problem == 14 )
    jac = p14_jac ( option, nvar, x );
  elseif ( problem == 15 )
    jac = p15_jac ( option, nvar, x );
  elseif ( problem == 16 )
    jac = p16_jac ( option, nvar, x );
  elseif ( problem == 17 )
    jac = p17_jac ( option, nvar, x );
  elseif ( problem == 18 )
    jac = p18_jac ( option, nvar, x );
  elseif ( problem == 19 )
    jac = p19_jac ( option, nvar, x );
  elseif ( problem == 20 )
    jac = p20_jac ( option, nvar, x );
  else
    jac = [];
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P00_JAC - Fatal error!\n' );
    fprintf ( 1, '  Unrecognized problem number = %d\n', problem );
    error ( 'P00_JAC - Fatal error!' );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function h = rxVisualizeRx(rx, loc);
%
% h = rxVisualizeRx(rx, [loc]);
%
% Display a visualization of the current prescription in a separate
% figure. Shows a montage of interpolated slices, and, if the
% Rx Figure is also open, adds a panel with sagittal, coronal, and
% axial views with the prescription on it.
%
% 'loc' is an optional argument for visualizing the rx, specifying
% the [axi cor sag] slices to show in the Rx views.
%
% ras, 01/2007.
if notDefined('rx'), rx = get(findobj('Tag','rxControlFig'), 'UserData'); end
if ishandle(rx),	rx = get(rx, 'UserData');				end

h = figure('Color', 'w');

nSlices = rx.rxDims(3);
for slice = 1:nSlices
	[img2D images{slice}] = rxInterpSlice(rx, slice);
end

montage = imageMontage(images);
imagesc(montage); axis image; axis off;
title('Prescribed Slices');

%% draw ROIs if selected
if ~isempty(rx.rois) & get(rx.ui.rxDrawRx, 'Value')==1
	for r = 1:length(rx.rois)
		% convert coords from vol -> rx
		rxCoords = vol2rx(rx, rx.rois(r).volCoords, 1);

		% now convert these coords to 2-D image coordinates
		pts = rxCoords2Montage(rxCoords, nSlices, size(img2D));

        if isempty(pts), continue; end

        % remove pts outside range
        outside = find(pts(1,:) < 1 | pts(1,:) > size(montage, 1) ...
                       | pts(2,:) < 1 | pts(2,:) > size(montage, 2));
        ok = setdiff(1:size(pts, 2), outside);
        pts = pts(:,ok);

		% draw the outline of this ROI
		% (rx.rois(r) doubles as a prefs struct, to set the color)
		outline(pts, rx.rois(r));
	end
end

%% add panel with views of the Rx on each plane
if ishandle(rx.ui.rxAxes)
	if notDefined('loc'),   loc = round( rx.volDims / 2 );  end
	if ischar(loc),         loc = str2num(loc);             end

	%% open a panel to the right of the main figure w/ the Rxs
	hp = mrvPanel('right', .3);

	%% put up the prescriptions
	for ori = 1:3
		vol = rx.vol;
		volSlice = uint8(loc(ori));

		%% orient the slice properly
		if ori~=3   % do nothing for sagittal view
			% allow for radiological L/R
			hRadiological = findobj('Tag','rxRadiologicalMenu');
			if isequal(get(hRadiological, 'Checked'), 'on')
				vol = flipdim(vol, 3);
			end

			% permute volume as needed
			if ori==1, vol = permute(vol,[2 3 1]); end   % axial
			if ori==2, vol = permute(vol,[1 3 2]); end   % coronal
		end

		%% plot the slice
		% get slice
		volImg = vol(:,:,volSlice);
		volImg = rxClip(volImg, [], rx.ui.volBright, rx.ui.volContrast);

		% make subplot
		pos = [0 1-(ori/3) 1 .3];
		hax(ori) = axes('Parent', hp, 'Units', 'norm', 'Position', pos);

		% put up image
		htmp = image(volImg); colormap gray; axis off; axis equal;

		% draw Rx, label directions
		rxDrawRx(rx, volSlice, ori);  hold on
		rxLabelVolAxes(rx, ori);

		%% draw ROIs on prescription, if selected
		if ~isempty(rx.rois) & get(rx.ui.rxDrawRx, 'Value')==1
			for r = 1:length(rx.rois)
				R  = rx.rois(r);  % ROI struct
				C  = R.volCoords; % coords (being terse b/c of indents)

				% map the [axi cor sag] positions in C to the
				% 2D [x y] positions in pts:
				inSlice = find( round(C(ori,:)) == volSlice );
				switch ori
					case 1, % axi
						pts = C([2 3],inSlice);
					case 2, % cor
						pts = C([1 3],inSlice);
					case 3, % sag
						pts = C([1 2],inSlice);
				end

				% Draw the ROI outline
				if ~isempty(pts)
					outline(pts, R);
				end
			end
		end

	end
end

return
% /------------------------------------------------------------/ %

% /------------------------------------------------------------/ %
function pts = rxCoords2Montage(coords, nSlices, dims);
% convert 3-D volume coordinates into 2-D points on a montage
% of prescription slices.  This assumes that imageMontage has
% been called on all slices (1:nSlices), with the default #
% of columns and rows.
ncols = ceil( sqrt( nSlices ) );
nrows = ceil( nSlices / ncols );

% initialize pts output:
pts = [];

for slice = 1:nSlices
	% find (row, col) of this slice in the montage
	row = ceil(slice / ncols);
	col = mod(slice-1, ncols) + 1;

	% find columns of coords in this slice
	I = find( round(coords(3,:)) == slice );

	% main part: compute (x,y) locations given the montage offset
	y = coords(1,I) + (row-1) * dims(1);
	x = coords(2,I) + (col-1) * dims(2);

	% add to points list:
	% we will shuffle the order of coordinates in coords,
	% but remove points outside the Rx.
	pts = [pts [y; x]];
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% COMPHEADS - plot multiple TOPOPLOT maps of ICA component topographies
%
% Usage:
%       >> compheads(winv,'spline_file',compnos,'title',rowscols,labels,view)
%
% Inputs:
%   winv       - Inverse weight matrix = EEG scalp maps. Each column is a
%                map; the rows correspond to the electrode positions
%                defined in the eloc_file. Normally, winv = inv(weights*sphere).
%  spline_file - Name of the eloctrode position file in BESA spherical coords.
%   compnos    - Vector telling which (order of) component maps to show
%                Indices <0 tell COMPHEADS to invert a map; = 0 leave blank subplot
%                Example [1 0 -2 3 0 -6] {default|0 -> 1:columns_in_winv}
%  'title'     - Title string for each page {default|0 -> 'ICA Component Maps'}
%   rowscols   - Vector of the form [m,n] where m is total vertical tiles and n
%                is horizontal tiles per page. If the number of maps exceeds m*n,
%                multiple figures will be produced {def|0 -> one near-square page}.
%   labels     - Vector of numbers or a matrix of strings to use as labels for
%                each map {default|0 -> 1:ncolumns_in_winv}
%   view       - TOPOPLOT view, either [az el] or keyword ('top',...)
%                See >> help TOPOPLOT for options.
%
% Note: Map scaling is to +/-max(abs(data); green = 0
%
% Author: Scott Makeig, SCCN/INC/UCSD, La Jolla, 4-28-1998
%
% See also: TOPOPLOT

% Copyright (C) 4-28-98 from compmap.m Scott Makeig, SCCN/INC/UCSD, scott@sccn.ucsd.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. 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 HOLDER 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.

% 01-25-02 reformated help & license, added links -ad

function compheads(Winv,eloc_file,compnos,titleval,pagesize,srclabels,view)

if nargin<1
   help compheads
   return
end

[chans, frames] = size (Winv);

DEFAULT_TITLE = 'ICA Component Maps';
DEFAULT_EFILE = 'chan_file';
NUMCONTOUR = 5;     % topoplot() style settings
OUTPUT = 'screen';  % default: 'screen' for screen colors,
                    %          'printer' for printer colors
STYLE = 'both';
INTERPLIMITS = 'head';
MAPLIMITS = 'absmax';
SQUARE    = 1; % 1/0 flag making topoplot() asex square -> round heads

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% check inputs and set defaults
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

printlabel = OUTPUT; % default set above
if nargin < 7
   view = [-127 30];
end

if nargin < 6
   srclabels = 0;
end
if nargin < 5
   pagesize = 0;
end
if nargin < 4
   titleval = 0;
end
if nargin < 3
   compnos = 0;
end
if nargin < 2
   eloc_file = 0;
end

if srclabels == 0
  srclabels = [];
end
if titleval == 0;
   titleval = DEFAULT_TITLE;
end
if compnos == 0
   compnos = (1:frames);
end
if pagesize == 0
   numsources = length(compnos);
   DEFAULT_PAGE_SIZE = ...
[floor(sqrt(numsources)) ceil(numsources/floor(sqrt(numsources)))];
   m = DEFAULT_PAGE_SIZE(1);
   n = DEFAULT_PAGE_SIZE(2);
elseif length(pagesize) ==1
   help compheads
   return
else
   m = pagesize(1);
   n = pagesize(2);
end
if eloc_file == 0
   eloc_file = DEFAULT_EFILE;
end

totalsources = length(compnos);
if ~isempty(srclabels)
  if ~ischar(srclabels(1,1)) % if numbers
    if size(srclabels,1) == 1
       srclabels = srclabels';
    end
  end
  if size(srclabels,1) ~= totalsources,
     fprintf('compheads(): numbers of components and component labels do not agree.\n');
     return
  end
end
pages = ceil(totalsources/(m*n));
if pages > 1
   fprintf('compheads(): will create %d figures of %d by %d maps: ',...
            pages,m,n);
end

pos = get(gcf,'Position');
off = [ 25 -25 0 0];  % position offsets for multiple figures

fid = fopen(eloc_file);
if fid<1,
  fprintf('compheads()^G: cannot open eloc_file (%s).\n',eloc_file);
  return
end
fclose(fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%% plot the maps %%%%%%%%%%%%%%%%%%%%%%%

for i = (1:pages)
  if i > 1
    figure('Position',pos+(i-1)*off); % place figures in right-downward stack
  end
  set(gcf,'Color','w') %CJH - set background color to white

  if (totalsources > i*m*n)
    sbreak = n*m;
  else
    sbreak = totalsources - (i-1)*m*n;
  end

  for j = (1:sbreak) % maps on this page
    comp = j+(i-1)*m*n; % compno index
    if compnos(comp)~=0
      if compnos(comp)>0
       source_var = Winv(:,compnos(comp))';       % plot map
      elseif compnos(comp)<0
       source_var = -1*Winv(:,-1*compnos(comp))'; % invert map
      end

      subplot(m,n,j)
      headplot(source_var,eloc_file,'electrodes','off','view',view);
      %topoplot(source_var,eloc_file,'style',STYLE,...
      % 'numcontour',NUMCONTOUR,'interplimits',INTERPLIMITS,...
      % 'maplimits',MAPLIMITS); % draw map
      if SQUARE,
         axis('square');
      end

      if isempty(srclabels)
          t=title(int2str(compnos(comp)));
          set(t,'FontSize',16);
      else
        if ischar(srclabels)
          t=title(srclabels(comp,:));
          set(t,'FontSize',16);
        else
	      t=title(num2str(srclabels(comp)));
          set(t,'FontSize',16);
        end
      end
      drawnow  % draw one map at a time
    end
  end

  ax = axes('Units','Normal','Position',[.5 .04 .32 .05],'Visible','Off');
  colorbar(ax)
  axval = axis;
  Xlim = get(ax,'Xlim');
  set(ax,'XTick',(Xlim(2)+Xlim(1))/2)
  set(gca,'XTickLabel','0')
  set(gca,'XTickMode','manual')

  axbig = axes('Units','Normalized','Position',[0 0 1 1],'Visible','off');
  t1 = text(.25,.070,titleval,'HorizontalAlignment','center','FontSize',14);
  if pages > 1
     fprintf('%d ',i);
  end
end
if pages > 1
   fprintf('\n');
end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function im = RegdoCheckboard(Im1, Im2, rows, cols)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

    [m n p] = size(Im1);
    classname = class(Im1);
    %Im1 = cast(Im1, classname);
    %Im2 = cast(Im2, classname);

    m1 = fix(m/rows); n1 = fix(n/cols);
    black = zeros(m1, n1, classname);
    white = ones(m1, n1, classname);
    tile = [black white; white black];
    I = repmat(tile, [ceil(m/(2*m1)) ceil(n/(2*n1)) p]);
    c1 = I(1:m, 1:n, 1:p);
    c2 = (1-c1);

    CA_Image = c1.*(Im2) + c2.*(Im1);
    CA_Image(:,:,2) = CA_Image(:,:,1);
    CA_Image(:,:,3) = CA_Image(:,:,1);
    CA_Image = single(CA_Image);

    c1 = single(c1);
    CA_Image = (CA_Image - min(CA_Image(:))) / (max(CA_Image(:))-min(CA_Image(:)));
    CA_Image(:,:,2) = CA_Image(:,:,2) + c1*0.136;
    CA_Image(:,:,3) = CA_Image(:,:,3) + c1*0.136;
    CA_Image = min(CA_Image,1);

%     r = 0.05;
%     frac1 = (max(Im1(:))-min(Im1(:)))*r;
%     frac2 = (max(Im2(:))-min(Im2(:)))*r;
%     CA_Image = c1.*(Im2+frac2) + c2.*(Im1-frac1);

    im = CA_Image;

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
% Compressed Sensing Examples
%
% Copyright 2016, All Rights Reserved
% Code by Steven L. Brunton (sbrunton@uw.edu)

% Download and install CVX to run this example: http://cvxr.com/cvx/download/

clear all, close all, clc

x = sort(4*(rand(25,1)-.5)); % Random data from [-2,2]
b = .9*x + .1*randn(size(x)); % Line y=.9x with noise
atrue = x\b;          % Least-squares slope (no outliers)

b(end) = -5.5;                  % Introduce outlier
acorrupt = x\b;                 % New slope

cvx_begin;       % L1 optimization to reject outlier
    variable aL1;     % aL1 is slope to be optimized
    minimize( norm(aL1*x-b,1) );     % aL1 is robust
cvx_end;

hold on
scatter(x(1:end-1),b(1:end-1),'bo') % Data
scatter(x(end),b(end),'ro')         % Outlier
xgrid = -2:.01:2;
plot(xgrid,xgrid*atrue,'k--')       % L2 fit (no outlier)
plot(xgrid,xgrid*acorrupt,'r--')    % L2 fit (outlier)
plot(xgrid,xgrid*aL1,'b--')         % L1 fit
axis([-2 2 -6 2])

%%
set(gcf,'Position',[100 100 600 400])
set(gcf,'PaperPositionMode','auto')
print('-depsc2', '-loose', '../figures/f_chCS_robustregression');

author: dynamicslab


text:
function title = p16_title ( )

%*****************************************************************************80
%
%% P16_TITLE returns the title for problem 16.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 August 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string TITLE, the title of the problem.
%
  title = 'cos(pi x / 2 ) / sqrt(x)';

  return
end


text:
function st_to_mm_test01 ( )

%*****************************************************************************80
%
%% ST_TO_MM_TEST01 tests ST_TO_MM.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 November 2008
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'ST_TO_MM_TEST01\n' );
  fprintf ( 1, '  Use ST_TO_MM to read an ST file,\n' );
  fprintf ( 1, '  and write an equivalent MM file.\n' );

  input_file = 'st_to_mm_test01.st';
  output_file = 'st_to_mm_test01.mm';

  st_to_mm ( input_file, output_file );

  return
end


text:
function y = perform_curvelet_transform(x,options)

% perform_curvelet_transform - a wrapper to curvlab
%
%     M = perform_curvelet_transform(MW,options);
%
%   Forward and backward curvelet transform
%   You must provide options.n (width of the image).
%
%   Visit www.curvelab.org for the full code.

options.null = 0;
finest = getoptions(options, 'finest',1)  ; % =1(curv)  =2(wav)
nbscales = getoptions(options, 'nbscales', 5); % log2(size(M,1))-2;
nbangles_coarse = getoptions(options, 'nbangles_coarse', 16);
is_real = getoptions(options, 'is_real', 0);
n = getoptions(options, 'n', 1,1);

if not(iscell(x))
    % fwd transform
    y = fdct_wrapping(x, is_real, finest, nbscales, nbangles_coarse);
else
    y = real( ifdct_wrapping(x, is_real, n,n ) );
end

%%

function C = fdct_wrapping(x, is_real, finest, nbscales, nbangles_coarse)

% fdct_wrapping.m - Fast Discrete Curvelet Transform via wedge wrapping - Version 1.0
%
% Inputs
%   x           M-by-N matrix
%
% Optional Inputs
%   is_real     Type of the transform
%                   0: complex-valued curvelets
%                   1: real-valued curvelets
%               [default set to 0]
%   finest      Chooses one of two possibilities for the coefficients at the
%               finest level:
%                   1: curvelets
%                   2: wavelets
%               [default set to 2]
%   nbscales    number of scales including the coarsest wavelet level
%               [default set to ceil(log2(min(M,N)) - 3)]
%   nbangles_coarse
%               number of angles at the 2nd coarsest level, minimum 8,
%               must be a multiple of 4. [default set to 16]
%
% Outputs
%   C           Cell array of curvelet coefficients.
%               C{j}{l}(k1,k2) is the coefficient at
%                   - scale j: integer, from finest to coarsest scale,
%                   - angle l: integer, starts at the top-left corner and
%                   increases clockwise,
%                   - position k1,k2: both integers, size varies with j
%                   and l.
%               If is_real is 1, there are two types of curvelets,
%               'cosine' and 'sine'. For a given scale j, the 'cosine'
%               coefficients are stored in the first two quadrants (low
%               values of l), the 'sine' coefficients in the last two
%               quadrants (high values of l).
%
% See also ifdct_wrapping.m, fdct_wrapping_param.m
%
% By Laurent Demanet, 2004

X = fftshift(fft2(ifftshift(x)))/sqrt(prod(size(x)));
[N1,N2] = size(X);
if nargin < 2, is_real = 0; end;
if nargin < 3, finest = 2; end;
if nargin < 4, nbscales = ceil(log2(min(N1,N2)) - 3); end;
if nargin < 5, nbangles_coarse = 16; end;

% Initialization: data structure
nbangles = [1, nbangles_coarse .* 2.^(ceil((nbscales-(nbscales:-1:2))/2))];
if finest == 2, nbangles(nbscales) = 1; end;
C = cell(1,nbscales);
for j = 1:nbscales
    C{j} = cell(1,nbangles(j));
end;

% Loop: pyramidal scale decomposition
M1 = N1/3;
M2 = N2/3;
if finest == 1,

    % Initialization: smooth periodic extension of high frequencies
    bigN1 = 2*floor(2*M1)+1;
    bigN2 = 2*floor(2*M2)+1;
    equiv_index_1 = 1+mod(floor(N1/2)-floor(2*M1)+(1:bigN1)-1,N1);
    equiv_index_2 = 1+mod(floor(N2/2)-floor(2*M2)+(1:bigN2)-1,N2);
    X = X(equiv_index_1,equiv_index_2);
        % Invariant: equiv_index_1(floor(2*M1)+1) == (N1 + 2 - mod(N1,2))/2
        % is the center in frequency. Same for M2, N2.
    window_length_1 = floor(2*M1) - floor(M1) - 1 - (mod(N1,3)==0);
    window_length_2 = floor(2*M2) - floor(M2) - 1 - (mod(N2,3)==0);
        % Invariant: floor(M1) + floor(2*M1) == N1 - (mod(M1,3)~=0)
        % Same for M2, N2.
    coord_1 = 0:(1/window_length_1):1;
    coord_2 = 0:(1/window_length_2):1;
    [wl_1,wr_1] = fdct_wrapping_window(coord_1);
    [wl_2,wr_2] = fdct_wrapping_window(coord_2);
    lowpass_1 = [wl_1, ones(1,2*floor(M1)+1), wr_1];
    if mod(N1,3)==0, lowpass_1 = [0, lowpass_1, 0]; end;
    lowpass_2 = [wl_2, ones(1,2*floor(M2)+1), wr_2];
    if mod(N2,3)==0, lowpass_2 = [0, lowpass_2, 0]; end;
    lowpass = lowpass_1'*lowpass_2;
    Xlow = X .* lowpass;

    scales = nbscales:-1:2;

else

    M1 = M1/2;
    M2 = M2/2;
    window_length_1 = floor(2*M1) - floor(M1) - 1;
    window_length_2 = floor(2*M2) - floor(M2) - 1;
    coord_1 = 0:(1/window_length_1):1;
    coord_2 = 0:(1/window_length_2):1;
    [wl_1,wr_1] = fdct_wrapping_window(coord_1);
    [wl_2,wr_2] = fdct_wrapping_window(coord_2);
    lowpass_1 = [wl_1, ones(1,2*floor(M1)+1), wr_1];
    lowpass_2 = [wl_2, ones(1,2*floor(M2)+1), wr_2];
    lowpass = lowpass_1'*lowpass_2;
    hipass = sqrt(1 - lowpass.^2);
    Xlow_index_1 = ((-floor(2*M1)):floor(2*M1)) + ceil((N1+1)/2);
    Xlow_index_2 = ((-floor(2*M2)):floor(2*M2)) + ceil((N2+1)/2);
    Xlow = X(Xlow_index_1, Xlow_index_2) .* lowpass;
    Xhi = X;
    Xhi(Xlow_index_1, Xlow_index_2) = Xhi(Xlow_index_1, Xlow_index_2) .* hipass;
    C{nbscales}{1} = fftshift(ifft2(ifftshift(Xhi)))*sqrt(prod(size(Xhi)));
    if is_real, C{nbscales}{1} = real(C{nbscales}{1}); end;

    scales = (nbscales-1):-1:2;

end;
for j = scales,

    M1 = M1/2;
    M2 = M2/2;
    window_length_1 = floor(2*M1) - floor(M1) - 1;
    window_length_2 = floor(2*M2) - floor(M2) - 1;
    coord_1 = 0:(1/window_length_1):1;
    coord_2 = 0:(1/window_length_2):1;
    [wl_1,wr_1] = fdct_wrapping_window(coord_1);
    [wl_2,wr_2] = fdct_wrapping_window(coord_2);
    lowpass_1 = [wl_1, ones(1,2*floor(M1)+1), wr_1];
    lowpass_2 = [wl_2, ones(1,2*floor(M2)+1), wr_2];
    lowpass = lowpass_1'*lowpass_2;
    hipass = sqrt(1 - lowpass.^2);
    Xhi = Xlow;                 % size is 2*floor(4*M1)+1 - by - 2*floor(4*M2)+1
    Xlow_index_1 = ((-floor(2*M1)):floor(2*M1)) + floor(4*M1) + 1;
    Xlow_index_2 = ((-floor(2*M2)):floor(2*M2)) + floor(4*M2) + 1;
    Xlow = Xlow(Xlow_index_1, Xlow_index_2);
    Xhi(Xlow_index_1, Xlow_index_2) = Xlow .* hipass;
    Xlow = Xlow .* lowpass;     % size is 2*floor(2*M1)+1 - by - 2*floor(2*M2)+1

    % Loop: angular decomposition
    l = 0;
    nbquadrants = 2 + 2*(~is_real);
    nbangles_perquad = nbangles(j)/4;
    for quadrant = 1:nbquadrants
        M_horiz = M2 * (mod(quadrant,2)==1) + M1 * (mod(quadrant,2)==0);
        M_vert = M1 * (mod(quadrant,2)==1) + M2 * (mod(quadrant,2)==0);
        if mod(nbangles_perquad,2),
            wedge_ticks_left = round((0:(1/(2*nbangles_perquad)):.5)*2*floor(4*M_horiz) + 1);
            wedge_ticks_right = 2*floor(4*M_horiz) + 2 - wedge_ticks_left;
            wedge_ticks = [wedge_ticks_left, wedge_ticks_right(end:-1:1)];
        else
            wedge_ticks_left = round((0:(1/(2*nbangles_perquad)):.5)*2*floor(4*M_horiz) + 1);
            wedge_ticks_right = 2*floor(4*M_horiz) + 2 - wedge_ticks_left;
            wedge_ticks = [wedge_ticks_left, wedge_ticks_right((end-1):-1:1)];
        end;
        wedge_endpoints = wedge_ticks(2:2:(end-1));         % integers
        wedge_midpoints = (wedge_endpoints(1:(end-1)) + wedge_endpoints(2:end))/2;
                % integers or half-integers

        % Left corner wedge
        l = l+1;
        first_wedge_endpoint_vert = round(2*floor(4*M_vert)/(2*nbangles_perquad) + 1);
        length_corner_wedge = floor(4*M_vert) - floor(M_vert) + ceil(first_wedge_endpoint_vert/4);
        Y_corner = 1:length_corner_wedge;
        [XX,YY] = meshgrid(1:(2*floor(4*M_horiz)+1),Y_corner);
        width_wedge = wedge_endpoints(2) + wedge_endpoints(1) - 1;
        slope_wedge = (floor(4*M_horiz) + 1 - wedge_endpoints(1))/floor(4*M_vert);
        left_line = round(2 - wedge_endpoints(1) + slope_wedge*(Y_corner - 1));
                                                            % integers
        [wrapped_data, wrapped_XX, wrapped_YY] = deal(zeros(length_corner_wedge,width_wedge));
        first_row = floor(4*M_vert)+2-ceil((length_corner_wedge+1)/2)+...
            mod(length_corner_wedge+1,2)*(quadrant-2 == mod(quadrant-2,2));
        first_col = floor(4*M_horiz)+2-ceil((width_wedge+1)/2)+...
            mod(width_wedge+1,2)*(quadrant-3 == mod(quadrant-3,2));
                % Coordinates of the top-left corner of the wedge wrapped
                % around the origin. Some subtleties when the wedge is
                % even-sized because of the forthcoming 90 degrees rotation
        for row = Y_corner
            cols = left_line(row) + mod((0:(width_wedge-1))-(left_line(row)-first_col),width_wedge);
            admissible_cols = round(1/2*(cols+1+abs(cols-1)));
            new_row = 1 + mod(row - first_row, length_corner_wedge);
            wrapped_data(new_row,:) = Xhi(row,admissible_cols) .* (cols > 0);
            wrapped_XX(new_row,:) = XX(row,admissible_cols);
            wrapped_YY(new_row,:) = YY(row,admissible_cols);
        end;
        slope_wedge_right = (floor(4*M_horiz)+1 - wedge_midpoints(1))/floor(4*M_vert);
        mid_line_right = wedge_midpoints(1) + slope_wedge_right*(wrapped_YY - 1);
                % not integers in general
        coord_right = 1/2 + floor(4*M_vert)/(wedge_endpoints(2) - wedge_endpoints(1)) * ...
            (wrapped_XX - mid_line_right)./(floor(4*M_vert)+1 - wrapped_YY);
        C2 = 1/(1/(2*(floor(4*M_horiz))/(wedge_endpoints(1) - 1) - 1) + 1/(2*(floor(4*M_vert))/(first_wedge_endpoint_vert - 1) - 1));
        C1 = C2 / (2*(floor(4*M_vert))/(first_wedge_endpoint_vert - 1) - 1);
        wrapped_XX((wrapped_XX - 1)/floor(4*M_horiz) + (wrapped_YY-1)/floor(4*M_vert) == 2) = ...
            wrapped_XX((wrapped_XX - 1)/floor(4*M_horiz) + (wrapped_YY-1)/floor(4*M_vert) == 2) + 1;
        coord_corner = C1 + C2 * ((wrapped_XX - 1)/(floor(4*M_horiz)) - (wrapped_YY - 1)/(floor(4*M_vert))) ./ ...
            (2-((wrapped_XX - 1)/(floor(4*M_horiz)) + (wrapped_YY - 1)/(floor(4*M_vert))));
        wl_left = fdct_wrapping_window(coord_corner);
        [wl_right,wr_right] = fdct_wrapping_window(coord_right);
        wrapped_data = wrapped_data .* (wl_left .* wr_right);

        switch is_real
            case 0
                wrapped_data = rot90(wrapped_data,-(quadrant-1));
                C{j}{l} = fftshift(ifft2(ifftshift(wrapped_data)))*sqrt(prod(size(wrapped_data)));
            case 1
                wrapped_data = rot90(wrapped_data,-(quadrant-1));
                x = fftshift(ifft2(ifftshift(wrapped_data)))*sqrt(prod(size(wrapped_data)));
                C{j}{l} = sqrt(2)*real(x);
                C{j}{l+nbangles(j)/2} = sqrt(2)*imag(x);
        end;

        % Regular wedges
        length_wedge = floor(4*M_vert) - floor(M_vert);
        Y = 1:length_wedge;
        first_row = floor(4*M_vert)+2-ceil((length_wedge+1)/2)+...
            mod(length_wedge+1,2)*(quadrant-2 == mod(quadrant-2,2));
        for subl = 2:(nbangles_perquad-1);
            l = l+1;
            width_wedge = wedge_endpoints(subl+1) - wedge_endpoints(subl-1) + 1;
            slope_wedge = ((floor(4*M_horiz)+1) - wedge_endpoints(subl))/floor(4*M_vert);
            left_line = round(wedge_endpoints(subl-1) + slope_wedge*(Y - 1));
            [wrapped_data, wrapped_XX, wrapped_YY] = deal(zeros(length_wedge,width_wedge));
            first_col = floor(4*M_horiz)+2-ceil((width_wedge+1)/2)+...
                mod(width_wedge+1,2)*(quadrant-3 == mod(quadrant-3,2));
            for row = Y
                cols = left_line(row) + mod((0:(width_wedge-1))-(left_line(row)-first_col),width_wedge);
                new_row = 1 + mod(row - first_row, length_wedge);
                wrapped_data(new_row,:) = Xhi(row,cols);
                wrapped_XX(new_row,:) = XX(row,cols);
                wrapped_YY(new_row,:) = YY(row,cols);
            end;
            slope_wedge_left = ((floor(4*M_horiz)+1) - wedge_midpoints(subl-1))/floor(4*M_vert);
            mid_line_left = wedge_midpoints(subl-1) + slope_wedge_left*(wrapped_YY - 1);
            coord_left = 1/2 + floor(4*M_vert)/(wedge_endpoints(subl) - wedge_endpoints(subl-1)) * ...
                (wrapped_XX - mid_line_left)./(floor(4*M_vert)+1 - wrapped_YY);
            slope_wedge_right = ((floor(4*M_horiz)+1) - wedge_midpoints(subl))/floor(4*M_vert);
            mid_line_right = wedge_midpoints(subl) + slope_wedge_right*(wrapped_YY - 1);
            coord_right = 1/2 + floor(4*M_vert)/(wedge_endpoints(subl+1) - wedge_endpoints(subl)) * ...
                (wrapped_XX - mid_line_right)./(floor(4*M_vert)+1 - wrapped_YY);
            wl_left = fdct_wrapping_window(coord_left);
            [wl_right,wr_right] = fdct_wrapping_window(coord_right);
            wrapped_data = wrapped_data .* (wl_left .* wr_right);
            switch is_real
                case 0
                    wrapped_data = rot90(wrapped_data,-(quadrant-1));
                    C{j}{l} = fftshift(ifft2(ifftshift(wrapped_data)))*sqrt(prod(size(wrapped_data)));
                case 1
                    wrapped_data = rot90(wrapped_data,-(quadrant-1));
                    x = fftshift(ifft2(ifftshift(wrapped_data)))*sqrt(prod(size(wrapped_data)));
                    C{j}{l} = sqrt(2)*real(x);
                    C{j}{l+nbangles(j)/2} = sqrt(2)*imag(x);
            end;
        end;

        % Right corner wedge
        l = l+1;
        width_wedge = 4*floor(4*M_horiz) + 3 - wedge_endpoints(end) - wedge_endpoints(end-1);
        slope_wedge = ((floor(4*M_horiz)+1) - wedge_endpoints(end))/floor(4*M_vert);
        left_line = round(wedge_endpoints(end-1) + slope_wedge*(Y_corner - 1));
        [wrapped_data, wrapped_XX, wrapped_YY] = deal(zeros(length_corner_wedge,width_wedge));
        first_row = floor(4*M_vert)+2-ceil((length_corner_wedge+1)/2)+...
            mod(length_corner_wedge+1,2)*(quadrant-2 == mod(quadrant-2,2));
        first_col = floor(4*M_horiz)+2-ceil((width_wedge+1)/2)+...
            mod(width_wedge+1,2)*(quadrant-3 == mod(quadrant-3,2));
        for row = Y_corner
            cols = left_line(row) + mod((0:(width_wedge-1))-(left_line(row)-first_col),width_wedge);
            admissible_cols = round(1/2*(cols+2*floor(4*M_horiz)+1-abs(cols-(2*floor(4*M_horiz)+1))));
            new_row = 1 + mod(row - first_row, length_corner_wedge);
            wrapped_data(new_row,:) = Xhi(row,admissible_cols) .* (cols <= (2*floor(4*M_horiz)+1));
            wrapped_XX(new_row,:) = XX(row,admissible_cols);
            wrapped_YY(new_row,:) = YY(row,admissible_cols);
        end;
        slope_wedge_left = ((floor(4*M_horiz)+1) - wedge_midpoints(end))/floor(4*M_vert);
        mid_line_left = wedge_midpoints(end) + slope_wedge_left*(wrapped_YY - 1);
        coord_left = 1/2 + floor(4*M_vert)/(wedge_endpoints(end) - wedge_endpoints(end-1)) * ...
            (wrapped_XX - mid_line_left)./(floor(4*M_vert) + 1 - wrapped_YY);
        C2 = -1/(2*(floor(4*M_horiz))/(wedge_endpoints(end) - 1) - 1 + 1/(2*(floor(4*M_vert))/(first_wedge_endpoint_vert - 1) - 1));
        C1 = -C2 * (2*(floor(4*M_horiz))/(wedge_endpoints(end) - 1) - 1);
        wrapped_XX((wrapped_XX - 1)/floor(4*M_horiz) == (wrapped_YY - 1)/floor(4*M_vert)) = ...
            wrapped_XX((wrapped_XX - 1)/floor(4*M_horiz) == (wrapped_YY - 1)/floor(4*M_vert)) - 1;
        coord_corner = C1 + C2 * (2-((wrapped_XX - 1)/(floor(4*M_horiz)) + (wrapped_YY - 1)/(floor(4*M_vert)))) ./ ...
            ((wrapped_XX - 1)/(floor(4*M_horiz)) - (wrapped_YY - 1)/(floor(4*M_vert)));
        wl_left = fdct_wrapping_window(coord_left);
        [wl_right,wr_right] = fdct_wrapping_window(coord_corner);

        wrapped_data = wrapped_data .* (wl_left .* wr_right);
        switch is_real
            case 0
                wrapped_data = rot90(wrapped_data,-(quadrant-1));
                C{j}{l} = fftshift(ifft2(ifftshift(wrapped_data)))*sqrt(prod(size(wrapped_data)));
            case 1
                wrapped_data = rot90(wrapped_data,-(quadrant-1));
                x = fftshift(ifft2(ifftshift(wrapped_data)))*sqrt(prod(size(wrapped_data)));
                C{j}{l} = sqrt(2)*real(x);
                C{j}{l+nbangles(j)/2} = sqrt(2)*imag(x);
        end;

        if quadrant < nbquadrants, Xhi = rot90(Xhi); end;
    end;
end;

% Coarsest wavelet level
C{1}{1} = fftshift(ifft2(ifftshift(Xlow)))*sqrt(prod(size(Xlow)));
if is_real == 1,
    C{1}{1} = real(C{1}{1});
end;

function x = ifdct_wrapping(C, is_real, M, N)

% ifdct_wrapping.m - Inverse Fast Discrete Curvelet Transform via wedge wrapping - Version 1.0
% This is in fact the adjoint, also the pseudo-inverse
%
% Inputs
%   C           Cell array containing curvelet coefficients (see
%               description in fdct_wrapping.m)
%   is_real     As used in fdct_wrapping.m
%   M, N        Size of the image to be recovered (not necessary if finest
%               = 2)
%
% Outputs
%   x           M-by-N matrix
%
% See also fdct_wrapping.m
%
% By Laurent Demanet, 2004

% Initialization
nbscales = length(C);
nbangles_coarse = length(C{2});
nbangles = [1, nbangles_coarse .* 2.^(ceil((nbscales-(nbscales:-1:2))/2))];
if length(C{end}) == 1, finest = 2; else finest = 1; end;
if finest == 2, nbangles(nbscales) = 1; end;
if nargin < 2, is_real = 0; end;
if nargin < 4,
    if finest == 1, error('Syntax: IFCT_wrapping(C,M,N) where the matrix to be recovered is M-by-N'); end;
    [N1,N2] = size(C{end}{1});
else
    N1 = M;
    N2 = N;
end;

M1 = N1/3;
M2 = N2/3;

if finest == 1;

    bigN1 = 2*floor(2*M1)+1;
    bigN2 = 2*floor(2*M2)+1;
    X = zeros(bigN1,bigN2);

    % Initialization: preparing the lowpass filter at finest scale
    window_length_1 = floor(2*M1) - floor(M1) - 1 - (mod(N1,3)==0);
    window_length_2 = floor(2*M2) - floor(M2) - 1 - (mod(N2,3)==0);
    coord_1 = 0:(1/window_length_1):1;
    coord_2 = 0:(1/window_length_2):1;
    [wl_1,wr_1] = fdct_wrapping_window(coord_1);
    [wl_2,wr_2] = fdct_wrapping_window(coord_2);
    lowpass_1 = [wl_1, ones(1,2*floor(M1)+1), wr_1];
    if mod(N1,3)==0, lowpass_1 = [0, lowpass_1, 0]; end;
    lowpass_2 = [wl_2, ones(1,2*floor(M2)+1), wr_2];
    if mod(N2,3)==0, lowpass_2 = [0, lowpass_2, 0]; end;
    lowpass = lowpass_1'*lowpass_2;

    scales = nbscales:-1:2;

else

    M1 = M1/2;
    M2 = M2/2;

    bigN1 = 2*floor(2*M1)+1;
    bigN2 = 2*floor(2*M2)+1;
    X = zeros(bigN1,bigN2);

    window_length_1 = floor(2*M1) - floor(M1) - 1;
    window_length_2 = floor(2*M2) - floor(M2) - 1;
    coord_1 = 0:(1/window_length_1):1;
    coord_2 = 0:(1/window_length_2):1;
    [wl_1,wr_1] = fdct_wrapping_window(coord_1);
    [wl_2,wr_2] = fdct_wrapping_window(coord_2);
    lowpass_1 = [wl_1, ones(1,2*floor(M1)+1), wr_1];
    lowpass_2 = [wl_2, ones(1,2*floor(M2)+1), wr_2];
    lowpass = lowpass_1'*lowpass_2;
    hipass_finest = sqrt(1 - lowpass.^2);

    scales = (nbscales-1):-1:2;

end;

% Loop: pyramidal reconstruction

Xj_topleft_1 = 1;
Xj_topleft_2 = 1;
for j = scales,

    M1 = M1/2;
    M2 = M2/2;
    window_length_1 = floor(2*M1) - floor(M1) - 1;
    window_length_2 = floor(2*M2) - floor(M2) - 1;
    coord_1 = 0:(1/window_length_1):1;
    coord_2 = 0:(1/window_length_2):1;
    [wl_1,wr_1] = fdct_wrapping_window(coord_1);
    [wl_2,wr_2] = fdct_wrapping_window(coord_2);
    lowpass_1 = [wl_1, ones(1,2*floor(M1)+1), wr_1];
    lowpass_2 = [wl_2, ones(1,2*floor(M2)+1), wr_2];
    lowpass_next = lowpass_1'*lowpass_2;
    hipass = sqrt(1 - lowpass_next.^2);
    Xj = zeros(2*floor(4*M1)+1,2*floor(4*M2)+1);

    % Loop: angles
    l = 0;
    nbquadrants = 2 + 2*(~is_real);
    nbangles_perquad = nbangles(j)/4;
    for quadrant = 1:nbquadrants

        M_horiz = M2 * (mod(quadrant,2)==1) + M1 * (mod(quadrant,2)==0);
        M_vert = M1 * (mod(quadrant,2)==1) + M2 * (mod(quadrant,2)==0);
        if mod(nbangles_perquad,2),
            wedge_ticks_left = round((0:(1/(2*nbangles_perquad)):.5)*2*floor(4*M_horiz) + 1);
            wedge_ticks_right = 2*floor(4*M_horiz) + 2 - wedge_ticks_left;
            wedge_ticks = [wedge_ticks_left, wedge_ticks_right(end:-1:1)];
        else
            wedge_ticks_left = round((0:(1/(2*nbangles_perquad)):.5)*2*floor(4*M_horiz) + 1);
            wedge_ticks_right = 2*floor(4*M_horiz) + 2 - wedge_ticks_left;
            wedge_ticks = [wedge_ticks_left, wedge_ticks_right((end-1):-1:1)];
        end;
        wedge_endpoints = wedge_ticks(2:2:(end-1));         % integers
        wedge_midpoints = (wedge_endpoints(1:(end-1)) + wedge_endpoints(2:end))/2;

        % Left corner wedge

        l = l+1;
        first_wedge_endpoint_vert = round(2*floor(4*M_vert)/(2*nbangles_perquad) + 1);
        length_corner_wedge = floor(4*M_vert) - floor(M_vert) + ceil(first_wedge_endpoint_vert/4);
        Y_corner = 1:length_corner_wedge;
        [XX,YY] = meshgrid(1:(2*floor(4*M_horiz)+1),Y_corner);
        width_wedge = wedge_endpoints(2) + wedge_endpoints(1) - 1;
        slope_wedge = (floor(4*M_horiz) + 1 - wedge_endpoints(1))/floor(4*M_vert);
        left_line = round(2 - wedge_endpoints(1) + slope_wedge*(Y_corner - 1));
        [wrapped_XX, wrapped_YY] = deal(zeros(length_corner_wedge,width_wedge));
        first_row = floor(4*M_vert)+2-ceil((length_corner_wedge+1)/2)+...
            mod(length_corner_wedge+1,2)*(quadrant-2 == mod(quadrant-2,2));
        first_col = floor(4*M_horiz)+2-ceil((width_wedge+1)/2)+...
            mod(width_wedge+1,2)*(quadrant-3 == mod(quadrant-3,2));

        for row = Y_corner
            cols = left_line(row) + mod((0:(width_wedge-1))-(left_line(row)-first_col),width_wedge);
            new_row = 1 + mod(row - first_row, length_corner_wedge);
            admissible_cols = round(1/2*(cols+1+abs(cols-1)));
            wrapped_XX(new_row,:) = XX(row,admissible_cols);
            wrapped_YY(new_row,:) = YY(row,admissible_cols);
        end;

        slope_wedge_right = (floor(4*M_horiz)+1 - wedge_midpoints(1))/floor(4*M_vert);
        mid_line_right = wedge_midpoints(1) + slope_wedge_right*(wrapped_YY - 1);
                                                            % not integers
                                                            % in general
        coord_right = 1/2 + floor(4*M_vert)/(wedge_endpoints(2) - wedge_endpoints(1)) * ...
            (wrapped_XX - mid_line_right)./(floor(4*M_vert)+1 - wrapped_YY);
        C2 = 1/(1/(2*(floor(4*M_horiz))/(wedge_endpoints(1) - 1) - 1) + 1/(2*(floor(4*M_vert))/(first_wedge_endpoint_vert - 1) - 1));
        C1 = C2 / (2*(floor(4*M_vert))/(first_wedge_endpoint_vert - 1) - 1);
        wrapped_XX((wrapped_XX - 1)/floor(4*M_horiz) + (wrapped_YY-1)/floor(4*M_vert) == 2) = ...
            wrapped_XX((wrapped_XX - 1)/floor(4*M_horiz) + (wrapped_YY-1)/floor(4*M_vert) == 2) + 1;
        coord_corner = C1 + C2 * ((wrapped_XX - 1)/(floor(4*M_horiz)) - (wrapped_YY - 1)/(floor(4*M_vert))) ./ ...
            (2-((wrapped_XX - 1)/(floor(4*M_horiz)) + (wrapped_YY - 1)/(floor(4*M_vert))));
        wl_left = fdct_wrapping_window(coord_corner);
        [wl_right,wr_right] = fdct_wrapping_window(coord_right);
        switch is_real
         case 0
          wrapped_data = fftshift(fft2(ifftshift(C{j}{l})))/sqrt(prod(size(C{j}{l})));
          wrapped_data = rot90(wrapped_data,(quadrant-1));
         case 1
          x = C{j}{l} + sqrt(-1)*C{j}{l+nbangles(j)/2};
          wrapped_data = fftshift(fft2(ifftshift(x)))/sqrt(prod(size(x)))/sqrt(2);
          wrapped_data = rot90(wrapped_data,(quadrant-1));
        end;
        wrapped_data = wrapped_data .* (wl_left .* wr_right);

        % Unwrapping data
        for row = Y_corner
            cols = left_line(row) + mod((0:(width_wedge-1))-(left_line(row)-first_col),width_wedge);
            admissible_cols = round(1/2*(cols+1+abs(cols-1)));
            new_row = 1 + mod(row - first_row, length_corner_wedge);
            Xj(row,admissible_cols) = Xj(row,admissible_cols) + wrapped_data(new_row,:);
                                % We use the following property: in an assignment
                                % A(B) = C where B and C are vectors, if
                                % some value x repeats in B, then the
                                % last occurrence of x is the one
                                % corresponding to the eventual assignment.
        end;

        % Regular wedges
        length_wedge = floor(4*M_vert) - floor(M_vert);
        Y = 1:length_wedge;
        first_row = floor(4*M_vert)+2-ceil((length_wedge+1)/2)+...
            mod(length_wedge+1,2)*(quadrant-2 == mod(quadrant-2,2));
        for subl = 2:(nbangles_perquad-1);
            l = l+1;
            width_wedge = wedge_endpoints(subl+1) - wedge_endpoints(subl-1) + 1;
            slope_wedge = ((floor(4*M_horiz)+1) - wedge_endpoints(subl))/floor(4*M_vert);
            left_line = round(wedge_endpoints(subl-1) + slope_wedge*(Y - 1));
            [wrapped_XX, wrapped_YY] = deal(zeros(length_wedge,width_wedge));
            first_col = floor(4*M_horiz)+2-ceil((width_wedge+1)/2)+...
                mod(width_wedge+1,2)*(quadrant-3 == mod(quadrant-3,2));
            for row = Y
                cols = left_line(row) + mod((0:(width_wedge-1))-(left_line(row)-first_col),width_wedge);
                new_row = 1 + mod(row - first_row, length_wedge);
                wrapped_XX(new_row,:) = XX(row,cols);
                wrapped_YY(new_row,:) = YY(row,cols);
            end;
            slope_wedge_left = ((floor(4*M_horiz)+1) - wedge_midpoints(subl-1))/floor(4*M_vert);
            mid_line_left = wedge_midpoints(subl-1) + slope_wedge_left*(wrapped_YY - 1);
            coord_left = 1/2 + floor(4*M_vert)/(wedge_endpoints(subl) - wedge_endpoints(subl-1)) * ...
                (wrapped_XX - mid_line_left)./(floor(4*M_vert)+1 - wrapped_YY);
            slope_wedge_right = ((floor(4*M_horiz)+1) - wedge_midpoints(subl))/floor(4*M_vert);
            mid_line_right = wedge_midpoints(subl) + slope_wedge_right*(wrapped_YY - 1);
            coord_right = 1/2 + floor(4*M_vert)/(wedge_endpoints(subl+1) - wedge_endpoints(subl)) * ...
                (wrapped_XX - mid_line_right)./(floor(4*M_vert)+1 - wrapped_YY);
            wl_left = fdct_wrapping_window(coord_left);
            [wl_right,wr_right] = fdct_wrapping_window(coord_right);
            switch is_real
             case 0
              wrapped_data = fftshift(fft2(ifftshift(C{j}{l})))/sqrt(prod(size(C{j}{l})));
              wrapped_data = rot90(wrapped_data,(quadrant-1));
             case 1
              x = C{j}{l} + sqrt(-1)*C{j}{l+nbangles(j)/2};
              wrapped_data = fftshift(fft2(ifftshift(x)))/sqrt(prod(size(x)))/sqrt(2);
              wrapped_data = rot90(wrapped_data,(quadrant-1));
            end;
            wrapped_data = wrapped_data .* (wl_left .* wr_right);

            % Unwrapping data
            for row = Y
                cols = left_line(row) + mod((0:(width_wedge-1))-(left_line(row)-first_col),width_wedge);
                new_row = 1 + mod(row - first_row, length_wedge);
                Xj(row,cols) = Xj(row,cols) + wrapped_data(new_row,:);
            end;

        end;    % for subl

        % Right corner wedge
        l = l+1;
        width_wedge = 4*floor(4*M_horiz) + 3 - wedge_endpoints(end) - wedge_endpoints(end-1);
        slope_wedge = ((floor(4*M_horiz)+1) - wedge_endpoints(end))/floor(4*M_vert);
        left_line = round(wedge_endpoints(end-1) + slope_wedge*(Y_corner - 1));
        [wrapped_XX, wrapped_YY] = deal(zeros(length_corner_wedge,width_wedge));
        first_row = floor(4*M_vert)+2-ceil((length_corner_wedge+1)/2)+...
            mod(length_corner_wedge+1,2)*(quadrant-2 == mod(quadrant-2,2));
        first_col = floor(4*M_horiz)+2-ceil((width_wedge+1)/2)+...
            mod(width_wedge+1,2)*(quadrant-3 == mod(quadrant-3,2));

        for row = Y_corner
            cols = left_line(row) + mod((0:(width_wedge-1))-(left_line(row)-first_col),width_wedge);
            admissible_cols = round(1/2*(cols+2*floor(4*M_horiz)+1-abs(cols-(2*floor(4*M_horiz)+1))));
            new_row = 1 + mod(row - first_row, length_corner_wedge);
            wrapped_XX(new_row,:) = XX(row,admissible_cols);
            wrapped_YY(new_row,:) = YY(row,admissible_cols);
        end;
        YY = Y_corner'*ones(1,width_wedge);
        slope_wedge_left = ((floor(4*M_horiz)+1) - wedge_midpoints(end))/floor(4*M_vert);
        mid_line_left = wedge_midpoints(end) + slope_wedge_left*(wrapped_YY - 1);
        coord_left = 1/2 + floor(4*M_vert)/(wedge_endpoints(end) - wedge_endpoints(end-1)) * ...
            (wrapped_XX - mid_line_left)./(floor(4*M_vert)+1 - wrapped_YY);
        C2 = -1/(2*(floor(4*M_horiz))/(wedge_endpoints(end) - 1) - 1 + 1/(2*(floor(4*M_vert))/(first_wedge_endpoint_vert - 1) - 1));
        C1 = -C2 * (2*(floor(4*M_horiz))/(wedge_endpoints(end) - 1) - 1);
        wrapped_XX((wrapped_XX - 1)/floor(4*M_horiz) == (wrapped_YY-1)/floor(4*M_vert)) = ...
            wrapped_XX((wrapped_XX - 1)/floor(4*M_horiz) == (wrapped_YY-1)/floor(4*M_vert)) - 1;
        coord_corner = C1 + C2 * (2-((wrapped_XX - 1)/(floor(4*M_horiz)) + (wrapped_YY - 1)/(floor(4*M_vert)))) ./ ...
            ((wrapped_XX - 1)/(floor(4*M_horiz)) - (wrapped_YY - 1)/(floor(4*M_vert)));
        wl_left = fdct_wrapping_window(coord_left);
        [wl_right,wr_right] = fdct_wrapping_window(coord_corner);
        switch is_real
         case 0
          wrapped_data = fftshift(fft2(ifftshift(C{j}{l})))/sqrt(prod(size(C{j}{l})));
          wrapped_data = rot90(wrapped_data,(quadrant-1));
         case 1
          x = C{j}{l} + sqrt(-1)*C{j}{l+nbangles(j)/2};
          wrapped_data = fftshift(fft2(ifftshift(x)))/sqrt(prod(size(x)))/sqrt(2);
          wrapped_data = rot90(wrapped_data,(quadrant-1));
        end;
        wrapped_data = wrapped_data .* (wl_left .* wr_right);

         % Unwrapping data
        for row = Y_corner
            cols = left_line(row) + mod((0:(width_wedge-1))-(left_line(row)-first_col),width_wedge);
            admissible_cols = round(1/2*(cols+2*floor(4*M_horiz)+1-abs(cols-(2*floor(4*M_horiz)+1))));
            new_row = 1 + mod(row - first_row, length_corner_wedge);
            Xj(row,fliplr(admissible_cols)) = Xj(row,fliplr(admissible_cols)) + wrapped_data(new_row,end:-1:1);
                                % We use the following property: in an assignment
                                % A(B) = C where B and C are vectors, if
                                % some value x repeats in B, then the
                                % last occurrence of x is the one
                                % corresponding to the eventual assignment.
        end;

        Xj = rot90(Xj);

    end;    % for quadrant

    Xj = Xj .* lowpass;
    Xj_index1 = ((-floor(2*M1)):floor(2*M1)) + floor(4*M1) + 1;
    Xj_index2 = ((-floor(2*M2)):floor(2*M2)) + floor(4*M2) + 1;
    Xj(Xj_index1, Xj_index2) = Xj(Xj_index1, Xj_index2) .* hipass;

    loc_1 = Xj_topleft_1 + (0:(2*floor(4*M1)));
    loc_2 = Xj_topleft_2 + (0:(2*floor(4*M2)));
    X(loc_1,loc_2) = X(loc_1,loc_2) + Xj;

    % Preparing for loop reentry or exit

    Xj_topleft_1 = Xj_topleft_1 + floor(4*M1) - floor(2*M1);
    Xj_topleft_2 = Xj_topleft_2 + floor(4*M2) - floor(2*M2);

    lowpass = lowpass_next;

end;    % for j

if is_real
    Y = X;
    X = rot90(X,2);
    X = X + conj(Y);
end

% Coarsest wavelet level
M1 = M1/2;
M2 = M2/2;
Xj = fftshift(fft2(ifftshift(C{1}{1})))/sqrt(prod(size(C{1}{1})));
loc_1 = Xj_topleft_1 + (0:(2*floor(4*M1)));
loc_2 = Xj_topleft_2 + (0:(2*floor(4*M2)));
X(loc_1,loc_2) = X(loc_1,loc_2) + Xj .* lowpass;

% Finest level
M1 = N1/3;
M2 = N2/3;
if finest == 1,

    % Folding back onto N1-by-N2 matrix
    shift_1 = floor(2*M1)-floor(N1/2);
    shift_2 = floor(2*M2)-floor(N2/2);
    Y = X(:,(1:N2)+shift_2);
    Y(:,N2-shift_2+(1:shift_2)) = Y(:,N2-shift_2+(1:shift_2)) + X(:,1:shift_2);
    Y(:,1:shift_2) = Y(:,1:shift_2) + X(:,N2+shift_2+(1:shift_2));
    X = Y((1:N1)+shift_1,:);
    X(N1-shift_1+(1:shift_1),:) = X(N1-shift_1+(1:shift_1),:) + Y(1:shift_1,:);
    X(1:shift_1,:) = X(1:shift_1,:) + Y(N1+shift_1+(1:shift_1),:);

else

    % Extension to a N1-by-N2 matrix
    Y = fftshift(fft2(ifftshift(C{nbscales}{1})))/sqrt(prod(size(C{nbscales}{1})));
    X_topleft_1 = ceil((N1+1)/2) - floor(M1);
    X_topleft_2 = ceil((N2+1)/2) - floor(M2);
    loc_1 = X_topleft_1 + (0:(2*floor(M1)));
    loc_2 = X_topleft_2 + (0:(2*floor(M2)));
    Y(loc_1,loc_2) = Y(loc_1,loc_2) .* hipass_finest + X;
    X = Y;

end;

x = fftshift(ifft2(ifftshift(X)))*sqrt(prod(size(X)));
if is_real, x = real(x); end;

function [wl,wr] = fdct_wrapping_window(x)

% fdct_wrapping_window.m - Creates the two halves of a C^inf compactly supported window
%
% Inputs
%   x       vector or matrix of abscissae, the relevant ones from 0 to 1
%
% Outputs
%   wl,wr   vector or matrix containing samples of the left, resp. right
%           half of the window
%
% Used at least in fdct_wrapping.m and ifdct_wrapping.m
%
% By Laurent Demanet, 2004

wr = zeros(size(x));
wl = zeros(size(x));
x(abs(x) < 2^-52) = 0;
wr((x > 0) & (x < 1)) = exp(1-1./(1-exp(1-1./x((x > 0) & (x < 1)))));
wr(x <= 0) = 1;
wl((x > 0) & (x < 1)) = exp(1-1./(1-exp(1-1./(1-x((x > 0) & (x < 1))))));
wl(x >= 1) = 1;
normalization = sqrt(wl.^2 + wr.^2);
wr = wr ./ normalization;
wl = wl ./ normalization;

author: gpeyre

repo: matlab-toolboxes


text:
function gT = simKernDiagGradX(kern, t);

% SIMKERNDIAGGRADX Gradient of SIM kernel's diagonal with respect to the
% input times t.
% FORMAT
% DESC computes the gradient of the diagonal of the single input motif
% kernel matrix with respect to the elements of the design matrix given
% in t.
% ARG kern : the kernel structure for which gradients are being computed.
% ARG t : the input data in the form of a design matrix.
% RETURN gT : the gradients of the diagonal with respect to each element
% of t. The returned matrix has the same dimensions as t.
%
% SEEALSO : simKernParamInit, kernDiagGradX, simkernGradX
%
% COPYRIGHT : Neil D. Lawrence, 2006
%
% COPYRIGHT : David Luengo, 2009

% KERN

if size(t, 2) > 1
  error('Input can only have one column');
end

sigma = sqrt(2/kern.inverseWidth);
t = t - kern.delay;
halfSigmaD = 0.5*sigma*kern.decay;

gT = zeros(size(t));

if (kern.isStationary == false)
    lnPart1 = lnDiffErfs(halfSigmaD, halfSigmaD-t/sigma);
    gT = kern.variance * exp(halfSigmaD*halfSigmaD - 2*kern.decay*t + lnPart1);
else
    gT = zeros(size(t));
end

if ~isfield(kern, 'isNormalised') || (kern.isNormalised == false)
    gT = gT * sigma * sqrt(pi);
end

if isfield(kern, 'gaussianInitial') && kern.gaussianInitial,
  error('simKerDiagGradX not implemented for gaussianInitial')
end


text:
function test_findcluster

% MEM 8gb
% WALLTIME 00:20:00
% DEPENDENCY findcluster

[dum, ftpath] = ft_version;
cd(fullfile(ftpath, 'private'));

ft_hastoolbox('spm12', 1);

% create a data matrix and a spatial neighbourhoodmatrix
nchan = 5;
C     = diag(ones(1,4),-1)+diag(ones(1,4),1)>0;
dat   = zeros(5,4,3);
dat(1:2,1:2,1:2) = 1;
dat(3,2:3,2:3) = 1;
dat(4,1,1:3) = 1;
dat(5,2:4,3) = 1;
dat = dat>0;

% input is space x something x something
[clus, num] = findcluster(dat,C);
assert(num==3);
assert(sum(clus(:)==1)==12);
assert(sum(clus(:)==2)==3);
assert(sum(clus(:)==3)==3);

% channel 2 is disconnected from the rest
C2 = C;
C2(1,1:3) = [0 0 1];
C2(3,1:3) = [1 0 0];
C2(2,1:3) = 0;
[clus, num] = findcluster(dat,C2);
assert(num==4);
assert(sum(clus(:)==1)==8);
assert(sum(clus(:)==2)==4);
assert(sum(clus(:)==3)==3);
assert(sum(clus(:)==4)==3);

% input is as if it's a spatial 3D reshapeable volume
[clus, num] = findcluster(shiftdim(dat,-1),false);
assert(num==3);
assert(sum(clus(:)==1)==12);
assert(sum(clus(:)==2)==3);
assert(sum(clus(:)==3)==3);

dat4d = repmat(dat, [1 1 1 4]);
dat5d = repmat(dat, [1 1 1 4 5]);
dat6d = repmat(dat, [1 1 1 4 5 6]);
dat7d = repmat(dat, [1 1 1 4 5 6 7]);

[clus4d, num4d] = findcluster(dat4d,C);
assert(num4d==3);
assert(sum(clus4d(:)==1)==12*4);
assert(sum(clus4d(:)==2)==3*4);
assert(sum(clus4d(:)==3)==3*4);

[clus5d, num5d] = findcluster(dat5d,C);
assert(num5d==3);
assert(sum(clus5d(:)==1)==12*4*5);
assert(sum(clus5d(:)==2)==3*4*5);
assert(sum(clus5d(:)==3)==3*4*5);

[clus6d, num6d] = findcluster(dat6d,C);
assert(num6d==3);
assert(sum(clus6d(:)==1)==12*4*5*6);
assert(sum(clus6d(:)==2)==3*4*5*6);
assert(sum(clus6d(:)==3)==3*4*5*6);

[clus7d, num7d] = findcluster(dat7d,C);
assert(num7d==3);
assert(sum(clus7d(:)==1)==12*4*5*6*7);
assert(sum(clus7d(:)==2)==3*4*5*6*7);
assert(sum(clus7d(:)==3)==3*4*5*6*7);

[clus4d, num4d] = findcluster(shiftdim(dat4d,-1), false);
assert(num4d==3);
assert(sum(clus4d(:)==1)==12*4);
assert(sum(clus4d(:)==2)==3*4);
assert(sum(clus4d(:)==3)==3*4);

[clus5d, num5d] = findcluster(shiftdim(dat5d,-1), false);
assert(num5d==3);
assert(sum(clus5d(:)==1)==12*4*5);
assert(sum(clus5d(:)==2)==3*4*5);
assert(sum(clus5d(:)==3)==3*4*5);

[clus6d, num6d] = findcluster(shiftdim(dat6d,-1), false);
assert(num6d==3);
assert(sum(clus6d(:)==1)==12*4*5*6);
assert(sum(clus6d(:)==2)==3*4*5*6);
assert(sum(clus6d(:)==3)==3*4*5*6);

[clus7d, num7d] = findcluster(shiftdim(dat7d,-1), false);
assert(num7d==3);
assert(sum(clus7d(:)==1)==12*4*5*6*7);
assert(sum(clus7d(:)==2)==3*4*5*6*7);
assert(sum(clus7d(:)==3)==3*4*5*6*7);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function result = ball_unit_f3_nd ( func, n )

%*****************************************************************************80
%
%% BALL_UNIT_F3_ND approximates an integral inside the unit ball in ND.
%
%  Integration region:
%
%    Points X(1:N) such that:
%
%      Sum ( X(1:N)**2 ) <= 1.
%
%  Discussion:
%
%    A 2**(N+1)-1 point 5-th degree formula is used, Stroud number SN:5-4.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 May 2004
%
%  Reference:
%
%    Arthur H Stroud,
%    Approximate Calculation of Multiple Integrals,
%    Prentice Hall, 1971.
%
%  Parameters:
%
%    Input, external FUNC, the name of the user supplied
%    function which evaluates F at the N-vector X, of the form
%      function value = func ( n, x )
%
%    Input, integer N, the dimension of the space.
%
%    Output, real RESULT, the approximate integral of the function.
%
  quad = 0.0E+00;
%
%  The first point is the center of the ball.
%
  x(1:n) = 0.0E+00;

  weight = 4.0E+00 / ( n + 2 )^2;
  quad = quad + weight * feval ( func, n, x );

  s = 1.0E+00 / sqrt ( ( n + 4 ) );

  for i = 1 : n

    ri = sqrt ( ( i + 2 ) / ( n + 4 ) );
%
%  Set up the first point, with (I-1) zeroes, RI, and then N-I S's.
%
    for j = 1 : n

      if ( j < i )
        x(j) = 0.0E+00;
      elseif ( j == i )
        x(j) = ri;
      else
        x(j) = s;
      end

    end

    weight = 2.0E+00^( i - n ) * ( n + 4 ) / ( ( i + 1 ) * ( i + 2 ) * ( n + 2 ) );
%
%  Now go through all sign permutations of the basic point.
%
    for j = 1 : 2^(n+1-i)

      jtemp = j - 1;

      for k = i : n

        if ( mod ( jtemp, 2 ) == 1 )
          x(k) = - abs ( x(k) );
        else
          x(k) = abs ( x(k) );
        end

        jtemp = floor ( jtemp / 2 );

      end

      quad = quad + weight * feval ( func, n, x );

    end

  end

  volume = ball_unit_volume_nd ( n );
  result = quad * volume;

  return
end


text:
function planC = createStapleStruct(structNumV,probCutoff,structureName,planC)
% function planC = createStapleStruct(structNumV,probCutoff,structureName,planC)
%
% Function to create STAPLE agreement structure at the passed probability
% probCutoff.
%
% APA, 11/13/2017

if ~exist('planC','var')
    global planC
end
indexS = planC{end};

scanNum = getStructureAssociatedScan(structNumV(1),planC);
numObs = length(structNumV);
siz = getUniformScanSize(planC{indexS.scan}(scanNum));
rateMatM = false(prod(siz),numObs);
for i=1:numObs
    mask3M = getUniformStr(structNumV(i),planC);
    rateMatM(:,i) = mask3M(:);
end
indV = sum(rateMatM,2) > 0;
rateMatM = rateMatM(indV,:);

iterlim=100;
senstart=0.9999*ones(1,numObs);
specstart=0.9999*ones(1,numObs);
[stapleV, sen, spec, Sall] = staple(rateMatM,iterlim, single(senstart), single(specstart));

staple3M = zeros(siz);
staple3M(indV) = stapleV;
newStrMask3M = staple3M >= probCutoff;

planC = maskToCERRStructure(newStrMask3M, 1, scanNum, structureName, planC);


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Spectral-spatial RF pulse design script in Matlab, based on
%"Spectral-spatial RF pulse design for through-plane phase precompensatory
%slice selection for T2*-weighted functional MRI", Chun-yu Yip et al,
%Magnetic Resonance in Medicine, 2009.
%
%Written by Chun-yu Yip, University of Michigan, Ann Arbor, 4/1/2009
%Current affiliation: Athinoula A. Martinos Center for Biomedical Imaging,
%Massachusetts General Hospital, Harvard Medical School, Charlestown, MA,
%USA.
%Email address: chunyuy@nmr.mgh.harvard.edu
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%Running this script requires that the image reconstruction toolbox be
%installed first: http://www.eecs.umich.edu/~fessler/irt/fessler.tgz. In
%your Matlab workspace, please "addpath" the folders in the Fessler
%toolbox. For example:
%
%addpath('fesslertoolboxlocation/utilities/');
%addpath('fesslertoolboxlocation/systems/');
%addpath('fesslertoolboxlocation/nufft/');
%addpath('fesslertoolboxlocation/wls/');
%addpath('fesslertoolboxlocation/mex/');
%addpath('fesslertoolboxlocation/mex/v7/');
%etc.
%
%You can add those lines in your startup.m file so that you do not have to
%add the paths manually everytime. You can type "path" at the matlab prompt
%to check that the paths are successfully added.
%
%
%To design pulse, you have to first load pulse design parameter values. To
%start, you can adopt those we used in our paper in the
%example_parameter_files folder. To keep track of parameters, we recommend
%that you create one folder for each pulse design occasion. For example,
%you can create a "siemens_parameter_files" folder and a
%"ge_parameter_files", or "parameter_files_1Jan2010".
%

[status,homepath] = system('pwd');
homepath = strtrim(homepath);
cd([homepath '/example_parameter_files/']);

kp = kparametersSPSP;           %k-space trajectory parameters
rfp = rfparametersSPSP;         %Pulse design parameters
iop = ioparametersSPSP;         %Input-output parameters

cd(homepath);
%
%kp, rfp, iop are matlab structures, whose fields can be accessed by, e.g.,
%"kp.pointtime". Please see parameter files for details of each field.

%Design z-gradient waveform, based on parameters in kp.
[kp,gz,kz,kf] = compute_gz_spsp(kp);

%Design complex-valued RF waveform iteratively using conjugate gradient
[b] = compute_rf_spsp_mgh(kp,rfp,gz,kz,kf);

%Write computed waveforms to files for simulation and/or scanner.
write2files_spsp(kp,rfp,iop,gz,b);

%Perform Bloch simulation in SPSP space
[mresult] = dosim7_spsp(kp,rfp,iop);

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
%Author: Moeti Ncube
% This code uses the Schwartz-Smith model to calibrate and simulate multiple assets such that:
%
% 1. The Calibration and simulation is consistent with the observable forward curve, adjusted for seasonality, at any given date
% 2. The Calibration and simulation is consistent with the observable ATM volatility at a given date
% 3. The Calibration and simulation is consistent with the observable correlation structure between the forward curve vectors at
% a given date
%
%
% In this example, I use four assets: 5x16,2x16,7x8 PJM forward prices and ATM volatilities along with natural gas forward prices
% and ATM volatilities
%
% Calibration of the parameters is done in Excel. By inputing the vectors of Forward and ATM volatilites for each commodity, I
% can compute the theoretical Schwartz Smith Forward Prices and Standard Deviations for a given maturity. I then use Excel solver
% to minimize the difference between the observed market values and their theoretical values to obtain the Schwartz-Smith
% parameters for each commodity. Note this methodology is drastically different from the one described in "Short Term Variation
% and Long-Term Dynamics in Commodity Prices" in which Schwartz and Smith calibrate their model to historical futures prices.
% Here I calibrate the model to the current forward and volatility curve as well as adjust for seasonality. This procedure is
% much more practical pricing methodology.
%
% The more difficult step was insuring that the correlation between and  asset(i) and asset(j) at maturity (t) was consistent
% with the implied correlation between the forward price vectors. This was done by adjusting the correlation between the
% short-term factors of commodities at each maturity. The matlab code factors in this adjustment and simulates the 4 commodities
% in this example to show that the theoretical prices, volatilities, and correlations match up with the observed market data.
%
% There is not, to my knowledge, a commodities methodology that incoporates so many market factors across multible commodities
% into one simulation. The advantages of such a model allows for more accurate modeling of spark spreads and pricing of deals
% that are dependent on multiple commodities prices. I have included all files, including excel, associated with this calibration
% and simulation.

clear all; close all

filename = 'SchwartzSmithCalibration.xls';

sim=1000;
str={'Fit 5x16';'Fit 2x16';'Fit 7x8';'Fit Fuel'};
iter=length(str);

for i=1:iter
data{i} = xlsread(filename, str{i}, 'B3:M86');
fwd(:,i)=data{i}(:,1);

kappa(i)=data{i}(1,end);
sigmax(i)=data{i}(3,end);
sigmae(i)=data{i}(4,end);
pxe(i)=data{i}(7,end);
T(:,i)=data{i}(:,3);
x0(i)=data{i}(8,end);
e0(i)=data{i}(9,end);
for j=1:length(T(:,i))
ve(j,i)=sigmae(i)^2*T(j,i);
vx(j,i)=(1-exp(-2*kappa(i)*T(j,i)))*(.5*sigmax(i)^2)/kappa(i);
covxe(j,i)=(1-exp(-kappa(i)*T(j,i)))*pxe(i)*sigmax(i)*sigmae(i)/kappa(i);
end
corxet(:,i)=covxe(:,i)./sqrt(ve(:,i).*vx(:,i));
end

%Compute Correlation matrix from observed forward curves
cmatrix=corrcoef(log(fwd));

%Find correlation structures needed to keep observed correlation structure
%during simulation
rho(:,1:1)=ones(length(T(:,1)),1);
for i=2:iter
for j=1:length(T(:,i))
atop1(j,1)=(1-exp(-(kappa(1)+kappa(i))*T(j,i)))*sigmax(1)*sigmax(i)/(kappa(1)+kappa(i));
atop2(j,1)=(1-exp(-kappa(1)*T(j,i)))*pxe(i)*sigmax(1)*sigmae(i)/kappa(1);
atop3(j,1)=(1-exp(-kappa(i)*T(j,i)))*pxe(1)*sigmax(i)*sigmae(1)/kappa(i);
atop4(j,1)=pxe(1)*pxe(i)*sigmae(1)*sigmae(i)*T(j,i);
abot1(j,1)=sqrt(vx(j,1)+ve(j,1)+2*covxe(j,1));
abot2(j,1)=sqrt(vx(j,i)+ve(j,i)+2*covxe(j,i));
acorxy(j,i)=(atop1(j,1)+atop2(j,1)+atop3(j,1)+atop4(j,1))/(abot1(j,1)*abot2(j,1));
end
rho(:,i)=cmatrix(1,i)./acorxy(:,i);
rho=min(rho,1);
end

iter2=1;
rmatrix0=randn(sim,length(data{1}(:,1)));
for j=1:iter

for k=1:length(data{1}(:,1))
    rmatrix{j}(:,k)=rho(k,j)*rmatrix0(:,k)+sqrt(1-rho(k,j)^2)*randn(sim,1);
end

[tspath,tmpath,tstdpath,tefwd,testdfwd,xpath1,epath1,r1,r2]=schwartzsmithsim(data{j}(22,end),data{j}(1:9,end),data{j}(10:21,end),data{j}(:,3),rmatrix{j},sim);
spath{j}=tspath;
mpath{j}=tmpath;
stdpath{j}=tstdpath;
efwd{j}=tefwd;
estdfwd{j}=testdfwd;
epath{j}=epath1;
xpath{j}=xpath1;
lnpath{j}=epath1+xpath1;
r1path{j}=r1;
r2path{j}=r2;

subplot(length(str),2,iter2)
title('Market Fwd (blue) vs Sim Fwd (red)')
hold on
plot(fwd(:,j))
plot(mpath{j},'r')

subplot(length(str),2,iter2+1)
title('Market Vol (blue) vs Sim Vol (red)')
hold on
plot(estdfwd{j})
plot(stdpath{j},'r')
iter2=iter2+2;
end

%Estimate empirical correlation matrix
for j=2:iter
for i=1:length(data{1}(:,1))
c(i)=corr(log(spath{1}(:,i)),log(spath{j}(:,i)));
end
empc(j)=mean(c);
end

MarketCorrelations=cmatrix(1,2:end)
SimCorrelations=empc(2:end)