MATLABを使った離散コサイン基底と離散コサイン変換入門No.1

1次元データ篇

Version 1.4
　
　
早稲田大学教育・総合科学学術院　新井仁之
　
　
概要　離散コサイン基底と離散コサイン変換を線形代数の基礎的な話を仮定して，丁寧に解説した．No. 1 では 1次元データを扱い，No.2 では2次元データ（画像データ）を扱う．
内容目次
  1. 記号と用語
  2. 離散コサイン基底
  3. 正規直交基底
  4. 離散コサイン基底によるフーリエ展開
  5. 離散コサイン基底から離散コサイン変換（DCT）
  6. 離散コサイン変換と逆離散コサイン変換
  7. データの解析例（東京の潮位）
  8. データ圧縮事始め
＊予備知識：線形代数の初歩．内積，基底，線形写像の定義をおおよそ知っていることが望ましい．
＊本ノートは MATLAB R2024a と MATLAB ライブエディターにより作成した．
＊MATLABを動かしながら，数学理論を修得することを目標とした．そのためプログラムは数学理論に沿うように作成した．
＊フーリエ解析ではインデックスは 0 から始めることが多い．MATLABのインデックスは 1 から始まるので要注意．
clear %とりあえず本ノートに必要ないものを最初にワークスペースから消しておく．
　

１．記号と用語

ℝを実数全体からなる集合とする．実数から成るN次元縦ベクトル を と略記する．また，実数から成る N 次元縦ベクトル全体のなす集合を により表す．特にすべての が 0 であるような N 次元縦ベクトルを と表し，ゼロベクトルという．
上で定義される次の演算を線形演算という．
線形演算
に対して
.
また，内積が次のように定義される．
内積
に対して
を と の内積という．特に のとき， と は直交するという．
% 内積のプログラム
% わざわざ作らなくてもよいが，数式にあわせるため．
function Ip = InnerProd(x,y)
Ip = x'*conj(y); % 本稿で扱うのは実数のみなので，conj は無くても良い．
end
「離散コサイン基底」のうち，「基底」という用語の説明をしておく．
定義1 （基底）
が の一組の基底であるとは，次の (B1), (B2) が成り立つことである．
(B1) 任意の は，ある実数 を選んで，必ず
　　　
　　と表すことができる．
(B2) は線形独立である．すなわち， をみたす
　　実数 は 0 に限る．
コメント　なお， が の基底ならば， になることが証明される．すなわち の基底をなすベクトルの個数は N 個である必要がある．
基底としてどのようなものがあるかというと，最もシンプルなものとしては，次の標準基底がある． として， と定義する． は の一組の基底であり，これは標準基底と呼ばれる．これから述べる離散コサイン基底は，標準基底とは異なる有用な基底である．
　

２．離散コサイン基底

離散コサイン基底は，次のように定義されるN個のN次元縦ベクトル からなるの基底である．ただし基底であることは第3節で確認する．証明するまでは基底と呼ぶべきではないかもしれないが，便宜上，証明する前から離散コサイン基底と呼ぶことにしたい．
,
に対して
.
これが実際に の基底であることは後で示すことにする．
まず離散コサイン基底をプログラムする．
Nを自然数とし，V=DCB(N) により，次の形で N 次元離散コサイン基底をまとめた行列が得られる．
=
function V = DCB(N)
% N 次元離散コサイン基底 V を行列形式（上述）で生成する．
V = zeros(N,N);
V(:,1)=1/sqrt(N)*ones(N,1);
for k = 2:N
for n = 1:N
V(n,k) = sqrt(2/N)*cos(pi/(2*N)*(2*(n-1)+1)*(k-1));
end
end
end
ここでは，の場合の離散コサイン基底の各ベクトルの数値のグラフを記述する．なおはJPEG圧縮で使われるケースである．
N=8;
v = DCB(N);
for l=1:2
for k=1:4
subplot(2,4,k+4*(l-1))
x=v(:,k+4*(l-1));
plot(x, 'o', 'LineWidth', 0.5, 'MarkerSize', 3, 'MarkerFaceColor', 'r')
 
for i = 1:N
line([i, i], [0, x(i)], 'LineStyle', '--', 'Color', 'k')
end
title([' v(:,' num2str(k+4*(l-1)-1) ') のグラフ'], 'FontSize', 8, ...
'FontWeight', 'bold')
xticks(1:N)
xticklabels(0:N-1)
 
xlim([1 N])
ylim([-.6 .6])
 
end
end

　

3. 正規直交性基底

ここで正規直交基底の概念を導入しておく．正規直交基底は非常に重要な基底のクラスである．その有用性は次節で説明する．まずは定義から始める．
定義 2（正規直交基底）
を の一組の基底とする．これが
(ON1)
をみたすとき， の一組の正規直交基底という．
ここで，与えられた正規直交しているベクトルが正規直交基底であるための有用な条件を述べておく．
定理 1
とする．もしもこれが (ON1) を満たしていれば， は の一組の正規直交基底になっている．
つまり，空間の次元と同じ個数のベクトルは，(ON1)をみたしていれば基底になっているということである．証明は [1] を参照してほしい．
したがって，離散コサイン基底が基底になっていることは，
=
を示せば十分である． の場合に，これをMATLABで確認しておこう．一般の場合の証明は例えば文献[3]参照．
（これまで離散コサイン基底と読んできたが，ここにおいて，基底であることが証明されたことになる．）
正規直交性の確認（ の場合）
N=8;
InnerProduct = zeros(N,N);
for k=1:N
for l=1:N
InnerProduct(k,l)=InnerProd(v(:,k),v(:,l));
end
end
disp('(v^(k)とv^(l)の内積)=');disp(InnerProduct);
(v^(k)とv^(l)の内積)= 1.0000 0 -0.0000 0.0000 0.0000 -0.0000 -0.0000 0.0000 0 1.0000 0.0000 -0.0000 0.0000 0.0000 0.0000 -0.0000 -0.0000 0.0000 1.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 0.0000 -0.0000 0.0000 1.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 0.0000 -0.0000 0.0000 1.0000 0.0000 -0.0000 -0.0000 -0.0000 0.0000 0.0000 -0.0000 0.0000 1.0000 0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 1.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 -0.0000 0.0000 -0.0000 1.0000

　

4. 離散コサイン基底によるフーリエ展開

線形代数の一般論として，次の定理が知られている．
定理 2. X を N次元線形空間とし，内積 を持つものとする．を正規直交基底とする．このとき，任意の に対して，
が成り立つ．
したがって特に離散コサイン基底に対しても成り立つ．
このとき，縦ベクトル ( ' は転置) を の基底による表現という．特に，本稿では離散コサイン基底による表現を周波数表現と呼ぶ．
　

5. 離散コサイン基底から離散コサイン変換（DCT）

離散コサイン基底と離散コサイン変換の関係を明確にしておく．離散コサイン基底 を定理1に当てはめると，次の展開式が得られる．
.
ここでこの展開式を次の図のように二つのフェーズに分けて考え，前半の部分からDCT（離散コサイン変換），後半の部分から IDCT（逆離散コサイン変換）が定義される．
Decomp-Syntheesis_DCT_MATLAB.jpg
　　

6. 離散コサイン変換と逆離散コサイン変換

ここでは上記の数学的な考え方を元に離散コサイン変換と逆離散コサイン変換のプログラムを作成する．数学理論に忠実に計算していくので，高速計算は考慮に入れていない．
% 離散コサイン変換（数学理論に沿ってプログラムしたもので，高速計算ではない）
function DC_X = DCTrans(x)
x_size = length(x);
v = DCB(x_size);
DC_X = zeros(x_size,1);
for k = 1:x_size
DC_X(k) = InnerProd(x,v(:,k));
end
end
% 逆コサイン離散変換（数学理論に沿ってプログラムしたもので，高速計算ではない）
function IDC_X = IDCTrans(x)
x_size = length(x);
v = DCB(x_size);
IDC_X= zeros(x_size,1);
for k = 1:x_size
IDC_X=IDC_X+x(k)*v(:,k);
end
end
コメント　簡単な検算
簡単な例で，検算と MATLAB の dct との比較をしておく．
x=[1 2 3 4 5 6 7 8]';
まず本稿の数学理論重視に沿って作った（高速計算を考慮してない）プログラムの場合．
tic;DCX = DCTrans(x); toc
経過時間は 0.001521 秒です。
MATLAB の dct の場合．
tic; MATLAB_dct = dct(x); toc
経過時間は 0.064718 秒です。
MATLAB の dct と等しいことの確認
er1=max(abs(DCX-MATLAB_dct));
disp(['誤差 =',num2str(er1)])
誤差 =2.3981e-14
本稿での DCT の反転公式の確認
z = IDCTrans(DCX)
z = 8×1
1.0000 2.0000 3.0000 4.0000 5.0000 6.0000 7.0000 8.0000

　

7. データの解析例（東京の潮位）

ここでは2024年4月の東京の潮位のデータを離散コサイン変換を使って解析する．
データは気象庁のホームページからダウンロードできる．
https://www.data.jma.go.jp/kaiyou/db/tide/genbo/index.php
まずこの中から，「東京，2024年4月，観測基準面表示，毎時潮位」を選択してデータを取り込む．それをExcelファイルにしたものが「東京潮位202404.xlsx」である．このExcelファイルは，上記ホームページにある表を取り込んだもので，行が時間（1時～24時），列が日（1日～30日）を表す 行列，720個の数値からなっている．
filename = '東京潮位202404.xlsx';
data1 = readmatrix(filename);
%第1行は「時間」，第1列は「日」の表示なのでそれを省いておく．
data2 = data1(2:end,2:end);
%DCT変換するためにそれを1列のデータ（縦ベクトル）に並べ替える
data3 = reshape(data2',1,[]);
%このデータをExcelファイルに変換する．（今は必要ないが，後で使う場合を考えて．）
writematrix(data3,'data202404.xlsx');
 
%潮位のデータのグラフ
figure
plot(data3,'k','LineWidth', 1)
xticks([1 25 49 73 97 121 145 169 193 217 241 265 289 313 337 361 385 409 433 457 ...
481 505 529 553 577 601 625 649 673 697]);
xticklabels({'1','2','3','4','5','6','7','8','9','10','11','12','13','14','15','16', ...
'17','18','19','20','21','22','23','24','25','26','27','28','29','30',});
xlim([1 length(data3)])
xlabel('日', 'FontSize', 8)
ylabel('潮位 (cm)', 'FontSize', 8)
title('東京 2024年4月 （気象庁によるデータに基づく）', 'FontSize', 8, 'FontWeight', 'bold')
% 300dpiで保存
print('Tokyo_Tyoui_2024_4', '-dtiff', '-r300');
 
tic; DX = DCTrans(data3'); toc %本稿の数学理論重視に沿って作った（高速計算を考慮してない）DCT
経過時間は 0.021878 秒です。
tic;MATLAB_DX = dct(data3); toc %MATLAB の dct
経過時間は 0.002317 秒です。
両者の誤差
er2= max(abs(DX - MATLAB_DX'));
disp(['誤差 =',num2str(er2)])
誤差 =9.2312e-11
DX_delete = DX;
DX_delete(1)=0; %直流成分は大きいので，直流成分は 0 に置き換えて表示する．
%表示上の変更なので，DCTの計算への影響はなし．
plot(DX_delete,'k')
xlim([1 length(DX)])
title('東京 2024年4月（気象庁によるデータ）のDCTによる周波数表現', 'FontSize', 8, 'FontWeight', 'bold')
print('Tokyo_Tyoui_Freq_2024_4', '-dtiff', '-r300');
%潮位のデータのグラフ
figure
subplot(2,1,1)
plot(data3,'k','LineWidth', 1)
xticks([1 25 49 73 97 121 145 169 193 217 241 265 289 313 337 361 385 409 433 457 ...
481 505 529 553 577 601 625 649 673 697]);
xticklabels({'1','2','3','4','5','6','7','8','9','10','11','12','13','14','15','16', ...
'17','18','19','20','21','22','23','24','25','26','27','28','29','30',});
xlim([1 length(data3)])
xlabel('日', 'FontSize', 8)
ylabel('潮位 (cm)', 'FontSize', 8)
title('東京 2024年4月 （気象庁によるデータに基づく）', 'FontSize', 8, 'FontWeight', 'bold')
 
subplot(2,1,2)
DX_delete = DX;
DX_delete(1)=0; %直流成分は大きいので，直流成分は 0 に置き換えて表示する．
%表示上の変更なので，DCTの計算への影響はなし．
plot(DX_delete,'k')
xlim([1 length(DX)])
title('東京 2024年4月（気象庁によるデータ）のDCTによる周波数表現', 'FontSize', 8, 'FontWeight', 'bold')
print('Tokyo_Tyoui_TimeFreq_2024_4', '-dtiff', '-r300');
この中で，特に離散コサイン変換値の絶対値が大きい順に，その大きな値をとる位置を抜き出す．
[~,max_perm] = sort(abs(DX),'descend');
そのうち，上位12番目までを抜き出す．
max_perm(1:12)'
ans = 1×12
1 117 121 56 62 118 116 61 114 58 64 120
すると，120付近と60付近にDCTの絶対値が大きいところが集中していることがわかる．離散コサイン基底（変換）の定義の仕方から，これは潮位の波形の周波数成分が，
　　　 時間周期と， 時間周期の部分
が大きいことを意味している．ところで
潮汐には月や太陽の影響で，一日周期と半日周期があることが知られている．離散コサイン変換による解析と一致していることがわかる．
　

8. データ圧縮事始め

そこで，主要な周波数成分だけを使って，潮位のグラフがどの程度再現できるかを見てみよう．これはデータ圧縮の側面からの考察でもある．1番絶対値の大きなの周波数成分から作った波形が Approx.1，1番目からk番目までの絶対値の大きな周波数成分から作った波形を Approx.k とする．つまり，原信号を復元するには， の720個の周波数データが必要ところを Approx.k の波形を得るには k 個の周波数データで済むということである．
以下のプログラムでは，Approx.12 までを描いている．ただし，グラフを見やすくするため，表示の範囲は 0～360 までに限定している．
% ここでは理論に沿って離散コサイン基底を使ってプログラムしたが，MATLABの dct と idct を使っても可能．
L=60; %近似に使うデータの個数 後で使う
G=360; %グラフの表示範囲
N=length(DX);
V = DCB(N);
A=zeros(N,1);
Z=abs(DX);
ss=zeros(1,L);
figure
subplot(3,6,1)
plot(data3,'k')
axis([0 G 0 300])
title('original')
for k=1:L
s=find(Z==max(Z),1);
ss(k)=s;
A=A+DX(s)*V(:,s);
if k<=12
subplot(3,6,k+6)
plot(A,'k')
axis([0 G 0 300])
title(['Approx.' num2str(k)])
end
Z(s)=0;
if k==12
AA =A;
end
end
Approx.60 と元データのグラフを書いて比較してみる．60まで使うとかなり元データを近似できていることが視認できるであろう．本来ならば 720 個のデータが必要なところ，60個のデータのグラフの概形がわかることになる（ただし概形は似ていても，潮位の観測データとは異なるので注意）．
figure
subplot(3,1,1)
plot(AA,'k')
axis([0 720 0 300])
title(['Approx.' num2str(12)])
subplot(3,1,2)
plot(A,'k')
axis([0 720 0 300])
title(['Approx.' num2str(L)])
subplot(3,1,3)
plot(data3,'k')
axis([0 720 0 300])
title('Original Data')
これは「ロスありデータ圧縮」の観点からの考察で，潮位のデータとしての意味は考えてない．
ロスありデータ圧縮は，音声データや画像など人の感覚にとって有用なデータの圧縮に使われる．これについては No. 2で画像データ圧縮を題材に解説する．

　

参考文献

線形代数について．（本稿の定理1，2 の証明は[1]にある．潮位の「フーリエ基底」による解析は[2]参照．)
[1] 新井仁之，線形代数 基礎と応用，日本評論社，2006 (2017年 update版は Kindle）．
[2] 新井仁之，基底とその応用 - 配列の数学と信号処理 (1) -，数学のたのしみ（日本評論社）2006年秋号，pp.101-126.
離散コサイン基底について．
[3] 新井仁之，フーリエ解析学，朝倉書店，2003.
2024/07/03 追記
No.2 画像データ圧縮篇を公開しました。
http://www.araiweb.matrix.jp/MATLAB/DCT_Tutorial_using_MATLAB2.html
-----------------------------------------
本ノートの履歴
Version 1.1 公開日 2024/06/19
Version 1.2 公開日 2024/06/20 a minor change
Version 1.3 公開日 2024/06/21 a minor change
Version 1.4 公開日 2024/08/02 minor changes
------------------------------------------
Copyright © 2024 by Hitoshi Arai.