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の標準関数と完全に同じ動作をするわけではありませんが、多くの場合で実用的に使えます。特に以下の点に注意してください：
- エラーチェックや境界条件の処理は最小限です
- 計算効率は最適化されていません
- 一部の高度なオプションは実装していません
- 多次元配列への対応は限定的です
実際の使用では、これらの実装をベースに、必要に応じて機能を追加・改良することをお勧めします。また、可能であれば適切なツールボックスのライセンスを取得することも検討してください。