はじめに
MATLABで、例えばfor文などでplot関数を何度も呼び、異なったデータを重ねて1枚の図にすることがよくあるかと思います。その際、5個や10個程度のデータなら色で判別できるのですが、20、30、もっと…とグラフを重ねると色での判別はつきにくくなります。そうした場合でも色を良い感じに分けて、可視化しやすくするためのMATLAB関数を自作しました。
具体的には、下記の図を見てもらうと分かりやすいと思います。下記は1枚に40個のグラフが載ってますが、10個ごとに①乱数、②ほぼ一定値、③指数的に減衰する関数、④三角関数、の4種類を表示しています。当然色をごっちゃにした左はそれぞれが全く区別できませんが、本記事で紹介する自作MATLAB関数を使うことで10個ごとの色のパターンを指定し、識別が容易になっています。
何をもってして良い感じの色分けか…は難しい話ですが、本記事での注意事項を先に載せておきます。
- 人の色覚や個人差に配慮した色彩理論…とかそういった話には踏み込みません。HSV色空間での彩度・明度は固定で色相だけ変える、ということやっていきます
- 信号の特徴量を使って自動的にグループ分けする…とかではありません。描画対象のグループ分けがあらかじめ想定出来ている前提での内容です
作った関数や実行例のMATLABコードは下記におきました。
R2021a、R2023bあたりで動作確認しました。
間違いや、もっとこうするといいよ、などあればコメントいただけると助かります。
使い方
まず、適当に描画したいデータを用意します。
やり方なんでも良いですが例えばこのように
Fs = 100; % サンプリング周波数
T = 1/Fs; % サンプリング時間
L = 100; % 信号の長さ
t = (0:L-1)*T; % 時間ベクトル
data = zeros(40, L);
rng("default")
% 1. 乱数
for i = 1:10
data(i, :) = randn(1, L);
end
% 2. 一定値
for i = 11:20
data(i, :) = 1 + 0.05*randn(1, L);
end
% 3. 0へ減少する指数関数
decay_constants = linspace(0.5, 5, 10); % 10種類の減衰定数
for i = 21:30
data(i, :) = exp(-decay_constants(i-20)*t) + 0.05*randn(1, L);
end
% 4. 正弦波
frequencies = linspace(0.5, 5, 10); % 10種類の周波数
for i = 31:40
data(i, :) = sin(2*pi*frequencies(i-30)*t) + 0.05*randn(1, L);
end
データを用意したうえで、本記事の主題の自作関数setColorPtn()を下記のように使用します。この関数に、「色分け対象が10個ずつ4グループである」という指定を意味するcolor_listを引数として渡すと、色配列cptnを生成します。MATLABの標準関数colororder()にこの色配列を渡すことで、プロットの色をカスタマイズできます。
color_list = [10 10 10 10];
[cptn,cptn_idx]= setColorPtn(color_list);
% [cptn,cptn_idx]= setColorPtn(color_list,[],[],1/3); %冒頭の画像と同じ色にする場合
colororder(cptn);
for i = 1:40
plot(t,data(i,:),'LineWidth',2)
hold on;
end
title('色配列指定')
hold off;
関数の解説
↓に関数の中身を置いておきます。これをsetColorPtn.mとして使います。
引数のエラーチェックやコーナーケースの挙動は怪しいので大目に見ていただきつつ、自由に改変して使ってもらえればと思います。
function [cptn,cptn_idx] = setColorPtn(item_list,num_ptn,color_div,color_border)
arguments
item_list = [];
num_ptn = 1;
color_div = 1;
color_border = 2/3;
end
% item_listが一次元配列であることのチェック
if ~isempty(item_list)
if ~isvector(item_list)
error('item_list は一次元配列（行ベクトルまたは列ベクトル）でなければなりません。');
end
end
hsv_num = 256; % HSV色相256色を使用
% 色ベクトルの指定
hsv_vec = hsv(hsv_num);
elem = [];
if isempty(item_list) & isempty(num_ptn) & isempty(color_div)
% いずれの引数も指定されてない場合、エラー
error('少なくとも一つの入力を指定してください。item_listか、num_ptn、またはcolor_divのいずれかが必要です。');
elseif ~isempty(item_list) % 引数でリスト指定されている場合を想定
num_ptn = sum(item_list);
color_div = length(item_list);
elem = item_list;
else % 引数でプロット総数と色の分割数が指定されている場合を想定
if color_div==0 % ゼロ割保護
error('色の分割数 color_div はゼロに設定できません。');
end
base_elem_size = floor(num_ptn / color_div);
remainder = mod(num_ptn, color_div);
elem = ones(1, color_div) * base_elem_size;
elem(1:remainder) = elem(1:remainder) + 1;
end
color_idx = [];
color_ptn = [];
for i=1:color_div
idx_start = 1+floor(hsv_num/color_div)*(i-1);
idx_end = idx_start + floor( (hsv_num/color_div)*(1-color_border) );
color_idx{i} = floor(linspace(idx_start, idx_end, elem(i)));
color_ptn{i} = hsv_vec(color_idx{i},:);
end
cptn = vertcat(color_ptn{:});
cptn_idx = horzcat(color_idx{:})';
end
しつこく書いてあるコメント等まで見たい方はこちら
function [cptn,cptn_idx] = setColorPtn(item_list,num_ptn,color_div,color_border)
% SETCOLORPTN 色を適切に分割しつつグループ化することで、colororder で使用するカラーパターンを生成します。
%
% 構文：
% [cptn, cptn_idx] = setColorPtn(item_list, num_ptn, color_div, color_border) % [cptn, cptn_idx] = setColorPtn(item_list, num_ptn, color_div, color_border)
%
% 説明：
% この関数は、色配列を作成することにより、プロットで使用する色パターンを生成します。
% この関数は、2 種類の入力構成を取ることができます：
% 1. item_listでプロットとグループの数を直接指定します。
% 例えば、3つのグループに分け、最初のグループから1色、2番目から3色、3番目から3色を選択するには、[1 3 3]を使用します。
% 2. または、プロットの総数 (num_ptn) とグループの数 (color_div) を指定します。この場合、色はグループ間で均等に分割されます。
% 割り切れない数は、適当にいずれかのグループに分配します。コード内の記述を参考にしてください。
%
% 両方の構成が指定された場合、関数は'item_list'によって提供された構成を優先します。
%
% 入力引数：
% - item_list : 各グループの色の数を指定するベクトル。指定された場合、他の入力を上書きします。(デフォルト: [])
% - num_ptn : (オプション) 生成する色の総数。item_list' が指定されていない場合に使用されます。(デフォルト: 1)
% - color_div : (オプション) 色を分割する色グループの数。item_list' が指定されていない場合に使用されます。(デフォルト: 1)
% - color_border : (オプション) ボーダー領域とみなす色グループの割合 (0 から 1)。
% 0にするとボーダー領域なしを意味し、1に近い場合ほぼすべての色がボーダー領域で、ごく限られた領域から色を取得します。(デフォルト： 2/3)
% つまり'color_border'パラメータは、グループをどの程度区別して表示するかを制御します。値が小さいとグループ間の区別が少なくなり、値が大きいと同じグループ内の色の区別が少なくなります。
%
% 出力引数:
% - cptn : プロットの色を指定するために colororder で使用できる色配列。
% - cptn_idx : HSV 色空間から選択された色のインデックス。
%
% 例1. item_list を使用して色を指定します：
% [cptn, cptn_idx] = setColorPtn([1, 4, 5]); % 10個のプロットを3グループに分けたカラーパターン生成
% colororder(cptn);
%
% 例2. num_ptn と color_div を使用して色を指定します：
% [cptn, cptn_idx] = setColorPtn([], 10, 3); % 10個のプロットを3グループに分けたカラーパターン生成
% colororder(cptn);
arguments
item_list = [];
num_ptn = 1;
color_div = 1;
color_border = 2/3;
end
% item_listが一次元配列であることのチェック
if ~isempty(item_list)
if ~isvector(item_list)
error('item_list は一次元配列（行ベクトルまたは列ベクトル）でなければなりません。');
end
end
hsv_num = 256; % HSV色相256色を使用
% 色ベクトルの指定
hsv_vec = hsv(hsv_num);
elem = [];
if isempty(item_list) & isempty(num_ptn) & isempty(color_div)
% いずれの引数も指定されてない場合、エラー
error('少なくとも一つの入力を指定してください。item_listか、num_ptn、またはcolor_divのいずれかが必要です。');
elseif ~isempty(item_list) % 引数でリスト指定されている場合を想定
num_ptn = sum(item_list);
color_div = length(item_list);
elem = item_list;
else % 引数でプロット総数と色の分割数が指定されている場合を想定
if color_div==0 % ゼロ割保護
error('色の分割数 color_div はゼロに設定できません。');
end
base_elem_size = floor(num_ptn / color_div);
remainder = mod(num_ptn, color_div);
% 例えば
% num_ptn=9, color_div=3 なら、elem=[3, 3, 3]に分配
% num_ptn=10, color_div=3 なら、elem=[4, 3, 3]に分配
% num_ptn=11, color_div=3 なら、elem=[4, 4, 3]に分配
% num_ptn=12, color_div=3 なら、elem=[4, 4, 4]に分配
elem = ones(1, color_div) * base_elem_size;
elem(1:remainder) = elem(1:remainder) + 1;
end
color_idx = [];
color_ptn = [];
for i=1:color_div
idx_start = 1+floor(hsv_num/color_div)*(i-1);
idx_end = idx_start + floor( (hsv_num/color_div)*(1-color_border) );
color_idx{i} = floor(linspace(idx_start, idx_end, elem(i)));
color_ptn{i} = hsv_vec(color_idx{i},:);
end
cptn = vertcat(color_ptn{:});
cptn_idx = horzcat(color_idx{:})';
% % テスト用 HSVのどの部分を取ってきているか可視化
% scatter(cptn_idx/hsv_num,linspace(0,0,length(cptn)),100, cptn, 'filled');
% hold on;
% scatter(linspace(0,1,hsv_num),linspace(0.2,0.2,hsv_num),100, hsv_vec, 'filled');
% text(0.2,0.1,sprintf('↑から%dグループ分を抽出したもの↓',color_div),'FontSize',14)
% hold off;
end
やっていることのイメージとしては、下記のコードを実行すると分かりやすいかと思います。setColorPtn()は引数の与え方で2通りの使い方ができるようにしており、通常の使い方の場合は2番目・3番目の引数は無視されます(念のため「[ ]」で空の配列を渡しています)。
color_list = [5 4 3 2];
color_group = length(color_list);
[cptn,cptn_idx]= setColorPtn(color_list,[],[],2/3); % 色配列を指定。2/3をボーダー領域として指定
% 描画用カラーパターン生成のテスト
scatter(cptn_idx/256,linspace(0,0,length(cptn)),100, cptn, 'filled');
hold on;
scatter(linspace(0,1,256),linspace(0.2,0.2,256),100, hsv(256), 'filled');
text(0.2,0.1,sprintf('↑から%dグループ分を抽出したもの↓',color_group),'FontSize',12)
実行すると、下記の画像が表示されます。Y軸の高さに意味はなく、適当にずらしてあるだけです。
さらに、この色配列生成関数で大事な「ボーダー領域」という概念も図示すると下記です。
要するに、
- HSVカラーマップ配列をもとに、配列color_listのサイズ(この例だと4)で、色のグループ分けをする
- ボーダー領域(この例だと2/3)を除いたところから、配列color_listで指定した個数分の色を均等に取ってくる
ということをしています。これによって、同一グループの色系統は近く、異なるグループの色系統は異なるものにしています。
したがってボーダー領域とは、異なるグループ間での色の区別をはっきりさせるか、同一グループ内での色の区別をはっきりさせるか、のトレードオフを調整するパラメータとなっています。ボーダー領域を大きくすると前者寄りに、ボーダー領域を小さくすると後者寄りになります。
別の使い方
色配列を渡すのではなく、「プロットの総数」と「色のグループ分け数」を指定して色分けする使い方もできます。その場合、1番目の引数に空の配列「[ ]」を渡して使います。
plot_num = 10;
color_group = 3;
[cptn,cptn_idx]= setColorPtn([],plot_num,color_group,2/3); % プロット総数と色分け数を指定。2/3をボーダー領域として指定
% 描画用カラーパターン生成のテスト
scatter(cptn_idx/256,linspace(0,0,length(cptn)),100, cptn, 'filled');
hold on;
scatter(linspace(0,1,256),linspace(0.2,0.2,256),100, hsv(256), 'filled');
text(0.2,0.1,sprintf('↑から%dグループ分を抽出したもの↓',color_group),'FontSize',12)
結果は下図となります。プロット総数10は3グループで割り切れないので、適当に4,3,3 と割り当てるようにしています。(詳しくはコード内のコメントを確認ください）
色配列を指定する、というもとの使い方とわけている理由としては
- 1グループのみでいいが、数十～百ぐらいのグラフを順番がわかる形で描画したい
- 全体をざっくり5色に分けたい。端数はどっちに配色されても気にしない
などの場合がしばしばあるので作りました。
活用例
以上で話したかった内容はほぼ終わりですが、一つだけ具体的な活用例を紹介します。(この例では、Control System Toolboxを使います)
詳細は割愛しますが、伝達関数の極ゼロの位置でステップ応答の様子が変わる、という内容をプロットしてみます。極$p_1$の位置で4つのグループ分け(色分け)して描画すると、まさに極$p_1$の位置によって波形が大きく変わっている…という様子がわかるかと思います。こういうのがやりたかったわけです。
clear;
impos = -10; %虚部は適当に大きめにしておく
p1 = -1;
p2 = -2 + impos*[1i -1i];
z1 = -3;
z2 = -4 + impos*[1i -1i];
tf = zpk([z1 z2],[p1 p2],1);
pzmap(tf)
pos_pattern = [1, 2, 3, 4];
pos_perm = -perms(pos_pattern);
legends = cell(length(pos_perm),1);
color_group = 4;
[cptn,cptn_idx]= setColorPtn([],length(pos_perm),color_group,2/3);
colororder(cptn);
for pos = 1:length(pos_perm)
p1 = pos_perm(pos,1);
p2 = pos_perm(pos,2) + impos*[1i -1i];
z1 = pos_perm(pos,3);
z2 = pos_perm(pos,4) + impos*[1i -1i];
all_poles = [p1, p2];
all_zeros = [z1, z2];
% 定常ゲインを1に調整
numerator = prod(-all_zeros);
denominator = prod(-all_poles);
gain_adjusted = denominator / numerator;
tf = zpk([z1 z2],[p1 p2],gain_adjusted);
[yout, tout] = step(tf,10);
plt = plot(tout,yout);
legends{pos} = "p1=" + string(p1) + ", p2=" + string(real(p2(1))) + ", z1=" + string(z1) + ", z2=" + string(real(z2(1)));
% 最大値が1+εを超えなかったものを表示
if(max(yout) < 1+0.0001)
display(legends{pos});
end
% データヒントの追加
plt.DataTipTemplate.DataTipRows(end+1).Label = legends{pos};
hold on;
end
ylim([0 2])
lgd = legend(legends);
lgd.ItemHitFcn = @hitcallback_ex1; %凡例クリックしたら表示/非表示
title('極ゼロの位置の違いによるステップ応答の変化')
hold off;
% 凡例クリックでグラフ表示/非表示切り替え
function hitcallback_ex1(src,evnt)
if strcmp(evnt.Peer.Visible,'on')
evnt.Peer.Visible = 'off';
else
evnt.Peer.Visible = 'on';
end
end
その他補足
上記のグラフで使った、たくさんグラフを描く場合に知っておくと便利な内容を紹介します。
凡例をクリックしたときにグラフ表示/非表示の切り替え
グラフをたくさん描くと、「どの凡例がどのグラフと対応してるんだっけ？」となりがちです。下記のMATLAB Answersにて、凡例をクリックしたときに対象のグラフを表示⇔非表示と切り替えてくれる自作の関数を紹介されています。(公式で標準的に入れておいて欲しい機能ですね…)
下記のような関数を定義しておいて、
function hitcallback_ex1(src,evnt)
if strcmp(evnt.Peer.Visible,'on')
evnt.Peer.Visible = 'off';
else
evnt.Peer.Visible = 'on';
end
end
下記のようにグラフを描画した後、凡例の定義をした後にコールバック関数として登録します。
% 何らかのプロットをしておき、legendsに凡例の文字列を入れておく
lgd = legend(legends);
lgd.ItemHitFcn = @hitcallback_ex1; %凡例クリックしたら表示/非表示
こうしておくと、凡例の内容をクリックしたとき、対応するグラフの表示⇔非表示が切り替わります。これによって、凡例がどのグラフと対応しているかがインタラクティブにわかります。
あわせて下記のような関数も定義しておくと、コマンド ウインドウ上で「showAllPlots」、「hideAllPlots」と入力するだけでアクティブなFigureのグラフをまとめて表示/非表示にできて便利です。
function showAllPlots()
ax = gca; % 現在の軸を取得
children = ax.Children; % すべてのプロットオブジェクトを取得
for i = 1:length(children)
children(i).Visible = 'on'; % 表示に設定
end
end
function hideAllPlots()
ax = gca; % 現在の軸を取得
children = ax.Children; % すべてのプロットオブジェクトを取得
for i = 1:length(children)
children(i).Visible = 'off'; % 非表示に設定
end
end
データヒントの追加
逆に、グラフの波形をみて「このデータはどの設定で描画したものだっけ？」となりがちなので、データヒントに表示項目を追加しておくと便利です。コード内の該当箇所を抽出すると下記の部分です。
% 何らかのプロットをしておき、pltとしてLineプロパティを格納しておく。legendsに凡例の文字列を入れておく
plt.DataTipTemplate.DataTipRows(end+1).Label = legends{pos};
こうしておくと、グラフにデータヒントを表示したとき、付加的な情報も表示することができます。
詳しくは公式の↓のドキュメントも参照ください。
おわりに
この記事では、MATLABで描いたたくさんのグラフをグループ分けし、色でだいたい傾向を把握する…というのに役立つ自作関数を紹介しました。こうするともっと簡単にできるよ、とかあったらぜひ教えてください。
ちなみに、本記事で紹介した色分けはHSV色空間で行いましたが、やたらと緑は見分けがつきにくい気がします。特に5グループとか取ってくると、たいがい緑系から２色取ってくることになって区別がつきづらいです。自作関数setColorPtn()内の下記の記述を変えるとカラーマップが変更できるので、自由に改変して試してみてください。
% 色ベクトルの指定
hsv_vec = hsv(hsv_num);
例えばjetにしてみると下記のように変わります。どちらも良し悪しありますね。
