text:
function M = load_image(type, n, options)

% load_image - load benchmark images.
%
%   M = load_image(name, n, options);
%
%   name can be:
%   Synthetic images:
%       'chessboard1', 'chessboard', 'square', 'squareregular', 'disk', 'diskregular', 'quaterdisk', '3contours', 'line',
%       'line_vertical', 'line_horizontal', 'line_diagonal', 'line_circle',
%       'parabola', 'sin', 'phantom', 'circ_oscil',
%       'fnoise' (1/f^alpha noise).
%   Natural images:
%       'boat', 'lena', 'goldhill', 'mandrill', 'maurice', 'polygons_blurred', or your own.
%
%   Copyright (c) 2004 Gabriel Peyre

if nargin<2
    n = 512;
end
options.null = 0;

if iscell(type)
    for i=1:length(type)
        M{i} = load_image(type{i},n,options);
    end
    return;
end

type = lower(type);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% parameters for geometric objects
eta         = getoptions(options, 'eta', .1);
gamma       = getoptions(options, 'gamma', 1/sqrt(2));
radius      = getoptions(options, 'radius', 10);
center      = getoptions(options, 'center', [0 0]);
center1     = getoptions(options, 'center1', [0 0]);
w           = getoptions(options, 'tube_width', 0.06);
nb_points   = getoptions(options, 'nb_points', 9);
scaling     = getoptions(options, 'scaling', 1);
theta       = getoptions(options, 'theta', 30 * 2*pi/360);
eccentricity = getoptions(options, 'eccentricity', 1.3);
sigma = getoptions(options, 'sigma', 0);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% for the line, can be vertical / horizontal / diagonal / any
if strcmp(type, 'line_vertical')
    eta = 0.5;              % translation
    gamma = 0;      % slope
elseif strcmp(type, 'line_horizontal')
    eta = 0.5;              % translation
    gamma = Inf;      % slope
elseif strcmp(type, 'line_diagonal')
    eta = 0;              % translation
    gamma = 1;      % slope
end

if strcmp(type(1:min(12,end)), 'square-tube-')
    k = str2double(type(13:end));
    c1 = [.22 .5]; c2 = [1-c1(1) .5];
    eta = 1.5;
    r1 = [c1 c1] + .21*[-1 -eta 1 eta];
    r2 = [c2 c2] + .21*[-1 -eta 1 eta];
    M = double( draw_rectangle(r1,n) | draw_rectangle(r2,n) );
    if mod(k,2)==0
        sel = n/2-k/2+1:n/2+k/2;
    else
        sel = n/2-(k-1)/2:n/2+(k-1)/2;
    end
    M( round(.25*n:.75*n), sel ) = 1;
    return;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
switch lower(type)

    case 'constant'
        M = ones(n);

    case 'ramp'
        x = linspace(0,1,n);
        [Y,M] = meshgrid(x,x);

    case 'bump'

        s = getoptions(options, 'bump_size', .5);
        c = getoptions(options, 'center', [0 0]);
        if length(s)==1
            s = [s s];
        end
        x = linspace(-1,1,n);
        [Y,X] = meshgrid(x,x);
        X = (X-c(1))/s(1); Y = (Y-c(2))/s(2);
        M = exp( -(X.^2+Y.^2)/2 );

    case 'periodic'
        x = linspace(-pi,pi,n)/1.1;
        [Y,X] = meshgrid(x,x);
        f = getoptions(options, 'freq', 6);
        M = (1+cos(f*X)).*(1+cos(f*Y));

    case {'letter-x' 'letter-v' 'letter-z' 'letter-y'}
        M = create_letter(type(8), radius, n);

    case 'l'
        r1 = [.1 .1  .3 .9];
        r2 = [.1 .1 .9 .4];
        M = double( draw_rectangle(r1,n) | draw_rectangle(r2,n) );

    case 'ellipse'
        c1 = [0.15 0.5];
        c2 = [0.85 0.5];
        x = linspace(0,1,n);
        [Y,X] = meshgrid(x,x);
        d = sqrt((X-c1(1)).^2 + (Y-c1(2)).^2) + sqrt((X-c2(1)).^2 + (Y-c2(2)).^2);
        M = double( d<=eccentricity*sqrt( sum((c1-c2).^2) ) );
    case 'ellipse-thin'
        options.eccentricity = 1.06;
        M = load_image('ellipse', n, options);
    case 'ellipse-fat'
        options.eccentricity = 1.3;
        M = load_image('ellipse', n, options);

    case 'square-tube'
        c1 = [.25 .5];
        c2 = [.75 .5];
        r1 = [c1 c1] + .18*[-1 -1 1 1];
        r2 = [c2 c2] + .18*[-1 -1 1 1];
        r3 = [c1(1)-w c1(2)-w c2(1)+w c2(2)+w];
        M = double( draw_rectangle(r1,n) | draw_rectangle(r2,n) | draw_rectangle(r3,n) );
    case 'square-tube-1'
        options.tube_width = 0.03;
        M = load_image('square-tube', n, options);
    case 'square-tube-2'
        options.tube_width = 0.06;
        M = load_image('square-tube', n, options);
    case 'square-tube-3'
        options.im = 0.09;
        M = load_image('square-tube', n, options);
    case 'polygon'
        theta = sort( rand(nb_points,1)*2*pi );
        radius = scaling*rescale(rand(nb_points,1), 0.1, 0.93);
        points = [cos(theta) sin(theta)] .* repmat(radius, 1,2);
        points = (points+1)/2*(n-1)+1; points(end+1,:) = points(1,:);
        M = draw_polygons(zeros(n),0.8,{points'});
        [x,y] = ind2sub(size(M),find(M));
        p = 100; m = length(x);
        lambda = linspace(0,1,p);
        X = n/2 + repmat(x-n/2, [1 p]) .* repmat(lambda, [m 1]);
        Y = n/2 + repmat(y-n/2, [1 p]) .* repmat(lambda, [m 1]);
        I = round(X) + (round(Y)-1)*n;
        M = zeros(n); M(I) = 1;
    case 'polygon-8'
        options.nb_points = 8;
        M = load_image('polygon', n, options);
    case 'polygon-10'
        options.nb_points = 8;
        M = load_image('polygon', n, options);
    case 'polygon-12'
        options.nb_points = 8;
        M = load_image('polygon', n, options);
    case 'pacman'
        options.radius = 0.45;
        options.center = [.5 .5];
        M = load_image('disk', n, options);
        x = linspace(-1,1,n);
        [Y,X] = meshgrid(x,x);
        T =atan2(Y,X);
        M = M .* (1-(abs(T-pi/2)<theta/2));
    case 'square-hole'
        options.radius = 0.45;
        M = load_image('disk', n, options);
        options.scaling = 0.5;
        M = M - load_image('polygon-10', n, options);

    case 'grid-circles'
        if isempty(n)
            n = 256;
        end
        f = getoptions(options, 'frequency', 30);
        eta = getoptions(options, 'width', .3);
        x = linspace(-n/2,n/2,n) - round(n*0.03);
        y = linspace(0,n,n);
        [Y,X] = meshgrid(y,x);
        R = sqrt(X.^2+Y.^2);
        theta = 0.05*pi/2;
        X1 = cos(theta)*X+sin(theta)*Y;
        Y1 = -sin(theta)*X+cos(theta)*Y;
        A1 = abs(cos(2*pi*R/f))<eta;
        A2 = max( abs(cos(2*pi*X1/f))<eta, abs(cos(2*pi*Y1/f))<eta );

        M = A1;
        M(X1>0) = A2(X1>0);

    case 'chessboard1'
        x = -1:2/(n-1):1;
        [Y,X] = meshgrid(x,x);
        M = (2*(Y>=0)-1).*(2*(X>=0)-1);

    case 'chessboard'
        width = getoptions(options, 'width', round(n/16) );
        [Y,X] = meshgrid(0:n-1,0:n-1);
        M = mod( floor(X/width)+floor(Y/width), 2 ) == 0;

    case 'square'
        if ~isfield( options, 'radius' )
            radius = 0.6;
        end
        x = linspace(-1,1,n);
        [Y,X] = meshgrid(x,x);
        M = max( abs(X),abs(Y) )<radius;

    case 'squareregular'
        M = rescale(load_image('square',n,options));
        if not(isfield(options, 'alpha'))
            options.alpha = 3;
        end
        S = load_image('fnoise',n,options);
        M = M + rescale(S,-0.3,0.3);

    case 'regular1'
        options.alpha = 1;
        M = load_image('fnoise',n,options);
    case 'regular2'
        options.alpha = 2;
        M = load_image('fnoise',n,options);
    case 'regular3'
        options.alpha = 3;
        M = load_image('fnoise',n,options);

    case 'sparsecurves'
        options.alpha = 3;
        M = load_image('fnoise',n,options);
        M = rescale(M);
        ncurves = 3;
        M = cos(2*pi*ncurves);

    case 'geometrical'

        J = getoptions(options, 'Jgeometrical', 4);
        sgeom = 100*n/256;
        options.bound = 'per';
        A = ones(n);
        for j=0:J-1
            B = A;
            for k=1:2^j
                I = find(B==k);
                U = perform_blurring(randn(n),sgeom,options);
                s = median(U(I));
                I1 = find( (B==k) & (U>s) );
                I2 = find( (B==k) & (U<=s) );
                A(I1) = 2*k-1;
                A(I2) = 2*k;
            end
        end
        M = A;

    case 'lic-texture'

        disp('Computing random tensor field.');
        options.sigma_tensor = getoptions(options, 'lic_regularity', 50*n/256);
        T = compute_tensor_field_random(n,options);
        Flow = perform_tensor_decomp(T); % extract eigenfield.
        options.isoriented = 0; % no orientation in streamlines
        % initial texture
        lic_width = getoptions(options, 'lic_width', 0);
        M0 = perform_blurring(randn(n),lic_width);
        M0 = perform_histogram_equalization( M0, 'linear');
        options.histogram = 'linear';
        options.dt = 0.4;
        options.M0 = M0;
        options.verb = 1;
        options.flow_correction = 1;
        options.niter_lic = 3;
        w = 30;
        M = perform_lic(Flow, w, options);

    case 'square_texture'
        M = load_image('square',n);
        M = rescale(M);
        % make a texture patch
        x = linspace(0,1,n);
        [Y,X] = meshgrid(x,x);
        theta = pi/3;
        x = cos(theta)*X + sin(theta)*Y;
        c = [0.3,0.4]; r = 0.2;
        I = find( (X-c(1)).^2 + (Y-c(2)).^2 < r^2 );
        eta = 3/n; lambda = 0.3;
        M(I) = M(I) + lambda * sin( x(I) * 2*pi / eta );

    case 'tv-image'
        M = rand(n);
        tau = compute_total_variation(M);
        options.niter = 400;
        [M,err_tv,err_l2] = perform_tv_projection(M,tau/1000,options);
        M = perform_histogram_equalization(M,'linear');

    case 'oscillatory_texture'
        x = linspace(0,1,n);
        [Y,X] = meshgrid(x,x);
        theta = pi/3;
        x = cos(theta)*X + sin(theta)*Y;
        c = [0.3,0.4]; r = 0.2;
        I = find( (X-c(1)).^2 + (Y-c(2)).^2 < r^2 );
        eta = 3/n; lambda = 0.3;
        M = sin( x * 2*pi / eta );

    case {'line', 'line_vertical', 'line_horizontal', 'line_diagonal'}
        x = 0:1/(n-1):1;
        [Y,X] = meshgrid(x,x);
        if gamma~=Inf
            M = (X-eta) - gamma*Y < 0;
        else
            M = (Y-eta) < 0;
        end

    case 'line-windowed'
        x = 0:1/(n-1):1;
        [Y,X] = meshgrid(x,x);
        eta = .3;
        gamma = getoptions(options, 'gamma', pi/10);
        parabola = getoptions(options, 'parabola', 0);
        M = (X-eta) - gamma*Y - parabola*Y.^2 < 0;
        f = sin( pi*x ).^2;
        M = M .* ( f'*f );

    case 'grating'
        x = linspace(-1,1,n);
        [Y,X] = meshgrid(x,x);
        theta = getoptions(options, 'theta', .2);
        freq = getoptions(options, 'freq', .2);
        X = cos(theta)*X + sin(theta)*Y;
        M = sin(2*pi*X/freq);

    case 'disk'
        if ~isfield( options, 'radius' )
            radius = 0.35;
        end
        if ~isfield( options, 'center' )
            center = [0.5, 0.5];    % center of the circle
        end
        x = 0:1/(n-1):1;
        [Y,X] = meshgrid(x,x);
        M = (X-center(1)).^2 + (Y-center(2)).^2 < radius^2;

    case 'diskregular'
        M = rescale(load_image('disk',n,options));
        if not(isfield(options, 'alpha'))
            options.alpha = 3;
        end
        S = load_image('fnoise',n,options);
        M = M + rescale(S,-0.3,0.3);

    case 'quarterdisk'
        if ~isfield( options, 'radius' )
            radius = 0.95;
        end
        if ~isfield( options, 'center' )
            center = -[0.1, 0.1];    % center of the circle
        end
        x = 0:1/(n-1):1;
        [Y,X] = meshgrid(x,x);
        M = (X-center(1)).^2 + (Y-center(2)).^2 < radius^2;

    case 'fading_contour'
        if ~isfield( options, 'radius' )
            radius = 0.95;
        end
        if ~isfield( options, 'center' )
            center = -[0.1, 0.1];    % center of the circle
        end
        x = 0:1/(n-1):1;
        [Y,X] = meshgrid(x,x);
        M = (X-center(1)).^2 + (Y-center(2)).^2 < radius^2;
        theta = 2/pi*atan2(Y,X);
        h = 0.5;
        M = exp(-(1-theta).^2/h^2).*M;

    case '3contours'
        radius = 1.3;
        center = [-1, 1];
        radius1 = 0.8;
        center1 = [0, 0];
        x = 0:1/(n-1):1;
        [Y,X] = meshgrid(x,x);
        f1 = (X-center(1)).^2 + (Y-center(2)).^2 < radius^2;
        f2 = (X-center1(1)).^2 + (Y-center1(2)).^2 < radius1^2;
        M = f1 + 0.5*f2.*(1-f1);

    case 'line_circle'

        gamma = 1/sqrt(2);

        x = linspace(-1,1,n);
        [Y,X] = meshgrid(x,x);
        M1 = double( X>gamma*Y+0.25 );
        M2 = X.^2 + Y.^2 < 0.6^2;
        M = 20 + max(0.5*M1,M2) * 216;

    case 'fnoise'

        % generate an image M whose Fourier spectrum amplitude is
        %   |M^(omega)| = 1/f^{omega}
        alpha = getoptions(options, 'alpha', 1);
        M = gen_noisy_image(n,alpha);

    case 'gaussiannoise'
        % generate an image of filtered noise with gaussian
        sigma = getoptions(options, 'sigma', 10);
        M = randn(n);
        m = 51;
        h = compute_gaussian_filter([m m],sigma/(4*n),[n n]);
        M = perform_convolution(M,h);
        return;

    case {'bwhorizontal','bwvertical','bwcircle'}

        [Y,X] = meshgrid(0:n-1,0:n-1);
        if strcmp(type, 'bwhorizontal')
            d = X;
        elseif strcmp(type, 'bwvertical')
            d = Y;
        elseif strcmp(type, 'bwcircle')
            d = sqrt( (X-(n-1)/2).^2 + (Y-(n-1)/2).^2 );
        end
        if isfield(options, 'stripe_width')
            stripe_width = options.stripe_width;
        else
            stripe_width = 5;
        end
        if isfield(options, 'black_prop')
            black_prop = options.black_prop;
        else
            black_prop = 0.5;
        end
        M = double( mod( d/(2*stripe_width),1 )>=black_prop );

    case 'parabola'

        % curvature
        c = getoptions(c, 'c', .1);
        % angle
        theta = getoptions(options, 'theta',  pi/sqrt(2));
        x = -0.5:1/(n-1):0.5;
        [Y,X] = meshgrid(x,x);
        Xs = X*cos(theta) + Y*sin(theta);
        Y =-X*sin(theta) + Y*cos(theta); X = Xs;
        M = Y>c*X.^2;

    case 'sin'

        [Y,X] = meshgrid(-1:2/(n-1):1, -1:2/(n-1):1);
        M = Y >= 0.6*cos(pi*X);
        M = double(M);

    case 'circ_oscil'

        x = linspace(-1,1,n);
        [Y,X] = meshgrid(x,x);
        R = sqrt(X.^2+Y.^2);
        M = cos(R.^3*50);

    case 'phantom'

        M = phantom(n);

    case 'periodic_bumps'
        nbr_periods = getoptions(options, 'nbr_periods', 8);
        theta = getoptions(options, 'theta', 1/sqrt(2));
        skew = getoptions(options, 'skew', 1/sqrt(2) );
        A = [cos(theta), -sin(theta); sin(theta), cos(theta)];
        B = [1 skew; 0 1];
        T = B*A;
        x = (0:n-1)*2*pi*nbr_periods/(n-1);
        [Y,X] = meshgrid(x,x);
        pos = [X(:)'; Y(:)'];
        pos = T*pos;
        X = reshape(pos(1,:), n,n);
        Y = reshape(pos(2,:), n,n);
        M = cos(X).*sin(Y);

    case 'noise'
        sigma = getoptions(options, 'sigma', 1);
        M = randn(n) * sigma;

    case 'disk-corner'
        x = linspace(0,1,n);
        [Y,X] = meshgrid(x,x);
        rho = .3; eta = .1;
        M1 = rho*X+eta<Y;
        c = [0 .2]; r = .85;
        d = (X-c(1)).^2 + (Y-c(2)).^2;
        M2 = d<r^2;
        M = M1.*M2;

    otherwise
        ext = {'gif', 'png', 'jpg', 'bmp', 'tiff', 'tif', 'pgm', 'ppm'};
        for i=1:length(ext)
            name = [type '.' ext{i}];
            if( exist(name) )
                M = imread( name );
                M = double(M);
                if not(isempty(n)) && (n~=size(M, 1) || n~=size(M, 2)) && nargin>=2
                    M = image_resize(M,n,n);
                end
                if strcmp(type, 'peppers-bw')
                    M(:,1) = M(:,2);
                    M(1,:) = M(2,:);
                end
                if strcmp(type, 'cameraman')
                    M(1:4,:) = M(8:-1:5,:);
                end
                if sigma>0
                    M = perform_blurring(M,sigma);
                end
                return;
            end
        end
        error( ['Image ' type ' does not exists.'] );
end

M = double(M);

if sigma>0
    M = perform_blurring(M,sigma);
end

M = rescale(M);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function M = create_letter(a, r, n)

c = 0.2;
p1 = [c;c];
p2 = [c; 1-c];
p3 = [1-c; 1-c];
p4 = [1-c; c];
p4 = [1-c; c];
pc = [0.5;0.5];
pu = [0.5; c];

switch a
    case 'x'
        point_list = { [p1 p3] [p2 p4] };
    case 'z'
        point_list = { [p2 p3 p1 p4] };
    case 'v'
        point_list = { [p2 pu p3] };
    case 'y'
        point_list = { [p2 pc pu] [pc p3] };

end
% fit image
for i=1:length(point_list)
    a = point_list{i}(2:-1:1,:);
    a(1,:) = 1-a(1,:);
    point_list{i} = round( a*(n-1)+1 );
end
M = draw_polygons(zeros(n),r,point_list);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function sk = draw_polygons(mask,r,point_list)

sk = mask*0;
for i=1:length(point_list)
    pl = point_list{i};
    for k=2:length(pl)
        sk = draw_line(sk,pl(1,k-1),pl(2,k-1),pl(1,k),pl(2,k),r);
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function sk = draw_line(sk,x1,y1,x2,y2,r)

n = size(sk,1);
[Y,X] = meshgrid(1:n,1:n);
q = 100;
t = linspace(0,1,q);
x = x1*t+x2*(1-t); y = y1*t+y2*(1-t);
if r==0
    x = round( x ); y = round( y );
    sk( x+(y-1)*n ) = 1;
else
    for k=1:q
        I = find((X-x(k)).^2 + (Y-y(k)).^2 <= r^2 );
        sk(I) = 1;
    end
end

function M = gen_noisy_image(n,alpha)

% gen_noisy_image - generate a noisy cloud-like image.
%
%   M = gen_noisy_image(n,alpha);
%
% generate an image M whose Fourier spectrum amplitude is
%   |M^(omega)| = 1/f^{omega}
%
%   Copyright (c) 2004 Gabriel Peyr?

if nargin<1
    n = 128;
end
if nargin<2
    alpha = 1.5;
end

if mod(n(1),2)==0
    x = -n/2:n/2-1;
else
    x = -(n-1)/2:(n-1)/2;
end

[Y,X] = meshgrid(x,x);
d = sqrt(X.^2 + Y.^2) + 0.1;
f = rand(n)*2*pi;

M = (d.^(-alpha)) .* exp(f*1i);
% M = real(ifft2(fftshift(M)));

M = ifftshift(M);
M = real( ifft2(M) );

function y = gen_signal_2d(n,alpha)

% gen_signal_2d -  generate a 2D C^\alpha signal of length n x n.
%   gen_signal_2d(n,alpha) generate a 2D signal C^alpha.
%
%   The signal is scale in [0,1].
%
%   Copyright (c) 2003 Gabriel Peyr?

% new new method

[Y,X] = meshgrid(0:n-1, 0:n-1);

A = X+Y+1;
B = X-Y+n+1;

a = gen_signal(2*n+1, alpha);
b = gen_signal(2*n+1, alpha);
y = a(A).*b(B);
% M = a(1:n)*b(1:n)';

return;

% new method
h = (-n/2+1):(n/2); h(n/2)=1;
[X,Y] = meshgrid(h,h);
h = sqrt(X.^2+Y.^2+1).^(-alpha-1/2);
h = h .* exp( 2i*pi*rand(n,n) );
h = fftshift(h);
y = real( ifft2(h) );

m1 = min(min(y));
m2 = max(max(y));
y = (y-m1)/(m2-m1);

return;

%% old code

y = rand(n,n);
y = y - mean(mean(y));
for i=1:alpha
    y = cumsum(cumsum(y)')';
    y = y - mean(mean(y));
end
m1 = min(min(y));
m2 = max(max(y));
y = (y-m1)/(m2-m1);

function newimg = image_resize(img,p1,q1,r1)

% image_resize - resize an image using bicubic interpolation
%
%   newimg = image_resize(img,nx,ny,nz);
% or
%   newimg = image_resize(img,newsize);
%
%   Works for 2D, 2D 2 or 3 channels, 3D images.
%
%   Copyright (c) 2004 Gabriel Peyr?

if nargin==2
    % size specified as an array
    q1 = p1(2);
    if length(p1)>2
        r1 = p1(3);
    else
        r1 = size(img,3);
    end
    p1 = p1(1);
end

if nargin<4
    r1 = size(img,3);
end

if ndims(img)<2 || ndims(img)>3
    error('Works only for grayscale or color images');
end

if ndims(img)==3 && size(img,3)<4
    % RVB image
    newimg = zeros(p1,q1, size(img,3));
    for m=1:size(img,3)
        newimg(:,:,m) = image_resize(img(:,:,m), p1, q1);
    end
    return;
elseif ndims(img)==3
    p = size(img,1);
    q = size(img,2);
    r = size(img,3);
    [Y,X,Z] = meshgrid( (0:q-1)/(q-1), (0:p-1)/(p-1), (0:r-1)/(r-1)  );
    [YI,XI,ZI] = meshgrid( (0:q1-1)/(q1-1), (0:p1-1)/(p1-1), (0:r1-1)/(r1-1) );
    newimg = interp3( Y,X,Z, img, YI,XI,ZI ,'cubic');
    return;
end

p = size(img,1);
q = size(img,2);
[Y,X] = meshgrid( (0:q-1)/(q-1), (0:p-1)/(p-1) );
[YI,XI] = meshgrid( (0:q1-1)/(q1-1), (0:p1-1)/(p1-1) );
newimg = interp2( Y,X, img, YI,XI ,'cubic');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function M = draw_rectangle(r,n)

x = linspace(0,1,n);
[Y,X] = meshgrid(x,x);
M = double( (X>=r(1)) & (X<=r(3)) & (Y>=r(2)) & (Y<=r(4)) ) ;

author: gpeyre

repo: matlab-toolboxes


text:
%ISPARALLEL Test on parallel mapping
%
%  N = ISPARALLEL(W)
%      ISPARALLEL(W)
%
% INPUT
%  W    input mapping
%
% OUTPUT
%  N    logical value
%
% DESCRIPTION
% Returns true for parallel mappings. If no output is required,
% false outputs are turned into errors. This may be used for
% assertion.
%
% SEE ALSO (PRTools Guide)
% ISMAPPING, ISSTACKED

% $Id: isparallel.m,v 1.2 2006/03/08 22:06:58 duin Exp $

function n = isparallel(w)

	if isa(w,'prmapping') & strcmp(w.mapping_file,'parallel')
		n = 1;
	else
		n = 0;
	end

	% generate error if input is not a parallel mapping
	% AND no output is requested (assertion)

	if nargout == 0 & n == 0
		error([newline '---- Parallel mapping expected -----'])
	end

return

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function [ a, b ] = p01_lim ( )

%*****************************************************************************80
%
%% P01_LIM returns the integration limits for problem 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real A, B, the limits of integration.
%
  a = 0.0;
  b = 1.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function inp = regCorrMeanInt(inp);
% regCorrMeanInt - Detects and Corrects for oscillation in the mean intensity
%                  of the inplanes, often caused by having the INTERLEAVE
%                  button OFF when acquiring the anatomy inplanes.
%                  Detection: based on the energy of the 1/2 cycle/sample
%                             harmonic of the mean signal.
%                  Correction: adjusts the mean intensity of the odd
%                  inplanes to a 2nd order polynomial, and then scales
%                  the even inplanes to match the values predicted by the
%                  polynomial.
%
%    inp = regCorrMeanInt(inp);
%
% Oscar Nestares - 5/99
%

% mean intensity of the inplanes
m = squeeze(mean(mean(inp)));

% FFT of the mean signal (using the next closest power of two)
L = length(m);
N = 2^ceil(log(L)/log(2));
M = abs(fft(m, N));

% oscillation if last harmonic greater than the mean of the other harmonics
if mean(M(2:N/2)) < M(N/2+1)
   disp('Mean intensities of the inplanes are oscillating.')
   disp('Was INTERLEAVE button ON?')
   disp('Correcting oscillation and continuing...')

   % fiting a 2nd order polynomial to the odd samples
   x = [1:length(m)]';
   po = polyfit(x(1:2:L),m(1:2:L),2);
   mo = polyval(po,x);

   % re-scaling the even samples
   for k=2:2:L
      inp(:,:,k) = inp(:,:,k) * mo(k)/m(k);
   end
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function v = cvx_value( x )
v = x;

% Copyright 2005-2014 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function [ a, det ] = spodi ( a, lda, n, job )

%*****************************************************************************80
%
%% SPODI computes the determinant and inverse of a certain matrix.
%
%  Discussion:
%
%    The matrix is real symmetric positive definite.
%    SPODI uses the factors computed by SPOCO, SPOFA or SQRDC.
%
%    A division by zero will occur if the input factor contains
%    a zero on the diagonal and the inverse is requested.
%    It will not occur if the subroutines are called correctly
%    and if SPOCO or SPOFA has set INFO == 0.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 November 2006
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Dongarra, Moler, Bunch and Stewart,
%    LINPACK User's Guide,
%    SIAM, (Society for Industrial and Applied Mathematics),
%    3600 University City Science Center,
%    Philadelphia, PA, 19104-2688.
%    ISBN 0-89871-172-X
%
%  Parameters:
%
%    Input, real A(LDA,N), the output A from SPOCO or SPOFA, or the output
%    X from SQRDC.
%
%    Input, integer LDA, the leading dimension of the array A.
%
%    Input, integer N, the order of the matrix A.
%
%    Input, integer JOB, specifies the task.
%    11, both determinant and inverse.
%    01, inverse only.
%    10, determinant only.
%
%    Output, real A(LDA,N), if SPOCO or SPOFA was used to factor A then
%    SPODI produces the upper half of inverse(A).  If SQRDC was used to
%    decompose X then SPODI produces the upper half of inverse(X'*X)
%    where X' is the transpose.  Elements of A below the diagonal are
%    unchanged.  If the units digit of JOB is zero, A is unchanged.
%
%    Output, real DET(2), the determinant of A or of X'*X
%    if requested.
%      determinant = DET(1) * 10.0**DET(2)
%    with 1.0 <= DET(1) < 10.0 or DET(1) == 0.0.
%

%
%  Compute the determinant.
%
  if ( job / 10 ~= 0 )

    det(1) = 1.0;
    det(2) = 0.0;
    s = 10.0;

    for i = 1 : n

      det(1) = a(i,i) * a(i,i) * det(1);

      if ( det(1) == 0.0 )
        break
      end

      while ( det(1) < 1.0 )
        det(1) = s * det(1);
        det(2) = det(2) - 1.0;
      end

      while ( s <= det(1) )
        det(1) = det(1) / s;
        det(2) = det(2) + 1.0;
      end

    end

  end
%
%  Compute inverse(R).
%
  if ( mod ( job, 10 ) ~= 0 )

    for k = 1 : n

      a(k,k) = 1.0 / a(k,k);
      t = -a(k,k);
      a(1:k-1,k) = sscal ( k-1, t, a(1:k-1,k), 1 );

      for j = k+1 : n
        t = a(k,j);
        a(k,j) = 0.0;
        a(1:k,j) = saxpy ( k, t, a(1:k,k), 1, a(1:k,j), 1 );
      end

    end
%
%  Form inverse(R) * (inverse(R))'.
%
    for j = 1 : n
      for k = 1 : j-1
        t = a(k,j);
        a(1:k,k) = saxpy ( k, t, a(1:k,j), 1, a(1:k,k), 1 );
      end
      t = a(j,j);
      a(1:j,j) = sscal ( j, t, a(1:j,j), 1 );
    end

  end

  return
end


text:
function msh = meshDefault(host,id,mmPerVox,relaxIter,meshName,decimateReduction)
%
%    msh = meshDefault([host],[id],[mmPerVox],[relaxIter],[meshName])
%
%Author: Wandell/RFD
%Purpose:
%   Return a default mesh structure used by mrVista (and mrMesh).  The
%   various parameters were originally set by RFD.  The structures and
%   access routines were built by BW and are still in progress.
%

if ieNotDefined('host'),host = 'localhost'; end
if ieNotDefined('id'), id = 1; end
if ieNotDefined('mmPerVox'), mmPerVox = [1,1,1]; end
if ieNotDefined('relaxIter'), relaxIter = 0; end
if ieNotDefined('meshName'), meshName = 'untitled'; end
if ieNotDefined('decimateReduction'), decimateReduction = 0.1; end

% These can be passed in.
msh.name = meshName;
msh.fileName = '';
msh.relaxIterations = relaxIter;
msh.mmPerVox = mmPerVox;
msh.host = host;
msh.id = id;
msh.actor  = [];
msh.curvature = [];
% We always save the initial vertices so we can go back to the unsmoothed
% form and so that we can compute the vertex to gray map without
% distortion.
msh.initVertices = [];

% Empty specifications for data.  These are the data slots that mrMesh
% uses.
msh.data.camera_space = [];
msh.data.triangles = [];
msh.data.normals = [];
msh.data.vertices = [];
msh.data.rotation = [];
msh.data.colors = [];
msh.data.origin = [];

% These are hard-wired for now.
msh.grayLayers = 0;
msh.vertexGrayMap = [];

msh.decimate_reduction = decimateReduction;

msh.decimate_iterations = 0;
msh.decimate_subiterations = 0;
msh.decimate_preserve_edges = 0;
msh.decimate_preserve_topology = 1;
msh.decimate_boudary_vertex_deletion = 1;
msh.decimate_aspect_ratio = 25;
msh.decimate_degree = 20;
msh.smooth_sinc_method = 0;

% If you use sinc, smooth_relaxation is the passband parameter. Lower
% values will give more smoothing. Try .05 with 20-40 iterations for light
% smoothing, or .0001 with 150 iterations to fully relax a mesh. The main
% advantage of windowed sinc smoothing over the regular smoothing
% (Laplacian) is that it doesn't shrink the mesh with heavy smoothing.
msh.smooth_iterations = 50;
msh.smooth_relaxation = .25;
msh.smooth_feature_angle = 45;
msh.smooth_edge_angle = 15;
msh.smooth_boundary = 1;
msh.smooth_feature_angle_smoothing = 0;

% This smooth_pre boolean flag param is new. If true, an additional
% smoothing will be applied BEFORE decimation. Currently, it uses the same
% parameters as the post-decimation smoothing. Note that this param only
% applies to the build_mesh, open_gray and open_class commands (NOT smooth).
msh.smooth_pre = 0;
msh.curvature_mod_depth = 0.25;
msh.curvature_color = 1;

msh.conMat = [];
msh.lights = [];

return;


text:
function value = r8vec_amin ( n, a )

%*****************************************************************************80
%
%% R8VEC_AMIN returns the minimum absolute value in an R8VEC.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 October 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries in the array.
%
%    Input, real A(N), the array.
%
%    Output, real VALUE, the value of the entry
%    of smallest magnitude.
%
  value = min ( abs ( a(1:n) ) );

  return
end


text:
function ind = subsind(ebsd,subs)
% subindexing of EBSD data
%

% ordinary indexing by id
if (ischar(subs{1}) && subs{1}(1) == ':') || ...
    (isnumeric(subs{1}) && ~islogical(subs{1}))
  ind = reshape(1:length(ebsd),size(ebsd));

  ind = subsref(ind,struct('type','()','subs',{subs}));
  return

end

ind = subsind@EBSD(ebsd,subs);

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function model = gpsimMapTest

% GPSIMMAPTEST Tests the gradients of the GPSIMMAP model.
% FORMAT
% DESC tests the gradients and log likelihoods of the MAP
% approximation for the GPSIM model.
% RETURN model : the model that was used for testing.
%
% SEEALSO : gpsimMapCreate
%
% COPYRIGHT : Neil D. Lawrence, 2006

% SHEFFIELDML

colordef white
load demBarenco1;
B = model.comp{1}.B;
D = model.comp{1}.D;
S = model.comp{1}.S;
clear model

[y, yvar, gene, times, scale, rawExp, rawVar] = gpsimLoadBarencoData;

numGenes = size(gene, 1);
options = gpsimMapOptions(numGenes);
options.intPoints = 17;
options.B = B;
options.D = D;
options.S = S;

nonLinearities = {'linear', 'quadratic', 'negLogLogit', 'exp', 'sigmoid'};
nonLinearities = {'exp'};
for i = 1:length(nonLinearities)
  options.nonLinearity = nonLinearities{i};
  fprintf('Nonlinearity: %s\n', options.nonLinearity);
  model = gpsimMapCreate(numGenes, 1, times, y{1}, yvar{1}, options);
  params = gpsimMapExtractParam(model);

  model = gpsimMapExpandParam(model, params);

  % Give reasonable values to f.
  f = gsamp(zeros(1, size(model.f, 1)), model.K, 1);
  model = gpsimMapFunctionalExpandParam(model, f);

  gradientCheck(params, 'modelObjective', 'modelGradient', model);

  model.type = 'gpsimMapFunctional';
  f = gpsimMapFunctionalExtractParam(model);
  gradientCheck(f, 'modelObjective', 'modelGradient', model);
  hessianCheck(f, 'modelObjective', 'modelHessian', model);
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function handWritingTest
%%
clc
clear
close all
%% 获取目录下的所有txt文件名称
d = dir(['digits/trainingDigits/' '*.txt']); % struct 类型
dircell = struct2cell(d); %cell 类型
trainSetLen = size(dircell,2);
K = 4;
dataSize = 1024;
trainLabels = zeros(trainSetLen,1);
trainSet = [];
simpleTrainSet = zeros(1,dataSize);
simpleTestSet = zeros(1,dataSize);

%% 加载数据
fprintf('loading data...')
for i = 1:trainSetLen
    trainName =  dircell(1,i);
    trainFilename = cell2mat(trainName);
    trainLabels(i) = str2num(trainFilename(1));

    fid = fopen(['digits/trainingDigits/' trainFilename],'r');
    traindata = fscanf(fid,'%s');
    for j = 1:dataSize
        simpleTrainSet(j) =  str2num(traindata(j));
    end
    trainSet = [trainSet ; simpleTrainSet];
    fclose(fid);
end

d = dir(['digits/testDigits/' '*.txt']); % struct 类型
dircell = struct2cell(d); %cell 类型
testSetLen = size(dircell,2);
error = 0;
%% 测试数据
for k = 1:testSetLen
    testName =  dircell(1,k);
    testFilename = cell2mat(testName);
    testLabels = str2num(testFilename(1));

    fid = fopen(['digits/testDigits/' testFilename],'r');
    testdata = fscanf(fid,'%s');
    for j = 1:dataSize
        simpleTestSet(j) =  str2num(testdata(j));
    end
    classifyResult = KNN(simpleTestSet,trainSet,trainLabels,K);
    fprintf('识别数字为：%d  真实数字为：%d\n' , [classifyResult , testLabels])
    if(classifyResult~=testLabels)
        error = error+1;
    end
    fclose(fid);
end

fprintf('识别准确率为：%f\n',1-error/testSetLen)

end

repo: MachineLearning


text:
classdef SegmentationLossLogistic < dagnn.Loss
    methods
        function outputs = forward(obj, inputs, params)
            sz = size(inputs{2});
            mass = sz(1) * sz(2) + 1;

            outputs{1} = vl_nnloss(inputs{1}, inputs{2}, [], ...
                'loss', obj.loss, ...
                'instanceWeights', 1./mass) ;
            n = obj.numAveraged ;
            m = n + size(inputs{1},4) ;
            obj.average = (n * obj.average + double(gather(outputs{1}))) / m ;
            obj.numAveraged = m ;
        end

        function [derInputs, derParams] = backward(obj, inputs, params, derOutputs)
            sz = size(inputs{2});
            mass = sz(1) * sz(2) + 1;
            derInputs{1} = vl_nnloss(inputs{1}, inputs{2}, derOutputs{1}, ...
                'loss', obj.loss, ...
                'instanceWeights', 1./mass) ;
            derInputs{2} = [] ;
            derParams = {} ;
        end

        function obj = SegmentationLossLogistic(varargin)
            obj.load(varargin) ;
        end
    end
end

author: aimerykong


text:
function [output] = F_tanh(input_layer)
input = input_layer.a;
output = tanh(input);

end

author: singaxiong

repo: SignalGraph


text:
classdef mixture < iosr.dsp.audio
%MIXTURE Class of sound source separation mixture.
%
%   iosr.bss.mixture objects contain information about a mixture of sound
%   sources.
%
%   IOSR.BSS.MIXTURE is a subclass of IOSR.DSP.AUDIO.
%
%   IOSR.BSS.MIXTURE properties:
%       azi_sep         - The azimuthal separation of the widest sources
%                         (read-only)
%       decomp          - The result of the time-frequency decomposition
%                         (mixture) (read-only)
%       decomp_i        - The result of the time-frequency decomposition
%                         (interferer) (read-only)
%       decomp_t        - The result of the time-frequency decomposition
%                         (target) (read-only)
%       decomposition   - Set the time-frequency decomposition. The options
%                         are:
%                             'stft'        : short-time fourier transform
%                                             (default)
%                             'gammatone'   : gammatone filterbank
%       elevation       - The median elevation of the mixture (read-only)
%       filename_t      - Name of the target audio file (based on the
%                         mixture filename) (read-only)
%       filename_i      - Name of the interferer audio file (based on the
%                         Mixture filename) (read-only)
%       gammatone       - Settings for the gammatone filterbank
%                         decomposition. The property is a structure
%                         containing the following fields:
%                             'cfs'     : the centre frequencies of the
%                                         gammatone filterbank
%                             'frame'   : the frame length (in samples)
%       sofa_path       - Path to a SOFA file containing spatialisation data
%       ibm             - The ideal binary mask
%       irm             - The ideal ratio mask
%       int_fns         - A char array containing the filenames of all of
%                         the interfering sources (read-only)
%       interferers     - An array of interferer sources of type
%                         iosr.bss.source
%       signal_t        - The sampled data (target) (read-only)
%       signal_i        - The sampled data (interferer) (read-only)
%       stft            - Settings for the STFT decomposition. The property
%                         is a structure containing the following fields:
%                             'hop'     : the hop size of the STFT
%                             'win'     : the STFT window
%                         See iosr.dsp.stft for more information
%       target          - The target source of type iosr.bss.source
%       tfPower         - The (frame-based) time-frequency power (mixture)
%                         (read-only)
%       tfPower_i       - The (frame-based) time-frequency power
%                         (interferer) (read-only)
%       tfPower_t       - The (frame-based) time-frequency power (target)
%                         (read-only)
%       tir             - The target-to-interferer ratio (target or
%                         interferers are attenuated in order that their
%                         RMS amplitudes have this ratio)
%       wdo             - The w-disjoint orthogonality (read-only)
%       wdo_lw          - The loudness-weighted w-disjoint orthogonality
%                         (read-only)
%
%   IOSR.BSS.MIXTURE methods:
%       mixture         - Create the mixture
%       clearCache      - Clears the internally cached data, reducing
%                         file size.
%       copy            - Create an independent copy of the mixture, its
%                         sources, and any rendered files
%       applyMask       - Apply a time-frequency mask
%       deleteFiles     - Delete the mixture audio files
%       sound_t         - Replay the target
%       sound_i         - Replay the interferer
%       write           - Save the mixture, target, and interferers to audio files
%     Static methods:
%       maskCentroids   - Calculate the centroids of a time-frequency mask
%       mixWdo          - WDO of a mixture
%       mixWdo_lw       - Loudness-weighted WDO of a mixture
%       mixWdo_Stokes   - WDO calculated using Stokes's method
%
%   Note that target and interferer properties may be modified as
%   MIXTURE.TARGET.PROPERTY_NAME and MIXTURE.INTERFERERS(N).PROPERTY_NAME,
%   except for the sampling frequency FS, which cannot be overridden. This
%   ensures that the sampling frequencies are identical for the mixture and
%   its sources.
%
%   See also IOSR.DSP.AUDIO, IOSR.BSS.SOURCE, SOFALOAD, IOSR.DSP.STFT,
%   IOSR.AUDITORY.GAMMATONEFAST.

%   Copyright 2016 University of Surrey.

%TODO: Add more decompositions.

    properties (AbortSet)
        decomposition = 'stft'  % The time-frequency decomposition
        gammatone = struct      % Settings for the gammatone filterbank decomposition
        sofa_path               % Path to a SOFA file containing spatialisation data
        stft = struct           % Settings for the STFT decomposition
        tir = 0                 % The target to interferer ratio
        interferers             % An array of interferer sources of type iosr.bss.source
        target                  % The target source of type iosr.bss.source
    end

    properties (Dependent, SetAccess = protected)
        signal          % The sampled data (read-only)
    end

    properties (Dependent, SetAccess = private)
        azi_sep     % The azimuthal separation of the widest sources (read-only)
        decomp      % The result of the time-frequency decomposition (mixture) (read-only)
        decomp_i    % The result of the time-frequency decomposition (interferer) (read-only)
        decomp_t    % The result of the time-frequency decomposition (target) (read-only)
        elevation   % The median elevation of the mixture (read-only)
        filename_t  % Name of the target audio file (read-only)
        filename_i  % Name of the interferer audio file (read-only)
        ibm         % Ideal binary mask (read-only)
        irm         % Ideal ratio mask (read-only)
        int_fns     % Filenames of all of the interfering sources (read-only)
        numchans    % Number of audio channels in the mixture (read-only)
        signal_t    % Return target (read-only)
        signal_i    % Return interferer (read-only)
        tfPower     % The time-frequency power (mixture) (read-only)
        tfPower_i   % The time-frequency power (interferer) (read-only)
        tfPower_t   % The time-frequency power (target) (read-only)
        wdo         % Return the w-disjoint orthogonality metric (read-only)
        wdo_lw      % Return the loudness-weighted w-disjoint orthogonality metric (read-only)
        wdo_stokes  % Return the w-disjoint orthogonality metric using Stokes's method (read-only)
    end

    properties (Access = private)
        cached_decomp               % the cached mixture decomposition
        decomp_cached = false       % whether the mixture decomposition is cached
        cached_decomp_i             % the cached interferer decomposition
        decomp_i_cached = false     % whether the interferer decomposition is cached
        cached_decomp_t             % the cached target decomposition
        decomp_t_cached = false     % whether the target decomposition is cached
        cached_tfPower              % the cached mixture TF power
        tfPower_cached = false      % whether the mixture TF power is cached
        cached_tfPower_i            % the cached interferer TF power
        tfPower_i_cached = false    % whether the interferer TF power is cached
        cached_tfPower_t            % the cached target TF power
        tfPower_t_cached = false    % whether the target TF power is cached
    end

    methods

        % constructor
        function obj = mixture(target,interferers,varargin)
        %MIXTURE Create a mixture
        %
        %   OBJ = IOSR.BSS.MIXTURE(TARGET,INTERFERER) creates a mixture by
        %   summing together the target source and the interferer
        %   source(s). The sources are mixed together such that their RMS
        %   amplitudes are equal (target-to-interferer ratio is 0dB). The
        %   OBJ sampling rate is equal to the target sampling rate.
        %   Information about the sources' spatial location is ignored.
        %
        %   OBJ = IOSR.BSS.MIXTURE(...,'PARAMETER',VALUE) allows additional
        %   options to be specified. The options are ({} indicate
        %   defaults):
        %
        %       'decomposition' : {'stft'} | 'gammatone'
        %       'filename'      : {[]} | str
        %           A filename used when writing the file with the write()
        %           method. The filename may also be set when calling the
        %           write() method. Filenames for the target and interferer
        %           are determined automatically, by append the filename
        %           with '_target' and '_interferer' respectively.
        %       'fs'            : {obj.target.fs} | scalar
        %           The sampling frequency of the mixture. All HRTFs and/or
        %           sources will be resampled to this frequency each time
        %           the signal is requested.
        %       'gammatone'
        %           Settings for the gammatone filterbank
        %           decomposition. The property is a structure
        %           containing the following fields:
        %               'cfs'     : the centre frequencies of the gammatone
        %                           filterbank (default is 64 channels
        %                           equally spaced on the ERB scale between
        %                           20Hz and the Nyquist limit
        %               'frame'   : the frame length (in samples) (the
        %                           default is 20ms)
        %       'sofa_path'     : {[]} | str
        %           A path to a SOFA file containing HRTFs that are
        %           convolved with sources in order to generate the
        %           mixture.
        %       'stft'
        %           Settings for the STFT decomposition. The property
        %           is a structure containing the following fields:
        %               'hop'     : the hop size of the STFT (the default
        %                           is 512)
        %               'win'     : the STFT window (the default is 1024)
        %       'tir'           : {0} | scalar
        %           The RMS ratio of the target and interfer sources.
        %           Interferer sources are individually set to this level
        %           prior to their summation.
        %
        %   OBJ = IOSR.BSS.MIXTURE creates an empty mixture object with
        %   empty target and interferer sources.
        %
        %   To speed up time-frequency operations, the  IOSR.BSS.MIXTURE
        %   class caches time-frequency data internally, which can lead to
        %   large file sizes if the object is saved. Use the CLEARCACHE()
        %   method to empty the internal cache.
        %
        %   Note that this is a handle class, as is IOSR.BSS.SOURCE. Target
        %   and interferer(s) are hence passed by reference. Use the COPY()
        %   method to create an independent copy of the mixture and its
        %   sources.

            if nargin > 0

                assert(nargin>1, 'iosr:mixture:nargin', 'Not enough input arguments')

                propNames = {'filename','fs','sofa_path','tir','decomposition','stft','gammatone'};

                % set sources
                obj.target = target;
                obj.interferers = interferers;

                % defaults
                obj.fs = obj.target.fs;
                obj.sofa_path = [];
                obj.tir = 0;
                obj.rendered = false;

                % default decomposition settings
                obj.stft = struct('win',1024,'hop',512);
                obj.gammatone = struct(...
                    'cfs',iosr.auditory.makeErbCFs(20, obj.fs/2, 64),...
                    'frame', round(0.01*obj.fs)...
                );

                % read parameter/value inputs
                if nargin > 1 % if parameters are specified
                    obj.set_properties(varargin,propNames)
                end

                % set sample rate to SOFA file if set and fs not specified
                if ~isempty(obj.sofa_path) && all(~strcmp('fs',varargin))
                    SOFAobj = SOFAload(obj.sofa_path); % load SOFA object
                    obj.fs = SOFAobj.Data.SamplingRate;
                end

                % ensure fs for sources matches instance
                obj.target.fs = obj.fs;
                for n = 1:numel(obj.interferers)
                    obj.interferers(n).fs = obj.fs;
                end

                % set parent
                obj.target.parent = obj;
                for n = 1:numel(obj.interferers)
                    obj.interferers(n).parent = obj;
                end

            end

        end

        function z = applyMask(obj,m)
        %APPLYMASK Apply a time-frequency mask to the mixture
        %
        %   Z = IOSR.BSS.MIXTURE.APPLYMASK(M) applies the time-frequency
        %   mask M to the mixture. The mixture is transformed according to
        %   the objects decomposition settings. The transform is then
        %   multiplied with M. Lastly, the result is then inverse
        %   transformed (or reynthesised). The time-domain output Z is
        %   returned.

            switch lower(obj.decomposition)
                case 'stft'
                    s = obj.decomp;
                    % correct dims
                    d = obj.stft_shifted_dims(s);
                    s = ipermute(s,d);
                    m = ipermute(m,d);
                    for c = 1:obj.numchans
                        % apply mask
                        z(:,c) = iosr.bss.applyMask( ...
                            s(:,:,c), ...
                            m(:,:,min(c,size(m,3))), ...
                            obj.stft.win, ...
                            obj.stft.hop, ...
                            obj.fs ...
                        ); %#ok<AGROW>
                    end
                case 'gammatone'
                    x = obj.signal;
                    for c = 1:obj.numchans
                        z(:,c) =  iosr.bss.resynthesise( ...
                            x(:,c), ...
                            obj.fs, ...
                            iosr.bss.cfs2fcs(obj.gammatone.cfs, obj.fs), ...
                            m(:,:,c), ...
                            'frame_length', obj.gammatone.frame, ...
                            'filter', 'sinc', ...
                            'kernel', gausswin(obj.gammatone.frame) ...
                        ); %#ok<AGROW>
                    end
            end

            z = obj.setlength(z,length(obj.signal));

        end

        function sound_t(obj)
        %SOUND_T Replay the target signal
        %
        %   IOSR.BSS.MIXTURE.SOUND_T() replays the target signal.

            obj.replay(obj.signal_t)

        end

        function sound_i(obj)
        %SOUND_I Replay the interferer signal
        %
        %   IOSR.BSS.MIXTURE.SOUND_I() replays the interferer signal.

            obj.replay(obj.signal_i)

        end

        function write(obj,filename)
        %WRITE Save the mixture to an audio file
        %
        %   IOSR.BSS.MIXTURE.WRITE() writes the mixture to an audio file
        %   specified by MIXTURE.FILENAME, and also writes the target and
        %   interferer signals to automatically-determined file names
        %   (MIXTURE.FILENAME_T and MIXTURE.FILENAME_I respectively).
        %
        %   IOSR.BSS.MIXTURE.WRITE(FILENAME) uses the specified FILENAME
        %   and updates MIXTURE.FILENAME.

            if exist('filename','var')==1
                fn = filename;
            else
                fn = obj.filename;
            end

            if obj.rendered && exist(fn, 'file') == 2
                if strcmp(obj.filename,fn)
                    % nothing to do
                else
                    copyfile(obj.filename,fn);
                    copyfile(obj.filename_t,obj.make_target_filename(fn));
                    copyfile(obj.filename_i,obj.make_interferer_filename(fn));
                end
            else
                % check filename is valid
                obj.filename = fn;
                assert(ischar(obj.filename) && ~isempty(obj.filename),'iosr:mixture:invalidFilename','FILENAME must be a non-empty char array. Set filename as MIXTURE.FILENAME or MIXTURE.WRITE(FILENAME).')

                % normalize
                M = obj.signal;
                T = obj.signal_t;
                I = obj.signal_i;
                [~,gain] = obj.normalize([M T I]);
                M = M.*gain;
                T = T.*gain;
                I = I.*gain;

                % ensure path
                obj.ensure_path(obj.filename)

                % write audio files
                audiowrite(obj.filename,M,obj.fs);
                audiowrite(obj.filename_t,T,obj.fs);
                audiowrite(obj.filename_i,I,obj.fs);

                % set flag to indicate mixture has been rendered to an audio file
                obj.rendered = true;

            end

        end

        function deleteFiles(obj)
        %DELETEFILES Delete the mixture audio files

            if exist(obj.filename, 'file') == 2
                delete(obj.filename);
            end
            if exist(obj.filename_t, 'file') == 2
                delete(obj.filename_t);
            end
            if exist(obj.filename_i, 'file') == 2
                delete(obj.filename_i);
            end
            obj.rendered = false;

        end

        function clearCache(obj)
        %CLEARCACHE Clear the internal cache

            obj.updateCachedFalse();
            obj.cached_decomp = [];
            obj.cached_decomp_i = [];
            obj.cached_decomp_t = [];
            obj.cached_tfPower = [];
            obj.cached_tfPower_i = [];
            obj.cached_tfPower_t = [];
        end

        % set/validate properties

        % set decomposition
        function set.decomposition(obj,val)
            assert(any(strcmpi(val,{'stft','gammatone'})), 'iosr:mixture:invalidDecomposition', '''decomposition'' must be ''stft'' or ''gammatone''');
            obj.updateCachedFalse();
            obj.decomposition = val;
        end

        % set gammatone settings
        function set.gammatone(obj,val)
            assert(isstruct(val), 'iosr:mixture:invalidGammatoneStruct', '''gammatone'' must be a struct')
            assert(all(ismember(fieldnames(val)',{'cfs','frame'})), 'iosr:mixture:invalidGammatoneFields', 'gammatone structure must contain fields ''cfs'' and ''frame''')
            obj.updateCachedFalse;
            obj.gammatone = val;
        end

        % set stft settings
        function set.stft(obj,val)
            assert(isstruct(val), 'iosr:mixture:invalidStftStruct', '''stft'' must be a struct')
            assert(all(ismember(fieldnames(val)',{'win','hop'})), 'iosr:mixture:invalidStftFields', 'stft structure must contain fields ''win'' and ''hop''')
            obj.updateCachedFalse;
            obj.stft = val;
        end

        % set tir
        function set.tir(obj,val)
            obj.tir = val;
            obj.property_changed('tir',val);
        end

        % validate sofa_path
        function set.sofa_path(obj,val)
            assert(ischar(val) || isempty(val), 'iosr:mixture:invalidSofaPath', 'sofa_path must be a char array or an empty array')
            if ~isempty(val)
                assert(exist(val,'file')==2, 'iosr:mixture:invalidSofaPath', 'SOFA file does not exist')
            end
            obj.sofa_path = val;
            obj.property_changed('sofa_path',val);
        end

        % validate interferers
        function set.interferers(obj,val)
            assert(isa(val,'iosr.bss.source'), 'iosr:mixture:invalidInterferers', 'INTERFERERS must be of type source')
            obj.interferers = val;
            obj.property_changed('interferers',val);
        end

        % validate target
        function set.target(obj,val)
            assert(isa(val,'iosr.bss.source') && numel(val)==1, 'iosr:mixture:invalidTarget', 'TARGET must be a scalar of type source')
            obj.target = val;
            obj.property_changed('target',val);
        end

        % dependent properties

        % get decomposition result (mixture)
        function s = get.decomp(obj)
            if obj.decomp_cached
                s = obj.cached_decomp;
            else
                s = obj.decompose(obj.signal);
                obj.cached_decomp = s;
                obj.decomp_cached = true;
            end
        end

        % get decomposition result (interferer)
        function s = get.decomp_i(obj)
            if obj.decomp_i_cached
                s = obj.cached_decomp_i;
            else
                s = obj.decompose(obj.signal_i);
                obj.cached_decomp_i = s;
                obj.decomp_i_cached = true;
            end
        end

        % get decomposition result (target)
        function s = get.decomp_t(obj)
            if obj.decomp_t_cached
                s = obj.cached_decomp_t;
            else
                s = obj.decompose(obj.signal_t);
                obj.cached_decomp_t = s;
                obj.decomp_t_cached = true;
            end
        end

        % get azimuthal separation
        function s = get.azi_sep(obj)
            [s,~] = get_loc(obj);
        end

        % get median elevation
        function e = get.elevation(obj)
            [~,e] = get_loc(obj);
        end

        % target filename
        function fn = get.filename_t(obj)
            fn = obj.make_target_filename(obj.filename);
        end

        % interferer filename
        function fn = get.filename_i(obj)
            fn = obj.make_interferer_filename(obj.filename);
        end

        % ibm
        function m = get.ibm(obj)
            switch lower(obj.decomposition)
                case 'stft'
                    for c = 1:obj.numchans
                        [~,m(:,:,c)] = iosr.bss.idealMasks( ...
                            obj.decomp_t(:,:,c), ...
                            obj.decomp_i(:,:,c) ...
                        ); %#ok<AGROW>
                    end
                case 'gammatone'
                    T = obj.tfPower_t;
                    I = obj.tfPower_i;
                    m = +(T>I);
            end
        end

        % irm
        function m = get.irm(obj)
            switch lower(obj.decomposition)
                case 'stft'
                    for c = 1:obj.numchans
                        m(:,:,c) = iosr.bss.idealMasks( ...
                            obj.decomp_t(:,:,c), ...
                            obj.decomp_i(:,:,c) ...
                        ); %#ok<AGROW>
                    end
                case 'gammatone'
                    T = obj.tfPower_t;
                    I = obj.tfPower_i;
                    m = T./(T+I);
            end
            m(isnan(m) | isinf(m)) = 0;
        end

        % filenames of all interferers
        function fns = get.int_fns(obj)
            for n = 1:length(obj.interferers)
                if n==1
                    fns = obj.interferers(n).filename;
                else
                    fns = [fns ', ' obj.interferers(n).filename]; %#ok<AGROW>
                end
            end
        end

        % number of audio channels
        function n = get.numchans(obj)
            if obj.hrtf_is_set()
                s = SOFAload(obj.sofa_path);
                n = size(s.Data.IR,2);
            else
                iN = [obj.interferers.numchans];
                n = max([obj.target.numchans; iN(:)]);
            end
        end

        % return target signal
        function signal_t = get.signal_t(obj)
            if obj.rendered && exist(obj.filename_t,'file')==2 % don't bother calculating
                signal_t = audioread(obj.filename_t);
            else % calculate
                signal_t = return_source(obj,obj.target);
                if obj.tir<0
                    % attenuate according to TIR
                    Trms = iosr.dsp.rms(signal_t(:));
                    Irms = iosr.dsp.rms(obj.signal_i(:));
                    signal_t = signal_t./(Trms/Irms); % match to interferer
                    signal_t = signal_t.*(10^(obj.tir/20)); % attenuate
                end
            end
        end

        % return interferer signal
        function signal_i = get.signal_i(obj)
            if obj.rendered && exist(obj.filename_i,'file')==2 % don't bother calculating
                signal_i = audioread(obj.filename_i);
            else % calculate
                signal_i = zeros(1, obj.numchans); % initialise
                maxlength = 0; % initialise
                for n = 1:numel(obj.interferers) % step through each interferer
                    % source signal (ensure 2-channel)
                    s = obj.return_source(obj.interferers(n));
                    % ensure signals are same length, or zero-pad
                    maxlength = max([length(s) maxlength]);
                    s = obj.setlength(s,maxlength);
                    signal_i = obj.setlength(signal_i,maxlength);
                    % add source to interferer
                    signal_i = signal_i + s;
                end
                if obj.tir>=0
                    % attenuate according to TIR
                    Trms = iosr.dsp.rms(obj.signal_t(:));
                    Irms = iosr.dsp.rms(signal_i(:));
                    signal_i = signal_i./(Irms/Trms); % match to interferer
                    signal_i = signal_i./(10^(obj.tir/20)); % attenuate
                end
            end
        end

        % return mixture signal
        function signal = get.signal(obj)
            if obj.rendered && exist(obj.filename,'file')==2 % don't bother calculating
                signal = audioread(obj.filename);
            else % calculate
                % return target and interferers
                T = obj.signal_t;
                I = obj.signal_i;
                % mix, ensuring equal length
                maxlength = max([length(T) length(I)]);
                signal = obj.setlength(T,maxlength) + obj.setlength(I,maxlength);
            end
        end

        % return time-frequency power (mixture)
        function val = get.tfPower(obj)
            if obj.tfPower_cached
                val = obj.cached_tfPower;
            else
                val = obj.tf_power(obj.decomp);
                obj.cached_tfPower = val;
                obj.tfPower_cached = true;
            end
        end

        % return time-frequency power (interferer)
        function val = get.tfPower_i(obj)
            if obj.tfPower_i_cached
                val = obj.cached_tfPower_i;
            else
                val = obj.tf_power(obj.decomp_i);
                obj.cached_tfPower_i = val;
                obj.tfPower_i_cached = true;
            end
        end

        % return time-frequency power (target)
        function val = get.tfPower_t(obj)
            if obj.tfPower_t_cached
                val = obj.cached_tfPower_t;
            else
                val = obj.tf_power(obj.decomp_t);
                obj.cached_tfPower_t = val;
                obj.tfPower_t_cached = true;
            end
        end

        % return wdo
        function w = get.wdo(obj)
            w = obj.mixWdo(abs(obj.decomp_t), abs(obj.decomp_i));
        end

        % return loudness-weighted wdo
        function w = get.wdo_lw(obj)
            [S,f] = obj.decompose(obj.signal_t);
            w = obj.mixWdo_lw(abs(S), abs(obj.decomp_i), f);
        end

        % return Stokes's wdo
        function w = get.wdo_stokes(obj)
            w = obj.mixWdo_Stokes(obj.irm);
        end

    end

    methods (Static, Access = public)

        function [C, Ct, Cf] = maskCentroids(m, fs, decomposition, nfft, hop)
        %MASKCENTROIDS Return various time-frequency mask centroids
        %
        %   C = IOSR.BSS.MIXTURE.MASKCENTROIDS(M,FS,NFFT,HOP) returns the
        %   time-frequency mask centroid C calculated from mask M, sample
        %   rate FS, FFT-length NFFT, and and hop size HOP. The mask
        %   centroid is a unitless measure of the centroid of a
        %   time-frequency mask. The mask centroid is a measure of the
        %   frequency content of a time-frequency mask.
        %
        %   [C,Ct,Cf] = IOSR.BSS.MIXTURE.MASKCENTROIDS(...) returns the
        %   time-related spectral centroid (in Hz) Ct and the
        %   frequency-related spectral centroid (in seconds) Cf.

            assert(ischar(decomposition), 'iosr:mixture:invalidDecomposition', '''decomposition'' must be a char array')

            numchans = size(m,3);
            frame_frequency = fs/hop;
            switch lower(decomposition)
                case 'stft'
                    assert(isscalar(nfft), 'iosr:mixture:invalidNfft', '''nfft'' must be a scalar')
                    quefrency = fs/nfft;
                case 'gammatone'
                    assert(isvector(nfft), 'iosr:mixture:invalidNfft', '''cfs'' must be a scalar')
                    cfs = nfft;
                    quefrency = fs/min(diff(iosr.auditory.erbRate2hz(cfs)));
                otherwise
                    error('iosr:mixture:unknownDecomp','Unknown decomposition.')
            end

            f_dct = repmat((((0:size(m,1)-1)./size(m,1)).*frame_frequency)',1,size(m,2));
            c_dct = repmat(((0:size(m,2)-1)./size(m,2)).*quefrency,size(m,1),1);

            C = zeros(1, numchans);
            Ct = zeros(1, numchans);
            Cf = zeros(1, numchans);
            for c = 1:numchans
                X = abs(dct(m(:,:,c)));
                sumX = sum(sum(X));
                C(c) = sum(sum(X.*f_dct.*c_dct))./sumX;
                Ct(c) = sum(sum(X.*f_dct))./sumX;
                Cf(c) = sum(sum(X.*c_dct))./sumX;
            end

        end

        function w = mixWdo(st, si)
        %MIXWDO Calculate the mixture w-disjoint orthogonality
        %
        %   W = IOSR.BSS.MIXTURE.MIXWDO(St,Si) calculates the w-disjoint
        %   orthogonality of a target and interferer(s) from their
        %   time-frequency magnitudes St and Si. The w-disjoint
        %   orthogonality is estimated by the two-dimensional
        %   cross-correlation coefficient of the source magnitudes.

            numchans = size(st,3);
            w = zeros(1,numchans);
            for c = 1:numchans
                w(c) = corr2(st(:,:,c),si(:,:,c));
            end

        end

        function w = mixWdo_lw(st, si, f)
        %MIXWDO_LW Calculate the loudness-weighted w-disjoint orthogonality
        %
        %   W = IOSR.BSS.MIXTURE.MIXWDO_LW(St,Si,F) calculates the
        %   loudness-weighted w-disjoint orthogonality of a target and
        %   interferer(s) from their time-frequency magnitudes St and Si.
        %   The weighting is calculated from the ISO-226 65-phon equal
        %   loudness contour using the frequency vector F containing the
        %   frequency of each bin of St and Si. The magnitudes are
        %   multiplied with the loudness-weighting factor prior to
        %   calculating the WDO.
        %
        %   See also IOSR.AUDITORY.LOUDWEIGHT.

            lw = iosr.auditory.loudWeight(f(:), 65)';
            lw = repmat(lw,size(st,1),1,size(st,3));
            w = iosr.bss.mixture.mixWdo(st.*lw,si.*lw);

        end

        function w = mixWdo_Stokes(irm)
        %MIXWDO_LW Calculate the w-disjoint orthogonality using Stokes's method
        %
        %   W = IOSR.BSS.MIXTURE.MIXWDO_STOKES(IRM) calculates the
        %   w-disjoint orthogonality using Stokes's method [1]. The method
        %   utilizes the ideal ratio mask (IRM), which is 0.5 when the
        %   target and interferer have equal power, and tends towards 0 or
        %   1 when either source dominates.
        %
        %   References
        %
        %   [1] Stokes, Tobias W. (2015) Improving the perceptual quality
        %       of single-channel blind audio source separation. Doctoral
        %       thesis, University of Surrey.

            h = hist(irm(:),11)./numel(irm);
            hw = [0:0.2:1 0.8:-0.2:0];
            w = sum(h.*hw);

        end

    end

    methods (Static, Access = private)

        function y = setlength(x,signal_length)
        %SETLENGTH Crop or zero-pad signal to specified length

            d = size(x);
            if size(x,1)>signal_length % need to crop
                subsidx = [{1:signal_length} repmat({':'},1,ndims(x)-1)];
                y = x(subsidx{:});
            elseif size(x,1)<signal_length % need to zero-pad
                y = [x; zeros([signal_length-size(x,1),d(2:end)])];
            else % do nothing
                y = x;
            end

        end

        function fn = append_filename(filename,append)
        %APPEND_FILENAME Append strings to MIXTURE.FILENAME

            if isempty(filename) % do nothing if filename is empty
                fn = [];
            else % append
                [filepath,name,ext] = fileparts(filename); % break up filename
                newname = [name append ext]; % append
                if isempty(filepath) % only filename specified
                    fn = newname;
                else % path specified
                    fn = [filepath filesep newname];
                end
            end

        end

        function d = stft_shifted_dims(s)
        %STFT_SHIFTED_DIMS Dim order for STFT, swapping time and bins

            dim_order = 1:ndims(s);
            if ndims(s) > 2 %#ok<ISMAT>
                d = dim_order([2 1 dim_order(3:end)]);
            else
                d = [2 1];
            end
        end

    end

    methods (Access = private)

        function c = hrtf_is_set(obj)
        %HRTF_IS_SET Determine whether HRTFs are specified

            c = ~isempty(obj.sofa_path);

        end

        function s = return_source(obj,src)
        %RETURN_SOURCE Return signal from specified source

            % ensure correct channel count
            if ~src.precomposed && obj.hrtf_is_set()
                src.numchans = 1; % signal will be spatialised
            else
                src.numchans = obj.numchans; % signal will be mixed directly
            end
            x = src.signal;
            % return/convolve signal
            if ~src.precomposed && obj.hrtf_is_set() % convolve
                s = obj.spat(obj.sofa_path,x,src.azimuth,src.elevation,obj.fs);
            else % return directly
                s = x;
            end

        end

        function [s,em] = get_loc(obj)
        %GET_LOC Return location information

            % get angles from the sources
            a = zeros(numel(obj.interferers)+1,1);
            e = a;
            for n = 1:numel(obj.interferers)
                a(n) = obj.interferers(n).azimuth;
                e(n) = obj.interferers(n).elevation;
            end
            a(end) = obj.target.azimuth;
            e(end) = obj.target.elevation;

            % get separation
            if any(a<0) % assume angles are -179:180
                s = abs(max(a)-min(a));
            else % assume angles are 0:359
                s = (360-max(a))+min(a);
            end
            s = mod(s,180);

            % get elevation
            em = median(e);

        end

        function fn = make_target_filename(obj,filename)
        %MAKE_TARGET_FILENAME Return the automated target filename

            fn = obj.append_filename(filename,'_target');

        end

        function fn = make_interferer_filename(obj,filename)
        %MAKE_INTERFERER_FILENAME Return the automated interferer filename

            fn = obj.append_filename(filename,'_interferer');

        end

        function [s,f,t] = decompose(obj,x)
        %DECOMPOSE Transform a signal into the time-frequency domain

            switch lower(obj.decomposition)
                case 'stft'
                    for c = 1:obj.numchans
                        [s(:,:,c),f,t] = iosr.dsp.stft(x(:,c), obj.stft.win, obj.stft.hop, obj.fs); %#ok<AGROW>
                    end
                    % put time down column
                    s = permute(s,obj.stft_shifted_dims(s));
                    s = obj.setlength(s, obj.get_frame_count(length(obj.signal)));
                case 'gammatone'
                    f = obj.gammatone.cfs;
                    t = (0:length(x)-1)./obj.fs;
                    for c = 1:obj.numchans
                        s(:,:,c) = iosr.auditory.gammatoneFast(x(:,c), obj.gammatone.cfs, obj.fs); %#ok<AGROW>
                    end
            end

        end

        function y = tf_power(obj,x)
        %TF_POWER Calculate TF power.

            switch lower(obj.decomposition)
                case 'gammatone'
                    frame_count = obj.get_frame_count(size(x,1));
                    y = zeros(frame_count,length(obj.gammatone.cfs),obj.numchans);
                    for m = 1:frame_count
                        samples = (m-1)*obj.gammatone.frame+1:m*obj.gammatone.frame;
                        for N = 1:size(x,3)
                            y(m,:,N) = sum(x(samples,:,N).^2);
                        end
                    end
                case 'stft'
                    y = abs(x).^2;
            end

        end

        function frame_count = get_frame_count(obj, signal_length)
        %GET_FRAME_COUNT Calculate number of T-F frames.

            assert(isscalar(signal_length), 'iosr:mixture:invalidSignalLength', 'SIGNAL_LENGTH must be a scalar.')

            switch lower(obj.decomposition)
                case 'gammatone'
                    frame_count = floor(signal_length/obj.gammatone.frame);
                case 'stft'
                    frame_count = fix((signal_length-(obj.nfft-obj.stft.hop))/obj.stft.hop);
            end

        end

        function n = nfft(obj)
        %NFFT return the fft length

            if isscalar(obj.stft.win)
                n = obj.stft.win;
            else
                n = length(obj.stft.win);
            end
        end

        function updateCachedFalse(obj)
            obj.decomp_cached = false;
            obj.decomp_i_cached = false;
            obj.decomp_t_cached = false;
            obj.tfPower_cached = false;
            obj.tfPower_i_cached = false;
            obj.tfPower_t_cached = false;
        end

    end

    methods(Access = protected)

        function property_changed(obj,name,val)
        %PROPERTY_CHANGED handle property changes

            obj.rendered = false;
            obj.updateCachedFalse();

            switch lower(name)
                case 'fs'
                    obj.target.fs = val;
                    for n = 1:numel(obj.interferers)
                        obj.interferers(n).fs = val;
                    end
            end

        end

        function cpObj = copyElement(obj)
        %COPYELEMENT Overload copy method with additional functionality

            % Make a shallow copy of all properties
            cpObj = copyElement@iosr.dsp.audio(obj);

            % Changed rendered file name
            cpObj.filename = obj.append_filename(cpObj.filename,'_copy');

            % copy files
            if obj.rendered
                if exist(obj.filename,'file')==2
                    copyfile(obj.filename,cpObj.filename)
                end
                if exist(obj.filename_i,'file')==2
                    copyfile(obj.filename_i,cpObj.filename_i)
                end
                if exist(obj.filename_t,'file')==2
                    copyfile(obj.filename_t,cpObj.filename_t)
                end
            end

            % Make a deep copy of target and interferers
            cpObj.target = copy(obj.target);
            cpObj.interferers = copy(obj.interferers);

        end

    end

end

author: IoSR-Surrey

repo: MatlabToolbox


text:
function [ SparsityMatrix ] = sparsity_grid( x, side, width, width_end)
%SIMILARITYNEIGHBOR Summary of this function goes here
%   Detailed explanation goes here

    % this assumes that the patch is laid out with first column, then second
    % column, ... final column (column major)

    SimilarityMatrix = zeros(side*side);
    for i=1:width
        SimilarityMatrix = (similarity_neighbor_grid_further(x, side, [1,2,3,4], i) | SimilarityMatrix);
    end

    SimilarityMatrix_end = zeros(side*side);
    for i=1:width_end
        SimilarityMatrix_end = (similarity_neighbor_grid_further(x, side, [1,2,3,4], i) | SimilarityMatrix_end);
    end

    SparsityMatrix = double(SimilarityMatrix_end & (~SimilarityMatrix));

    assert(isequal(SparsityMatrix, SparsityMatrix'));
end

author: TadasBaltrusaitis

repo: OpenFace


text:
function grid = getGrid(~, N, dom)
%GETGRID   Returns a grid correspoding to a SPINOP object.
%
% See also SPINOP.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

xx = trigpts(N, dom(1:2));
grid = {xx};

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function new_X = extract(X,idx)
%EXTRACT Creates a new ktensor with only the specified components.
%
%   Y = EXTRACT(X,S) selected the subset of components in X as defined by
%   S. It should be the case that S is a subset of [1,...,NCOMPONENTS(X)].
%
%   See also KTENSOR, NCOMPONENTS.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2012) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

%% Set-up
N = ndims(X);
%% Extract
new_lambda = X.lambda(idx);
new_U = cell(N,1);
for i = 1 : N
    new_U{i} = X.u{i}(:,idx);
end
new_X = ktensor(new_lambda, new_U);

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
classdef prtClassMatchedSubspace < prtClass

    properties (SetAccess=private)

        name = 'Matched Subspace'
        nameAbbreviation = 'MatchedSubspace'
        isNativeMary = true;

    end

    properties
        inferH1subspace = true;
        inferH0subspace = false;

        nH1components = 2;
        nH0components = 2;

        proj1
        proj0
    end

    methods
        function self = prtClassKnn(varargin)
            self = prtUtilAssignStringValuePairs(self,varargin{:});
            self.verboseStorage = true;
        end
    end

    methods (Access=protected, Hidden = true)

        function self = trainAction(self,ds)

            if self.inferH1subspace
                ds1 = ds.retainClassesByInd(2);
                x = ds1.X;
                [u,s] = svds(x',self.nH1components);
                p = u*(u'*u)^(-1/2)*u';
            else
                p = eye(ds.nFeatures);
            end
            self.proj1 = p;

            if self.inferH0subspace
                ds0 = ds.retainClassesByInd(1);
                x = ds0.X;
                [u,s] = svds(x,self.nH0components);
                p = u*(u'*u)^(-1/2)*u';
            else
                p = eye(ds.nFeatures);
            end
            self.proj0 = p;

        end

        function yOut = runAction(self,ds)

            x = ds.X;
            h1 = diag(x*self.proj1*x');
            h0 = diag(x*self.proj0*x');
            yOut = ds;
            %             yOut.X = log(h1)-log(h0);
            h1 = abs(h1);
            h0 = abs(h0);
            yOut.X = h1./h0;

        end
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function varargout = removeMeshFaces(v, f, fI)
%REMOVEMESHFACES Remove faces from a mesh by face indices.
%   [V2, F2] = removeMeshFaces(V, F, FI) removes faces from the mesh by
%   the face indices FI into faces F of the mesh. The mesh is represented
%   by the vertex array V and the face array F. The result is the new set
%   of vertices V2 and faces F2 without the faces indexed by FI. FI can be
%   either a linear or a logical index.
%
%   [V2, F2] = removeMeshFaces(MESH, FI) with the struct MESH containing
%   the fields "vertices" (V) and "faces" (F)
%
%   MESH2 = removeMeshFaces(V, F, FI) with the struct MESH2 containing the
%   fields "vertices" (V2) and "faces" (F2)
%
%   MESH2 = removeMeshFaces(MESH, FI) with the structs MESH and MESH2
%   containing the fields "vertices" (V, V2) and "faces" (F, F2)
%
%   Example
%     [v, f] = createSoccerBall;
%     f = triangulateFaces(f);
%     fI = true(length(f),1);
%     fI(1:length(f)/2) = false;
%     [v2, f2] = removeMeshFaces(v, f, fI);
%     drawMesh(v, f, 'faceColor', 'none', 'faceAlpha', .2);
%     drawMesh(v2, f2, 'faceAlpha', .7);
%     view(3); axis equal
%
%   See also
%   meshes3d, drawMesh
%
% ---------
% Authors: oqilipo, David Legland
% Created: 2017-07-04

% parse inputs
narginchk(2,3)
nargoutchk(1,2)

if nargin == 2
    fI = f;
    [v, f] = parseMeshData(v);
end

p = inputParser;
isIndexToFaces = @(x) ...
    (islogical(x) && isequal(length(fI), size(f,1))) || ...
    (all(floor(x)==x) && min(x)>=1 && max(x)<=size(f,1));
addRequired(p,'fI',isIndexToFaces)
parse(p, fI);
if ~islogical(p.Results.fI)
    fI=false(size(f,1),1);
    fI(p.Results.fI)=true;
else
    fI=p.Results.fI;
end

% algorithm
f2 = f(~fI,:);
[unqVertIds, ~, newVertIndices] = unique(f2);
v2 = v(unqVertIds,:);
f2 = reshape(newVertIndices,size(f2));

% parse outputs
if nargout == 1
    mesh2.vertices=v2;
    mesh2.faces=f2;
    varargout{1}=mesh2;
else
    varargout{1}=v2;
    varargout{2}=f2;
end

end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function res = isreal(Y)
%ISREAL (overloaded)

res = isreal(Y.basis);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function value = p05_f ( n, x )

%*****************************************************************************80
%
%% P05_F evaluates the function for problem p05.
%
%  Discussion:
%
%    Step function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 August 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of evaluation points.
%
%    Input, real X(N,1), the evaluation points.
%
%    Output, real VALUE(N,1), the function values.
%
  value =       cos (  7.0 * x ) ...
        + 5.0 * cos ( 11.2 * x ) ...
        - 2.0 * cos ( 14.0 * x ) ...
        + 5.0 * cos ( 31.5 * x ) ...
        + 7.0 * cos ( 63.0 * x );

  return
end