text:
function [parameters, LL, ht, VCVrobust, VCV, scores, diagnostics] = tarch(epsilon, p, o, q, error_type, tarch_type, startingvals, options)
% TARCH(P,O,Q) parameter estimation with different error distributions:
% Normal, Students-T, Generalized Error Distribution, Skewed T
% Estimation of ARCH or GARCH models if o=0 and tarch_type=2
% Estimation of TARCH or GJR asymmetric models if o>0 and tarch_type=1 or 2
%
% USAGE:
%   [PARAMETERS] = tarch(EPSILON,P,O,Q)
%   [PARAMETERS,LL,HT,VCVROBUST,VCV,SCORES,DIAGNOSTICS] =
%                                   tarch(EPSILON,P,O,Q,ERROR_TYPE,TARCH_TYPE,STARTINGVALS,OPTIONS)
%
% INPUTS:
%   EPSILON      - A column of mean zero data
%   P            - Positive, scalar integer representing the number of symmetric innovations
%   O            - Non-negative scalar integer representing the number of asymmetric innovations (0
%                    for symmetric processes)
%   Q            - Non-negative, scalar integer representing the number of lags of conditional
%                    variance (0 for ARCH)
%   ERROR_TYPE   - [OPTIONAL] The error distribution used, valid types are:
%                    'NORMAL'    - Gaussian Innovations [DEFAULT]
%                    'STUDENTST' - T distributed errors
%                    'GED'       - Generalized Error Distribution
%                    'SKEWT'     - Skewed T distribution
%   TARCH_TYPE   - [OPTIONAL] The type of variance process, either
%                    1 - Model evolves in absolute values
%                    2 - Model evolves in squares [DEFAULT]
%   STARTINGVALS - [OPTIONAL] A (1+p+o+q), plus 1 for STUDENTST OR GED (nu), plus 2 for SKEWT
%                    (nu,lambda), vector of starting values.
%                     [omega alpha(1) ... alpha(p) gamma(1) ... gamma(o) beta(1) ... beta(q) [nu lambda]]'.
%   OPTIONS      - [OPTIONAL] A user provided options structure. Default options are below.
%
% OUTPUTS:
%   PARAMETERS   - A 1+p+o+q column vector of parameters with
%                  [omega alpha(1) ... alpha(p) gamma(1) ... gamma(o) beta(1) ... beta(q) [nu lambda]]'.
%   LL           - The log likelihood at the optimum
%   HT           - The estimated conditional variances
%   VCVROBUST    - Robust parameter covariance matrix
%   VCV          - Non-robust standard errors (inverse Hessian)
%   SCORES       - Matrix of scores (# of params by t)
%   DIAGNOSTICS  - Structure of optimization outputs and other values useful for functions calling TARCH.
%
% COMMENTS:
% The following (generally wrong) constraints are used:
%   (1) omega > 0
%   (2) alpha(i) >= 0 for i = 1,2,...,p
%   (3) gamma(i) + alpha(i) > 0 for i=1,...,o
%   (3) beta(i)  >= 0 for i = 1,2,...,q
%   (4) sum(alpha(i) + 0.5*gamma(j) + beta(k)) < 1 for i = 1,2,...p and
%   j = 1,2,...o, k=1,2,...,q
%   (5) nu>2 of Students T and nu>1 for GED
%   (6) -.99<lambda<.99 for Skewed T
%
%   The conditional variance, h(t), of a TARCH(P,O,Q) process is modeled as follows:
%
%    g(h(t)) = omega
%            + alpha(1)*f(r_{t-1}) + ... + alpha(p)*f(r_{t-p})+...
%            + gamma(1)*I(t-1)*f(r_{t-1}) +...+ gamma(o)*I(t-o)*f(r_{t-o})+...
%            beta(1)*g(h(t-1)) +...+ beta(q)*g(h(t-q))
%
%     where f(x) = abs(x)  if tarch_type=1
%          g(x) = sqrt(x) if tarch_type=1
%          f(x) = x^2     if tarch_type=2
%          g(x) = x       if tarch_type=2
%
%   Default Options
%    options  =  optimset('fminunc');
%    options  =  optimset(options , 'TolFun'      , 1e-005);
%    options  =  optimset(options , 'TolX'        , 1e-005);
%    options  =  optimset(options , 'Display'     , 'iter');
%    options  =  optimset(options , 'Diagnostics' , 'on');
%    options  =  optimset(options , 'LargeScale'  , 'off');
%    options  =  optimset(options , 'MaxFunEvals' , '400*numberOfVariables');
%
%  See also TARCH_LIKELIHOOD, TARCH_CORE, TARCH_PARAMETER_CHECK, TARCH_STARTING_VALUES,
%  TARCH_TRANSFORM, TARCH_ITRANSFORM
%
%  You should use the MEX files (or compile if not using Win64 Matlab) as they provide speed ups of
%  approx 100 times relative to the m file.

% Copyright: Kevin Sheppard
% kevin.sheppard@economics.ox.ac.uk
% Revision: 3    Date: 9/1/2005

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
switch nargin
    case 4
        [p,o,q,error_type,tarch_type,startingvals,options]=tarch_parameter_check(epsilon, p, o, q);
    case 5
        [p,o,q,error_type,tarch_type,startingvals,options]=tarch_parameter_check(epsilon, p, o, q, error_type);
    case 6
        [p,o,q,error_type,tarch_type,startingvals,options]=tarch_parameter_check(epsilon, p, o, q, error_type, tarch_type);
    case 7
        [p,o,q,error_type,tarch_type,startingvals,options]=tarch_parameter_check(epsilon, p, o, q, error_type, tarch_type, startingvals);
    case 8
        [p,o,q,error_type,tarch_type,startingvals,options]=tarch_parameter_check(epsilon, p, o, q, error_type, tarch_type, startingvals, options);
    otherwise
        error('Number of inputs must be between 4 and 8');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Initial setup
m  =  max([p o q]);

%Augment the data with back casts to avoid costly memory allocations
if tarch_type==1
    %fepsilon is f(epsilon), as above
    fepsilon   =  [mean(abs(epsilon))*ones(m,1) ; abs(epsilon)];
    %fIepsilon is fepsilon*(epsilon<0)
    fIepsilon   =  [0.5*mean(abs(epsilon))*ones(m,1) ; abs(epsilon).*(epsilon<0)];

    % Local back casting
    back_cast_length = max(floor(length(epsilon)^(1/2)),1);
    back_cast_weights = .05*(.9.^(0:back_cast_length ));
    back_cast_weights = back_cast_weights/sum(back_cast_weights);
    back_cast = back_cast_weights*(abs(epsilon(1:back_cast_length+1)));
    if back_cast==0
        back_cast=mean(abs(epsilon));
    end
else
    %fepsilon is f(epsilon), as above
    fepsilon   =  [mean(epsilon.^2)*ones(m,1) ; epsilon.^2];
    %fIepsilon is fepsilon*(epsilon<0)
    fIepsilon   =  [0.5*mean(epsilon.^2)*ones(m,1) ; epsilon.^2.*(epsilon<0)];
    % Local back casting
    back_cast_length = max(floor(length(epsilon)^(1/2)),1);
    back_cast_weights = .05*(.9.^(0:back_cast_length ));
    back_cast_weights = back_cast_weights/sum(back_cast_weights);
    back_cast = back_cast_weights*((epsilon(1:back_cast_length+1)).^2);
    if back_cast==0
        back_cast=mean(epsilon.^2);
    end
end
epsilon_augmented=[zeros(m,1);epsilon];
%Compute the length of the augmented epsilon
T     = size(fepsilon,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Starting values
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%This flag is for the robustness check below, if the user supplies starting
%values, there will be no robustness check
if isempty(startingvals)
    startingflag=0;
else
    startingflag=1;
end
%Grid search for starting values.
[startingvals,nu,lambda,~,ordered_parameters]=tarch_starting_values(startingvals,epsilon_augmented,fepsilon,fIepsilon,p,o,q,T,error_type,tarch_type);
%Finally, initialize the starting values
startingvals = [startingvals; nu; lambda];
%Transform the starting vals
[garch_params_transformed,nu_transformed,lambda_transformed]=tarch_transform(startingvals,p,o,q,error_type);
%Re-append nu, lambda
startingvals_transformed = [garch_params_transformed; nu_transformed; lambda_transformed];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Starting values
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Estimate the parameters. Note the 1 in the last argument to indicate it
% is a constrained optimization

%LL0 is used to make sure the log likelihood improves
LL0=tarch_likelihood(startingvals_transformed,epsilon_augmented,fepsilon,fIepsilon,p,o,q,error_type,tarch_type,back_cast,T,1);
%Parameter estimation
[parameters,LL,exitflag,output]=fminunc('tarch_likelihood',startingvals_transformed,options,epsilon_augmented,fepsilon,fIepsilon,p,o,q,error_type,tarch_type,back_cast,T,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Estimation Robustness
%This portion of the code is to make sure that the optimization converged
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%This is the case where the optimization did not converge, but improved on
%the initial log likelihood
if  exitflag<=0 && LL<LL0
    %Try more iterations, only do more iterations if the final likelihood is
    %actually better than the initial

    %Increase the max iterations and max fun evals
    %Also switch to steepest descent
    if ischar(options.MaxFunEvals)
        options.MaxIter=2*100*length(parameters);
    else
        options.MaxIter=2*options.MaxIter;
    end
    if ischar(options.MaxFunEvals)
        options.MaxFunEvals=4*100*length(parameters);
    else
        options.MaxFunEvals=2*options.MaxFunEvals;
    end
    options.HessUpdate='steepdesc';
    % Estimate the parameters.
    [parameters,LL,exitflag,output]=fminunc('tarch_likelihood',parameters,options,epsilon_augmented,fepsilon,fIepsilon,p,o,q,error_type,tarch_type,back_cast,T,1);
end

%If the optimization still hasn't converged, try other starting values
if startingflag==0 && exitflag<=0
    %Keep track of the final estimates, if nothing converges, we will
    %return the
    robust_parameters(1,:)=parameters';
    %Also keep the LL
    robust_LL = zeros(2,1);
    robust_LL(1) = LL;
    %Keep track of the iteration
    index=2;
    while exitflag<=0
        %This condition checks that we haven't converged
        %OR that the best objective is worse than best grid search
        %Sort the original grid search log likelihoods and parameters

        startingvals=[ordered_parameters(index,:)' ; nu; lambda];
        %Transform the starting vals
        [garch_params_transformed,nu_transformed,lambda_transformed]=tarch_transform(startingvals,p,o,q,error_type);
        %Reappend nu
        startingvals_transformed = [garch_params_transformed; nu_transformed; lambda_transformed];

        LL0=tarch_likelihood(startingvals_transformed,epsilon_augmented,fepsilon,fIepsilon,p,o,q,error_type,tarch_type,back_cast,T,1);
        options.HessUpdate='bfgs';
        %Try the second set of starting values
        [parameters,LL,exitflag,output]=fminunc('tarch_likelihood',startingvals_transformed,options,epsilon_augmented,fepsilon,fIepsilon,p,o,q,error_type,tarch_type,back_cast,T,1);
        if  exitflag<=0 && LL<LL0
            %Again, if the LL improved, try more iterations
            %Increase the max iterations and max fun evals
            %Also switch to steepest descent
            options.MaxIter=2*options.MaxIter;
            options.MaxFunEvals=2*options.MaxFunEvals;
            options.HessUpdate='steepdesc';
            % Estimate the parameters.
            [parameters,LL,exitflag,output]=fminunc('tarch_likelihood',parameters,options,epsilon_augmented,fepsilon,fIepsilon,p,o,q,error_type,tarch_type,back_cast,T,1);
        end
        %Save the parameter estimates
        robust_parameters(index,:)=parameters';
        robust_LL(index)=LL;
        %Increment the index
        index=index+1;

        if index>size(ordered_parameters,1);
            %save the best LL and parameters and break
            warning('MFEToolbox:Convergence','Convergence not achieved.  Use results with caution');
            [LL,index]=min(robust_LL);
            parameters=robust_parameters(index,:)';
            break
        end
    end
end

%Transform the parameters from the real line to the restricted space
[parameters,nu,lambda]=tarch_itransform(parameters,p,o,q,error_type);
parameters=[parameters;nu;lambda];
%Compute the log likelihood if needed
if nargout>1
    [LL, ~, ht]=tarch_likelihood(parameters,epsilon_augmented,fepsilon,fIepsilon,p,o,q,error_type,tarch_type,back_cast,T);
    LL=-LL;
end

%Compute standard errors using RobustVCV if needed.
if nargout>3
    nw=0; %No newey west on scores
    [VCVrobust,A,~,scores,hess]=robustvcv('tarch_likelihood',parameters,nw,epsilon_augmented,fepsilon,fIepsilon,p,o,q,error_type,tarch_type,back_cast,T);
    VCV=hess^(-1)/(T-m);
    diagnostics.A = A;
end

%Report diagnostics in case requested
diagnostics.EXITFLAG=exitflag;
diagnostics.ITERATIONS=output.iterations;
diagnostics.FUNCCOUNT=output.funcCount;
diagnostics.MESSAGE=output.message;
diagnostics.m = m;
diagnostics.T = T;
diagnostics.fdata = fepsilon;
diagnostics.fIdata = fIepsilon;
diagnostics.back_cast = back_cast;

author: bashtage

repo: mfe-toolbox


text:
% test for reducing the TV norm
name = 'lena';
n = 256;
M = load_image(name);
M = rescale(crop(M, n));

t = compute_total_variation(M);

tvtgt = t*.1;
niter = 400;

options.tvtgt = tvtgt;
options.niter = niter;
options.nrefresh = 200;
[Mtv,err,tv,lalist] = perform_tv_denoising(M,options);

L = 1/2*err.^2+lalist(end)*tv;
clf;
subplot(2,1,1);
plot(tv-tvtgt, '.-'); title('TV-TV_0'); axis tight;
subplot(2,1,2);
plot(L, '.-'); title('Lagrangian'); axis tight;

author: gpeyre

repo: matlab-toolboxes


text:
function plot_detector_wiggle_voc07()

  colormap = [ ...
    83, 136, 173 ; ...
    219, 135, 45 ; ...
    218, 71, 56 ; ...
    135, 130, 174 ; ...
    142, 195, 129 ; ...
    138, 180, 66 ; ...
    223, 200, 51 ; ...
    92, 172, 158 ; ...
    ] ./ 256;

  config = [];
  config(1).name = 'R-CNN';
  config(1).func = @detector_wiggle_rcnn;
  config(1).filename = '/BS/candidate_detections/work/v1/pascal_voc07_test_rcnn_scoremaps-3d.mat';
  config(2).name = 'LM-LLDA';
  config(2).func = @detector_wiggle_rcnn;
  config(2).filename = '/BS/candidate_detections/work/v1/pascal_voc07_test_dpm_scoremaps-3d.mat';
  config(3).name = 'LM-LLDA bboxpred';
  config(3).func = @detector_wiggle_rcnn;
  config(3).filename = '/BS/candidate_detections/work/v1/pascal_voc07_test_dpm_bboxpred_scoremaps-3d.mat';
  config(4).name = 'Fast R-CNN';
  config(4).func = @detector_wiggle_from_iou_score;
  config(4).filename = '/BS/candidate_detections/work/v1/pascal_voc07_test_frcn_iou_score.mat';
  config(5).name = 'Fast R-CNN bboxpred';
  config(5).func = @detector_wiggle_from_iou_score;
  config(5).filename = '/BS/candidate_detections/work/v1/pascal_voc07_test_frcn_regr_iou_score.mat';

  autoscale = true;

  bin_edges = 0.2:0.05:1;
  bin_centers = (bin_edges(1:end-1) + bin_edges(2:end)) / 2;

  hists = build_histograms(config, bin_edges);
  if autoscale
    hists = autoscale_per_class(hists);
  end

  plot_histograms(config, bin_centers, bin_edges, hists.avg_scores, hists.class_avg_scores, colormap);
%   plot_histograms(config, bin_centers, bin_edges, hists.median_scores, hists.class_median_scores, colormap);

%   figure;
%   bar(bin_centers', hists.bin_counts');
%   legend({config.name}, 'Location', 'SouthEast');

  if false
    scale_factors = [0.8182    0.8485    0.8864    0.9242    0.9621    1.0000    1.0455    1.0833    1.1288  1.1742    1.2273];
    x = (-5:5) .* 4*227/128;
    y = x;
    ticks = [-5 -2 0 2 5] .* 4*227/128;
    tick_labels = arrayfun(@(t) sprintf('%.0f', t), ticks(:), 'UniformOutput', false);
    method_i = 1;
    n_steps = size(hists.avg_cubes{method_i},1);
    class_cubes = hists.method_class_avg_cubes(method_i,:);
    n_classes = numel(class_cubes);
    % normalize per class
    for c = 1:n_classes
      clims = [min(class_cubes{c}(:)), max(class_cubes{c}(:))];
      class_cubes{c} = (class_cubes{c} - clims(1)) ./ (clims(2) - clims(1));
    end
    cube = nanmean(cat(4, class_cubes{:}), 4);

    clims = [min(cube(:)), max(cube(:))];

    figure;
    scale_idxs = [1 3 6 9 11];
    for j = 1:numel(scale_idxs), i = scale_idxs(j); % 1:n_steps
      subplot(1,6,j);
      imagesc(x, y, cube(:,:,i), clims);
      set(gca, 'XTick', ticks, 'XTickLabels', tick_labels);
      set(gca, 'YTick', ticks, 'YTickLabels', tick_labels);
      axis square
      title(sprintf('scale %.2f', scale_factors(i)));
    end
    subplot(1,6,6);
    axis off
    cbh = colorbar();
    P = get(cbh, 'Position');
    set(cbh, 'Position', [P(1)-0.01, P(2)+0.1, P(3)+0.01, P(4)-0.2]);
    hei = 3.5; wid = 25;
    set(gcf, 'Units','centimeters', 'Position',[0 0 wid hei]);
    set(gcf, 'PaperPositionMode','auto');
    printpdf('figures/detector_scoremaps_pos_scale_voc07.pdf');
    printpng('figures/detector_scoremaps_pos_scale_voc07.png');

%     figure;
%     for i = 1:n_steps
%       subplot(4,3,i);
%       imagesc(reshape(hists.avg_cubes{method_i}(:,i,:), [n_steps n_steps]), clims);
%       axis equal
%       title(sprintf('x jitter %d', i - floor((n_steps+1)/2)));
%     end
%
%     figure;
%     for i = 1:n_steps
%       subplot(4,3,i);
%       imagesc(reshape(hists.avg_cubes{method_i}(i,:,:), [n_steps n_steps]), clims);
%       axis equal
%       title(sprintf('y jitter %d', i - floor((n_steps+1)/2)));
%     end
  end
end

function plot_histograms(config, bin_centers, bin_edges, avg_scores, class_avg_scores, colormap)
  ld = load('data/classes.mat');
  classes = ld.classes;

  fh_all_classes = figure; hold on;
  for i = 1:numel(config)
    plot(bin_centers, avg_scores(i,:), '.-', 'LineWidth', 1.5, 'Color', colormap(i,:), 'MarkerSize', 20);
  end
  legend({config.name}, 'Location', 'SouthEast');
  xlabel('IoU with GT'); ylabel('detector score');
  xlim([min(bin_edges), max(bin_edges)]);
  hei = 8; wid = 12;
  set(gcf, 'Units','centimeters', 'Position',[0 0 wid hei]);
  set(gcf, 'PaperPositionMode','auto');
  printpdf('figures/detector_scoremaps_voc07.pdf');

  fh_per_class = figure;
  for c = 1:20
    for i = 1:numel(config)
      figure(fh_per_class); subplot(4,5,c);
      plot(bin_centers, class_avg_scores{c}(i,:), '-o', 'LineWidth', 2, 'Color', colormap(i,:));
      hold on;
      xlim([min(bin_edges), max(bin_edges)]);
      title(classes{c});
    end
  end
end

% function [avg_scores, bin_counts, class_avg_scores, class_bin_counts, avg_cubes] = build_histograms(config, bin_edges)
function hists = build_histograms(config, bin_edges)
  histogram_cache_file = 'data/score_histogram_cache.mat';
  try
    hists = load(histogram_cache_file);
  catch
    n_methods = numel(config);
    n_classes = 20;
    n_bins = numel(bin_edges)-1;

    % allocate
    avg_scores = zeros(n_methods, n_bins);
    median_scores = zeros(n_methods, n_bins);
    bin_counts = zeros(n_methods, n_bins);
    class_avg_scores = cell(1, n_classes);
    class_median_scores = cell(1, n_classes);
    class_bin_counts = cell(1, n_classes);
    for c = 1:n_classes
      class_avg_scores{c} = zeros(n_methods, n_bins);
      class_median_scores{c} = zeros(n_methods, n_bins);
      class_bin_counts{c} = zeros(n_methods, n_bins);
    end
    avg_cubes = cell(n_methods,1);
    method_class_avg_cubes = cell(n_methods,n_classes);

    % compute stuff
    for i = 1:n_methods, conf = config(i);
      fprintf('computing %d/%d %s\n', i, numel(config), conf.name);
      [overlap, score, class_ids, avg_cubes{i}, class_avg_cubes] = conf.func(conf.filename);
      if ~isempty(class_avg_cubes)
        method_class_avg_cubes(i,:) = class_avg_cubes;
      end

      [bin_counts(i,:), avg_scores(i,:), median_scores(i,:)] = build_histogram(overlap, score, bin_edges);

      fprintf('per class: ');
      for c = 1:n_classes
        fprintf('%d ', c);
        class_mask = class_ids == c;
        class_overlap = overlap(class_mask);
        class_score = score(class_mask);
        [class_bin_counts{c}(i,:), class_avg_scores{c}(i,:), class_median_scores{c}(i,:)] = build_histogram(class_overlap, class_score, bin_edges);
      end
      fprintf('\n');
    end

    save(histogram_cache_file, 'avg_scores', 'median_scores', 'bin_counts', ...
      'class_avg_scores', 'class_median_scores', 'class_bin_counts', ...
      'avg_cubes', 'method_class_avg_cubes');
    hists = load(histogram_cache_file);
  end
end

function [hists] = autoscale_per_class(hists)
  n_classes = 20;
  n_methods = size(hists.avg_scores,1);
  n_bins = size(hists.avg_scores,2);

  denom = zeros(n_methods,n_bins);
  hists.avg_scores = zeros(n_methods,n_bins);
  hists.median_scores = zeros(n_methods,n_bins);
  for c = 1:n_classes
    hists.class_avg_scores{c} = autoscale_scores(hists.class_avg_scores{c});
    hists.class_median_scores{c} = autoscale_scores(hists.class_median_scores{c});

    hists.avg_scores = hists.avg_scores + hists.class_avg_scores{c} .* hists.class_bin_counts{c};
    hists.median_scores = hists.median_scores + hists.class_median_scores{c} .* hists.class_bin_counts{c};
    denom = denom + hists.class_bin_counts{c};
  end
  hists.avg_scores = hists.avg_scores ./ denom;
  hists.median_scores = hists.median_scores ./ denom;
  assert(all(denom(:) == hists.bin_counts(:)));
end

function scaled_scores = autoscale_scores(avg_scores)
  minscore = min(avg_scores, [], 2);
  maxscore = max(avg_scores, [], 2);
  minscore = repmat(minscore, [1, size(avg_scores, 2)]);
  maxscore = repmat(maxscore, [1, size(avg_scores, 2)]);
  scaled_scores = (avg_scores - minscore) ./ (maxscore - minscore);
end

function [counts, avg_scores, median_scores] = build_histogram(overlap, score, bin_edges)
%     [counts,~,bin_idxs] = histcounts(overlap, bin_edges);
  [counts,bin_idxs] = histc(overlap, bin_edges);
  counts = counts(1:end-1);

  n_centers = numel(bin_edges)-1;
  avg_scores = zeros(1, n_centers);
  median_scores = zeros(1, n_centers);
  % average over the elements in each bin
  for j = 1:n_centers
    bin_scores = score(bin_idxs == j);
    if ~isempty(bin_scores)
      avg_scores(j) = mean(bin_scores);
      median_scores(j) = median(bin_scores);
    end
  end
end

author: hosang

repo: detection-proposals


text:
function [P,P2,d] = get_filename2(dwcard,varargin)
% :Usage:
% ::
%
%     function [P,P2,d] = get_filename2(search string (as with ls command),[verbose])
%
% Start in directory above individual subject directories
%
% :Inputs:
%
%   **dwcard:**
%        Enter dwcard for wildcard of directories to look in; e.g. '02*'
%
%   **wcard:**
%        Enter wcard for image or file to get - e.g., 'beta_0001.img'
%        This can also include subdirectories (no *'s) 'anatomy/beta*img'
%
% :Outputs:
%
%     Returns list of all files in individual directories  in string matrix
%
%   **P:**
%        file names with full paths
%
%   **P2:**
%        file names only
%
%   **d:**
%        list of directories searched for files
%
% Missing files, or entries in directory that do not contain files, are
% removed from the list.
% NOT entering a * in dwcard seems to produce an error.
%
% :Examples:
% ::
%
%    P = get_filename2('02*/beta_0001*')
%
% one * is allowed in the directory structure right now,
% multiple *'s in the filename.
%
% ..
%    Tor Wager 5/2/03
% ..

if length(varargin) > 0, verb = varargin{1};, else, verb = 0;,end

clear d, clear D

P = []; P2 = [];

[dd,fstr,ext] = fileparts(dwcard);

% get list of directories to search

if any(dd == '*')
    d = dir(dd);

    [CWD,newd] = fileparts(dd);

    if ~isempty(d),
        for i = 1:length(d), d(i).name = fullfile(CWD,d(i).name);,end
    end

    pdname = newd;
    while any(CWD == '*')
        d = [d dir(CWD)];

        if ~isempty(d),
            for i = 1:length(d), d(i).name = fullfile(fileparts(CWD),d(i).name,pdname);,end
        end

        [CWD,newd] = fileparts(newd);
        pdname = fullfile(pdname,newd);
        if verb, disp(['Looking in ' CWD]);,end
    end

    if verb, disp('List of dirs to search:'),str2mat(d.name),disp(' '),end
else
    clear d,
    d(1).name = dd; if isdir(d(1).name), d(1).isdir = 1;, else, d(1).isdir = 0;,end
end

% search for files in dirs

for i = 1:length(d),
	if d(i).isdir,
		D = fullfile(d(i).name,[fstr ext]);,

        if verb, str = ['dir(''' D ''')'];, disp(str),end
        f = dir(D);

        % file name only
        if ~isempty(f), P2 = str2mat(P,str2mat(f(~cat(1,f.isdir)).name));,end

        % full path name
        for j = 1:length(f)
            if ~f(j).isdir,
                if isempty(P)
                    P = fullfile(fileparts(D),f(j).name);
                else
                    P = str2mat(P,fullfile(fileparts(D),f(j).name));,
                end
            end
        end
    end
end

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
%GET Get datafile parameter fields
%
%   [VALUE1,VALUE2,...] = GET(A,FIELD1,FIELD2,...)
%
% INPUT
%   A       Datafile
%   FIELDx  Field names (strings)
%
% OUTPUT
%   VALUEx  Field values
%
% SEE ALSO (PRTools Guide)
% PRDATAFILE, SET

author: marianux

repo: ecg-kit

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


text:
function varargout = parseMeshData(varargin)
%PARSEMESHDATA Conversion of data representation for meshes.
%
%   MESH = parseMeshData(VERTICES, FACES)
%   MESH = parseMeshData(VERTICES, EDGES, FACES)
%   Returns the mesh info into a single structure with fields "vertices",
%   "edges" and "faces".
%
%   [VERTICES, FACES] = parseMeshData(MESH)
%   Returns the mesh info into two output variables containing coordinates
%   of vertices (as a Nv_by_3 array) and the list of vertex indices for
%   each face (either as a Nf-by-3 or Nf-by-4 int array, or as a cell
%   array).
%
%   [VERTICES, EDGES, FACES] = parseMeshData(MESH)
%   Also returns the vertex indices of each edge, as a Ne-by-2 array of
%   vertex indices.
%
%   See also
%   meshes3d, formatMeshOutput
%

% ------
% Author: David Legland
% E-mail: david.legland@inra.fr
% Created: 2010-12-06, using Matlab 7.9.0.529 (R2009b)
% Copyright 2010-2022 INRA - Cepia Software Platform

% initialize edges to empty variable
edges = [];

% Process input arguments
switch nargin
    case 1
        % input is a data structure
        mesh = varargin{1};
        vertices = mesh.vertices;
        faces = mesh.faces;
        if isfield(mesh, 'edges')
            edges = mesh.edges;
        end

    case 2
        % input are vertices and faces
        vertices = varargin{1};
        faces = varargin{2};

    case  3
        % input are vertices, edges and faces
        vertices = varargin{1};
        edges = varargin{2};
        faces = varargin{3};

    otherwise
        error('Wrong number of arguments');
end

% returns either a struct or several variables
varargout = formatMeshOutput(nargout, vertices, edges, faces);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function pass = test_linearScalarODEs(pref)
% A linear CHEBOP test. This test tests a scalar ODE, both with and without
% breakpoints, as well as with discontinuous coefficients. It solves the
% problems using chebcolloc1, chebcolloc2 and ultraS discretizations.

%% Setup
dom = [0 pi];
if ( nargin == 0 )
    pref = cheboppref;
end
%pref.bvpTol = 1e-11;

%% Simple scalar problem
N = chebop(@(x,u) diff(u,2) + x.*u, dom);
N.lbc = 2;
N.rbc = 3;

x = chebfun(@(x) x, dom);
rhs = sin(x);

%% Try different discretizations
% Start with chebcolloc2
pref.discretization = @chebcolloc2;
u1 = solvebvp(N, rhs, pref);

%% Change to ultraS
pref.discretization = @ultraS;
u2 = solvebvp(N, rhs, pref);

%% Change to chebcolloc1
pref.discretization = @chebcolloc1;
u3 = solvebvp(N, rhs, pref);

%% Did we pass?
% To pass, both residuals have to be small, but we should not expect u1 and u2
% to be identical!
tol = 1e3*pref.bvpTol;
pass(1) = norm(N(u1)-rhs) < tol && ( u1(0) - 2 < tol) && ( u1(pi) - 3 < tol);
pass(2) = norm(N(u2)-rhs) < tol && ( u2(0) - 2 < tol) && ( u2(pi) - 3 < tol);
pass(3) = norm(N(u3)-rhs) < tol && ( u3(0) - 2 < tol) && ( u3(pi) - 3 < tol);
pass(4) = ( (norm(u1 - u2) ~= 0) && (norm(u2 - u3) ~= 0) && ...
    (norm(u1 - u3) ~= 0));

%% Problem with breakpoints
dom = [-1 0 pi];
N = chebop(@(x,u) diff(u,2) + cos(x).*u, dom);
N.lbc = 2;
N.rbc = -1;

x = chebfun(@(x) x, dom);
rhs = sin(x);

%% Try different discretizations
% Start with chebcolloc2
pref.discretization = @chebcolloc2;
u4 = solvebvp(N, rhs, pref);

%% Change to ultraS
pref.discretization = @ultraS;
u5 = solvebvp(N, rhs, pref);

%% Change to chebcolloc1
pref.discretization = @chebcolloc1;
u6 = solvebvp(N, rhs, pref);

%% Did we pass?
% To pass, both residuals have to be small, but we should not expect u3 and u4
% to be identical!
tol = 1e3*pref.bvpTol;
pass(5) = norm(N(u4)-rhs) < tol && ( u4(-1) - 2 < tol) && ( u4(pi) + 1 < tol);
pass(6) = norm(N(u5)-rhs) < tol && ( u5(-1) - 2 < tol) && ( u5(pi) + 1 < tol);
pass(7) = norm(N(u6)-rhs) < tol && ( u6(-1) - 2 < tol) && ( u6(pi) + 1 < tol);
pass(8) = norm(jump(u4, 0)) < tol && norm(jump(u5, 0)) < tol && ...
    norm(jump(u6, 0)) < tol;
pass(9) = ( (norm(u4 - u5) ~= 0) && (norm(u4 - u6) ~= 0) && ...
    (norm(u5 - u6) ~= 0));

%% Problem with discontinuous coefficients
dom = [-1 0 pi];
N = chebop(@(x,u) diff(u,2) + abs(x-1).*u, dom);
N.lbc = 2;
N.rbc = -1;

x = chebfun(@(x) x, dom);
rhs = sin(x);

%% Try different discretizations
% Start with collocation
pref.discretization = @chebcolloc2;
u7 = solvebvp(N, rhs, pref);

%% Change to ultraS
pref.discretization = @ultraS;
u8 = solvebvp(N, rhs, pref);

%% Change to chebcolloc1
pref.discretization = @chebcolloc1;
u9 = solvebvp(N, rhs, pref);

%% Did we pass?
% To pass, both residuals have to be small, but we should not expect u3 and u4
% to be identical!
tol = 1e1*pref.bvpTol;
pass(10) = norm(N(u7)-rhs) < 50*tol && ...
    ( u7(-1) - 2 < tol) && ( u7(pi) + 1 < tol);
pass(11) = norm(N(u8)-rhs) < 10*tol && ...
    ( u8(-1) - 2 < tol) && ( u8(pi) + 1 < tol);
pass(12) = norm(N(u9)-rhs) < 10*tol && ...
    ( u9(-1) - 2 < tol) && ( u9(pi) + 1 < tol);
pass(13) = norm(jump(u7, 0)) < tol && norm(jump(u8, 0)) < tol && ...
    norm(jump(u9, 0)) < tol;
pass(14) = ( norm(u5-u6) ~= 0 );

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function t = tvec_even_bracket2 ( nt, theta1, theta2 )

%*****************************************************************************80
%
%% TVEC_EVEN_BRACKET2 computes an evenly spaced set of angles between THETA1 and THETA2.
%
%  Example:
%
%    NT = 5
%    THETA1 = 30
%    THETA2 = 90
%
%    T = ( 40, 50, 60, 70, 80 )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 May 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NT, the number of values to compute.
%
%    Input, real THETA1, THETA2, the limiting angles.
%
%    Output, real TVEC(NT), the evenly spaced angles.
%
  for i = 1 : nt
    t(i) = ( ( nt + 1 - i ) * theta1   ...
           + (          i ) * theta2 ) ...
           / ( nt + 1     );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function directories_added=moxunit_set_path()
% sets the search path for MOxUnit
%
% directories_added=moxunit_set_path()
%
% Output:
%    directories_added          Cell string with directories that were
%                               added to the search path. If no
%                               directories were added to the search path
%                               (i.e., all functions for MOxUnit were
%                               already in the search path), then this is
%                               an empty cell array.
% Note:
%   - This adds directory in which this function resides, as well as the
%     'util' sub-directory to the search path

    return_directories_added=nargout>=1;

    root_dir=fileparts(mfilename('fullpath'));

    sub_dirs={'','util'};

    % Cannot use moxunit_util_strjoin because it may not be in the path.
    % Instead construct the path to be added by adding the pathsep() after
    % each subdirectory
    make_full_path_with_sep=@(sub_dir)[fullfile(root_dir,sub_dir) ...
                                        pathsep()];

    full_dirs=cellfun(make_full_path_with_sep,sub_dirs,...
                        'UniformOutput',false);

    if return_directories_added
        % store original path elements
        orig_path_cell=get_path_cell();
    end

    addpath(cat(2,full_dirs{:}));

    if return_directories_added
        new_path_cell=get_path_cell();
        directories_added=setdiff(new_path_cell, orig_path_cell);
    end

function pc=get_path_cell()
    pc=regexp(path(), pathsep(), 'split');

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function [X,rho,eta] = prrgmres(A,L,N,b,k)
%PRRGMRES Preconditioned RRGMRES algorithm for square inconsistent systems
%
% [X,rho,eta] = rrgmres(A,L,N,b,k)
%
% PRRGMRES applies smoothing-norm preconditioning to the RRGMRES method,
% which is a variant of RRGMRES for square linear systems A x = b, with
% starting vector A*b (instead of b as in GMRES).  This function returns
% all k iterates, stored as columns of the matrix X.  The solution norm
% and residual norm are returned in eta and rho, resp.
%
% The preconditioner uses two matrices: the matrix L that defines the
% smoothing norm, and the matrix N whose columns span the null space
% of L.  It is assumed that L is p-times-n with p < n.
%
% For symmetric matrices, use the function pmr2 instead.

% Reference: P. C. Hansen and T. K. Jensen, "Smoothing-norm preconditioning
% for regularizing minimum-residual methods", SIAM J. Matrix Anal. Appl.
% 29 (2006), 1-14.

% Per Christian Hansen, IMM, September 21, 2007.

% Check input arguments.
if (k < 1), error('Number of steps k must be positive'), end
[m,n] = size(A);
p = size(L,1);
if (m ~= n), error('A must be square'), end
if nargin < 5, error('Too few input arguments'), end

% Allocate space.
X = zeros(n,k); % Matrix of solutions.
V = zeros(p,k); % Orthonormal vectors spanning Krylov subspace
h = zeros(k,1); % New column of Hessenberg matrix
Q = zeros(k+1); % H = Q*T, Q orthogonal
T = zeros(k);   %          T upper triangular
W = zeros(p,k); % W = V*inv(T)
if (nargout>1), rho = zeros(k,1); end
if (nargout>2), eta = zeros(k,1); end

% Initialization for working with pseudoinverses of L.
[Q0,R0] = qr(A*N,0);  % Compate QR factorization of A*N.
T0 = N'*Q0;
[QL,RL] = qr(L',0);   % Compact QR factgorization of L'.
TN = pinit(N,A);      % Prepare for A-weighted pseudoinverse computations.
bb = RL\( QL'*( b - Q0*(T0\(N'*b)) ) );

% Initialize variables.
v1     = A*( QL*(RL'\bb) );
v2     = Q0*( T0\(N'*v1) );
r      = RL\( QL'*(v1-v2) );
alpha  = norm(r);
V(:,1) = r/alpha;  % Initial vector of Krylov subspace is A*b.
Q(1,1) = 1;
x0     = N*( R0\(Q0'*b) );
xi     = zeros(p,1);
beta   = V(:,1)'*bb;

% Begin iterations.
for i=1:k

  v1 = A*( QL*(RL'\V(:,i)) );
  v2 = Q0*( T0\(N'*v1) );
  r  = RL\( QL'*(v1-v2) );

  % Modified Gram-Schmidt on the new vector.
  for j=1:k
    h(j) = V(:,j)'*r;
    r = r - V(:,j)*h(j);
  end
  alpha = norm(r);

  % Store new Arnoldi vector and update projected rhs.
  V(:,i+1) = r/alpha;
  beta = [beta;V(:,i+1)'*bb];

  % Apply previous rotations to h.
  T(1:i,i) = Q(1:i,1:i)'*h(1:i);

  % Compute Givens rotation parameters.
  rc = T(i,i);
  if alpha == 0
    c = 1; s = 0;
  elseif abs(alpha) > abs(rc)
    tau = -rc/alpha;
    s = 1 / sqrt(1 + abs(tau)^2);
    c = s*tau;
  else
    tau = -alpha/rc;
    c = 1 / sqrt(1 + abs(tau)^2);
    s = c*tau;
  end

  % Apply givens rotations.
  T(i,i) = c'*rc - s'*alpha;
  Q(1:i,[i,i+1]) = Q(1:i,i)*[c s];
  Q(i+1,[i,i+1]) = [-s c];

  if abs(T(i,i)) <= eps
    disp('Hession matrix is (numerically) singular')
  end

  % Update W = V*inv(T);
  W(:,i) = (V(:,i) - W(:,1:i-1)*T(1:i-1,i))/T(i,i);

  % Update solution.
  xi = xi + (Q(1:i+1,i)'*beta)*W(:,i);

  % Update output variables.
  X(:,i) = lsolve(L,xi,N,TN) + x0;
  if nargout>1, rho(i) = norm(A*X(:,i)-b); end
  if nargout>2, eta(i) = norm(xi); end

end

repo: REST


text:
function t = tvec_even2 ( nt )

%*****************************************************************************80
%
%% TVEC_EVEN2 computes an evenly spaced set of angles between 0 and 2*PI.
%
%  Discussion:
%
%    The computation realizes that 0 = 2 * PI.  The values are equally
%    spaced in the circle, do not include 0, and are symmetric about 0.
%
%  Example:
%
%    NT = 4
%
%    T = ( PI/4, 3*PI/4, 5*PI/4, 7*PI/4 )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 May 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NT, the number of values to compute.
%
%    Output, real TVEC(NT), the evenly spaced angles, in radians.
%
  t(1:nt) = ( 1:2:(2*nt-1) ) * pi / nt;

  return
end


text:
function Xptm_run = fast_polytrendmtx(run,ntrs,nruns,order)
% Xptm = fast_polytrendmtx(run,ntrs,nruns,order)
%

if(nargin ~= 4)
  msg = 'USAGE: Xptm = fast_polytrendmtx(run,ntrs,nruns,order)';
  qoe(msg);error(msg);
end

Xptm = ones(ntrs,1);
t = [0:ntrs-1]'; %'
for n = 1:order
  r0 = t.^n;
  M = eye(ntrs) - Xptm*inv(Xptm'*Xptm)*Xptm';
  r = M*r0;
  r = r/std(r);
  Xptm = [Xptm r];
end

Xptm_run = zeros(ntrs,nruns*(order+1));
n1 = (run-1)*(order+1) + 1;
n2 = n1 + order;
Xptm_run(:,n1:n2) = Xptm;

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function y = choosvd( n, d)

if n <= 100
    if d / n <= 0.02
        y = 1;
    else
        y = 0;
    end
elseif n <= 200
    if d / n <= 0.06
        y = 1;
    else
        y = 0;
    end
elseif n <= 300
    if d / n <= 0.26
        y = 1;
    else
        y = 0;
    end
elseif n <= 400
    if d / n <= 0.28
        y = 1;
    else
        y = 0;
    end
elseif n <= 500
    if d / n <= 0.34
        y = 1;
    else
        y = 0;
    end
else
    if d / n <= 0.38
        y = 1;
    else
        y = 0;
    end
end

author: andrewssobral

repo: lrslibrary


text:
function K = knGauss(X, Y, s)
% Gaussian (RBF) kernel K = exp(-|x-y|/(2s));
% Input:
%   X: d x nx data matrix
%   Y: d x ny data matrix
%   s: sigma of gaussian
% Ouput:
%   K: nx x ny kernel matrix
% Written by Mo Chen (sth4nth@gmail.com).
if nargin < 3
    s = 0.4;
end

if nargin < 2 || isempty(Y)
    K = ones(1,size(X,2));            % norm in kernel space
else
    D = bsxfun(@plus,dot(X,X,1)',dot(Y,Y,1))-2*(X'*Y);
    K = exp(D/(-2*s^2));
end

author: PRML

repo: PRMLT

save_path: github-repos/MATLAB/PRML-PRMLT


text:
%KDOG Difference of Gaussian kernel
%
% K = KDOG(SIGMA1) is a 2-dimensional difference of Gaussian kernel equal
% to KGAUSS(SIGMA1) - KGAUSS(SIGMA2), where SIGMA1 > SIGMA2.  By default
% SIGMA2 = 1.6*SIGMA1.  The kernel is centred within the matrix K whose
% half-width H = 3xSIGMA and W=2xH+1.
%
% K = KDOG(SIGMA1, SIGMA2) as above but SIGMA2 is specified directly.
%
% K = KDOG(SIGMA1, SIGMA2, H) as above but the kernel half-width is specified.
%
% Notes::
% - This kernel is similar to the Laplacian of Gaussian and is often used
%   as an efficient approximation.
%
% See also KGAUSS, KDGAUSS, KLOG, ICONV.

% Copyright (C) 1993-2011, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% MVTB 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with MVTB.  If not, see <http://www.gnu.org/licenses/>.

function m = kdog(sigma1, sigma2, w)

    % sigma1 > sigma2
    if nargin == 1
        sigma2 = 1.6*sigma1;
        w = ceil(2*sigma2);
    elseif nargin == 2
        sigma2 = 1.6*sigma1;
        w = ceil(2*sigma2);
    elseif nargin == 3
        if sigma2 < sigma1
            t = sigma1;
            sigma1 = sigma2;
            sigma2 = t;
        end
    end
	if nargin < 3
        w = ceil(3*sigma1);
	end

    % sigma2 > sigma1
    m1 = kgauss(sigma1, w);     % thin kernel
    m2 = kgauss(sigma2, w);     % wide kernel

    m = m2 - m1;

author: petercorke

repo: machinevision-toolbox-matlab


text:
function outdata = luisFilter(data,  TR, cutoff, verbose)
% This is a low pass FIR filter using the Parks Mclellan design algorithm
% The phase introduced by the filter is linear and is un-done by
% filtering the data again, backwards
% if you want to see what it does to the data, use verbose=1
% if just want to filter the data, don't use the argument at all.
%
% :Usage:
% ::
%
%     function outdata = luisFilter(data,  TR, cutoff [, verbose])
%
% ..
%    Luis Hernandez. University of Michigan.  Last Edit 9/13/01
% ..

% ..
% close all
%
%%%%%%  These are the important lines of the code   %%%%%%%%%%%%%
% ..
    cutoff = cutoff / (2*TR)
    b = remez(10, [0 cutoff-0.05  cutoff+0.05  1], [1 1 0 0]);
    outdata = filtfilt(b,1, data);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if (nargin==4)
    freqz(b)
    figure

    % Read the data from file
%    data = load(data_file);
%    data = data(:,2);
    len = size(data, 1);
    time = [1:len] * TR;
    subplot(2,1,1) , plot(time, abs(data));
    xlabel('sec.')

    % Fourier Transform the data
    fdata = fftshift(fft(data));
    scalefactor = max(abs(fdata));

    %frequency range:
    f = [1:len]';
    f = (f - len/2 ) * 1/(TR * len);
    whos
    % Plotting the data before the filter:

    subplot(2,1,1) , hold on,plot(time, abs(data))
    subplot(2,1,2), plot(f, abs(fdata));
    subplot(2,1,2), axis([0  1/(2*TR)  0  scalefactor]);
    xlabel('Hz.');

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Make  Fourier Domain Filter by hand
    %ffilt = ones(len, 1);
    %co  = len/2  +  cutoff * (TR * len);
    %ffilt(co:end) = exp(-((f(co:end) - cutoff).^2) / 0.001 );
    %co  = len/2  - cutoff * (TR * len);
    %ffilt(1:co)   = exp( -((f(1:co) + cutoff).^2) / 0.001 );

    %subplot(2,1,2) , hold on , plot(f, ffilt*scalefactor/2 ,'g');
    %whos

    % Apply the filter to the data and IFT it
    %fdata  = fdata .* ffilt;
    %data = ifft(fdata);

     % Fourier Transform the data
    outfdata = fftshift(fft(outdata));

    % Let's look at the impulse response of the filter:
    impulse= zeros(size(data));

    impulse(size(data,1)/2) = 0.1 *scalefactor;
    impulse_response = filtfilt(b,1,impulse);

    freq_response = fftshift(fft(impulse_response));
    freq_impulse = fftshift(fft(impulse));

    freq_gain = scalefactor/2 + 20*log10( abs(freq_response).^2 ./ abs(freq_impulse) ) ;

    %Plotting the reponse of the filter
    subplot(2,1,1) , hold on,plot(time, abs(outdata), 'r'), title ('Time Domain'), legend('Before', 'After')
    subplot(2,1,2), hold on, plot(f, abs(outfdata), 'r'), title ('Frequency Domain');

    subplot(2,1,2), hold on, plot(f, freq_gain, 'g'),legend('Before','After','Filter Response');

    %figure;
    %plot(f, freq_gain, 'g');
end

return


text:
function w = sign(w)
   %SIGN Signum function for waveforms.
   % WAVEFORM = SIGN(WAVEFORM)
   % see sign

   % AUTHOR: Celso Reyes, Geophysical Institute, Univ. of Alaska Fairbanks
   % $Date$
   % $Revision$

   for n=1:numel(w)
      % w(n) = set(w(n),'data',sign(double(get(w(n),'data'))));
      w(n).data = sign(w(n).data);
   end
end

author: geoscience-community-codes

repo: GISMO


text:
function RoiFileName=dtiExportRoiToNifti(roi, refImage, outFile)
%
% dtiExportRoiToNifti(roi, ref, [outFile='roi'])
% Note:  refImage and outFile filenames: prefixes only
% Save a mrDiffusion ROI as a binary NIFTI image. Useful for
% loading the ROI into another program, like FSL or SPM.
%
% HISTORY:
% 2007.10.06 Kaustubh Supekar wrote it.

% Allow roi to be either a fileName or the roi struct.
if(ischar(roi))
   roi = dtiReadRoi(roi);
end

if(~exist('outFile','var')||isempty(outFile))
    outFile = 'roi';
end
l = length(outFile);
if(l<5||(~strcmpi(outFile(l-3:l),'.nii')&&~strcmpi(outFile(l-2:l),'.gz')))
    outFile = strcat(outFile,'.nii.gz');
end
% ref can either be a NIFTI filename or a NIFTI struct
if(ischar(refImage))
   l = length(refImage);
   if(l<5||(~strcmpi(refImage(l-3:l),'.nii')&&~strcmpi(refImage(l-2:l),'.gz')))
      refImage = [refImage '.nii.gz'];
   end
   ref = niftiRead(refImage);
end
c = mrAnatXformCoords(ref.qto_ijk, roi.coords);
c = round(c);
roiIm = ref;
roiIm.data = zeros(size(roiIm.data),'uint8');
roiIm.data(sub2ind(size(roiIm.data), c(:,1), c(:,2), c(:,3))) = 1;
roiIm.fname = outFile;
writeFileNifti(roiIm);
RoiFileName=roiIm.fname;
return;


text:
% Simple 1st order RC transient response:
% File:  rctime.m
% 9/19/02
%
clear;clc
%
K=1e3;uF=1e-6;us=1e-6;ms=1e-3;
%
R1=20*K;R2=40*K;C1=0.5*uF;Ein=1;E1=1;
%
U=1;N=1;M=1;Y=1;
%
A1=[1/R1+1/R2 1;1 0];
%
B2=[0 Ein/R1;E1 0];
%
P=C1;
%
% Template matrix equations:
%
V=A1\B2;H=V(U+1:U+N,1:N+M);AB=inv(P)*H;I=eye(N);

A=AB(1:N,1:N);B=AB(1:N,N+1:N+M);
D=V(Y:Y,1:N);E=V(Y:Y,N+1:N+M);
%
% Display A, B, D, & E
A
B
D
E
%
Rp=R1*R2/(R1+R2);
dt=500*us;kmax=100;Ein=4.5; % Ein can be changed after A, B, D, E have been computed.
%dt=100*us;kmax=500;Ein=4.5; % Try this dt after using the dt above.
Vc=zeros(1,kmax);Vo=zeros(1,kmax); % Vo(k) = Vc(k);
%
for k=2:kmax
   Vc(k)=(A*Vc(k-1)+B*Ein)*dt+Vc(k-1);
   Vo(k)=D*Vc(k)+E*Ein;
   T(k)=k*dt;
   F(k)=(Ein*Rp/R1)*(1-exp(-T(k)/(Rp*C1))); % F(k) = inverse LaPlace tranform
end
%
h=plot(T/ms,Vo,'k',T/ms,F,'r');
set(h,'LineWidth',2);
grid on
xlabel('Time (ms)');
ylabel('Volts');
title('RC Time Response');
legend('Vo(k)','F(k)');
figure(1);
%
% Note that difference between F(k) and the time iteration Vo(k).  This can be reduced by
% decreasing dt with a corresponding increase in kmax to retain the same sweep time.
% For example, if dt = 250us, kmax should be set to 200.

author: Sable

repo: mcbench-benchmarks


text:
function g = localthresh(f, nhood, a, b, meantype)
%LOCALTHRESH Local thresholding.
%	G = LOCALTHRESH(F, NHOOD, A, B, MEANTYPE) thresholds image F by
%	computing a local threshold at the center,(x, y), of every
%	neighborhood in F. The size of the neighborhoods is defined by
%	NHOOD, an array of zeros and ones in which the nonzero elements
%	specify the neighbors used in the computation of the local mean
%	and standard deviation. The size of NHOOD must be odd in both
%	dimensions.
%
%	The segmented image is given by
%
%           1  if (F > A*SIG) AND (F > B*MEAN)
%       G =
%           0  otherwise
%
%	where SIG is an array of the same size as F containing the local
%	standard deviations. If MEANTYPE = 'local' (the default), then
%	MEAN is an array of local means.  If MEANTYPE = 'global', then
%	MEAN is the global (image) mean, a scalar.  Constants A and B
%	are nonnegative scalars.
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt

% Intialize.
f = tofloat(f);

% Compute the local standard deviations.
SIG = stdfilt(f, nhood);
% Compute MEAN.
if nargin == 5 && strcmp(meantype,'global')
   MEAN = mean2(f);
else
   MEAN = localmean(f, nhood); % This is a custom function.
end

% Obtain the segmented image.
g = (f > a*SIG) & (f > b*MEAN);

author: dipum

repo: dipum-toolbox

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


text:
function plotexport(prefix)
saveas(gcf,prefix,'fig');
print( gcf, '-depsc2', prefix );
print( gcf, '-dpng', prefix );
switch isunix
    case 0
    print( gcf, '-dmeta', prefix );
end


text:
function [x, y, z] = rfPlot(params, RF, parent, peak)
% rfPlot - script to visualize cropped RF
%
% [x, y, z] = rfPlot(params, RF, [parent=axes in new figure]);
%
% Will produce a plot illustrating the estimated location of a
% 2D Gaussian receptive field, with a grid.
%
% RF is a
% the 'parent' argument directs where to display the plot. Default
% is to create a new figure with its own axes.
%
% 2006/02 SOD: wrote it.
% 2006/09 RAS: added optional 'parent' argument, so you can
% direct the plot to a subplot axes.
% 2008/06 RAS: updated calculation of RF grid to use the X, Y sample points
% in params.analysis.X (and .Y). This replaces a previous method using the
% sample rate; I've found that when you recompute the stimulus (e.g.
% rmRecomputeParams), X/Y sample points for which no stimulus was presented
% are omitted from the analysis. We use this sampling grid to be
% consistent, and prevent bugs in code like rmPlotGUI.
% 2008/07 SOD: reverted back to original. See comments below at the
% relevant code. Must validate rmRecomputeParams.
if ~exist('parent','var') || isempty(parent), figure; parent = gca;      end;
if ~exist('peak','var'), peak = [];      end;

[x,y] = prfSamplingGrid(params);
z    = NaN(size(x));
z(params.stim(1).instimwindow) = RF;
z    = reshape(z,size(x));

% plot
axes(parent);
cla;
hold on;
surf(x,y,z,'LineStyle','none');

%draw lines at every degree
mylines = 0; %[-floor(params.analysis.fieldSize):floor(params.analysis.fieldSize)];
for ll = 1:numel(mylines),
    for n=1:2,
        if n==1,
            ii = find(x==mylines(ll) & isfinite(z));
            [xs, is] = sort(y(ii));
        else
            ii = find(y==mylines(ll) & isfinite(z));
            [xs, is] = sort(x(ii));
        end;
        if ~isempty(is),
            ii = ii(is);
            h  = line(x(ii),y(ii),z(ii));
            if mylines(ll)==0,
                set(h,'LineWidth',1,'Color',[0 0 0]);
            else
                set(h,'LineWidth',0.5,'Color',[0 0 0]);
            end;
        end;
    end;
end;

minz = min(z(:));
maxz = peak;
if isempty(maxz), maxz = max(z(:)); end;
if isnan(maxz), maxz = 0.1; minz = -0.1; end;
if minz==maxz,
    minz = minz - 0.1;
    maxz = maxz + 0.1;
end;

% scale axis
axis([min(x(:)) max(x(:)) min(y(:)) max(y(:)) minz maxz]);
axis image;

% scale colorbar to be centered on zero
caxis([-1 1].*max(abs(minz),abs(maxz)));

% axis labels
xlabel('x (deg)');
ylabel('y (deg)');
zlabel('BOLD amplitude (%/deg^2/sec)');
title('pRF profile');
hold off;

return


text:
function [varargout]=triThick(varargin)

% function [E,V,Fp1,Fp2]=triThick(Fp,Vp,dirSet,layerThickness,numSteps)
% ------------------------------------------------------------------------
%
% Kevin Mattheus Moerman
% gibbon.toolbox@gmail.com
%
%
% 2019/04/18 Created based on quadThick
% ------------------------------------------------------------------------

%%

[E,V,Fp1,Fp2]=patchThick(varargin{:});

%Collect output
varargout{1}=E;
varargout{2}=V;
varargout{3}=Fp1;
varargout{4}=Fp2;

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program 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/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
classdef MimZoomPlugin < MimGuiPlugin
    % MimZoomPlugin. Gui Plugin for using a preset bone window/level
    %
    %     You should not use this class within your own code. It is intended to
    %     be used by the gui of the TD MIM Toolkit.
    %
    %     MimZoomPlugin is a Gui Plugin for the MIM Toolkit.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties
        ButtonText = 'Zoom'
        SelectedText = 'Zoom'
        ToolTip = 'Enables Matlab''s zoom tool'
        Category = 'Tools'
        Visibility = 'Dataset'
        Mode = 'Toolbar'

        HidePluginInDisplay = false
        PTKVersion = '1'
        ButtonWidth = 6
        ButtonHeight = 1
        Icon = 'zoom.png'
        Location = 15
    end

    methods (Static)
        function RunGuiPlugin(gui_app)
            gui_app.ImagePanel.SetControl('Zoom');
        end

        function enabled = IsEnabled(gui_app)
            enabled = gui_app.IsDatasetLoaded;
        end

        function is_selected = IsSelected(gui_app)
            is_selected = strcmp(gui_app.ImagePanel.SelectedControl, 'Zoom');
        end

    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
% Copyright (C) 2013 Quan Wang <wangq10@rpi.edu>,
% Signal Analysis and Machine Perception Laboratory,
% Department of Electrical, Computer, and Systems Engineering,
% Rensselaer Polytechnic Institute, Troy, NY 12180, USA

% this is a demo showing the use of our dynamic time warping package
% we provide both Matlab version and C/MEX version
% the C/MEX version is much faster and highly recommended

clear;clc;close all;

mex dtw_c.c;

a=rand(500,3);
b=rand(520,3);
w=50;

tic;
d1=DTWComputer.mddtw(a,b);
t1=toc;

tic;
d2=dtw_c(a,b,w);
t2=toc;

fprintf('Using Max version: distance=%f, running time=%f\n',d1,t1);
fprintf('Using C/MEX version: distance=%f, running time=%f\n',d2,t2);

author: avenix

repo: WDK

save_path: github-repos/MATLAB/avenix-WDK


text:
function varargout = dtmfdet(x,fs,varargin)
% DTMFDET   DTMF Tone Detection Using Goertzel Algorithm.
%
%  Y = DTMFDET(X,FS) returns the DTMF coded sequence
%  in present vector X. The sound sample frequency
%  is given by FS.
%
%  Y = DTMFDET(X,FS,'strong') uses a more robust algorithm
%  to detect two equall DTMF signals in a sequence. Some times
%  not using this strong and longer method can result in
%  only one DTMF signal detected, instead of the two.
%

LIM = .7; % Limit for detecting tone (if the intensity of the freq from goertzel is less the LIM* the max intensity, the tone wont be detected)
T = ['1' '2' '3' 'A' '4' '5' '6' 'B' '7' '8' '9' 'C' '*' '0' '#' 'D'];
FL = [697 770 852 941];
FH = [1209 1336 1477 1633];
OL = .0; % Overlapping block in goertzel.

pl = 0; % Plotting flag

for i = 1 : length(varargin)
    if strcmp(varargin{i},'strong')
        OL = .4;
    elseif strcmp(varargin{i},'plot')
        pl = 1;
    else
        error('Unknown parameter.')
    end
end

time = 0.02;
N = time*fs;

ADV = max(round((1-OL)*N),1);
YL = zeros(4,floor(length(x)/ADV));
YH = YL;

for i = 1 : 4
    YL(i,:) = goertzel(x,fs,FL(i),N,'mag','overlap',OL);
    YH(i,:) = goertzel(x,fs,FH(i),N,'mag','overlap',OL);
end

[a,a] = sort(YL,1);
[b,b] = sort(YH,1);

YL(YL<LIM*max(YL(:))) = 0;
YH(YH<LIM*max(YH(:))) = 0;

[l1,l2] = max(YL);
[h1,h2] = max(YH);

r = T((l2-1)*4+h2);
r(l1==0) = '-';
r(h1==0) = '-';

if pl == 1
    plot(YL')
    figure(gcf+1)
    plot(YH')
end

ct = 1;
for i = 2 : length(r)
    if r(i)~='-' & r(i)~=r(i-1)
        y(ct) = r(i);
        ct = ct+1;
    end
end

varargout{1} = y;


text:
% Check the configuration section to make the necessary adjustments to
% run this script.

% @
% Copyright (C) 2017 Jonas Ruesch
%
% This program 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/>.
% @

% Configuration
% ===================================================================================
tornadoDirectory = 'F:\svn\dev\matlab\tornado\T135_export'; % where is Tornado installed?
nameTornadoAircraftGeometry = 'ExperimentalCarrier'; % which Tornado model to use?
estimatedVelocity = 12; % 12; % the velocity at which aerodynamic coefficients are computed.
velocityDeviationTol = 0.5; % how much the found trimmed velocity can differ from the estimated velocity to accept the solution.
airDensity = 1.225;
elevatorFlapIndex = 1; % which is the elevator flap in the Tornado model?
rudderFlapIndex = 2; % which is the rudder flap in the Tornado model?
centerOfGravity = 0.092; % in meters from zero of Tornado aircraft model reference frame. Tornado x-axis extends aft (!) -> positive.
mass = 1.56; % aircraft mass in kg (overwrites the value set in the Tornado model).
InertiaTensor = [1 0 0; 0 1 0; 0 0 1]; % inertia tensor (TODO: compute a realistic one :)
% The alpha-sweep configuration
alphaStart = 4; %2;
alphaEnd = 6; % 6;
numAlphas = 5;
% The beta-sweep configuration
betaStart = -1;
betaEnd = 1;
numBetas = 5;
computeStaticMargin = 0; % [0/1]
% ===================================================================================

nameLongitudinalModel = 'ExperimentalCarrier_longitudinal';
nameLateralModel = 'ExperimentalCarrier_lateral';

disp(' ');

% Compute longitudinal and lateral LTIs for trimmed gliding.
[ltiLongitudinal, trimmedStateVariablesLongitudinal, ltiLateral, trimmedStateVariablesLateral, Cm_longitudinal, Cn_lateral] = computeLTIs(...
    tornadoDirectory, ... % absolute path to the Tornado root directory
    nameTornadoAircraftGeometry, ... % the name of the Tornado geometry file
    nameLongitudinalModel, ... % the name of the Simulink long. model.
    nameLateralModel,... % the name of the Simulink lateral model.
    alphaStart, alphaEnd, ... % range in degrees over which coefficients are computed for the longitudinal model
    numAlphas, ... % the number of alphas to sample
    betaStart, betaEnd, ... range in degrees over which coefficients are computed for the lateral model
    numBetas, ... the number of betas to sample
    estimatedVelocity, ... % velocity in m/s for which aerodyn. coeffs. are computed
    velocityDeviationTol, ... % factor [0 ... 1].
    airDensity, ... % kg/m^3 -> 1.225 at sea level
    elevatorFlapIndex, ...
    rudderFlapIndex, ...
    mass, ...
    centerOfGravity, ...
    InertiaTensor); % inertia tensor
trimmedAlphaRad = atan(trimmedStateVariablesLongitudinal(2)/trimmedStateVariablesLongitudinal(1)); % atan(w/u);
trimmedBetaRad = atan(trimmedStateVariablesLateral(1)/trimmedStateVariablesLongitudinal(1)); % atan(w/v);

% Compute static margin of the aircraft.
if(computeStaticMargin)
    disp('Computing static margin of longitudinal dynamics...');
    [ac, staticMarginMAC, staticMarginInMeters] = computeStaticMargin(nameTornadoAircraftGeometry, ...
        estimatedVelocity, airDensity, tornadoDirectory, centerOfGravity);
    disp(['Static margin as MAC factor  : ' num2str(staticMarginMAC)]);
    disp(['Static margin (sm)      (m)  : ' num2str(staticMarginInMeters)]);
    disp(['Aerodynamic center (ac) (m)  : ' num2str(ac)]);
    disp(['ac - sm = center of gravity  : ' num2str(ac - staticMarginInMeters)]);
end

% Longitudinal
% ========================================================================
% disp('The longitudinal LTI:');
% ltiLongitudinal

% sort complex numbers
polesLongitudinal = cplxpair(pole(ltiLongitudinal));  % sort first complex pairs (each with negative imag. part first, then pure real (ascending?)).

% Poles Map
figure();
h = subplot(2,2,1);
hold on;
plot(h, real(polesLongitudinal(1)), imag(polesLongitudinal(1)), 'rx');
plot(h, real(polesLongitudinal(2)), imag(polesLongitudinal(2)), 'bx');
plot(h, real(polesLongitudinal(3)), imag(polesLongitudinal(3)), 'ro');
plot(h, real(polesLongitudinal(4)), imag(polesLongitudinal(4)), 'bo');
hold off;
legend('phugoid mode 1', 'phugoid mode 2', 'short period mode 1', 'short period mode 2');
title('Longitudinal Poles Map');
xlabel('re');
ylabel('im');

% Step response
% Remember: for a trimmed glide, LTI state variables are all 0. Thus, the plot
% shows the resulting deviation of state variables from the trimmed state.
% Also note: with the current model (very little air resistance -> very
% good gliding ratio, the non-linear matlab model shows that it is sensible
% that w (the velocity in body-z) is becoming negative (moving 'upwards' in
% body-z).
h = subplot(2,2,2);
elevatorDeflectionStep = 2*pi/180;
stepOpt = stepDataOptions('StepAmplitude', elevatorDeflectionStep); % step response for elevatorDeflectionStep radian elevator deflection
stepplot(h, ltiLongitudinal, stepOpt);
title('Longitudinal Step Response');
hAllAxes = findobj(gcf,'type','axes');
hLeg = findobj(hAllAxes,'tag','legend');
hAxes = setdiff(hAllAxes,hLeg); % All axes which are not
for k=1:length(hAxes)
   if ~isempty(strfind(hAxes(k).YLabel.String, '(1)'))
       hAxes(k).Title.String = ['State change for an elevator step of ' num2str(elevatorDeflectionStep*180/pi) ' degrees'];
       hAxes(k).YLabel.String = 'u (m/s)'; % velocity in body x-axis
   elseif ~isempty(strfind(hAxes(k).YLabel.String, '(2)'))
       hAxes(k).YLabel.String = 'w (m/s)'; % velocity in body z-axis
   elseif ~isempty(strfind(hAxes(k).YLabel.String, '(3)'))
       hAxes(k).YLabel.String = 'q (rad/s)'; % rotation velocity around y-axis
   elseif ~isempty(strfind(hAxes(k).YLabel.String, '(4)'))
       hAxes(k).YLabel.String = 'theta (rad)'; % angle body-x vs earth-x
   end
end

% Impulse response
h = subplot(2,2,3);
impulseplot(h, ltiLongitudinal);
title('Longitudinal Impulse Response');
hAllAxes = findobj(gcf,'type','axes');
hLeg = findobj(hAllAxes,'tag','legend');
hAxes = setdiff(hAllAxes,hLeg); % All axes which are not
for k=1:length(hAxes)
   if ~isempty(strfind(hAxes(k).YLabel.String, '(1)'))
       hAxes(k).Title.String = 'State change for an elevator dirac impulse';
       hAxes(k).YLabel.String = 'u (m/s)'; % velocity in body x-axis
   elseif ~isempty(strfind(hAxes(k).YLabel.String, '(2)'))
       hAxes(k).YLabel.String = 'w (m/s)'; % velocity in body z-axis
   elseif ~isempty(strfind(hAxes(k).YLabel.String, '(3)'))
       hAxes(k).YLabel.String = 'q (rad/s)'; % rotation velocity around y-axis
   elseif ~isempty(strfind(hAxes(k).YLabel.String, '(4)'))
       hAxes(k).YLabel.String = 'theta (rad)'; % angle body-x vs earth-x
   end
end

% Plot pitch moment curve
h = subplot(2,2,4);
hold on;
plot(h, [alphaStart:(alphaEnd-alphaStart)/(numAlphas-1):alphaEnd], Cm_longitudinal, 'r-');
plot(h, trimmedAlphaRad*180/pi, 0, 'bx');
hold off;
legend('Cm', ['trimmed alpha = ' num2str(trimmedAlphaRad*180/pi) ' (deg)']);
title('Pitch moment');
xlabel('alpha');
ylabel('Cm');

% Compute damping, frequency, period from poles. See also p83 in Caughey.

% Short period poles
polesShortPeriod = polesLongitudinal(3:4); % choose last ones from sorted complex numbers
polesShortPeriodReal = real(polesShortPeriod);
polesShortPeriodImag = imag(polesShortPeriod);

Damping_sp = sqrt(1./(1+(polesShortPeriodImag./polesShortPeriodReal).^2)); % Damping ratio
Frequency_sp = -polesShortPeriodReal./Damping_sp; % Undampted, natural frequency per second
Period_sp = 2*pi./(Frequency_sp.*sqrt(1-Damping_sp.^2)); % seconds
Cycles_half_sp = ((log(2))/(2*pi))*sqrt(1-Damping_sp.^2)./Damping_sp; % Number of cycles to damp to half the amplitude

% Output result for longitudinal trimmed state
disp('Short Period Mode Properties (longitudinal): ');
disp('                                                 Pole 1      Pole 2');
disp('==============================================================================');
disp(sprintf(['Damping ratio                                 : ' num2str(Damping_sp(:)')]));
disp(sprintf(['Undampted, natural frequency             (1/s): ' num2str(Frequency_sp(:)')]));
disp(sprintf(['Period                                     (s): ' num2str(Period_sp(:)')]));
disp(sprintf(['Num. cycles to damp to half the amplitude     : ' num2str(Cycles_half_sp(:)')]));
disp(' ');

% Phugoid poles
polesPhugoid = polesLongitudinal(1:2); % choose first ones from sorted complex numbers
polesPhugoidReal = real(polesPhugoid);
polesPhugoidImag = imag(polesPhugoid);

Damping_ph = sqrt(1./(1+(polesPhugoidImag./polesPhugoidReal).^2)); % Damping ratio
Frequency_ph = -polesPhugoidReal./Damping_ph; % Undampted, natural frequency per second
Period_ph = 2*pi./(Frequency_ph.*sqrt(1-Damping_ph.^2)); % seconds
Cycles_half_ph = ((log(2))/(2*pi))*sqrt(1-Damping_ph.^2)./Damping_ph; % Number of cycles to damp to half the amplitude

% Output result for longitudinal trimmed state
disp('Phugoid Mode Properties (longitudinal): ');
disp('                                                 Pole 1      Pole 2');
disp('==============================================================================');
disp(sprintf(['Damping ratio                                 : ' num2str(Damping_ph(:)')]));
disp(sprintf(['Undampted, natural frequency             (1/s): ' num2str(Frequency_ph(:)')]));
disp(sprintf(['Period                                     (s): ' num2str(Period_ph(:)')]));
disp(sprintf(['Num. cycles to damp to half the amplitude     : ' num2str(Cycles_half_ph(:)')]));
disp(' ');

% Lateral
% ========================================================================

%disp('The lateral LTI:');
%ltiLateral

% sort complex numbers
polesLateral = cplxpair(pole(ltiLateral)); % sort first complex pairs (each with negative imag. part first, then pure real (ascending?)).

% Poles Map
figure();
h = subplot(2,2,1);
hold on;
plot(h, real(polesLateral(1)), imag(polesLateral(1)), 'rx');
plot(h, real(polesLateral(2)), imag(polesLateral(2)), 'bx');
plot(h, real(polesLateral(3)), imag(polesLateral(3)), 'ro');
plot(h, real(polesLateral(4)), imag(polesLateral(4)), 'bo');
hold off;
legend('dutch roll mode 1', 'dutch roll mode 2', 'rolling mode', 'spiral mode');
title('Lateral Poles Map');
xlabel('re');
ylabel('im');

% Step response
% Remember: for a trimmed glide, LTI state variables are all 0. Thus, the plot
% shows the resulting deviation of state variables from the trimmed state.
h = subplot(2,2,2);
rudderDeflectionStep = 2*pi/180;
stepOpt = stepDataOptions('StepAmplitude', rudderDeflectionStep); % step response for elevatorDeflectionStep radian elevator deflection
stepplot(h, ltiLateral, stepOpt);
title('Lateral Step Response');
hAllAxes = findobj(gcf,'type','axes');
hLeg = findobj(hAllAxes,'tag','legend');
hAxes = setdiff(hAllAxes,hLeg); % All axes which are not
for k=1:length(hAxes)
   if ~isempty(strfind(hAxes(k).YLabel.String, '(1)'))
       hAxes(k).Title.String = ['State change for a rudder step of ' num2str(rudderDeflectionStep*180/pi) ' degrees'];
       hAxes(k).YLabel.String = 'v (m/s)'; % velocity in body y-axis
   elseif ~isempty(strfind(hAxes(k).YLabel.String, '(2)'))
       hAxes(k).YLabel.String = 'p (rad/s)'; % rotation velocity around x-axis
   elseif ~isempty(strfind(hAxes(k).YLabel.String, '(3)'))
       hAxes(k).YLabel.String = 'phi (rad)'; % roll angle (euler phi)
   elseif ~isempty(strfind(hAxes(k).YLabel.String, '(4)'))
       hAxes(k).YLabel.String = 'r (rad/s)'; % rotation velocity around z-axis
   end
end

% Impulse response
h = subplot(2,2,3);
impulseplot(h, ltiLateral);
title('Lateral Impulse Response');
hAllAxes = findobj(gcf,'type','axes');
hLeg = findobj(hAllAxes,'tag','legend');
hAxes = setdiff(hAllAxes,hLeg); % All axes which are not
for k=1:length(hAxes)
   if ~isempty(strfind(hAxes(k).YLabel.String, '(1)'))
       hAxes(k).Title.String = 'State change for a rudder dirac impulse';
      hAxes(k).YLabel.String = 'v (m/s)'; % velocity in body y-axis
   elseif ~isempty(strfind(hAxes(k).YLabel.String, '(2)'))
       hAxes(k).YLabel.String = 'p (rad/s)'; % rotation velocity around x-axis
   elseif ~isempty(strfind(hAxes(k).YLabel.String, '(3)'))
       hAxes(k).YLabel.String = 'phi (rad)'; % roll angle (euler phi)
   elseif ~isempty(strfind(hAxes(k).YLabel.String, '(4)'))
       hAxes(k).YLabel.String = 'r (rad/s)'; % rotation velocity around z-axis
   end
end

% Plot yaw moment curve
h = subplot(2,2,4);
hold on;
plot(h, [betaStart:(betaEnd-betaStart)/(numBetas-1):betaEnd], Cn_lateral, 'r-');
plot(h, trimmedBetaRad*180/pi, 0, 'bx');
hold off;
legend('Cn', ['trimmed beta = ' num2str(trimmedBetaRad*180/pi) ' (deg)']);
title('Yaw moment');
xlabel('beta');
ylabel('Cn');

% Compute damping, frequency, period from poles. See also p83 in Caughey.

% Roll pole
poleRoll = polesLateral(3);
poleRollReal = real(poleRoll);
poleRollImag = imag(poleRoll);

Damping_roll = sqrt(1./(1+(poleRollImag./poleRollReal).^2)); % Damping ratio
Frequency_roll = -poleRollReal./Damping_roll; % Undampted, natural frequency per second
Period_roll = 2*pi./(Frequency_roll.*sqrt(1-Damping_roll.^2)); % seconds
Cycles_half_roll = ((log(2))/(2*pi))*sqrt(1-Damping_roll.^2)./Damping_roll; % Number of cycles to damp to half the amplitude

disp('Roll Mode Properties (lateral): ');
disp('                                                 Pole 1');
disp('==============================================================================');
disp(sprintf(['Damping ratio                                 : ' num2str(Damping_roll(:)')]));
disp(sprintf(['Undampted, natural frequency             (1/s): ' num2str(Frequency_roll(:)')]));
disp(sprintf(['Period                                     (s): ' num2str(Period_roll(:)')]));
disp(sprintf(['Num. cycles to damp to half the amplitude     : ' num2str(Cycles_half_roll(:)')]));
disp(' ');

% Spiral pole
poleSpiral = polesLateral(4);
poleSpiralReal = real(poleSpiral);
poleSpiralImag = imag(poleSpiral);

Damping_spiral = sqrt(1./(1+(poleSpiralImag./poleSpiralReal).^2)); % Damping ratio
Frequency_spiral = -poleSpiralReal./Damping_spiral; % Undampted, natural frequency per second
Period_spiral = 2*pi./(Frequency_spiral.*sqrt(1-Damping_spiral.^2)); % seconds
Cycles_half_spiral = ((log(2))/(2*pi))*sqrt(1-Damping_spiral.^2)./Damping_spiral; % Number of cycles to damp to half the amplitude

disp('Spiral Mode Properties (lateral): ');
disp('                                                 Pole 1');
disp('==============================================================================');
disp(sprintf(['Damping ratio                                 : ' num2str(Damping_spiral(:)')]));
disp(sprintf(['Undampted, natural frequency             (1/s): ' num2str(Frequency_spiral(:)')]));
disp(sprintf(['Period                                     (s): ' num2str(Period_spiral(:)')]));
disp(sprintf(['Num. cycles to damp to half the amplitude     : ' num2str(Cycles_half_spiral(:)')]));
disp(' ');

% Dutch Roll poles
polesDutchRoll = polesLateral(1:2); % choose first ones from sorted complex numbers
polesDutchRollReal = real(polesDutchRoll);
polesDutchRollImag = imag(polesDutchRoll);

Damping_dutch = sqrt(1./(1+(polesDutchRollImag./polesDutchRollReal).^2)); % Damping ratio
Frequency_dutch = -polesDutchRollReal./Damping_dutch; % Undampted, natural frequency per second
Period_dutch = 2*pi./(Frequency_dutch.*sqrt(1-Damping_dutch.^2)); % seconds
Cycles_half_dutch = ((log(2))/(2*pi))*sqrt(1-Damping_dutch.^2)./Damping_dutch; % Number of cycles to damp to half the amplitude

% Output result for longitudinal trimmed state
disp('Dutch Roll Mode Properties (lateral): ');
disp('                                                 Pole 1      Pole 2');
disp('==============================================================================');
disp(sprintf(['Damping ratio                                 : ' num2str(Damping_dutch(:)')]));
disp(sprintf(['Undampted, natural frequency             (1/s): ' num2str(Frequency_dutch(:)')]));
disp(sprintf(['Period                                     (s): ' num2str(Period_dutch(:)')]));
disp(sprintf(['Num. cycles to damp to half the amplitude     : ' num2str(Cycles_half_dutch(:)')]));
disp(' ');

% =========================================================================

% Observability
disp(['Number of unobservable states (longitudinal)  : ' num2str(length(ltiLongitudinal.a) - rank(obsv(ltiLongitudinal)))]);
disp(['Number of unobservable states (lateral)       : ' num2str(length(ltiLateral.a) - rank(obsv(ltiLateral)))]);
disp(' ');

% Controllability
disp(['Number of uncontrollable states (longitudinal): ' num2str(length(ltiLongitudinal.a) - rank(ctrb(ltiLongitudinal)))]);
disp(['Number of uncontrollable states (lateral)     : ' num2str(length(ltiLateral.a) - rank(ctrb(ltiLateral)))]);
disp(' ');

% TODO
% ========================================================================

% 1) Run this whole thing for different center of gravity settings and plot
% Cm curve and poles.

author: jrgenerative

repo: fixed-wing-sim


text:
function start_num = p06_start_num ( )

%*****************************************************************************80
%
%% P06_START_NUM returns the number of starting points for problem 6.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 May 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer START_NUM, the number of starting points.
%
  start_num = 2;

  return
end


text:
function [meg] = read_ctf_meg4(fname, hdr, begsample, endsample, chanindx)

% READ_CTF_MEG4 reads specified samples from a CTF continous datafile
% It neglects all trial boundaries as if the data was acquired in
% non-continous mode.
%
% Use as
%   [meg] = read_ctf_meg4(filename, hdr, begsample, endsample, chanindx)
% where
%   filename    name of the datafile, including the .meg4 extension
%   header      with all data information (from read_ctf_meg4)
%   begsample   index of the first sample to read
%   endsample   index of the last sample to read
%   chanindx    index of channels to read (optional, default is all)
%
% See also READ_CTF_MEG4

% "VSM MedTech Ltd. authorizes the release into public domain under the
% GPL licence of the Matlab source code files "read_ctf_res4.m" and
% "read_ctf_meg4.m" by the authors of said files from the F.C. Donders
% Centre, Nijmegen, The Netherlands."

% Author(s): Jim McKay November 1999
% Last revision: Jim McKay
% Copyright (c) 1999-2000 CTF Systems Inc. All Rights Reserved.
%
% modifications Copyright (C) 2002, Ole Jensen
% modifications Copyright (C) 2003, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.ru.nl/neuroimaging/fieldtrip
% 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: read_ctf_meg4.m 945 2010-04-21 17:41:20Z roboos $

% use global flag for feedback
global fb
if isempty(fb)
  fb = 0;
end

nsmp = hdr.nSamples;
ntrl = hdr.nTrials;
nchn = hdr.nChans;

if begsample<1,              error('cannot read before the start of the data');        end
if endsample>nsmp*ntrl*nchn, error('cannot read beyond the end of the data');          end
if begsample>endsample,      error('cannot read a negative number of samples');        end
if nargin<5,                 chanindx = 1:nchn;                                        end
if isempty(chanindx),        error('no channels were specified for reading CTF data'); end

%open the .meg4 file
fid = fopen(fname,'r','ieee-be');
if fid == -1,
  error('could not open datafile');
end

%check whether it is a known format
CTFformat=char(fread(fid, 8, 'uint8'))';
% This function was written for MEG41RS, but also seems to work for some other formats
if ~strcmp(CTFformat(1,1:7),'MEG41CP') && ~strcmp(CTFformat(1,1:7),'MEG4CPT')
  warning('meg4 format (%s) is not supported for file %s, trying anyway...', CTFformat(1,1:7), fname);
end

%determine size of .meg4 file
fseek(fid, 0, 'eof');
nbytes   = ftell(fid);

%number of trials per 2GB file FIXME assumes constancy across the .meg4 files
ntrlfile = round((nbytes-8)/(4*nchn*nsmp));
%ntrlfile = (nbytes-8)/(4*nchn*nsmp);
openfile = 0;

%determine which trials have to be read
begtrial = ceil(begsample/nsmp);
endtrial = ceil(endsample/nsmp);
trials   = begtrial:endtrial;

%to ensure correct sample handling in the case of multiple trials
sumsmp = [(begtrial-1):(endtrial-1)]*nsmp;
rawbeg = 1;

minchn  = min(chanindx); %to ensure correct channel handling in the case of blockwise reading
maxchn  = max(chanindx);
loopchn = length(trials)==1 || (maxchn-minchn+1)<=nchn; %decide whether to read in blockwise, or the specified samples per channel

raw     = zeros(endsample-begsample+1, length(chanindx)); %allocate memory
for trllop = 1:length(trials)
  trlnr  = trials(trllop);
  filenr = floor(trlnr/(ntrlfile+0.1));

  %ensure that the correct .meg4 file is open
  if filenr~=openfile && filenr>0,
    fclose(fid);
    nextname = sprintf('%s.%d_meg4', fname(1:(end-5)), filenr);
    if fb
      fprintf('data goes beyond 2GB file boundary, continuing with %s\n', nextname);
    end
    fid = fopen(nextname,'r','ieee-be');
    fseek(fid, 0, 'eof');
    openfile = filenr;
  end

  %this is relative to the current datafile
  rawtrl = mod(trlnr-1, ntrlfile) + 1;
  offset = 8 + 4*(rawtrl-1)*nsmp*nchn;

  %begin and endsamples expressed as samples with respect to the current trial
  tmpbeg = max(begsample-sumsmp(trllop), 1);
  tmpend = min(endsample-sumsmp(trllop), nsmp);
  rawend = rawbeg+tmpend-tmpbeg;

  %either read per channel or read entire trialblock and postselect the channels
  if loopchn,
    for chnlop = 1:length(chanindx)
      %this is relative to the current trial
      chanoffset = 4*(chanindx(chnlop)-1)*nsmp;
      sampoffset = 4*(tmpbeg-1);
      fseek(fid, offset+chanoffset+sampoffset, 'bof');
      [tmp, count]               = fread(fid,[tmpend-tmpbeg+1,1],'int32');
      raw(rawbeg:rawend, chnlop) = tmp;
    end
  else
    %this is relative to the current trial
    chanoffset = 4*(minchn-1)*nsmp;
    fseek(fid, offset+chanoffset, 'bof');
    ntmpchn               = maxchn - minchn + 1;
    [tmp, count]          = fread(fid,[nsmp,length(ntmpchn)],'int32');
    selchn                = chanindx - minchn + 1; %relative to the first channel to be read
    raw(rawbeg:rawend, :) = tmp(tmpbeg:tmpend, selchn);
  end
  rawbeg = rawend+1;
end
fclose(fid);

% multiply the dimensionless values with the calibration value
gain = hdr.gainV(chanindx); % only for selected channels
meg  = raw';            % transpose the raw data
for i=1:size(meg,1)
  meg(i,:) = gain(i)*meg(i,:);
end

author: PatternRecognition

repo: OpenBMI


text:
function varargout = daqwaterfall(varargin)
% daqwaterfall
%
%  Displays a waterfall plot with data streamed in from the Data
%  Acquisition Toolbox.  This demo was created using Guide in MATLAB 6.5.
%  This demo uses a lot of the code from demoai_fft.
%
%  The input source is hardcoded to channel one of the system soundcard.
%  You should be able to change this if necessary by modifying the
%  parameter in the I N P U T S section of the file daqwaterfall.m
%
% Last Modified by GUIDE v2.5 06-Jan-2003 23:34:34
% Author: Daniel Lee (dlee@mathworks.com)

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @daqwaterfall_OpeningFcn, ...
                   'gui_OutputFcn',  @daqwaterfall_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin & isstr(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before daqwaterfall is made visible.
function daqwaterfall_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to daqwaterfall (see VARARGIN)

% Choose default command line output for daqwaterfall
handles.output = hObject;

%%%%%%%%%%%%%%%%%%%%%%%
% ++++ I N P U T S ++++
%%%%%%%%%%%%%%%%%%%%%%%
% Change these settings to select a different source.
handles.adaptor = 'winsound';
handles.id      = 0;
handles.chan    = 1;

handles.samplesPerTrigger = 1024;
handles.sampleRate = 44100;
handles.numTraces = 10;   % number of traces to show in the waterfall.
handles.cycleTime = .9;   % Proportional to the amount of time spent per
                          %  visualization on CycleAll setting.
%%%%%%%%%%%%%%%%%%%%%%%
% ---- I N P U T S ----
%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ++++ S E T U P   T H E   F I G U R E ++++
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*%%%%%%
set(handles.figure1,'Color',get(handles.tTitle,'BackgroundColor'));

axes(handles.axes1);
handles.hLine1 = plot(zeros(1,handles.samplesPerTrigger)');
set(handles.hLine1,'Color', [.1 .1 0.5]);
set(handles.axes1,'Color',[235/255 255/255 235/255])
set(handles.axes1,'XGrid','on','YGrid','on')
t=title('Time Domain Signal','Color',[.05 .05 .25],'FontWeight','Bold','FontSize',9);
xlabel('Time (s)','FontSize',8);
ylabel('Voltage (V)','FontSize',8);

axes(handles.axes2);
handles.hLine2 = plot(zeros(1,handles.samplesPerTrigger/2)');
set(handles.hLine2,'Color', [.1 0.5 .1]);
set(handles.axes2,'Color',[235/255 255/255 235/255])
set(handles.axes2,'XGrid','on','YGrid','on')
t=title('Frequency Domain Signal','Color',[.05 0.25 .05],'FontWeight','Bold','FontSize',9);
xlabel('Frequency (Hz)','FontSize',8);
ylabel('Magnitude (dB)','FontSize',8);

axes(handles.axes3);
set(handles.axes3,'View',[103 10]);
set(handles.axes3,'Color',[234/255 234/255 255/255]);
grid(handles.axes3,'on');
h = get(handles.axes3,'title');
set(h,'string','Waterfall Plot','FontWeight','Bold','Color',[.25 .05 .05],'FontSize',9);
h = get(handles.axes3,'ylabel');
set(h,'string','Frequency (Hz)','FontSize',8);
h = get(handles.axes3,'zlabel');
set(h,'string','Magnitude (dB)','FontSize',8);

set(hObject,'RendererMode','Manual')  %  If you don't do this, the surface plot
set(hObject,'Renderer','OpenGL')      %    will draw VERY slowly.

set(handles.tSource,'String',sprintf('%s:%d',handles.adaptor,handles.id));
set(handles.tChannel,'String',num2str(handles.chan));
set(handles.poSampleRate,'String',[{'44100'},{'22000'},{'8000'}]);
set(handles.poPlotType,'String',[{'CycleAll'},{'Classic'},{'Classic(Top)'},{'Mosaic'},{'Waterfall'},{'Rotate'}]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ---- S E T U P   T H E   F I G U R E ----
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ++++ D R A W  T H E  L O G O ++++
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
L = 40*membrane(1,25);
axes(handles.axes4);
set(handles.axes4,...
    'CameraPosition', [-193.4013 -265.1546  220.4819],...
    'XLim',[1 51], ...
    'YLim',[1 51], ...
    'Visible','off', ...
    'ZLim',[-13 40]);

s = surface(L, ...
    'EdgeColor','none', ...
    'FaceColor',[0.9 0.2 0.2], ...
    'FaceLighting','phong', ...
    'AmbientStrength',0.3, ...
    'DiffuseStrength',0.6, ...
    'Clipping','off',...
    'BackFaceLighting','lit', ...
    'SpecularStrength',1.1, ...
    'SpecularColorReflectance',1, ...
    'SpecularExponent',7);
l1 = light('Position',[40 100 20], ...
    'Style','local', ...
    'Color',[0 0.8 0.8]);
l2 = light('Position',[.5 -1 .4], ...
    'Color',[0.8 0.8 0]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ---- D R A W  T H E  L O G O ----
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ai=localSetupAI(handles);
handles.ai = ai;

% Update handles structure
guidata(hObject, handles);

localStartAI(ai);

% UIWAIT makes daqwaterfall wait for user response (see UIRESUME)
% uiwait(handles.figure1);

function localStartAI(ai)
%%%%%%%%%%%%%%%%%%%%%%%%%%
% ++++ S T A R T  A I ++++
%%%%%%%%%%%%%%%%%%%%%%%%%%
start(ai);
trigger(ai);
%%%%%%%%%%%%%%%%%%%%%%%%%%
% ---- S T A R T  A I ----
%%%%%%%%%%%%%%%%%%%%%%%%%%

function localStopAI(ai)
%%%%%%%%%%%%%%%%%%%%%%%%
% ++++ S T O P  A I ++++
%%%%%%%%%%%%%%%%%%%%%%%%
stop(ai);
delete(ai);
%%%%%%%%%%%%%%%%%%%%%%%%
% ---- S T O P  A I ----
%%%%%%%%%%%%%%%%%%%%%%%%

function ai=localSetupAI(handles)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ++++ S E T U P   T H E   A N A L O G   I N P U T ++++
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Object Configuration.
% Create an analog input object with one channel.
ai = analoginput(handles.adaptor, handles.id);
addchannel(ai, handles.chan);

% Configure the callback to update the display.
set(ai, 'TimerFcn', @localfftShowData);

% Configure the analog input object.
set(ai, 'SampleRate', handles.sampleRate);

% Configure the analog input object to trigger manually twice.
%  We do this because we are using peekdata to acquire the data in
%   a timer callback function.
%  The first trigger will fill the buffer with handles.samplesPerTrigger
%   number of samples.  We'll know we have enough samples to start
%   processing data when the analog input object's SamplesAvailable property
%   is equal to handles.samplesPerTrigger.
%  The analog input object will then wait for
%   another manual trigger, and while it is waiting the object will still be
%   in its running state, which means the timer event will run. To keep the
%   object in the running state, we need only never manually trigger this
%   second trigger.
%  Had we set the TriggerRepeat to 0, the analog input object would stop
%   after the first trigger and the timer functions would stop running.
%
set(ai, 'SamplesPerTrigger', handles.samplesPerTrigger);
set(ai, 'TriggerRepeat', 1);
set(ai, 'TriggerType', 'manual');

% Initialize callback parameters.  The TimerAction is initialized
% after figure has been created.
set(ai, 'TimerPeriod', 0.01);
set(ai, 'BufferingConfig',[handles.samplesPerTrigger*2,20]);

% Initialize time and frequency plots with lines of y=0
d=zeros(1,handles.samplesPerTrigger);
time = 1:handles.samplesPerTrigger;
f=1:handles.samplesPerTrigger/2;
mag=zeros(1,handles.samplesPerTrigger/2);

% Store state information in the analog input objects UserData area.
data.storedFFTsIndex = 1;
data.plotSurf        = 0;
data.ai              = ai;
data.getdata         = [d time];
data.daqfft          = [f mag];
data.handle          = [];
data.figureHandles   = handles;
data.view            = [103 10];
data.rotateStep      = 4;
data.counter         = 0;

% Set the object's UserData to data.
set(data.ai, 'UserData', data);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ---- S E T U P   T H E   A N A L O G   I N P U T ----
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% --- Outputs from this function are returned to the command line.
function varargout = daqwaterfall_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

% --- Executes on button press in pbExit or when you press
%     the figure close 'X' button (I set this function to
%     the figures CloseRequestFcn in GUIDE).
function pbExit_Callback(hObject, eventdata, handles)
% hObject    handle to pbExit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
localStopAI(handles.ai);
closereq;

% --- Executes during object creation, after setting all properties.
function poSampleRate_CreateFcn(hObject, eventdata, handles)
% hObject    handle to poSampleRate (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ++++ C H A N G E   T H E   S A M P L E   R A T E  ++++
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% --- Executes on selection change in poSampleRate.
function poSampleRate_Callback(hObject, eventdata, handles)

% hObject    handle to poSampleRate (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = get(hObject,'String') returns poSampleRate contents as cell array
%        contents{get(hObject,'Value')} returns selected item from poSampleRate

% First, stop and delete the current analog input object
localStopAI(handles.ai);

% Extract the new samplerate.
v=get(handles.poSampleRate,'Value');
s=get(handles.poSampleRate,'String');
handles.sampleRate = str2num(s{v});

% Create a new analog input with the new sample rate.
handles.ai = localSetupAI(handles);

% Update handles structure
guidata(hObject, handles);

% Restart the analog input
localStartAI(handles.ai);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ---- C H A N G E   T H E   S A M P L E   R A T E  ----
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% ***********************************************************************
% Calculate the fft of the data.  (Copied from demoai_fft.m)
function [f, mag] = localDaqfft(data,Fs,blockSize)

% Calculate the fft of the data.
xFFT = fft(data);
xfft = abs(xFFT);

% Avoid taking the log of 0.
index = find(xfft == 0);
xfft(index) = 1e-17;

mag = 20*log10(xfft);
mag = mag(1:blockSize/2);

f = (0:length(mag)-1)*Fs/blockSize;
f = f(:);

% ***********************************************************************
% Update the plot. This routine is a Timer callback, it is called
%  automatically at a preset time interval. See line 184 for where
%  this routine is assigned as a callback
function localfftShowData(obj,event)

if (get(obj,'SamplesAvailable') >= obj.SamplesPerTrigger)

	% Get the handles.
	data = obj.UserData;

	handles = data.figureHandles;

	% Execute a peekdata.
	x = peekdata(obj, obj.SamplesPerTrigger);

	% FFT calculation.
	Fs = obj.SampleRate;
	blockSize = obj.SamplesPerTrigger;
	[f,mag] = localDaqfft(x,Fs,blockSize);

	% Dynamically modify Analog axis as we go.
	maxX=max(x);
	minX=min(x);
	yax1=get(handles.axes1,'YLim');
	yax1(1)=minX - .0001; % need to subtract a value to make sure yax(1) never equals yax(2)
	yax1(2)=maxX + .0001;
	set(handles.axes1,'YLim',yax1)
	set(handles.axes1,'XLim',[0 (obj.SamplesPerTrigger-1)/obj.SampleRate])

	% Dynamically modify Frequency axis as we go.
	maxF=max(f);
	minF=min(f);
	xax=get(handles.axes2,'XLim');
    xax(1)=minF;
    xax(2)=maxF;
	set(handles.axes2,'XLim',xax)

	% Dynamically modify Magnitude axis as we go.
	maxM=max(mag);
	minM=min(mag);
	yax2=get(handles.axes2,'YLim');
	yax2(1)=minM - .0001;
	yax2(2)=maxM + .0001;
	set(handles.axes2,'YLim',yax2)

	% Update the line plots.
	set(handles.hLine1, 'XData', [0:(obj.SamplesPerTrigger-1)]/obj.SampleRate, 'YData', x(:,1));
	set(handles.hLine2, 'XData', f(:,1), 'YData', mag(:,1));

    % Find the frequency at which the max signal strength is at.
    [ymax,maxindex] = max(mag);
    set(handles.tFreq,'String',sprintf('%4.1d Hz',f(maxindex)));

    % Store the current FFT into the array of FFTs used for the waterfall.
	data.storedFFTs(data.storedFFTsIndex,:) = mag';

    % This circular shift is used so that when we display the 3D plot, the
    %  newest FFT will appear in 'front' and the oldest in 'back'.
    % To understand this, note how the plotting routines are using this fftOrder
    %  array to reorder the FFTs stored in data.storedFFTs and also note
    %  how data.storedFFTsIndex is used to store FFTs in data.storedFFTs.
    %
	fftOrder = 1:handles.numTraces;
	fftOrder = circshift(fftOrder,[ 1 -data.storedFFTsIndex ]);

	data.storedFFTsIndex = data.storedFFTsIndex + 1;
	if (data.storedFFTsIndex > handles.numTraces)
        data.storedFFTsIndex = 1;
        data.plotSurf        = 1; % Indicates a full history is stored.
	end

	% Update the surface plot if we have a full history.
	if (data.plotSurf)
        cla(handles.axes3);

        v=get(handles.poPlotType,'Value');
        s=get(handles.poPlotType,'String');
        switch s{v}
            case 'Classic'
                data.view  = [103 30];
                data=localClassic(handles,data,f,fftOrder);
            case 'Classic(Top)'
                data.view  = [90 -90];
                data=localClassic(handles,data,f,fftOrder);
            case 'Mosaic'
                data.view  = [90 -90];
                data=localMosaic(handles,data,f,fftOrder);
            case 'Waterfall'
                data.view = [103 30];
                data=localWaterfall(handles,data,f,fftOrder,yax2);
            case 'Rotate'
                data=localRotate(handles,data,f,fftOrder);
            case 'CycleAll'
                data=localCycleAll(handles,data,f,fftOrder);
        end

    end

	set(data.ai, 'UserData', data);

	drawnow;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ++++ V I S U A L I Z A T I O N ++++
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function data=localClassic(handles,data,f,fftOrder)
    [X,Y] = meshgrid(1:handles.numTraces,f(1:end));
    surf(X,Y,data.storedFFTs(fftOrder,:)','parent',handles.axes3);
  	set(handles.axes3,'XLim',[1 handles.numTraces],'YLim',[0 f(end)])
    shading(handles.axes3,'interp');
    set(handles.axes3,'View',data.view)

function data=localMosaic(handles,data,f,fftOrder)
    [X,Y] = meshgrid(1:handles.numTraces,f(1:10:end));
    surf(X,Y,data.storedFFTs(fftOrder,(1:10:end))','parent',handles.axes3);
  	set(handles.axes3,'XLim',[1 handles.numTraces],'YLim',[0 f(end)])
    set(handles.axes3,'View',data.view)

function data=localWaterfall(handles,data,f,fftOrder,yax2)
    [X,Y] = meshgrid(1:handles.numTraces,f(1:end));
    p=plot3(X,Y,data.storedFFTs(fftOrder,:)','parent',handles.axes3);

    % rotate the color map of the lines in the plot3
    map  = linspace(0,1,handles.numTraces);
    map2 = linspace(1,0,handles.numTraces);
    rotatemap  = map(fftOrder);
    rotatemap2 = map2(fftOrder);
    for k=1:handles.numTraces;
        set(p(k),'Color',[rotatemap(k) .1 rotatemap2(k)]);
    end

  	set(handles.axes3,'XLim',[1 handles.numTraces],'YLim',[0 f(end)]);
    shading(handles.axes3,'interp');
    set(handles.axes3,'View',data.view)

function data=localRotate(handles,data,f,fftOrder)
    [X,Y] = meshgrid(1:handles.numTraces,f(1:8:end));
    surf(X,Y,data.storedFFTs(fftOrder,(1:8:end))','parent',handles.axes3);
  	set(handles.axes3,'XLim',[1 handles.numTraces],'YLim',[0 f(end)])
    set(handles.axes3,'View',data.view)

    % Rotate the view point.
    data.view(1) = 90;
    if data.view(2) >= 90-data.rotateStep
        data.rotateStep = -4;
    elseif data.view(2) <= -90-data.rotateStep
        data.rotateStep = 4;
    end
    data.view(2) = data.view(2)+data.rotateStep;

function data=localCycleAll(handles,data,f,fftOrder)
    data.counter = data.counter + get(data.ai,'TimerPeriod');
    if data.counter > 5*handles.cycleTime
        data.counter = 0;
    elseif data.counter > 4*handles.cycleTime
        data=localRotate(handles,data,f,fftOrder);
    elseif data.counter > 3*handles.cycleTime
        data.view  = [103 30];
        data=localWaterfall(handles,data,f,fftOrder);
    elseif data.counter > 2*handles.cycleTime
        data.view  = [90 -90];
        data=localMosaic(handles,data,f,fftOrder);
    elseif data.counter > 1*handles.cycleTime
        data.view  = [90 -90];
        data=localClassic(handles,data,f,fftOrder);
    else
        data.view  = [103 30];
        data=localClassic(handles,data,f,fftOrder);
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ---- V I S U A L I Z A T I O N ----
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% --- Executes during object creation, after setting all properties.
function poPlotType_CreateFcn(hObject, eventdata, handles)
% hObject    handle to poPlotType (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

% --- Executes on selection change in poPlotType.
function poPlotType_Callback(hObject, eventdata, handles)
% hObject    handle to poPlotType (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = get(hObject,'String') returns poPlotType contents as cell array
%        contents{get(hObject,'Value')} returns selected item from poPlotType


text:
function cU = chebTcoeffs2chebUcoeffs(cT)
%CHEBTCOEFFS2CHEBUCOEFFS   Convert Chebyshev-T coeffs. to Chebyshev-U coeffs.
%   CU = CHEBTCOEFFS2CHEBUCOEFFS(CT) takes a column vector of coefficients of a
%   polynomial represented as an expansion in the Chebyshev polynomials of the
%   first kind T_k(x) (ordered starting with the coefficient of T_0(x)) and
%   returns the column vector of coefficients of the same polynomial expanded
%   in Chebyshev polynomials of the second kind U_k(x) (ordered starting with
%   the coefficient of U_0(x)).
%
%   If CT is a matrix, the conversion is performed column-wise.

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

if ( isempty(cT) )
    cU = [];
    return;
end

% We compute 2nd-kind coefficients from the 1st-kind coefficients
% using the recurrence T_n(x) = (1/2)*(U_n(x) - U_{n-2}(x)).  The
% coefficient of U_n requires the coefficients of T_n and T_{n + 2}, so we
% pad with two zero-coefficients initially.
nCols = size(cT, 2);
cU = [cT ; zeros(2,nCols)];

% Run the recurrence.
cU(1,:) = 2*cT(1,:);
cU = 0.5*[cU(1:end-2,:) - cU(3:end,:) ; cU(end-1:end,:)];
cU = cU(1:end-2,:);

end