離散コサイン変換について

質問

Webにあるサンプルソースや書籍のソースを見て、
感じの違うソースなのですが、
以下のソースはDCTになっているのでしょうか？
書籍には、
Data[N]=1.2.3[.3.2.1<-付け足]
して、2NでDFTすればよいとありました。
ですが、式を見てもそういう記述にはなっていません。
よろしくお願いします。
#include<stdio.h>
#define _USE_MATH_DEFINES
#include<math.h>

int main()
{
double *Data;
double *DataAfter;
double *DataConvert;
int DataLen;
DataLen=16;
Data=new double[DataLen*2];
DataAfter=new double[DataLen*2];
DataConvert=new double[DataLen*2];
for(int i=0;i<DataLen;i++)
{
Data[i]=1.5*i;
}
for(int i=0;i<DataLen;i++)
{
Data[DataLen+i]=Data[DataLen-1-i];
}
for(int i=0;i<DataLen*2;i++)
{
DataAfter[i]=0;
for(int j=0;j<DataLen*2;j++)
{
DataAfter[i]+=1.0*Data[i]*cos(M_PI*i*j/(DataLen*2));
}
DataAfter[i]/=DataLen;
}
for(int i=0;i<DataLen*2;i++)
{
DataConvert[i]=0;
for(int j=0;j<DataLen*2;j++)
{
DataConvert[i]+=1.0*DataAfter[i]*cos(-M_PI*i*j/(DataLen*2));
}
DataConvert[i]*=DataLen;
}
for(int i=0;i<DataLen*2;i++)
{
printf("%d %f\t%d %f\t%d %f\n",i,Data[i],i,DataAfter[i],i,DataConvert[i]);
}
delete[] Data;
delete[] DataAfter;
delete[] DataConvert;
getchar();
return 0;
}

回答No.1

bulgaris · ベストアンサー

元データから生成した偶対称な信号に対するDFT。
これがDCTの正体(の1つ)です(だったと思います^^;)。
この事は偶対称な信号に対するフーリエ変換式はsin項が消えるという性質が係わってます。

元データ(N個)から生成した偶対称な信号(2N個)に対してDFTを行った結果(2N個)は、偶対称な結果になります。
そのため片方半分の結果(N個)を保持すれば、もう片方の結果が得られるため、変換前<->変換後の相互の変換が可能になります。
この半分の結果は元データに対してDCTを行った結果(N個)と等価なものになります。
なぜならDCTの式は、偶対称な信号に対するフーリエ変換式はsin項が消えるという性質を基に作ることが出来るからです。

質問には答えていませんが、2NのDFTとDCTの関係について述べてみました。
自分自身は大雑把な理解しかしていないので、真実は下記の資料を参考にして下さい。

http://momonga.t.u-tokyo.ac.jp/~ooura/fftman/index.html
http://laputa.cs.shinshu-u.ac.jp/~yizawa/InfSys1/basic/index.htm
http://laputa.cs.shinshu-u.ac.jp/~yizawa/InfSys1/advanced/index.htm
http://www.amazon.co.jp/exec/obidos/ASIN/4789836797/249-4541087-5530767

お礼

補足の内容が2*i+1にすると、
偶関数（データを倍にする）の意味がわからなくなりました。
忘れてください。
ありがとうございました。

補足

回答ありがとうございます。
マルチメディア技術の基礎は途中までですが読みました。
やり直しのための信号数学も一通り読みました。

参考URLも一通り知っています。ですが、少し難解です。
後、
http://www5.airnet.ne.jp/tomy/cpro/etc19.htm
http://d.hatena.ne.jp/ryoko_komachi/20000601
このあたりでしょうか。
もう少し質問させてもらってよいですか？
若干省略しますが、
以下のようにひとつおきにデータが0となるので、
2*i+1とし、0項目の周波数をスキップしました。
IDCTも同様です。
それとは別に、
for(int j=0;j<DataLen*2;j++)のところで、
*2→*1になりそうなのですが、なりそうでしょうか？
半分の計算で済ませるためには、
重み係数？みたいなsqrt(2.0/N)をかけてやるのでしょうか？たぶん、変換式の導出と関係ありますか？
>>そのため片方半分の結果(N個)を保持すれば、
>>もう片方の結果が得られるため、変換前<->変換後の
>>相互の変換が可能になります。
と関係ありそうですが、、、。
何度もすみませんがよろしくお願いします。

//////////////////////////////////////////////////////////////////////
// DCT？
for(int i=0;i<DataLen;i++)
{
DataAfter[i]=0;
for(int j=0;j<DataLen*2;j++)
{
DataAfter[i]+=1.0*Data[i]*cos(M_PI*(2*i+1)*j/(DataLen*2));
}
DataAfter[i]/=DataLen;
}
//////////////////////////////////////////////////////////////////////
// IDCT？
for(int i=0;i<DataLen;i++)
{
DataConvert[i]=0;
for(int j=0;j<DataLen*2;j++)
{
DataConvert[i]+=1.0*DataAfter[i]*cos(-M_PI*(2*i+1)*j/(DataLen*2));
}
DataConvert[i]*=DataLen;
}
//////////////////////////////////////////////////////////////////////