入力ベクトルをxとした時の,離散コサイン変換は
で定義される.また逆離散コサイン変換(IDCT)は次の式で定義される.

これだけではプログラムは組めないので.行列計算ができる形になるように考える.
DCTは(そしてIDCTも)有限個(N個)のデータに対しての変換なので,入力データはN点ごとに分割する必要がある.通常は,2のべき乗の長さをとり,一般的には8個が良く用いられる.以下では8個の場合で話を進める.
![]()
![]()

一方逆変換は

で示される.変換行列Dは,逆行列と転地行列が同一の直交行列である.ということで,一旦行列Dが生成できれば離散コサイン変換に関する実験を行うことができる.
の数式でXは,DCT係数とも呼ばれる.DCTでは相関の強い信号は,低周波数領域に信号の電力が集中する傾向がある.従ってDCT係数でも電力が集中する係数と集中しない係数が存在する.

DCTを用いたデータ圧縮とは,この電力集中を応用して,電力の集中する領域に多くのビット数を,電力が集中しない領域では少ないビット量を割り当てることで行う.
その最も単純な方式が,高周波領域の一部が0に近いことを利用して,0値近似してしまう(すなわちこの部分のデータを切り捨てる)ことで行う方法である.復元する場合は,切り捨てた部分に0の値を挿入してIDCTを行う.本来DCTは可逆変換であるが,この場合は切り捨てたデータを近似的に補完することになるので不可逆変換(オリジナルデータが復元されない変換)となる.今回のプログラミングでは,この0値近似による復元データへの影響を観察する.
今回のデータ処理の流れ
実際の(JPEG/MPEGなど)のデータ圧縮では,元データにDCTを行った後で,量子化及び符号化を行っている点は,承知しておいてほしい.
このページでは,画像などで圧縮率を高くすると発生するノイズとして,下記が挙げられている.
1) MATLAB用関数のスケルトン
% function dctt1() % function for discrite cosine transfer traning1 % Input: none % Output: Graph out % By your_name, your class, build date % % function y=dctt1()2) 特徴的なパラメータの設定,ここではDCT変換係数行列中で何個を0にするか,と変換行列のサイズを変更可能にするための定数を定義しておく.
znum=0; %この値が0のままだと0値近似の個数が0なのでデータ圧縮はされないことになる3) 変換する対象のデータ.この例では整数倍の波長の正弦波を幾つか合成したものを用いている.このデータは,高周波成分が少ないので,高周波成分を切り捨てる圧縮方式の影響を受けにくいことが予想される.高周波成分が多い別の種類のデータも後で生成させて実験する.
clear x; for(n=1:1:32) x(n)=sin(pi*n/64)+sin(pi*n/32)+sin(pi*n/16)+sin(pi*n/8); end4) 入力されるデータxを,nsize行の行列に変換する.reshape(x,m,n)は第1引数をm行n列の行列に変形する.この形式では端数のデータは切り捨てられる.
x1=reshape(x,nsize,round(length(x)/nsize));5) ここまで作成したら実行させて,x及びx1の値を確認すること.
6) 変換行列の生成
for(n=0:1:nsize-1) for(k=0:1:nsize-1) d(k+1,n+1)=cos(pi*k*(2.0*n+1.0)/(nsize*2)); if(k==0) d(k+1,n+1)=d(k+1,n+1)/sqrt(2.0); end end end d=d*sqrt(2.0/nsize);7) x1の各列に対してDCTを行うDCT係数を得る.
y1=d*x1; y2=y1;8) 各DCT係数の中で,高周波成分の中のnsizeで指定した個数だけ0にする.
y2((nsize-znum+1):nsize,:)=zerosy2の各列の下側の行がnsizeの個数だけ0になっていることを確認すること.
9)次にy2に対してDCT変換行列の転地行列をかけてIDCT(逆離散コサイン変換)を実行する.
x2=d'*y210) これだけでは,結果の判定や圧縮の効果を理解しにくいので,グラフ化を行う.先に示したデー手処理の流れに従って,元データ,DCT係数,0値を入れたDCT係数,復元したデータからなるグラフを図と同じ形式で一括して表示する.このためにsubplot関数を用いる.subplot(m,n,k) はグラフ画面をm行n列に区切り,その中のk番目を使用することを指定する.順番は左上から1,右隣が2という形で指定される.
※ x(:)は複数次元の配列から,1列目,2列目と順に取り出して1次元の列ベクトルを作り出す.
subplot(2,2,1); plot(1:1:length(x1(:)), x1(:)); subplot(2,2,2); plot(1:1:length(x1(:)), y1(:)); subplot(2,2,4); plot(1:1:length(x1(:)), y2(:)); subplot(2,2,3); plot(1:1:length(x1(:)), x2(:));
左下の復元データと実データを重ねて表示するためには,subplot(2,2,3)での描画部分を下のplot関数に差し替える.左下のグラフ中でオリジナルデータが○で表示される.
plot(1:1:length(x1(:)),x1(:),'o',1:1:length(x2(:)),x2(:)); grid onプログラムのパラメータの値をグラフのタイトル中に示すためには,例えば
strtitle1=['Higher ',num2str(znum/nsize*100),'% are set to 0']; strtitle2=[num2str(100-znum/nsize*100),'% compressed result'];という形で,タイトル用文字列を生成しておき,グラフを描画する際のタイトル指定に変数として入力することで指定できる.
subplot(2,2,4); plot(1:1:length(x1(:)), y2(:));grid on;title(strtitle1); subplot(2,2,3); plot(1:1:length(x1(:)),x1(:),'o',1:1:length(x2(:)),x2(:)); grid on;title(strtitle2);

このように高周波成分が大きい信号(画像)で圧縮率を上げると,細かなノイズが発生することになる.画像として再現すると輪郭線の周りに細かなノイズが出現することになる.これをモスキートノイズと呼んでいる.
DCT変換するデータの種類を変更して観察する.方形波,三角波は,適当なアルゴリズムで生成してみてほしい.
圧縮率を上げるにつれて元データと復元データの差が大きくなる様子を,subplotを使ってグラフを並べて表示せよ.