前回の記事では、Pythonを利用して、作成した正弦波の音声をFFT処理し、周波数のグラフ化までを行いました。
今回の記事では、周波数データに変換したピー音の周波数を変更し、元の波形に戻すとどうなるか実験してみます。
その前に、前回の記事で出力したFFTのデータをもう一度良く見て見るために、少しプログラムを変えてグラフ出力します。具体的には20行目に下記変更を加えて、データ全体をグラフに表します。
【変更前】plt.plot(freq[1:int(N/2)], Amp[1:int(N/2)])
【変更後】plt.plot(F.imag[1:N])
実行して得られた周波数領域のグラフが以下になります。
続いて、7行目に下記変更を加えて周波数を変更します。
【変更前】f0 = 440
【変更後】f0 = 880
実行して得られた周波数領域のグラフは以下になります。
図1 440Hzの周波数グラフと図2 880Hzの周波数グラフの違いを確認してみましょう。
図2のほうが縦棒が中央に寄っていることがわかると思います。
これは周波数が上がれば、周波数領域のデータが中央に寄ることを意味しています。
この結果を元に、440Hzの音声データをFFTし、周波数データを中央に寄せて、逆FFTで元の音声データに戻してみます。
プログラム全体は下記となります。
| import wave | |
| import numpy as np | |
| import matplotlib.pyplot as plt | |
| import struct | |
| # 周波数シフト関数 | |
| def shift_freq(F, s_freq_hz, fs, sec=1): | |
| N = fs * sec | |
| s_freq = s_freq_hz * sec | |
| # プラス方向へのシフト | |
| if s_freq > 0: | |
| # 前半部分のシフト | |
| for i in reversed(range(0,int(N/2))): | |
| si = i - s_freq | |
| if si >= 0: | |
| F.real[i] = F.real[si] | |
| F.imag[i] = F.imag[si] | |
| elif si < 0: | |
| F.imag[i] = 0 | |
| # 後半部分のシフト | |
| for i in range(int(N/2)+1,N): | |
| si = i + s_freq | |
| if si < N: | |
| F.real[i] = F.real[si] | |
| F.imag[i] = F.imag[si] | |
| elif si >= N: | |
| F.imag[i] = 0 | |
| # マイナス方向へのシフト | |
| elif s_freq < 0: | |
| # 前半部分のシフト | |
| for i in range(0,int(N/2)): | |
| si = i - s_freq | |
| if si < int(N/2): | |
| F.real[i] = F.real[si] | |
| F.imag[i] = F.imag[si] | |
| elif si >= int(N/2): | |
| F.imag[i] = 0 | |
| # 後半部分のシフト | |
| for i in reversed(range(int(N/2)+1,N)): | |
| si = i + s_freq | |
| if si > int(N/2)+1: | |
| F.real[i] = F.real[si] | |
| F.imag[i] = F.imag[si] | |
| elif si < int(N/2)+1: | |
| F.imag[i] = 0 | |
| return F | |
| if __name__ == '__main__': | |
| a = 1 #振幅 | |
| fs = 8192 #サンプリング周波数 | |
| f0 = 440 #周波数 | |
| sec = 5 #秒 | |
| N = fs * sec | |
| swav=[] | |
| for n in np.arange(fs * sec): | |
| #サイン波を生成 | |
| s = a * np.sin(2.0 * np.pi * f0 * n / fs) | |
| swav.append(s) | |
| # 正弦波プロット | |
| plt.plot(swav[0:100]) | |
| plt.savefig("sin_440.png") | |
| plt.gca().clear() | |
| # FFT | |
| F = np.fft.fft(swav) # 変換結果 | |
| # Fを弄って周波数を変更 | |
| F = shift_freq(F, 440, fs, sec=5) | |
| # 逆FFT | |
| swav_ifft = np.fft.ifft(F) | |
| plt.plot(swav_ifft.real[0:100]) | |
| plt.savefig("sin_ifft.png") | |
| plt.gca().clear() | |
| #サイン波を-32768から32767の整数値に変換(signed 16bit pcmへ) | |
| swav = [int(x * 32767.0) for x in swav_ifft.real] | |
| #バイナリ化 | |
| binwave = struct.pack("h" * len(swav), *swav) | |
| #サイン波をwavファイルとして書き出し | |
| w = wave.Wave_write("output_ifft.wav") | |
| p = (1, 2, 8000, len(binwave), 'NONE', 'not compressed') | |
| w.setparams(p) | |
| w.writeframes(binwave) | |
| w.close() |
少し長くなりましたので、周波数変更処理の部分は関数として定義(7行目から46行目)しています。
48行目から86行目がメイン処理になります。
メイン処理はやっていることは単純なので、流れのみ説明します。
１．440Hzの正弦波を作成し、グラフ化しています。
２．440HzのデータをFFTし、周波数データに変換しています。
３．周波数データを新たに定義したshift_freqという関数を使って440Hzプラスします。
４．変更した周波数データに対して、逆FFT処理をかけて元の波形に戻し、グラフ化します。
５．元に戻したデータをWaveファイルとして出力しています。
続いて、shift_freqの処理について解説していきます。
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
def shift_freq(F, s_freq_hz, fs, sec=1):
N = fs * sec
s_freq = s_freq_hz * sec
7行目で、「shift_freq」という関数を定義しています。引数として、FFTの実行結果、シフトする周波数、サンプリング周波数、音声の長さ(秒)を定義しています。シフトする周波数はプラスでもマイナスでも構いません。今回のメイン処理ではわかりやすいように440Hzプラスし、880Hzにしています。
8行目で、サンプリング周波数に音声の長さ(秒)を乗算して、信号データの総数Nを算出しています。
9行目では、シフトする周波数に音声の長さ(秒)を乗算して、実際にシフトする量を算出しています。
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# プラス方向へのシフト
if s_freq > 0:
# 前半部分のシフト
for i in reversed(range(0,int(N/2))):
si = i - s_freq
if si >= 0:
F.real[i] = F.real[si]
F.imag[i] = F.imag[si]
elif si < 0:
F.imag[i] = 0
# 後半部分のシフト
for i in range(int(N/2)+1,N):
si = i + s_freq
if si < N:
F.real[i] = F.real[si]
F.imag[i] = F.imag[si]
elif si >= N:
F.imag[i] = 0
10行目から27行目で、周波数をプラス方向にシフトする場合の処理を定義しています。
13行目から19行目のループでは、データの中央部から左方向に向かって値を見ていき、シフト分先(si)にある値を現在(i)の位置にシフトします。シフト分先(si)がマイナスになったら、現在(i)の位置に0を入れます。
21行目から27行目のループでは、データの中央部から右方向に向かって値を見ていき、シフト分先(si)にある値を現在(i)の位置にシフトします。こちらもシフト分先(si)が総数Nを超えたら、現在(i)の位置に0を入れます。
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# マイナス方向へのシフト
elif s_freq < 0:
# 前半部分のシフト
for i in range(0,int(N/2)):
si = i - s_freq
if si < int(N/2):
F.real[i] = F.real[si]
F.imag[i] = F.imag[si]
elif si >= int(N/2):
F.imag[i] = 0
# 後半部分のシフト
for i in reversed(range(int(N/2)+1,N)):
si = i + s_freq
if si > int(N/2)+1:
F.real[i] = F.real[si]
F.imag[i] = F.imag[si]
elif si < int(N/2)+1:
F.imag[i] = 0
28行目から45行目の処理は、周波数をマイナス方向にシフトする場合の処理を定義しています。
処理内容の説明はプラス方向と同様のため割愛します。
最後の46行目でシフトしたFを戻り値としてメイン処理に渡して関数処理は完了します。
本プログラムを実行して得られた波形画像は以下になります。
図3が作成した440Hzの波形で、図4が、440Hzの波形をFFTして440Hzシフトし880Hzとして、波形データに戻したものです。周波数を上げたので波形が細かくなっているのがわかります。
音声でも確認してみましょう。音量にご注意ください。
プー音からピー音になりましたね。
周波数領域のデータをシフトすることで、周波数を変更できることが確認できました。
次の記事では、作成した波形ではなく、実際の音声を利用して処理してみたいと思います。
コメント