周波数分布の出力過程コードKl[数値]の意味

Question

前提・実現したいこと

jupyter notebookでwavファイルの離散型フーリエ変換を実施し、周波数分布図の出力を試みております。
補足情報の参照ページを参考に実施しておりますが、下記に記載の2点が分からない状態です。

意味が分かれば進められるかと思われますので、ご教授願います。

発生している問題

~~1. "def fourier (x, n, w):"　の"(x,n,w)"の意味
~~
2. "amp = [np.sqrt(c.real ** 2 + c.imag ** 2) for c in Kl[2500]]"
"amp = [np.sqrt(c.real ** 2 + c.imag ** 2) for c in Kr[2500]]"
のkl[],kr[]内の数値"2500"の意味

Kl = fourier(left, N, span)
Kr = fourier(right, N, span)
freqlist = np.fft.fftfreq(N, d=1/fr) #周波数リスト
amp = [np.sqrt(c.real ** 2 + c.imag ** 2) for c in Kl[1]] #振幅スペクトル #実部と虚部を取り出すには、".real" と ".imag" を使用
plot(freqlist, amp, marker= 'o', linestyle='-') #周波数リスト、振幅スペクトル、点、線スタイル
axis([0, fr / 2 , 0, 100000])

amp = [np.sqrt(c.real ** 2 + c.imag ** 2) for c in Kr[1]]
plot(freqlist, amp, marker= 'o', linestyle='-')

ValueError Traceback (most recent call last)
<ipython-input-30-de4a7fd136d4> in <module>
----> 1 Kl = fourier(left, N, span)
2 Kr = fourier(right, N, span)
3 freqlist = np.fft.fftfreq(N, d=1/fr) #周波数リスト
4 amp = [np.sqrt(c.real ** 2 + c.imag ** 2) for c in Kl[1]] #振幅スペクトル #実部と虚部を取り出すには、".real" と ".imag" を使用
5 plot(freqlist, amp, marker= 'o', linestyle='-') #周波数リスト、振幅スペクトル、点、線スタイル

<ipython-input-27-ad18b91aff2f> in fourier(x, n, w)
4 for i in range(0, w-2): #2番目おきに要素を得ているため、0～N-2範囲となる
5 sample = x[i * n:( i + 1) * n] #i～（i+1）番目の要素を得る
----> 6 partial = np.fft.fft(sample) #"sample"をフーリエ変換
7 K.append(partial) #"K"に"partial"を追加
8

D:\ProgramData\Anaconda3\lib\site-packages\mkl_fft_numpy_fft.py in fft(a, n, axis, norm)
158 """
159 x = _float_utils.__downcast_float128_array(a)
--> 160 output = mkl_fft.fft(x, n, axis)
161 if _unitary(norm):
162 output *= 1 / sqrt(output.shape[axis])

mkl_fft_pydfti.pyx in mkl_fft._pydfti.fft()

mkl_fft_pydfti.pyx in mkl_fft._pydfti._fft1d_impl()

mkl_fft_pydfti.pyx in mkl_fft._pydfti.__process_arguments()

ValueError: Dimension n should be a positive integer not larger than the shape of the array along the chosen axis

該当のソースコード

現在進んでいる状態です。

**第一工程**
import wave
import struct
from scipy import fromstring, int32
import numpy as np
from pylab import *
%matplotlib inline

wavfile = '440Hz.wav' #サンプルwavファイル
wr = wave.open(wavfile, "rb")
ch = wr.getnchannels() # モノラルなら1，ステレオなら2
width = wr.getsampwidth() # バイト数 (1byte=8bit)
fr = wr.getframerate() #サンプリンググレート（サンプリング周波数）
fn = wr.getnframes() # 総フレーム数（データ分割数）⇒サンプリング周波数で割れば時間

#wavfileから(N*span)のデータを先頭から切り出す「先頭から"N"個ずつ、"span"回フーリエ変換」
#連続的なデータを離散的に表すためには、データの成分は"fr"の1/2未満(ナイキスト周波数)でなければならない。
#高速フーリエ変換 ( FFT ) は、一度に変換するサンプル数が"2**n"になっていると、最も効率がいい
N = 22050 #サンプリングレート"fr"の半分の値
span =220500 #Nの整数倍　#周波数分解能：分解できる周波数の細かさ

print('チャンネル', ch)
print('バイト数', width)
print('サンプリンググレート', fr)
print('総フレーム数', fn)
print('サンプル時間', 1.0 * N * span / fr, '秒')

origin = wr.readframes(wr.getnframes()) #メソッドreadframes(n)でnフレーム分のデータを読み込む、ここでは総フレーム数分の読み込み
data = origin[:N * span * ch * width] #"origin"から要素の範囲[ ]を指定
wr.close()

print('現配列長', len(origin)) #"origin"の要素数
print('サンプル配列長: ', len(data)) #"data"の要素数

チャンネル 1
バイト数 2
サンプリンググレート 44100
総フレーム数 441000
サンプル時間 110250.0 秒
現配列長 882000
サンプル配列長: 882000

**第二工程**
X = np.frombuffer(data, dtype="int16")#"data"をバイナリ表記から16bitsの整数数列に変換

# ステレオ前提、左右音に分ける　※モノラルは単に１つおきにデータを読みこむため、必要ない工程 
left = X[::2] #"0から2番目おき"に要素を得る
right = X[1::2] #"1から２番目おき"に要素を得る

print(X)
print(left)
print(right)

[ 0 2052 4097 ... -6126 -4097 -2052]
[ 0 4097 8130 ... -12036 -8130 -4097]
[ 2052 6126 10103 ... -10103 -6126 -2052]

**第三工程**
#各サンプル区間ごとの周波数分布を配列で返してきます
def fourier (x, n, w):
    K = []
    for i in range(0, w-2): #2番目おきに要素を得ているため、0～N-2範囲となる
        sample = x[i * n:( i + 1) * n] #i～（i+1）番目の要素を得る
        partial = np.fft.fft(sample) #"sample"をフーリエ変換
        K.append(partial) #"K"に"partial"を追加

    return K

#周波数分布をもとに、実空間での波形を生成しています
def inverse_fourier (k):
    ret = []
    for sample in k:
        inv = np.fft.ifft(sample) #"sample"を逆フーリエ変換
        ret.extend(inv.real) #"inv.real"を"ret"に追加

    print (len(sample))
    return ret

試したこと

wavファイルはモノラルですが、フーリエ変換の理解のため、まずは参照ページの通りに実行してみました。

補足情報（FW/ツールのバージョンなど）

windows10 Python3.7.4
参照：https://qiita.com/niisan-tokyo/items/764acfeec77d8092eb73

Comments

[t_obara]

1の問いですが、参考ページ（のコメント）に補足があるようですが、そちらもご覧になった上で、わからないのでしょうか？

[退会済みユーザー]

ご回答ありがとうございます。すみません、こちらのページを見過ごしておりました。確認いたします。

[退会済みユーザー]

N個の複素数列x (0) ,⋯,x (N−1) に対してDFTすることで、N個の複素数列X (0) ,⋯,X (N−1)となるので、xはDFTする前のN個の数列成分、ｎは数列成分の個数、ｗは範囲を指定しているようでした。
第二工程で2番目おきに要素をえているため、"for i in range(0, w-2):"となっていることがりかいできました。

Answer 1

補足情報にあるウェブサイトの内容を確認しました。少し間違いがあるので正しながら説明したいと思います。
wavファイルを読み込む冒頭の部分については下記のようにするのが良いでしょう。

wr = wave.open(wavfile, "rb")
ch = wr.getnchannels()
width = wr.getsampwidth()
fr = wr.getframerate()
fn = wr.getnframes()
buf = wr.readframes(fn)
wr.close()

print('チャンネル', ch)
print('サンプル長(bytes)', width)
print('サンプル周波数', fr)
print('オーディオフレーム数', fn)
print('サンプル時間', fn / fr, '秒')
print('オーディオフレームのバイト数', len(buf))

全体のオーディオフレーム数はfnです。フレーム数というと別のイメージを持たれる可能性があるため、単にデータ点数と思えば良いでしょう。これをN個ずつspan回取り出します。
オリジナルのウェブサイトではNとspanの両方が指定されていますが、オーディオフレーム全体に適用する場合、spanは計算で求まります。

span = fn // N  # int(fn/N) と同じ

絵で書くと下記のような状態です。
|<---N個--->|<---N個--->|<---N個--->| ... |<---N個--->| ← Nがspan個ある

関数fourier()で行っていることは、|<---N個--->|を先頭から順々にフーリエ変換し、結果をリストに格納しています。FFTのアルゴリズム上、Nは2のn乗(オリジナルのサイトはN=2の8乗=256)である必要があります。

関数fourier()でrange(0, w-2)とされている箇所は間違っています。正しい関数は以下になります。

def fourier(x, n, w):
    K = []
    for i in range(w):
        sample = x[i * n:(i + 1) * n]
        partial = np.fft.fft(sample)
        K.append(partial)
    return K

まずbufに格納されたデータをデコードして整数配列にします。モノラルなので変数名monoとでもします。

mono = np.frombuffer(buf, dtype='int16')

先に修正した関数fourier()に、このmonoを渡してフーリエ変換を行うと、span個のフーリエ変換済みのリストが返ってきます。

K = fourier(mono, N, span)

print(len(K))とやると結果を格納した変数Kがspan個の要素を持つことがわかります。そうすると振幅を計算する部分

amp = [np.sqrt(c.real ** 2 + c.imag ** 2) for c in K]

は何をやっているかというと、python特有の内包表記というもので、分解すると次のことと等価です。

amp[0] = np.sqrt(K[0].real ** 2 + K[0].imag ** 2)
amp[1] = np.sqrt(K[1].real ** 2 + K[1].imag ** 2)
...
amp[span-1] = np.sqrt(K[span-1].real ** 2 + K[span-1].imag ** 2)

ここで注意する必要があるのはK[0]は単純な複素数ではなく、フーリエ変換済みの結果であり、要素をN個を持つ配列であることです(numpyのndarray)。numpyの様々な関数は引数に配列を渡しても演算可能で、結果を配列で返してくれます。すなわち

amp[0] = np.sqrt(K[0].real ** 2 + K[0].imag ** 2)

の一行は下記と同じことになります。

amp[0][0] = np.sqrt(K[0][0].real ** 2 + K[0][0].imag ** 2)
amp[0][1] = np.sqrt(K[0][1].real ** 2 + K[0][1].imag ** 2)
...
amp[0][N-1] = np.sqrt(K[0][N-1].real ** 2 + K[0][N-1].imag ** 2)

参考までに、ここまでをまとめた全体のソースコードを示します。

import wave
import numpy as np
import matplotlib.pyplot as plt


def fourier(x, n, w):
    K = []
    for i in range(w):
        sample = x[i * n : (i + 1) * n]
        partial = np.fft.fft(sample)
        K.append(partial)
    return K


def main():
    wavfile = "440Hz.wav"

    wr = wave.open(wavfile, "rb")
    ch = wr.getnchannels()
    width = wr.getsampwidth()
    fr = wr.getframerate()
    fn = wr.getnframes()
    buf = wr.readframes(fn)
    wr.close()

    print("チャンネル", ch)
    print("サンプル長(bytes)", width)
    print("サンプル周波数", fr)
    print("オーディオフレーム数", fn)
    print("サンプル時間", fn / fr, "秒")
    print("オーディオフレームのバイト数", len(buf))

    mono = np.frombuffer(buf, dtype="int16")
    N = 256
    span = mono.size // N  # int(fn/N) と同じ

    K = fourier(mono, N, span)
    freqlist = np.fft.fftfreq(N, d=1 / fr)
    amp = [np.sqrt(c.real ** 2 + c.imag ** 2) for c in K]

    # 先頭の結果のみ表示(場合に応じてamp[1]やamp[2]等を指定すること)
    plt.plot(freqlist, amp[0], marker="o", linestyle="-")
    plt.xlim(0, fr / 2)
    plt.ylim(0, 1e5)
    plt.show()


if __name__ == "__main__":
    main()

Comments

[退会済みユーザー]

ありがとうございます。
頂いたコードから、spanが1722ですので、その分フーリエ変換が行われていると理解できました。
spanから誤った理解をしていたのがエラーの原因だったと思われます。
また、モノラルをステレオとして左音右音にわけるといういらない作業もエラーの原因だったのでしょうか？

[yymmt]

モノラルをステレオとして分ける作業は、ソースコードを見る限りデータ数が半分になるだけで直接的なエラーの原因ではないと思われます。同じエラーを再現できていないので正確なところは分かりませんが、spanの設定が不適切なためにデータを上手く切り出せず、そのデータをFFTに掛けたのが原因ではないかと思います。

[退会済みユーザー]

ご回答ありがとうございます。調べてみます。

Answer 2

"def fourier (x, n, w):"　の"(x,n,w)"の意味

その関数では引き数を３つ取るってことですが。

Comments

[退会済みユーザー]

ご回答ありがとうございます。
x,n,wはそれぞれ何を意味しているのでしょうか。X＝ｘではないと思っているのですがいかがでしょうか。nは任意の値かと思われますが、wが分からない状態です。