オーディオ信号のMDCTからFFTへの変換方法

このQ&Aのポイント
  • オーディオ信号の圧縮アルゴリズムにおいて、MDCT(Modified Discrete Cosine Transform)の式や計算方法について調査中です。
  • MDCTの式であるy[n] = (2/N)*Σ{x[k]*cos((2π/N)*(n+a)*(k+0.5))}をFFTで計算する方法やcos(x)の計算回数について、ご教示いただけないでしょうか。
  • また、cos(x)の値をテーブルで持つ際、何種類の値を用意する必要があるのかも知りたいです。
回答を見る
  • ベストアンサー

MDCT→FFTへの変換?(オーディオ信号)

いまオーディオ圧縮アルゴリズムの勉強をしています。 IMDCTの式が    y[n]= (2/N)*Σ{x[k]*cos((2π/N)*(n+a)*(k+0.5))}         k ただし0≦n<N, aは定数, Σはk=0,1,…N/2-1 とあるのですが、これをFFTで計算するには どのような処理をすればいいのでしょうか? また、cos(x)の計算は何回すればいいのでしょうか? #cos(x)の値をテーブルで持っておくには、 #何種類の値をもっていればいいのでしょうか? 本やネットで調べたら、DCTの式はあるのですが、 MDCTについてはあまり触れられていないので… よろしくお願いします。

  • oddo
  • お礼率86% (174/201)

質問者が選んだベストアンサー

  • ベストアンサー
  • stomachman
  • ベストアンサー率57% (1014/1775)
回答No.1

DCTをFFTでやっちゃおうという訳ですか。 FFTは F[n] = Σ(f[k] exp(-2πikn/M))  (Σはk=0,1,…M-1) を計算するんでしたよね、確か。 一方、求めたいのは Σ(x[k]cos((2π(n+a)(k+0.5)/N))) (Σはk=0,1,…N/2-1 ) で、これは Q[n]=Σ(x[k] exp(-2πi(n+a)(k+0.5)/N)) (Σはk=0,1,…N/2-1 ) の実部です。そして、 exp(-2πi(n+a)(k+0.5)/N)  =exp(-πia(2k+1)/N) exp(-2πin(2k+1)/(2N)) 以上から、k=0,1,...,N/2-1について  f[2k]=0, f[2k+1]=x[k] exp(-πia(2k+1)/N) k=N/2, ...., N-1について  f[2k]=0, f[2k+1]=0 とし、 M=2N とすれば、 F[n] = Σ(f[k] exp(-2πikn/M)) (Σはk=0,1,…M-1) = Σ(f[k] exp(-2πikn/2N)) (Σはk=0,1,…2N-1) = Σ(f[k] exp(-2πikn/2N)) (Σはk=0,1,…N-1) = Σ(f[2k] exp(-2πikn/2N)+f[2k+1]exp(-2πi(2k+1)n/2N))(Σはk=0,1,…N/2-1) = Σ(f[2k+1]exp(-2πi(2k+1)n/2N))(Σはk=0,1,…N/2-1) = Σ(x[k] exp(-πia(2k+1)/N)exp(-2πi(2k+1)n/2N))(Σはk=0,1,…N/2-1) = Σ(x[k] exp(-2πi(2k+1)(n+a)/2N))(Σはk=0,1,…N/2-1) =Q[n] ですから、その実部を取れば良いことになります。 もうちょっとましな方法もありそうですが、お急ぎのようなので。

oddo
質問者

お礼

アドバイスどうもありがとうございました。 半分あきらめかけていたので、感激しています。 いまから式をじっくり見ていきたいと思います。 教えていただいた式を見ると、cos(x)は2N個あれば 対応できそうですね。 本当にありがとうございました。

関連するQ&A

  • FFT(高速フーリエ変換)の考え方

    ///////////質問(1)//////////////////// 離散フーリエ変換が以下の式で表せることは理解できます。      N-1       nk          -2πi / N A =  Σ  a   W ,   W = e    n   k=0   k                 n:基本周波数のn倍 N:データ数 ここで、実際のプログラム上では、 (1)akに サンプリングしたデータとサンプリング周期を掛けたもの を入れ、 (2)それをWと掛けあわせて行くと思います。 では、Wには何を代入すればいいのでしょうか? Wを実部・虚部にわける、もしくは、 絶対値・位相にわけて計算する方法で あっているのでしょうか? Anは周波数n成分の面積を表していると思うのですが、面積に虚数概念がで てくるのも変な話なような気がして、質問させていただきました。 ///////////質問(2)////////////////////    N/2-1    2nk   n N/2-1     2nk A = Σ a  W + W  Σ  a  W  n  k=0  2k        k=0  2k+1 Cooley-Tukey 型 FFTは質問(1)の式を上記のように2項に分解することで、 計算を減らすことができるとのことです。 しかし、私には理解できません。 左側の項でN/2回の掛け算を行い、右側の項でもN/2回の掛け算を行うのでは、 結局式を分解しただけで、何も変わっていないように思えます。 どのように考えればいいのでしょうか? アドバイスよろしくお願いします。

  • フーリエ級数展開の問題

    フーリエ級数展開の問題 cを実数の定数とし、fは周期関数2πの関数で区間[-π,π)において f(x)= (c-2)(x+π/2) :-π<=x<0 (2c-3)(x-π/2):0<=x<π であるとする。この時のフーリエ級数展開 a_(0)/2+Σ[n=1,∞]{a_(n)cos(nx) + b_(n)sin(nx)} について各問に答えよ (1)関数fが偶関数になるような定数cの値を求め、その時のフーリエ係数a_(1)の値を求めよ。 切片が同じで、傾きが逆になればいいので、 (c-2)=-(2c-3)と式を立てて c=5/3 a_(n)=1/π∫[-π,π]f(x)cos(nx) dx -π<=x<0の時と、0<=x<πの時とを分けて積分 a_(n)=1/π{∫[-π,0](-x/3-π/6)cos(nx) dx + ∫[0,π](x/3-π/6)cos(nx) dx} n=1の時を求めればいいだけなのでn=1を代入して a_(1)=1/π{∫[-π,0](-x/3-π/6)cos(x) dx + ∫[0,π](x/3x-π/6)cos(x) dx} 式の∫[-π,0](-x/3-π/6)cos(x) dx の部分を計算 部分積分で計算し、 ∫[-π,0](-x/3-π/6)cos(x) dx=[(-x/3-π/6)sin(x)-cos(x)/3][-π,0] =-1/3-1/3==-2/3 ∫[0,π](x/3-π/6)cos(x) dx の部分を同じく計算 ∫[0,π](x/3-π/6)cos(x) dx=2/3 よって a_(1)=1/π{-2/3+2/3}=0 となってしまいました。0となり不安です間違っている気がすごくします。これで合っているんでしょうか? あと、この次の小問(2)で (2)関数fが奇関数になるような定数cの値を求め、その時のフーリエ係数b_(1)の値を求めよ。 という問題があるのですが、これはcの求め方からして分かりません。 存在しない気すらします。どのように求めればいいんでしょうか?

  • 窓関数の処理の仕方によるFFT結果の違いについて

    ある1周期分のデータX(n)={X0,X1,・・・,Xn-1}があるとき(nは2の乗数)、ハミング窓関数K=0.54-0.46cos(2πn/N-1)を以下のように処理 (1)Y(n)=K*X(n) (2)Y(n)=K*X(n)+A (3)Y(n)=K*{X(n)+A} ただし、AはX(n)の平均値。 ここで、Y(n)をFFTした場合の違いについてシミュレーションしましたが、振幅スペクトルの結果が次のようになりました。  (1)と(2)は0次のみ異なる(1次以降は同じ)  (2)と(3)((1)と(3))は0次と1次が異なる(2次以降は同じ) 【質問1】 (2)と(3)の2次以降が同じになる理由はどうしてでしょうか。

  • FFTをc言語でプログラミング

    助けてください↓c言語でFFTをプログラミングしているのですがバグが発見できずに困り果て居ます。 プログラミングは以下の通りです。画像はカラーではなく、グレースケールのMANDRILLに対して行っています。 for(n=0;n<M;n++){    for(i=0;i<N;i++){ for(j=0;j<N;j++){ Imbatx[n][i][j]=0; } } } //行FFT //まずはソート for(i=0;i<N;i++){ int x=1; for(j=0;j<N;j++){ if(j<N/2){        //左半分 if(j%2!=0) //jが奇数 batx[0][i][j]=imagein[i][j+N/2-1]; if(j%2==0) //jが偶数 batx[0][i][j]=imagein[i][j]; } if(N/2<=j && j<N){ //右半分 if(j%2!=0) //jが奇数 batx[0][i][j]=imagein[i][x+N/2-1];      if(j%2==0) //jが偶数 batx[0][i][j]=imagein[i][x]; x++; } } } for(i=0;i<N;i++){ for(n=1;n<M+1;n++){ k=0; for(j=0;j<N;j++){ mul=(int)pow(2.0,(double)n); if(j%mul < mul/2){ //偶数列 batx[n][i][j] = batx[n-1][i][j]+cos(2*PI*k/mul)*batx[n-1][i][j+mul/2]- sin(2*PI*k/mul)*Imbatx[n-1][i][j+mul/2]; Imbatx[n][i][j]=Imbatx[n-1][i][j]+sin(2*PI*k/mul)*batx[n-1][i][j+mul/2]+ cos(2*PI*k/mul)*Imbatx[n-1][i][j+mul/2]; k++; } else{ //奇数列 batx[n][i][j] = cos(2*PI*k/mul)*batx[n-1][i][j] + batx[n-1][i][j-mul/2]- sin(2*PI*k/mul)*Imbatx[n-1][i][j]; Imbatx[n][i][j]=Imbatx[n-1][i][j-mul/2]+sin(2*PI*k/mul)*batx[n-1][i][j]+cos(2*PI*k/mul)*Imbatx[n-1][i][j]; k++; } Re_countx[i][j] = batx[M][i][j]; Im_countx[i][j] = Imbatx[M][i][j]; } } } for(i=0;i<N;i++){ for(j=0;j<N;j++){ FFT[i][j]=pow(pow(Re_countx[i][j],2)+pow(Im_countx[i][j],2),0.5); FFT[i][j]=log10(FFT[i][j]); if(FFT[i][j]>max) max=FFT[i][j]; imageout[i][j]=(unsigned char)255*FFT[i][j]/max; } } 結果的に出力される画像は以下の画像になります。 ごらんのように縞模様が出てしまいます。 分かる方どうかお願いいたします。

  • 離散フーリエ変換について

    double DiscreteFourierTransformClass::DiscreteFourierTransformProcess(int n) { double e=2.7182818; double Pai=3.1415926535; int Tau=10; int N=DataLen/Tau; double Sum=0; for(int k=0;k<N;k++) { Sum+=Tau*Data[k*Tau]*pow(e,i*2*Pai*k*n/N); } return Sum; } Tau*Data[k*Tau]=帯 pow(e,i*2*Pai*k*n/N)=ほしい周波数の乗算 ということで、離散フーリエ展開と同じなのでしょうか? パソコン上で計算するには、離散フーリエ展開で計算して、 FFTというのは、pow(e,i*2*Pai*k*n/N)この回転子が3,4個程度のcos,sinで決まるから その値で計算すれば高速になるよということでしょうか? 離散フーリエ変換 G(n/N)=Σ[k=0::N-1]{τ*f(kτ)e^-i2πkn/N} k:何番目の値か τ:値を読んだ間隔 N:値を読んだ数 n:基本周波数の何倍か G(n/N)->ほしいnの成分を得る

  • フーリェ変換(特に一次元のFFT)

    フーリェ変換について勉強しているものです。このページ(http://www.kurims.kyoto-u.ac.jp/~ooura/fftman/ftmndl.html)で配布している、fft.tgz (71 KB)というファイルのdfst()という関数を使おうと思っているんですが、その関数の返り値(返り配列と言った正しい?)の詳しい意味がよくわからないのです。 --------------------------------------------------------- 以下の変換は : S[k] = sum_j=1^n-1 a[j]*sin(pi*j*k/n), 0<k<n このように分割される : S[2*k] = sum_j=1^n/2-1 (a[j]-a[n-j])*sin(pi*j*k/(n/2)), ◎ S[2*k+1] = sum'_j=1^n/2 (a[j]+a[n-j])*sin(pi*j*(k+1/2)/(n/2)) 注意:sum_■^▲はΣを表しており、■から▲までの部分和となっています。 nは要素数です。 ------------------------------------------------------ この◎の関数はどのような意味を持つのでしょうか?

  • 三角関数の問題です。

    三角関数の問題です。 2次方程式 5x^2-7x+k=0 の2つの解が、sinΘ、cosΘであるとき、  定数k の値と sin^3Θ+cos^3Θの値を求めよ。 です。 「sinΘ+cosΘ=7/5」 「sinΘcosΘ=k/5」 を使って計算するらしいのですが、 この2つの式はどうやって求めたのでしょうか?

  • FFTの基底の選び方について教えてください。

    FFTの基底の選び方について教えてください。 NETLIB-NAPACKから任意個数のFFTをダウンロードし、個数が素数の計算をしたところ計算時間がものすごくかかります。例えば、N=16411で複素数の場合、そのままだと約15秒、後にゼロをつけて2のべき乗(N=32768)にすると0.02秒、N=17496(=2^3*3^7、試行で求めた)にすると0.03秒となります。そこで後につけるゼロの個数をできるだけ減らし、かつ元個数に最も近いような個数(例えば今回の17496)を、2,3,5,7のべき乗で表す個数(今回の17496)を求めるアルゴリズム、もしくはプログラムは無いものでしょうか。

  • cos(x/2)*cos(x/2^2)*・・・・・cos(x/2^n)

    実数x及び自然数nに対して a_n=cos(x/2)*cos(x/2^2)*・・・・・cos(x/2^n) とする。 (1)2^n*a_n*sin(x/2^n)の値はnと無関係に一定であることを証明せよ。 (2)log|a_n|をxで微分することにより、 Σ(n=2~∞)1/2^n *tan(π/2^n)=1/π であることを証明せよ この問題に取り組んでいます。 (1)で2^n*a_n*sin(x/2^n)の計算を行っていて、いろいろな三角関数の公式を利用してみたのですが全然うまくいきません。「nと無関係」ということはnが消えればいいということだと思うのですが・・・。 (2)はloga_nを微分したところ -1/2 tan(x/2) - 1/2^2 tan(x/2^2) -・・・となったのですがここから証明すべき式に変形するにはどうしたらいいのでしょうか? 回答いただければありがたいです。よろしくお願いします

  • FFTについて。

    今、FFTについて勉強しています。 参考に以下のURLから参考プログラムをよんで、試しに 自分のサンプル値を入力して出力を試みました。 http://www.cn.kagawa-nct.ac.jp/~kusama/study/fdtd/fdtd_3d/fft.f90 しかし、このURLでの出力(output)のn, n/(s*t), cabs(x(n))で周波数の値がでてきていないことに 疑問を持ちました。 もしかしたら、周波数の値の意味を持っているのかもしれませんが、 まだ勉強不足で理解ができていません。 お手数ですが、どなたかご指摘くださる方がいらっしゃったらよろしくお願いします。