【MATLAB】gnuplotのカラーマップをmatlabで再現する

はじめに

12月はエンジニアの皆様が1日1人記事を投稿してクリスマスを指折り数える文化（Advent Calendar）があるが，そういうのは無視して投稿することにする．(スキルに自信がないのでそういうのに参加するのは躊躇してしまうため)

今回はgnuplotカラーマップをmatlabにて再現し，スペクトログラムを描画したいと思う．

このカラーマップが欲しい

仕事でmatlabを用いてスペクトログラムを描画することがたまにあり，カラーマップはいつもjetを使っている．インターネットでスペクトログラムなどで検索すると，下図のような紫～黄色のカラーマップで作成された画像をよく見るが，どうやらこのカラーマップはmatlabには無いようなのである．

現在のカラーマップの表示と設定 - MATLAB colormap

調べてみたところ，どうやらこのカラーマップはmatplotlibではgnuplotという名前で入っているらしい．（gnuplotとは昔の描画ソフトのことで，このソフトで用いられてきたカラーマップなのだと思うが，コンピュータの歴史は詳しくないのでその辺は分からないです）
Choosing Colormaps in Matplotlib

カッコいいカラーマップがmatplotlibでは使えるのにmatlabで使えないのは悔しいので，自分で実装してみることにする．

画像を拝借

上記のサイトよりカラーマップの画像をダウンロードして，MSペイントを駆使してgnuplotのカラーバー部分をトリミングした．

こいつをmatlabで読み取ればいいんじゃねえのかしら．
倫理的にアアアではあるが・・・

clc
clear all
close all

% 画像からカラーマップを生成する
IM = imread('gnuplot_colormap.png');% IMの型はuint8かuint16である
IM = squeeze(IM(1,:,:));% IMを1行分だけ抽出する
% colormapは0~1の範囲なのでuintの最大値で割って正規化する
myColormap = double(IM) / double(intmax(class(IM))); 


% 分析する信号
[y, Fs] = audioread('異国のシグナル.mp3');
y = y(:,1);% 1chだけ用いる

% STFTする
Nw = 512;                   % 時間窓の長さ
Nol = round(Nw/2);          % オーバーラップ
NFFT = Nw;                  % FFT点数
spectrogram(y,Nw,Nol,NFFT,Fs,'yaxis');  % スペクトログラムを描画
caxis([-120 -40]);          % レベルの範囲を指定
colormap(myColormap);       % カラーマップを設定

結果

うむ，期待通りの画像が得られている．

カラーマップを解析

前の節で終わってもいいのだが，せっかくなので，読み取ったカラーマップ画像のRGBのグラフを描いてみる．

clc
clear all
close all

IM = imread('gnuplot_colormap.png');
IM = squeeze(IM(1,:,:));
myColormap = double(IM) / double(intmax(class(IM)));
% ここまでは上と同じ

% RGBのグラフを描く
figure
c = ['r'; 'g'; 'b'];
L = length(myColormap);
n = ((1:L)-1)/L;% 横軸を[0 1]で正規化
for ind = 1:3
    plot(n, myColormap(:,ind), c(ind));
    hold on
end
legend(c);
title('カラーマップ（画像より生成）のRGB値');
hold off

結果

結果は割と分かりやすい感じになった．
これなら手で実装できそうである．多分こういうことでしょ．

{\rm RED} = \sqrt{x}\\
{\rm GREEN} = {x}^3\\
{\rm BLUE} = \sin({2 \pi x}) \ \ (0 \le x \le 0.5), \ \ 0 \ \ ({\rm else})\\

カラーマップを手実装

r = sqrt(n);
g = n.^3;
b = sin(2*pi*1*n);
b(round(length(b)/2):end) = 0;
myColormap2 = [r', g', b'];
figure
n = ((1:L)-1)/L;
for ind = 1:3
    plot(n,myColormap2(:,ind), c(ind));
    hold on
end
legend(c);
title('カラーマップ（手実装）のRGB値');
hold off

結果

うむ，いい感じ．
これを重ね書きしてみる．

figure
for ind = 1:3
    plot(n,myColormap(:,ind), c(ind));
    hold on
end
for ind = 1:3
    plot(n,myColormap2(:,ind), c(ind), 'LineStyle', ':');
    hold on
end
legend({'R(画像)', 'G(画像)', 'B(画像)', 'R(手実装)', 'G(手実装)', 'B(手実装)'});

完全に一致した．ノルムの比較とかはしなくてもいいでしょ．
ヨシ！

まとめ

最後に，手実装したカラーマップでスペクトログラムを描いてみる．

% 分析する信号
[y, Fs] = audioread('異国のシグナル.mp3');
y = y(:,1);

% STFTする
Nw = 512;                   % 時間窓の長さ
Nol = round(Nw/2);          % オーバーラップ
NFFT = Nw;                  % FFT点数
figure
spectrogram(y,Nw,Nol,NFFT,Fs,'yaxis');  % スペクトログラムを描画
caxis([-120 -40]);            % レベルの範囲を指定
% ここまでは同じ

colormap(myColormap2);               % 手実装したcolormapを指定する

所望の結果が得られた．

ちなみにこちらはいつも使っているjetである．
見やすさは正直jetのほうが良いと個人的には思うが，気分で使い分けていきたいと思う．