ぼんさい塾

ぼんさいノートと補遺に関する素材や注釈です.ミスが多いので初稿から1週間を経た重要な修正のみ最終更新日を残しています.

離散フーリエ変換 (9)

2013-11-11 20:06:33 | 暮らし

sys.pdf
sys-s.pdf
記事一覧

差し替え: 2013-11-20


                       アナログモデルの公式

武部先生のゼミで再度 DFT でのゼロ詰めが話題になったので,Web で検索してみました.孤立波形の場合は本当のゼロなので,ゼロ詰めすれば当然周波数分解能は向上します.問題は,窓で切り出した波形を FFT を適用するために行った場合です.ゼロ詰めする前と後のスペクトルの関係を,アナログモデルで考えて見ましょう.

周期1で標本化された信号 x(t)comb(t) を長さ N の方形窓で切り取った信号を

    y(t) = x(t)comb(t)rect((t - 2-1(N-1)) / N)

とします.y(t) のフーリエ変換は

    (F1y)(ν) = (F1x)(ν) * comb(ν) * α sinc(Nν)

であり(αは定数),(F1y)(ν) = (F1y)(ν-1) なので

    (F1y)(ν) = Y(ν) * comb(ν)  ( Y(ν) = rect(ν)Y(ν) )

と表現できます(Y(ν) は (F1y)(ν) を周期1で折り返した関数).Y(ν) を標本化した Y(ν)comb(Mν) を求めるということは y(t) * comb(t / M) のスペクトルを求めることに相当します(定数倍は無視).

    x(t) = Σk ck exp(i2πkt)

のスペクトルは exp(i2πkt) のスペクトルの加重加算なので,簡単のため x(t) = exp(i2πkt) とし,N が十分に大きいとして

    Y(ν) = (δ(ν - k) * comb(ν) * α sinc(Nν))rect(ν)
          ≒ α((sinc(N(ν - k)rect(ν - k)) * comb(ν))rect(ν)

と近似した場合について考えます.sinc(N(ν - k))comb(Nν) = δ(ν - k) なので M = N であれば窓の影響はありませんが,M > N のときは Y(ν)comb(Mν) に ν = k の近傍に高調波が表れます.連続関数としての Y(ν) でなく,N 点 DFT の標本値がほしいときは,M > 4 N 程度にして Y(ν) を折れ線近似したグラフを作り,1次補間で所望の周波数での値を近似できます--- M が大きいと Y(ν) は緩やかに変化します.
※ Y(ν) の近似式は α sinc(Nν)rect(ν) を巡回シフトしたもので α sinc(Nνk)rect(ν) (νk は (ν - k/N),(ν - k/N + 1) のいずれか) と表現できます. 

補足:(0) 周期に等しくない窓をかけると基本周波数近傍のスペクトルは原信号のスペクトルとは別物になっています.このことは窓より低い周波数の正弦波が基本波成分と多数の高調波成分に分解されてしまうことからも明らかです --- 非周期信号には通常(方形でない)窓をかけます.
(1)  x(t) = exp(i2πkt) の例から,ゼロ詰めしたときのスペクトルをゼロ詰めしないときのスペクトルに変換する行列の sparsity を予想できます.
(2) sinc(Nν) との畳み込みによって輝線スペクトルが拡がる幅 Δν は ν に依存しないので周波数特性を log ν で考えるような分野では中高域の部分に窓の影響はあまりありません.
(3) 振幅変調された信号の場合,被変調波の帯域が 1/2 以下であればそのまま周期1で標本化して得られたスペクトルを搬送周波数に応じて巡回シフトすれば被変調波のスペクトルになります.
(4) rect(t) の定義を 旧sys-deleted.pdf から変えました.この結果,ほとんどいたるところでなく,厳密に (F1sinc)(ν) = rect(ν) が成立しますが,rect(t / N)comb(t) は N + 1 個の非 0 の標本値を持つことになり,周期が N でない関数の標本値を切り取るには上記のようなシフトが必要になります --- 煩わしいのですが,離散コサイン変換も視野に入れて,こちらに変更.

[1] FFT Zero Padding | BitWeenie
  http://www.bitweenie.com/listings/fft-zero-padding/
  If you apply a windowing function to your waveform, the windowing function needs to be
  applied before zero padding the data.
[2] Zero Padding Does Not Buy Spectral Resolution - National ...
  http://www.ni.com/white-paper/4880/en/
  3. Padded Records Not Always a Good Fit
[3] 振幅推定とゼロ パディング - MATLAB & Simulink - MathWorks 日本
  http://www.mathworks.co.jp/jp/help/signal/ug/amplitude-estimation-and-zero-padding.html
  DFT を求める前にデータにゼロ パディングを行えば、振幅推定を改善する上で役立つことがよくあります。
[4] Padding with Zeros
  http://www.cs.hmc.edu/~kperdue/zeroPadding.html
  However, we must make sure that we are still getting a reasonable approximation of ・・・
[5] Increasing FFT resolution using zero padding using ...
  http://forums.ni.com/t5/High-Speed-Digitizers/Increasing-FFT-resolution-using-zero-padding-using/td-p/2246228
[6] EE 4078: Zero-Padding the FFT - ECE User Home Pages
  http://users.ece.gatech.edu/~mcclella/courses/ee4078/FFT-interpol.pdf
  The zero padding must bee done "in the middle."
[7] FFTにおけるゼロ追加、補間や分解能について | 数学のQ&A【OKWave】
  http://okwave.jp/qa/q3825994.html
  http://soudan1.biglobe.ne.jp/qa3825994.html
[8] フーリエ変換について。 - Yahoo!知恵袋
  http://detail.chiebukuro.yahoo.co.jp/qa/question_detail/q1331890093
[9] フーリエ変換のデータの補間について - 数学 - 教えて!goo
  http://oshiete.goo.ne.jp/qa/5012398.html
[10] fft - Why should I zero-pad a signal before taking the Fourier ...
  http://dsp.stackexchange.com/questions/741/why-should-i-zero-pad-a-signal-before-taking-the-fourier-transform
  But, essentially, zero padding before a DFT/FFT is a computationally efficient method of ・・・
[11] Converting between zero padded and non zero padded FFT ...
  http://www.mathworks.com/matlabcentral/answers/67326