MATLABの拡張ツールボックス関数を自前で実装する

はじめに

MATLABを使っていると、便利な関数が追加のツールボックスに含まれていて、ライセンスの関係で使えないことがよくあります。この記事では、そんな頻出関数を自前で実装する方法をまとめました。

完全に同じ機能を再現することは難しいですが、基本的な用途では十分使える実装を目指しています。

1. 標準化（zscore関数の代替）

Statistics and Machine Learning Toolboxに含まれるzscore関数は、データの標準化（Zスコア変換）を行います。

function [z, mu, sigma] = my_zscore(x, flag, dim)
    % データを標準化する関数（Zスコア変換）
    % x: 入力データ
    % flag: 0または省略で標本標準偏差、1で母標準偏差
    % dim: 演算を行う次元
    
    % デフォルト引数の処理
    if nargin < 2 || isempty(flag)
        flag = 0;
    end
    if nargin < 3 || isempty(dim)
        % 最初の非シングルトン次元を見つける
        dim = find(size(x) ~= 1, 1);
        if isempty(dim)
            dim = 1;
        end
    end
    
    % 平均と標準偏差を計算
    mu = mean(x, dim);
    if flag == 0
        sigma = std(x, 0, dim);  % 標本標準偏差
    else
        sigma = std(x, 1, dim);  % 母標準偏差
    end
    
    % 標準化
    z = (x - mu) ./ sigma;
    
    % 標準偏差が0の場合の処理
    z(sigma == 0) = 0;
end

使用例：

% テストデータ
data = randn(100, 3) * 10 + 50;

% 標準化
[z_data, mu, sigma] = my_zscore(data);

% 確認
disp(['平均: ', num2str(mean(z_data))]);  % ほぼ0
disp(['標準偏差: ', num2str(std(z_data))]);  % ほぼ1

2. 移動平均フィルタ（movmean関数の代替）

Signal Processing Toolboxに含まれるmovmean関数は、時系列データの移動平均を計算する際によく使います。

function y = my_movmean(x, k, varargin)
    % 移動平均を計算する関数
    % x: 入力データ
    % k: ウィンドウサイズ（スカラーまたは[kb kf]の形式）
    % varargin: 'Endpoints'オプション（'shrink', 'fill', 'discard'）
    
    % デフォルト設定
    endpoints = 'shrink';
    
    % オプション引数の解析
    if nargin > 2
        for i = 1:2:length(varargin)
            if strcmpi(varargin{i}, 'Endpoints')
                endpoints = varargin{i+1};
            end
        end
    end
    
    % ウィンドウサイズの処理
    if isscalar(k)
        kb = floor(k/2);
        kf = floor((k-1)/2);
    else
        kb = k(1);
        kf = k(2);
    end
    
    n = length(x);
    y = zeros(size(x));
    
    for i = 1:n
        % ウィンドウの範囲を決定
        start_idx = max(1, i - kb);
        end_idx = min(n, i + kf);
        
        % 平均を計算
        y(i) = mean(x(start_idx:end_idx));
    end
    
    % エンドポイントの処理
    switch lower(endpoints)
        case 'discard'
            y = y((kb+1):(n-kf));
        case 'fill'
            % 端の値で埋める
            y(1:kb) = y(kb+1);
            y((n-kf+1):n) = y(n-kf);
    end
end

3. 正規化（normalize関数の代替）

データの正規化を行うnormalize関数の実装です。

function [y, C, S] = my_normalize(x, method, varargin)
    % データを正規化する関数
    % x: 入力データ
    % method: 'zscore', 'norm', 'range', 'scale', 'center'
    % varargin: メソッド固有のオプション
    
    if nargin < 2
        method = 'zscore';
    end
    
    switch lower(method)
        case 'zscore'
            % Zスコア正規化
            [y, C, S] = my_zscore(x);
            
        case 'norm'
            % ノルム正規化（デフォルトは2-ノルム）
            p = 2;
            if ~isempty(varargin)
                p = varargin{1};
            end
            
            if isvector(x)
                norm_val = norm(x, p);
                y = x / norm_val;
                C = 0;
                S = norm_val;
            else
                % 列ごとに正規化
                norm_vals = sqrt(sum(x.^p, 1));
                y = x ./ norm_vals;
                C = zeros(1, size(x, 2));
                S = norm_vals;
            end
            
        case 'range'
            % 範囲正規化 [0, 1]
            min_val = min(x);
            max_val = max(x);
            range_val = max_val - min_val;
            
            % ゼロ除算を避ける
            range_val(range_val == 0) = 1;
            
            y = (x - min_val) ./ range_val;
            C = min_val;
            S = range_val;
            
        case 'scale'
            % 標準偏差でスケーリング
            S = std(x);
            S(S == 0) = 1;  % ゼロ除算を避ける
            y = x ./ S;
            C = 0;
            
        case 'center'
            % 中心化（平均を引く）
            C = mean(x);
            y = x - C;
            S = 1;
            
        otherwise
            error('Unknown normalization method');
    end
end

4. ピーク検出（findpeaks関数の代替）

Signal Processing Toolboxのfindpeaks関数は、信号のピークを検出する際に重宝します。

function [pks, locs] = my_findpeaks(x, varargin)
    % 信号のピークを検出する関数
    % x: 入力信号
    % varargin: オプション引数
    
    % デフォルト値の設定
    minPeakHeight = -inf;
    minPeakDistance = 1;
    minPeakProminence = 0;
    threshold = 0;
    
    % オプション引数の解析
    for i = 1:2:length(varargin)
        switch lower(varargin{i})
            case 'minpeakheight'
                minPeakHeight = varargin{i+1};
            case 'minpeakdistance'
                minPeakDistance = varargin{i+1};
            case 'minpeakprominence'
                minPeakProminence = varargin{i+1};
            case 'threshold'
                threshold = varargin{i+1};
        end
    end
    
    % ピーク候補の検出
    n = length(x);
    peak_indices = [];
    
    for i = 2:n-1
        % 局所最大値かチェック
        if x(i) > x(i-1) && x(i) > x(i+1)
            % 高さの条件
            if x(i) >= minPeakHeight
                % しきい値の条件
                if x(i) - max(x(i-1), x(i+1)) >= threshold
                    peak_indices = [peak_indices, i];
                end
            end
        end
    end
    
    % プロミネンスの計算と適用
    if minPeakProminence > 0 && ~isempty(peak_indices)
        prominences = zeros(size(peak_indices));
        for i = 1:length(peak_indices)
            idx = peak_indices(i);
            % 左側の最小値
            left_min = min(x(1:idx));
            % 右側の最小値
            right_min = min(x(idx:end));
            % プロミネンス
            prominences(i) = x(idx) - max(left_min, right_min);
        end
        keep = prominences >= minPeakProminence;
        peak_indices = peak_indices(keep);
    end
    
    % 最小ピーク間隔の適用
    if ~isempty(peak_indices) && minPeakDistance > 1
        % ピークを高さでソート
        [~, sort_idx] = sort(x(peak_indices), 'descend');
        sorted_indices = peak_indices(sort_idx);
        
        % 間隔をチェックしながら選択
        selected = false(size(sorted_indices));
        for i = 1:length(sorted_indices)
            if ~any(abs(sorted_indices(i) - sorted_indices(selected)) < minPeakDistance)
                selected(i) = true;
            end
        end
        
        peak_indices = sort(sorted_indices(selected));
    end
    
    locs = peak_indices;
    pks = x(peak_indices);
end

5. 移動標準偏差（movstd関数の代替）

移動窓での標準偏差を計算する関数です。

function y = my_movstd(x, k, flag)
    % 移動標準偏差を計算する関数
    % x: 入力データ
    % k: ウィンドウサイズ
    % flag: 0で標本標準偏差（デフォルト）、1で母標準偏差
    
    if nargin < 3
        flag = 0;
    end
    
    n = length(x);
    y = zeros(size(x));
    
    % ウィンドウの半分のサイズ
    half_k = floor(k/2);
    
    for i = 1:n
        % ウィンドウの範囲を決定
        start_idx = max(1, i - half_k);
        end_idx = min(n, i + half_k);
        
        % 標準偏差を計算
        window_data = x(start_idx:end_idx);
        y(i) = std(window_data, flag);
    end
end

6. 外れ値除去（rmoutliers関数の代替）

統計的手法で外れ値を検出・除去する関数です。

function [y, TF, lower_bound, upper_bound] = my_rmoutliers(x, method, threshold)
    % 外れ値を除去する関数
    % x: 入力データ
    % method: 'median'（デフォルト）, 'mean', 'quartiles'
    % threshold: しきい値（デフォルトは3）
    
    if nargin < 2
        method = 'median';
    end
    if nargin < 3
        threshold = 3;
    end
    
    switch lower(method)
        case 'median'
            % 中央値絶対偏差（MAD）法
            med = median(x);
            mad = median(abs(x - med));
            % スケール因子（正規分布の場合）
            c = -1 / (sqrt(2) * erfcinv(3/2));
            lower_bound = med - threshold * c * mad;
            upper_bound = med + threshold * c * mad;
            
        case 'mean'
            % 平均値と標準偏差を使用
            mu = mean(x);
            sigma = std(x);
            lower_bound = mu - threshold * sigma;
            upper_bound = mu + threshold * sigma;
            
        case 'quartiles'
            % 四分位範囲（IQR）法
            q1 = prctile(x, 25);
            q3 = prctile(x, 75);
            iqr = q3 - q1;
            lower_bound = q1 - threshold * iqr;
            upper_bound = q3 + threshold * iqr;
    end
    
    % 外れ値の判定
    TF = x < lower_bound | x > upper_bound;
    
    % 外れ値を除去
    y = x(~TF);
end

7. バターワースフィルタ（butter関数の簡易版）

Signal Processing Toolboxのbutter関数は、バターワースフィルタの設計に使われます。

function [b, a] = my_butter(n, Wn, type)
    % バターワースフィルタの設計（簡易版）
    % n: フィルタ次数（1または2のみサポート）
    % Wn: 正規化カットオフ周波数（0-1）
    % type: 'low'（デフォルト）, 'high'
    
    if nargin < 3
        type = 'low';
    end
    
    % プレワーピング
    w = 2 * tan(pi * Wn / 2);
    
    switch n
        case 1
            % 1次フィルタ
            if strcmpi(type, 'low')
                k = w / (1 + w);
                b = k * [1, 1];
                a = [1, (w - 1)/(w + 1)];
            else  % high
                k = 1 / (1 + w);
                b = k * [1, -1];
                a = [1, (w - 1)/(w + 1)];
            end
            
        case 2
            % 2次フィルタ
            if strcmpi(type, 'low')
                k = w^2 / (1 + sqrt(2)*w + w^2);
                b = k * [1, 2, 1];
                a = [1, ...
                     2*(w^2 - 1)/(1 + sqrt(2)*w + w^2), ...
                     (1 - sqrt(2)*w + w^2)/(1 + sqrt(2)*w + w^2)];
            else  % high
                k = 1 / (1 + sqrt(2)*w + w^2);
                b = k * [1, -2, 1];
                a = [1, ...
                     2*(w^2 - 1)/(1 + sqrt(2)*w + w^2), ...
                     (1 - sqrt(2)*w + w^2)/(1 + sqrt(2)*w + w^2)];
            end
            
        otherwise
            error('この実装では1次または2次のみサポートしています');
    end
end

8. 百分位数（prctile関数の代替）

Statistics and Machine Learning Toolboxのprctile関数の基本的な実装です。

function y = my_prctile(x, p)
    % データの百分位数を計算する関数
    % x: 入力データ（ベクトル）
    % p: 百分位数（0-100）
    
    % ベクトルに変換
    x = x(:);
    
    % NaNを除去
    x = x(~isnan(x));
    
    % ソート
    x_sorted = sort(x);
    n = length(x_sorted);
    
    % 百分位数の計算
    y = zeros(size(p));
    
    for i = 1:length(p)
        if p(i) == 0
            y(i) = x_sorted(1);
        elseif p(i) == 100
            y(i) = x_sorted(end);
        else
            % 線形補間
            k = (n - 1) * p(i) / 100 + 1;
            k_floor = floor(k);
            k_ceil = ceil(k);
            
            if k_floor == k_ceil
                y(i) = x_sorted(k);
            else
                % 線形補間
                alpha = k - k_floor;
                y(i) = (1 - alpha) * x_sorted(k_floor) + alpha * x_sorted(k_ceil);
            end
        end
    end
end

まとめ

これらの実装は、MATLABの標準関数と完全に同じ動作をするわけではありませんが、多くの場合で実用的に使えます。特に以下の点に注意してください：

エラーチェックや境界条件の処理は最小限です
計算効率は最適化されていません
一部の高度なオプションは実装していません
多次元配列への対応は限定的です

実際の使用では、これらの実装をベースに、必要に応じて機能を追加・改良することをお勧めします。また、可能であれば適切なツールボックスのライセンスを取得することも検討してください。

参考リンク