🔦
【離散コサイン変換】Part2. 欠損データをはっきりくっきり?復元する

2025/07/01に公開
 はじめに前回の記事で、離散コサイン変換(DCT)の仕組みと特徴を説明しました。

DCTをゆっくりじっくり理解する
今回はその続きです。

内容としてはDCTの特徴を応用して、計測とか制御にどう使えるのか、ということを考えていきます。

具体的にはタイトルにも書いていますが、スパース性を利用して欠損データからの信号復元、いわゆる圧縮センシングについて見ていこうと思います。
※記事中で使用したPythonとMATLABスクリプトの全体コードは下記に置いています。

https://github.com/otter1602/DCT_and_Min_L1norm_Sample

 DCTの特徴おさらいですが、DCTの定義式と特徴を振り返っておきます。
DCT(タイプ2)の定義式は以下でした。

X[k]=\sum_{n=0}^{N-1}x[n]cos(\frac{π}{N}k(n+\frac{1}{2}))
詳しくは前回の記事にて解説していますが、DCTによって周波数領域でcos波に分解された信号は下記の特徴を持っています。　※どんな信号でも必ず、ではありません。少なくとも傾向としては、です。
特徴(1)：エネルギーが低周波に集中する

特徴(2)：同じ信号を離散フーリエ変換した場合より少ない非ゼロのDCT係数で正確に表現できる
今回はこの2つの特徴を信号復元に利用していきます。

 信号復元の考え方DCTが持つ2つの特徴から、次のように言うことができます。

「高周波成分があまり含まれない滑らかな信号をDCTした場合、DCT係数はスパースになりやすい」
DCTはDFTと同様に逆変換(IDCT)が存在します。なので、少し乱暴な表現ですが上記を逆に言うと、

「DCT係数がスパースでも元信号が滑らかならば、IDCTによって元信号を比較的正確に復元できる」
と言うことができます。
IDFTの場合、フーリエ係数が欠損していると基本的に元の信号を正確に復元することはできません。

IDCTでも同様に、DCT係数が欠損していると厳密な復元はできませんが、元信号が滑らかでDCT係数がスパースになる性質を持っていれば、IDCTによって元信号の良い近似が得られる可能性があります。
なので今回のテーマである「DCTを応用して信号復元を行う」を言い換えると、次の通りになります。

すなわち、「良さそうなDCT係数を求め、IDCTによって元信号を近似する」

 問題設定早速、良さそうなDCT係数を求める方法について考えていきましょう。

必要な作業としては、第1に問題設定、次にそれを解く方法を考えることですね。
まず問題設定のために前提条件をいくつか出してみます。
＜前提条件＞

(1)：元信号はDCTによって係数がスパースになることが期待されるような信号である

　　　⇒つまり、滑らかでゆっくりな時間変化を持つような信号です。
(2)：観測した信号には欠損がある(大体30%ぐらいしか観測できていない)

　　　⇒例えばロギング機器の設定ミスや外乱ノイズ等、ありがちな状況ですね。
(3)：欠損した部分の値は補間が必要となるため、信号に対する何らかの構造的な仮定が必要である

　　　⇒単純な補間ではなく構造の仮定とそれに対する結果を見たいということです。
(4)：信号に対する構造的な仮定はできるだけ少ない変数にまとめたい

　　　⇒構造の仮定といっても、オッカムの剃刀よろしく不要な複雑さは避けるべきですね。
(5)：復元によって得たいのは「精度の良いデータ」よりも「信号の構造に対する知見」である

　　　⇒一番強調したいところです。構造が分かればデータの精度も上げられるかもしれません。
ここで、(1),(4)は信号の性質に、(2),(3)は問題設定の方向性に、そして(5)は問題の解法にそれぞれ関わってくるところです。
結論から言うと、(1),(4)によってL1ノルムを使うこと、(2),(3),(5)によって最適化問題として定式化することが有効だと言えます。
L1ノルムについて次節で説明しますが、ここで最適化問題として定式化することの意味を説明します。
(2),(3)より、今回取得した観測データから直接に信号の構造を考えるのは難しいです。

なので、観測データを使用して、「(何らかの意味で)自然と思われる値を取る」 という操作が必要になります。
そして、「自然である」ということは言い換えると、少なくともある観点からは 「最適である」 と考えることができます。これが今回最適化問題として定式化することの意味です。
まとめると、今回の問題設定は次のようになります。

「L1ノルムを最小化するような最適化問題を解くことで、スパース性のあるDCT係数を求めよう！」
続いて、なぜL1ノルムなのか、最適化問題として定式化して解けるのかについて説明していきます。

 L1ノルム最適化L1ノルムを使う理由を説明する前に、信号処理で頻出のノルムを確認しておきます。
・L0ノルム (L0疑似ノルム)　※厳密にはノルムではないが、慣例的にノルム扱い

　⇒ベクトルの非ゼロ要素の個数を数えた結果を表す。

　　L0ノルムの定義式：||x||_0 = \#\{i \mid x_i \ne 0\}
・L1ノルム

　⇒ベクトルの各要素の絶対値の総和を表す。

　　L1ノルムの定義式：||x||_1 = \sum_{i=1}^n|x_i|
・L2ノルム

　⇒ベクトルの各要素の二乗和平方根を表す。

　　L2ノルムの定義式：||x||_2 = (\sum_{i=1}^nx_i^2)^\frac{1}{2}
・L∞ノルム

　⇒ベクトルの各要素の絶対値中の最大値を表す。

　　L∞ノルムの定義式：||x||_∞ = \max_{1≦i≦n} |x_i|
一応、一般化してLpノルムとした場合は次のように書くことができます。
・一般化(Lp)ノルム

　⇒Lpノルムの定義式：||x||_p = (\sum_{i=1}^n|x_i|^p)^\frac{1}{p}
さて、上記L0~L∞までのノルムのなかで、ベクトルのスパース性を表す指標として最も適しているのはどのノルムでしょうか？
定義から考えると、答えはL0ノルムですね。一応、数値上からも確かめてみましょう。

次のような2種類のベクトルを用意します。

(1)：{2, 0, 0}

(2)：{1, 1, 1}
この場合、よりスパースなベクトルは非ゼロ要素が1つしか存在しない(1)の方ですね。

つまり 「(1)に対して計算した値が(2)に対して計算した値よりも小さくなるようなノルム」がベクトルのスパース性を表す指標として、より優れていると言えます。
では、それぞれのベクトルについてL0~L∞ノルムを求めてみます。

(1)：||x||_0 = 1, 　||x||_1 = 2, 　||x||_2 = \sqrt{4} = 2, 　　　||x||_∞ = 2

(2)：||x||_0 = 3, 　||x||_1 = 3, 　||x||_2 = \sqrt{3} \approx 1.732, 　||x||_∞ = 1
L0ノルムとL1ノルムは(1)に対して計算した方がより小さい値になりますが、L2ノルムとL∞ノルムは(2)に対して計算した方がより値が小さくなっています。
したがって、直感では最も直接的にスパース性を表すことができるL0ノルムが、実際に指標として適していることが数値上からも分かります。
じゃあL0ノルムを最小化するような問題を設定して、それを解けば良いんじゃないの？と思ってしまいますが、実際はL0ノルムは計算上の課題を抱えています。
本節の始めに述べた通り、L0ノルムは個数を数えるという組み合わせ論的な値です。そのため次元が大きくなると組み合わせ爆発が起こり、問題のサイズが急激に大きくなってしまうため、現実的な時間で最適化問題を解くことが困難になります。
※100次元のベクトルから50個の要素を非ゼロにする問題を解く場合、組み合わせは10^{29}通りです！
そこでL0ノルムとならんで(1)に対して計算した方がより小さい値になった、L1ノルムの出番です。
L1ノルムはベクトルのスパース性を表す指標としてL0ノルムと同様に適しており、組み合わせ論的な値でもありません。なので、L1ノルムはL0ノルムを数学的に扱いやすく、緩やかに近似したものと考えることができ、実際の定式化にはL0ノルムの近似としてL1ノルムを使用するのが定石です。
これが前節で出た疑問その1、「なぜL1ノルムを使うのか？」に対する答えです。
加えてさらに重要な事実として、L1ノルムは凸関数になります。

凸関数についての説明は省略しますが、簡単に言うと初期値がどうあれ、大域的な最適解に収束することが保証されている関数です。要するに凸関数の最適化問題は解けるに決まっているのです。
これが前節で出た疑問その2、「最適化問題として定式化して解けるのか？」に対する答えです。
それではここから、「L1ノルムを最小化する最適化問題」を定式化していきましょう！

 最適化問題の定式化まず各記号を下記の通り定義します。

x：元信号の列ベクトル

y：観測信号の列ベクトル

M：観測信号に欠損を作るためのマスク行列(0,1の対角行列)

D：DCT変換行列(つまり、DCTで使用する基底ベクトル(cos~)を並べた行列)

c：DCT係数

その他：行列Aの転置をA^Tと書く
これらを使って、観測された信号を表してみましょう。
観測信号yは欠損マスクMを元信号xにかけたものとして、y=Mxですね。
そして、元信号xを行列DによってDCTした結果がDCT係数であるcなので、cを逆DCTするとxになると考えると、x=D^Tcになりますね。

すなわち、観測信号yはy=Mxとx=D^Tcを使って、y=MD^Tcとcを使った形で表せます。
今回の問題設定を改めて振り返ると、

「L1ノルムを最小化するような最適化問題を解くことで、スパース性のあるDCT係数を求めよう！」

ということなので、目的関数はDCT係数であるcのL1ノルムであり、制約条件にcを含んだ形の式であるy=MD^Tcを使って定式化します。
つまり最終的な最適化問題は、

\mathbf {min　||c||_1\  \ subject\ to\ \ y=MD^Tc}

と書くことができます。
定式化自体はこれで完了ですが、せっかくなので実際に数値を入れて行列形式でDCTを計算し、IDCTで元に戻るところを見てみましょう。※元信号の長さN=4とします。
まず、 元信号x = \left(\begin{matrix}1 \\2 \\3 \\4 \end{matrix}\right)とします。これに対して、
DCTタイプ2の定義式：X_k=\sum_{n=0}^{N-1}x[n]\cdot cos(\frac{π}{N}k(n+\frac{1}{2}))から
DCTタイプ2の変換行列　D_{k,n}=\begin{cases}\sqrt\frac{1}{N}\text{ 　　　　　　　　if}　 k = 0 \\ \sqrt \frac{2}{N} \cdot cos(\frac{π(2n+1)k}{2N})\text{　　if}　 k > 0  \end{cases}　となり、計算すると、
D \approx \begin{bmatrix} 0.5 & 0.5 & 0.5 & 0.5 \\ 0.653 & 0.271 & -0.271 & -0.653 \\ 0.5 & -0.5 & -0.5 & 0.5 \\ 0.271 & -0.653 & 0.653 & -0.271 \\ \end{bmatrix}　が出てきます。これがDCT変換行列ですね。
※D_{k,n}のkが行方向のインデックス, nが列方向のインデックスで考えています。
ここで、x=D^Tcからc=D \cdot xなので、
c = D \cdot x \approx \begin{bmatrix} 5.0 \\ -2.230 \\ 0.0 \\ -0.158 \end{bmatrix}
でDCT係数cが分かります。これがDCT変換行列Dを元信号xにかけてDCT係数cが出てくる形ですね。
IDCTの場合は先ほどとは逆にx=D^Tcを計算すると、
x = D^T \cdot c\approx \begin{bmatrix} 1 \\ 2 \\ 3 \\ 4 \end{bmatrix}
まさしく、逆変換なのでDCT係数cから元信号xが計算できます。
※一連の計算は下記のコードで確認できます。　PythonとMATLABでDCT処理の仕様が違うため、

　PythonコードにはMATLAB版と同じ仕様になるコードも載せています。

　ScipyのDCTは行優先で処理するようで、そのままだとDが転置された形で出てきます。

　そのため、ここまでの説明通りの式にする場合はDを転置する必要があります。
＜Pythonコード＞
import numpy as np
from scipy.fft import dct

x = np.array([1, 2, 3, 4])  #元信号

D = dct(np.eye(4), type=2, norm='ortho')    #DCT変換行列
D2 = D.T                                    #DCT変換行列 (MATLAB形式)

c = D @ x   #DCT係数
c2 = D2 @ x #DCT係数 (MATLAB形式)

x_rec = D.T @ c     # IDCT
x_rec2 = D2.T @ c2  # IDCT (MATLAB形式)

print("DCT basis:", D)                      #基底ベクトル
print("DCT basis(MATLAB):", D2)             #基底ベクトル (MATLAB形式)
print("DCT coefficient:", c)                #DCT係数
print("DCT coefficient(MATLAB):", c2)       #DCT係数 (MATLAB形式)
print("Original signal:", x_rec)            #IDCT結果
print("Original signal(MATLAB):", x_rec2)   #IDCT結果 (MATLAB形式)
＜MATLABコード＞
x = [1; 2; 3; 4]; %元信号
D = dct(eye(4), type=2); %DCT変換行列
c = D*x; %DCT係数
x_rec = D'*c; %IDCT

disp("DCT basis:")
disp(D)
disp("DCT coefficient:")
disp(c)
disp("Original signal:")
disp(x_rec)

 最適化問題を解く長々説明しましたが、後は解くだけですね。定式化した最適化問題を再掲します。

\mathbf {min　||c||_1\  \ subject\ to\ \ y=MD^Tc}
この問題を解くわけですが、今回はMATLABの場合はcvx, Pythonの場合はcvxのPython版であるcvxpyを使用します。
cvx及びcvxpyは凸最適化問題に特化したモデリングツールで、数式を数式通りにプログラムに落とし込めばそのまま問題を解いてくれる優れものです。
一応軽く説明すると、実際に最適化問題を解くのはcvxやcvxpyのようなモデリングツールではなく、ソルバーと呼ばれるプログラムで、色々な種類があります。しかし直接ソルバーへ問題を渡して解いてもらおうとすると、問題に対して適切なソルバーを選んだり、ソルバーが理解できる形に式変形したりとなかなか面倒です。
そこでモデリングツールの出番なわけですが、元の数式に近い形のままモデリングツールへ問題を渡すと、問題に応じて適切なソルバーを選びそのソルバーにふさわしい形に変形して渡してくれます。

つまりユーザーとソルバーの橋渡しをしてくれるのがモデリングツールということです。
コード全文は別途あげますが、定式化の記述部分だけ抜粋すると次のようになります。

最適化問題の式　\mathbf {min　||c||_1\  \ subject\ to\ \ y=MD^Tc}　をコードで再現しているのが見て取れるかと思います。
＜Pythonコード＞
# 観測点数
N = 128

# 観測マスク（30%観測）
mask_bool = np.random.rand(N) < 0.30
    M = np.diag(mask_bool.astype(float))#マスクするための対角行列

# DCT行列
D = dct(np.eye(N), type=2, norm='ortho')

#定式化ブロック#
def reconstruct(x_true):
    #変数の定義#
    y_observed = M @ x_true     #y = M・x
    c = cp.Variable(N)          #cを未知の変数とする
    x_reconstructed = D.T @ c   #x = D^T・c
    
    #制約条件と目的関数を設定し、それらを使って定式化#
    constraints = [M @ x_reconstructed == y_observed] #制約条件：y = M・x 
    objective = cp.Minimize(cp.norm1(c))         #目的関数：変数cのL1ノルムを最小化する
    prob = cp.Problem(objective, constraints) #設定した目的関数と制約条件で定式化
    
    #凸最適化問題を解く#
    prob.solve()
    return x_reconstructed.value, c.value   #xとcを返す
＜MATLABコード＞　※やってることはPythonコードと同じ
% 観測点数
N = 128;

% 観測マスク（30%観測）
mask_bool = rand(1, N) < 0.3;
M = diag(double(mask_bool));  % 対角行列

% DCT行列
D = dct(eye(N));

% 定式化ブロック
function [x_rec, c] = reconstruct(x_true, M, D)
    y_obs = M * x_true(:);
    N = length(x_true);
    %定式化開始
    cvx_begin quiet 
        variable c(N)
        x_est = D' * c;
        minimize(norm(c, 1))
        subject to
            M * x_est == y_obs
    cvx_end

    x_rec = x_est;
end
今回は単調なsin波を元信号とした場合と、ランダムなノイズを元信号とした場合のそれぞれに対して70%の欠損、すなわち30%しかデータが取得できなかった観測を行ったと仮定して、30%のデータから元信号を復元できるか試してみましょう。

 信号復元結果Pythonで実行してみた結果を下記に示します。sin波の方は\frac{1}{128}[rad/sample]と\frac{2}{128}[rad/sample]の重ね合わせで、ノイズ信号の方はガウスノイズです。

上段がsin波、下段がノイズ信号として、それぞれ元信号と復元信号、観測点を図中に記載しています。

どうでしょうか？sin波の方は30%しかデータを取得できていないにもかかわらず、かなりの精度で復元できています。一方でノイズ信号の方はさっぱり復元できていないですね。
この結果はDCT係数がスパース性を持つか否かが復元結果を左右するということを示しています。

すなわち、低周波成分しか持たないsin波の方はDCT係数がスパース性を持つためかなりの精度で復元可能ですが、ガウスノイズの方は全帯域に渡って成分が存在するためDCT係数がスパースにならず、復元精度が低下する、ということですね。
もう少し観測点を減らして(=欠損数を多くして)様子を見てみましょう。

結果を下記に示します。
<25%観測>



<20%観測>



<15%観測>

25%までは良い精度で復元できていますが、20%以下の観測点になると精度はかなり低下してしまっています。DCT係数がスパースであるならば精度良く復元できるとすると、20%以下の観測点の場合はDCT係数もスパースにできていないのでしょうか？確認してみましょう。
<元信号のDCT係数(つまり、100%観測の場合)>

上段がsin波の方のDCT係数で下段がノイズ信号のDCT係数になります。元信号からはsin波のDCT係数はスパース性を持っており、ノイズ信号は広い帯域に係数がばらけて存在することが分かります。
<30%観測>



<20%観測>



<15%観測>

予想通り、sin波のDCT係数がスパースで無くなるにつれて復元精度が低下するのが見て取れます。

 まとめ今回はDCT係数のスパース性を利用してL1ノルム最小化による、欠損データからの信号復元の基本的な考え方について見てきました。
実際半分以上のデータが欠損している状態から、それなりの精度で意味のある信号が復元できるというのは非常に面白く、また応用の利く内容だと思います。
しかし、私の周囲でもそうですが取得したデータに欠損があった場合やサンプリング周波数が足りない場合等、真っ先にやろうとするのが線形補間やスプライン補間だと思います。
実際それらの手法は現場で非常に役に立つものですが、それでは補間と信号復元の違いは一体どこにあるのでしょうか？

わざわざ最適化問題の定式化->求解という割と面倒な手順を踏んでまで信号復元を行うメリットとは？
次回の記事ではこの疑問について考えるとともに、信号のスパース性を活かした別の応用例について、もう少し見ていこうと思います。