いぬおさんのおもしろ数学実験室

おいしい紅茶でも飲みながら数学、物理、工学、プログラミング、そして読書を楽しみましょう

離散フーリエ係数の値を変更するときの注意点

 

 この記事に直接関連した内容の書籍をブログ管理人が出版しました。記事をリファインして、極力詳しく解説しています。見ていただけるとうれしいです(2022年5月15日(日))↓

 

 これまで、高音を抑える実験などを試してきました。

www.omoshiro-suugaku.com

そのとき、離散フーリエ係数の値を変えるのですが、守るべきルールがあります。上の記事の通りです。今回は、これを守らないとどうなるか……という話です。

 

 実データをFFTを使って周波数分析し、例えばデータから高い周波数の成分をカットしたいとしましょう。実験なので、データは8点としておきます(N=8)。FFTを実行すると次の8個の離散フーリエ係数が得られます。

f:id:Inuosann:20220321094318p:plain

これらに対し、FFTの性質から次の式が成立します。

f:id:Inuosann:20220321094043p:plain

高周波成分の係数を0にすることを考えるので、次のようにします。

f:id:Inuosann:20220321094521p:plain

これで高い周波数の成分はカットされます。しかし、離散フーリエ係数の成分の調整をこれだけにしてIFFTを実行してもとのデータがどうなるのか確認すると、まずいことになっています。具体的には元が実データだったのに、虚数が混ざってしまうのです。そもそも実データを相手にFFTを施すと

f:id:Inuosann:20220321094043p:plain

たちが成立する(共役複素数になる)はずなのですから、これらが不成立になっていると元が実データではなくなってしまうのですね。

 

 ルールを守らない例と守った例をお見せしましょう。まず、守らない例。

#===============================

import sys

import numpy as np

import scipy as sp

from scipy.fftpack import fft

y = [10, 20, 30, 40, 50, 60, 70, 80]  #実データ

yf = sp.fftpack.fft(y) #yf[k]には第k次高調波のフーリエ係数が入る

for i in range(8):

    print("■yf", i, yf[i])

yf[1] = 3-4j

yf[7] = 5  #yf[1]の共役複素数でない

yf[2] = 6j

yf[6] = 1+3j  #yf[2]の共役複素数でない

yif = sp.fftpack.ifft(yf) #ifftで元の波形に戻してみる

for i in range(8):

    print("■■yif", i, yif[i])

sys.exit()

#========================================

実行結果は次のようになります。元が実データなのに、IFFTの結果が実データに戻っていません。周波数成分を変えたのだから元のデータも変わるのは当たり前なのですが、これはまずいでしょう……。

 

■yf 0 (360+0j)

■yf 1 (-40+96.5685424949238j)

■yf 2 (-40+40j)

■yf 3 (-40+16.568542494923804j)

■yf 4 (-40+0j)

■yf 5 (-40-16.568542494923804j)

■yf 6 (-40-40j)

■yf 7 (-40-96.5685424949238j)

■■yif 0 (31.125+0.625j)

■■yif 1 (54.827795795510774-0.6553300858899105j)

■■yif 2 (44.51713562373095-1.375j)

■■yif 3 (40.02144660940672+0.3017766952966374j)

■■yif 4 (49.125+1.625j)

■■yif 5 (44.422204204489226+0.4053300858899105j)

■■yif 6 (35.23286437626905-0.875j)

■■yif 7 (60.72855339059328-0.05177669529663745j)

 

 jは虚数単位です(数学ではiを用いますが、工学の分野では電流にiを使うので虚数単位にはjを使う習慣だそうで、その関係かも知れません)。

 

 今度は離散フーリエ係数を、ルールを守って書き換える例です。

#===============================

import sys

import numpy as np

import scipy as sp

from scipy.fftpack import fft

y = [10, 20, 30, 40, 50, 60, 70, 80]

yf = sp.fftpack.fft(y) #yf[k]には第k次高調波のフーリエ係数が入る

for i in range(8):

    print("■yf", i, yf[i])

yf[1] = 3-4j

yf[7] = 3+4j  #yf[1]の共役複素数

yf[2] = 6j

yf[6] = -6j  #yf[2]の共役複素数

yif = sp.fftpack.ifft(yf) #ifftで元の波形に戻してみる

for i in range(8):

    print("■■yif", i, yif[i])

sys.exit()

#========================================

実行結果は以下です。

■yf 0 (360+0j)

■yf 1 (-40+96.5685424949238j)

■yf 2 (-40+40j)

■yf 3 (-40+16.568542494923804j)

■yf 4 (-40+0j)

■yf 5 (-40-16.568542494923804j)

■yf 6 (-40-40j)

■yf 7 (-40-96.5685424949238j)

■■yif 0 (30.75+0j)

■■yif 1 (53.879572490807405+9.18454765366783e-17j)

■■yif 2 (45.14213562373095+0j)

■■yif 3 (41.67677669529664+7.963329431634469e-16j)

■■yif 4 (49.25+0j)

■■yif 5 (43.120427509192595+9.18454765366783e-17j)

■■yif 6 (34.85786437626905+0j)

■■yif 7 (61.32322330470336-9.800238962368035e-16j)

 

離散フーリエ係数を変えたのですから元のデータも変わるのは当たり前ですが、すべて実データに戻っていることを確認してください。端数処理の関係でしょうから、

f:id:Inuosann:20220321094834p:plain

は実数とみてよいでしょう。

 

 離散フーリエ係数を変更するとき、例えばyf[3]を変更するならyf[7]もルールを守りながら変更(共役複素数にする)しなければならない、ということです。ちなみに……最初の例はどうすればよいのでしょうか。

f:id:Inuosann:20220321095418p:plain

 

 なぜか「離散フーリエ係数を変更するときは、こうしなければなりません」と教えてくれている本が見当たりません。間違えているのか……?と、ややビクビクしていますが実験結果なのだし、これが正しいのでしょう!! 「当たり前だから書いていないだけ」とか言われそうな気もしますが、初心者は困るんじゃないかな……。