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

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

Python基本サウンドプログラミング～周波数分析，高速フーリエ変換(FFT)入門～　専用ページ

Python基本サウンドプログラミング～周波数分析，高速フーリエ変換(FFT)入門～

Python基本サウンドプログラミング～周波数分析，高速フーリエ変換(FFT)入門～

Amazon

コードをコピーし、読者の開発環境で貼り付けてご使用ください。

目次

2－4　2音を鳴らすwavファイルを作る

#==================
#waveファイルを作る。左は440Hz、右は700Hzで10秒間。
#
import sys
import wave
import numpy as np
import datetime 
dt = datetime.datetime.now()
print('▼▼▼▼ {}:{}:{}'.format(dt.hour, dt.minute, dt.second))
pi = np.pi
fs = 44100      #サンプリング周波数
s  = 10         #10秒の長さのwaveファイルを作る
#s秒分のサウンドデータを作るには、fs*sフレームが必要。
N  = fs * s     #フレーム数（左右で1フレーム）
t = np.arange(0, s, 1/fs)   #横軸
#sin(2π*□＊t) の周波数は□の部分。16ビットで録音するので、
#－32768～32767がデータの範囲。だから－15000～15000の範囲の
#大きさの音を作ることにした。
ldata = 15000 * np.sin(2*pi*440*t)  #440Hz
rdata = 15000 * np.sin(2*pi*700*t)  #700Hz
sounddata = np.zeros(N * 2, dtype= "int16") #左右のデータ分が必要
sounddata[0::2] = ldata
sounddata[1::2] = rdata
writewave = wave.Wave_write("keshitene.wav") #保存先ファイル
#setparamの引数は
#(nchannels, sampwidth, framerate, nframes, comptype, compname)
#先頭から順にステレオ、サンプルのサイズ2バイト、
#サンプリング周波数fs、フレーム数N。
writewave.setparams((2, 2, fs, N, 'NONE', 'NONE'))
writewave.writeframes(sounddata) #編集したサウンドデータを保存
writewave.close()
dt = datetime.datetime.now()
print('▲▲▲▲ {}:{}:{}'.format(dt.hour, dt.minute, dt.second))
sys.exit()
#==================

2－6　平均律、純正律のドミソの和音を聞き比べる

#==================
#wavファイルを作る。ドミソの和音。
#
import sys
import numpy as np
import wave
pi = np.pi
fs = 22050  #サンプリング周波数
s  = 20  #20秒の長さのwaveファイルを作る
#s秒分のサウンドデータを作るには、fs*sフレームが必要。
N  = fs * s  #フレーム数
t = np.arange(0, s, 1/fs)  #横軸
#sin(2π*□＊t) の周波数は□の部分。16ビットで録音するので、
#－32768～32767がデータの範囲であることに注意。
#この範囲の外に出ると音が潰れる。
do = 523.3 #ドの音の周波数
r12 = 2 ** (1/12) #2の12乗根（≒1.059463094）
#--------------
# A
#
#純正律でドミソ
#ldata = 10000 * np.sin(2*pi*do*t) + 10000 * np.sin(2*pi*do*5/4*t)
#rdata = 10000 * np.sin(2*pi*do*3/2*t)
#--------------
# B
#
#平均律でドミソ
ldata = 10000 * np.sin(2*pi*do*t) + 10000 * np.sin(2*pi*do*(r12**4)*t)
rdata = 10000 * np.sin(2*pi*do*(r12**7)*t)
#--------------
sounddata = np.zeros(N * 2, dtype= "int16") #左右のデータ分が必要
sounddata[0::2] = ldata
sounddata[1::2] = rdata
writewave = wave.Wave_write("keshitene.wav") #保存先ファイル
#setparamの引数は
#(nchannels, sampwidth, framerate, nframes, comptype, compname)
#先頭から順にステレオ、サンプルのサイズ2バイト、
#サンプリング周波数fs、フレーム数N。
writewave.setparams((2, 2, fs, N, 'NONE', 'NONE'))
writewave.writeframes(sounddata) #編集したサウンドデータを保存
writewave.close()
sys.exit()
#==================

2－7　周波数の近い2音を鳴らしてうなりを確認する

#==================
import sys
import numpy as np
import wave
pi = np.pi
fs = 22050      #サンプリング周波数
s  = 20         #20秒の長さのwavファイルを作る
#s秒分のサウンドデータを作るには、fs*sフレームが必要。
N  = fs * s     #フレーム数
t = np.arange(0, s, 1/fs)   #横軸
#sin(2π*□＊t) の周波数は□の部分。16ビットで録音するので、
#－32768～32767がデータの範囲であることに注意。
#この範囲の外に出ると音が潰れる。
do = 523.3 #ドの音の周波数
r12 = 2 ** (1/12) #2の12乗根
#--------------
ldata = 15000 * np.sin(2*pi*do*5/4*t) #左スピーカー。純正律でミ
rdata = 15000 * np.sin(2*pi*do*(r12**4)*t) #右スピーカー。平均律でミ
#--------------
sounddata = np.zeros(N * 2, dtype= "int16") #左右のデータ分が必要
sounddata[0::2] = ldata
sounddata[1::2] = rdata
writewave = wave.Wave_write("keshitene.wav") #保存先ファイル
#setparamの引数は
#(nchannels, sampwidth, framerate, nframes, comptype, compname)
#先頭から順にステレオ、サンプルのサイズ2バイト、
#サンプリング周波数fs、フレーム数N。
writewave.setparams((2, 2, fs, N, 'NONE', 'NONE'))
writewave.writeframes(sounddata) #編集したサウンドデータを保存
writewave.close()
sys.exit()
#==================

2－8　プッシュボタンを使わずに電話をかけられるか

#==================
#wavファイルを作る。1をプッシュ
#
import sys
import numpy as np
import wave
pi = np.pi
fs = 22050  #サンプリング周波数
s  = 20  #20秒の長さのwavファイルを作る
#s秒分のサウンドデータを作るには、fs*sフレームが必要。
N  = fs * s  #フレーム数
t = np.arange(0, s, 1/fs)  #横軸
#sin(2π*□＊t) の周波数は□の部分。16ビットで録音するので、
#－32768～32767がデータの範囲であることに注意。
#この範囲の外に出ると音が潰れる。
#----------------------
#1……697Ｈｚ、1209Ｈｚ
#7……852Ｈｚ、1209Ｈｚ
#
#697の部分を852に書き換えると「7」になる
ldata = 15000 * np.sin(2*pi*697*t)
rdata = 15000 * np.sin(2*pi*1209*t)  #1209は1，7で共通
#----------------------
sounddata = np.zeros(N * 2, dtype= "int16")  #左右のデータ分が必要
sounddata[0::2] = ldata
sounddata[1::2] = rdata
writewave = wave.Wave_write("keshitene.wav")  #保存先ファイル
#setparamの引数は
#(nchannels, sampwidth, framerate, nframes, comptype, compname)
#先頭から順にステレオ、サンプルのサイズ2バイト、
#サンプリング周波数fs、フレーム数N。
writewave.setparams((2, 2, fs, N, 'NONE', 'NONE'))
writewave.writeframes(sounddata)  #編集したサウンドデータを保存
writewave.close()
sys.exit()
#==================

2－10　wavファイルのパラメータを読み取る

#==================
import sys
import datetime
import wave
dt = datetime.datetime.now()
print('▼▼▼▼ {}:{}:{}'.format(dt.hour, dt.minute, dt.second))
wavefile= wave.open("C:\\Users\\xxxx\\Desktop\\sound\\'yyyy.wav","rb")
print('■parameter = ', wavefile.getparams())
print('■framerate = ', wavefile.getframerate())
wavefile.close()
dt = datetime.datetime.now()
print('▲▲▲▲ {}:{}:{}'.format(dt.hour, dt.minute, dt.second))
sys.exit()
#==================

2－11　matplotlibモジュールでグラフを描く

1本目

#==================
#グラフを描画する
import numpy as np
import matplotlib.pyplot as plt
t = np.arange(0, 1, 0.2)
y = t * t
plt.plot(y)
plt.show()
#==================

2本目

#==================
import numpy as np
import matplotlib.pyplot as plt
t = np.arange(0, 1, 0.2)
y = t * t
plt.plot(t, y)
plt.show()
#==================

3本目

#==================
import numpy as np
import matplotlib.pyplot as plt
pi = np.pi #円周率、π
N = 100 #1秒でデータ100個
t = np.arange(0, 1, 1/N) #時刻0以上1秒未満、1/N秒刻みの時刻の配列
#A
y = np.sin(2*pi*4*t) #周波数は4Hz
#B
#y = np.sin(2*pi*4*t)+2*np.sin(2*pi*10*t) #周波数は4Hz，10Hz
plt.plot(t, y) #波形を描画
plt.show()
#==================

3－1　wavファイルで、サウンドデータの符号を一斉に反転する

#==================
import sys
import numpy as np
import wave
wavefile = wave.open("BGM.wav", "rb")
buf = wavefile.readframes(wavefile.getnframes())#サウンドデータ部分の読み取り
buf = np.frombuffer(buf, dtype= "int16")#16bit符号付き整数に変換
params = wavefile.getparams()
wavefile.close()
print(params)
if params.sampwidth != 2:
    sys.exit() #16ビットで記録していなければ終了
sounddata = buf.copy() #サウンドデータの操作用バッファ
print(sounddata)
sounddata = -sounddata  #サウンドデータの符号を一斉に反転
print(sounddata) #サウンドデータを見てみる
writewave = wave.Wave_write("keshitene.wav") #保存先ファイル
writewave.setparams(wavefile.getparams()) #保存先にパラメータをセット
writewave.writeframes(sounddata) #編集したサウンドデータを保存
writewave.close()
sys.exit()
#==================

3－2－1　クリッピングされたデータでグラフを描く

#==================
import numpy as np
import matplotlib.pyplot as plt
pi = np.pi #円周率、π
N = 1000 #1秒でデータ1000個
t = np.arange(0, 1, 1/N) #0以上1秒未満、1/N秒刻みの時刻の配列
y = np.sin(2*pi*2*t) #周波数は2Hz
y[y < -0.8] = -0.8
y[y > 0.8] = 0.8
plt.ylim(-1, 1)
plt.plot(t, y) #波形を描画
plt.show()
#==================

3ー2－2　wavファイルでクリッピングする

#==================
#クリッピング
#
import sys
import wave
import numpy as np
import datetime
#-------------------------------------------------
dt = datetime.datetime.now()
print('▼▼▼▼ {}:{}:{}'.format(dt.hour, dt.minute, dt.second))
wavefile = wave.open("C:\\Users\\xxxx\\Desktop\\BGM.wav","rb")
params = wavefile.getparams()
if params.sampwidth != 2:
    wavefile.close()
    sys.exit() #16ビットで記録していなければ終了
nframes = params.nframes #フレーム数の取得
buf = wavefile.readframes(nframes) #サウンドデータの読み取り
buf = np.frombuffer(buf, dtype= "int16") #16bit符号付き整数に変換
wavefile.close()
bufcopy = buf.copy()
lim = 1000
bufcopy[bufcopy > lim] = lim
bufcopy[bufcopy < -lim] = -lim
writewave = wave.Wave_write("C:\\Users\\xxxx\\Desktop\\keshitene.wav")
writewave.setparams(params) #元のwavファイルと同じパラメータをセット
writewave.writeframes(bufcopy) #編集したサウンドデータを保存
writewave.close()
dt = datetime.datetime.now()
print('▲▲▲▲ {}:{}:{}'.format(dt.hour, dt.minute, dt.second))
sys.exit()
#==================

3－3　逆再生する、ドレミファソラシドを鳴らす

逆再生する

#==================
#逆再生
#
import sys
import wave
import numpy as np
import datetime
dt = datetime.datetime.now()
print('▼▼▼▼ {}:{}:{}'.format(dt.hour, dt.minute, dt.second))
wavefile = wave.open("C:\\Users\\xxxx\\Desktop\\keshitene.wav","rb")
params = wavefile.getparams()
nframes = params.nframes #フレーム数の取得
if (params.sampwidth != 2):
    wavefile.close()
    sys.exit() #16ビットで記録していなければ終了
buf = wavefile.readframes(nframes) #サウンドデータの読み取り
buf = np.frombuffer(buf, dtype= "int16") #16bit符号付き整数に変換
wavefile.close()
buf1 = np.zeros(buf.size, dtype = "int16")
for i in range(buf.size): #bufの内容を逆転してbuf1へ格納
    buf1[i] = buf[buf.size - i - 1]
writewave= wave.Wave_write("C:\\Users\\xxxx\\Desktop\\keshitene.wav")
writewave.setparams(params)
writewave.writeframes(buf1) #編集したサウンドデータを保存
writewave.close()
#
dt = datetime.datetime.now()
print('▲▲▲▲ {}:{}:{}'.format(dt.hour, dt.minute, dt.second))
sys.exit()
#==================

ドレミファソラシドを鳴らす

#==================
#wavファイルを作る。ドレミファソラシド
#
import sys
import numpy as np
import wave
import datetime
dt = datetime.datetime.now()
print('▼▼▼▼ {}:{}:{}'.format(dt.hour, dt.minute, dt.second))
pi = np.pi
fs = 22050 #サンプリング周波数
s  = 0.5
#s秒分のサウンドデータを作るには、fs*sフレームが必要。
N  = int(fs * s) #フレーム数
t = np.arange(0, s, 1/fs) #横軸
do = 523.3 #ドの音の周波数
x = do
r12 = 2 ** (1/12) #2の12乗根（≒1.059463094）
#--------------
#平均律でドレミファソラシド
ldata = 15000 * np.sin(2*pi*x*t) #ド
x = x * r12 * r12
ldata = np.concatenate([ldata, 15000 * np.sin(2*pi*x*t)]) #レ
x = x * r12 * r12
ldata = np.concatenate([ldata, 15000 * np.sin(2*pi*x*t)]) #ミ
x = x * r12
ldata = np.concatenate([ldata, 15000 * np.sin(2*pi*x*t)]) #ファ
x = x * r12 * r12
ldata = np.concatenate([ldata, 15000 * np.sin(2*pi*x*t)]) #ソ
x = x * r12 * r12
ldata = np.concatenate([ldata, 15000 * np.sin(2*pi*x*t)]) #ラ
x = x * r12 * r12
ldata = np.concatenate([ldata, 15000 * np.sin(2*pi*x*t)]) #シ
x = x * r12
ldata = np.concatenate([ldata, 15000 * np.sin(2*pi*x*t)]) #ド
rdata = ldata
sounddata = np.zeros(int(N*2*8), dtype= "int16")
sounddata[0::2] = ldata
sounddata[1::2] = rdata
writewave= wave.Wave_write("C:\\Users\\xxxx\\Desktop\\keshitene.wav") #保存先ファイル
#setparamの引数は
#(nchannels, sampwidth, framerate, nframes, comptype, compname)
#先頭から順にステレオ、サンプルのサイズ2バイト、
#サンプリング周波数fs、フレーム数N。
writewave.setparams((2, 2, fs, N, 'NONE', 'NONE'))
writewave.writeframes(sounddata) #編集したサウンドデータを保存
writewave.close()
dt = datetime.datetime.now()
print('▲▲▲▲ {}:{}:{}'.format(dt.hour, dt.minute, dt.second))
sys.exit()
#==================

3－4　変声機の原理を確認する

#==================
#左のスピーカーからA-Bが流れ、右から Bが流れるwavファイル
import sys
import numpy as np
import scipy as sp
from scipy.fftpack import fft
import wave
import datetime
dt = datetime.datetime.now()
print('▼▼▼▼ {}:{}:{}'.format(dt.hour, dt.minute, dt.second))

#音楽A
wavefile = wave.open("C:\\Users\\xxxx\\Desktop\\音楽A.wav","rb") #こちらの方が長い曲
params = wavefile.getparams()
if (params.sampwidth != 2) or (params.nchannels != 2):
    wavefile.close()
    sys.exit() #16ビットで記録していないか、ステレオでなければ終了
#サウンドデータ部分の読み取り
buf = wavefile.readframes(wavefile.getnframes())
buf = np.frombuffer(buf, dtype= "int16")#16bit符号付き整数に変換
print(wavefile.getparams())
wavefile.close()
sdata = buf.copy()
ldata = sdata[0::2] #左のデータ
rdata = sdata[1::2] #右のデータ

#音楽B
#こちらの方が短い曲
wavefile1 = wave.open("C:\\Users\\xxxx\\Desktop\\音楽B.wav","rb")
params = wavefile1.getparams()
if (params.sampwidth != 2) or (params.nchannels != 2):
    wavefile1.close()
    sys.exit() #16ビットで記録していないか、ステレオでなければ終了
#サウンドデータ部分の読み取り
buf = wavefile1.readframes(wavefile1.getnframes())
buf = np.frombuffer(buf, dtype= "int16")#16bit符号付き整数に変換
print(wavefile1.getparams())
wavefile1.close()
sdata1 = buf.copy()
ldata1 = sdata1[0::2] #左のデータ
rdata1 = sdata1[1::2] #右のデータ
soundLen = rdata1.size #短い方の曲の長さ

#2つの曲(A,B)ともステレオの右のデータだけを使う
#rdataにはA-B、ldataにはBを入れる
rdata[:soundLen] = rdata[:soundLen]//2 - rdata1//2  #A-B を入れる
ldata[:soundLen] = rdata1//2  #B を入れる
sdata[0::2] = ldata
sdata[1::2] = rdata

writewave = wave.Wave_write("C:\\Users\\xxxx\\Desktop\\keshitene.wav")
writewave.setparams(wavefile.getparams()) #保存先にパラメータをセット
writewave.writeframes(sdata) #編集したサウンドデータを保存
writewave.close()
dt = datetime.datetime.now()
print('▲▲▲▲ {}:{}:{}'.format(dt.hour, dt.minute, dt.second))
sys.exit()
#==================

4－5　SciPyのFFT、IFFTを使ったわかりやすい例（1）

#==================
import numpy as np
import scipy as sp
from scipy.fftpack import fft
import matplotlib.pyplot as plt

pi = np.pi
fs = 600 #サンプリング周波数
N = 200 #サンプル数

t = np.arange(0, N/fs, 1/fs) #1/fs秒刻みの時刻の配列
y = np.sin(2*pi*30*t)+2*np.sin(2*pi*45*t) #周波数は30Hz，45Hz
plt.ylim(-4, 4) #描画するｙ軸の範囲をセット
plt.plot(t, y) #元の波形を描画
plt.show()

#yf[k]には第ｋ次高調波のフーリエ係数が入る
yf = sp.fftpack.fft(y[:N])
plt.plot(np.abs(yf)) #含まれる周波数成分を描画
plt.show()
yf[10] = 0 #第10次高調波をカット
yf[N-10] = 0 #第10次高調波をカット
yif = sp.fftpack.ifft(yf[:N]).real
plt.ylim(-4, 4)
plt.plot(t, yif)
plt.show()
#==================

4－6　SciPyのFFT、IFFTを使ったわかりやすい例（2）

#==================
import numpy as np
import scipy as sp
from scipy.fftpack import fft
import matplotlib.pyplot as plt
#
fs = 600
t = np.arange(0, 6, 1/fs)
y= np.sin(2*np.pi*50*t)+2*np.sin(2*np.pi*80*t)+3*np.sin(2*np.pi*100*t)
yf = sp.fftpack.fft(y[:300])
plt.plot(np.abs(yf))
plt.show()
#==================

4－7　SciPyのFFT、IFFTを使ったわかりやすい例（3）

#==================
import numpy as np
import scipy as sp
from scipy.fftpack import fft
import matplotlib.pyplot as plt

pi = np.pi
fs = 100 #サンプリング周波数
N = fs #サンプル数はfsとする

t = np.arange(0, 1, 1/fs) #1/fs秒刻みの時刻の配列。基本周期T = 1。
y = np.sin(2*pi*4*t)+2*np.sin(2*pi*10*t) #周波数は4Hz，10Hz
plt.plot(t, y) #元の波形を描画
plt.show()

#yf[k]には第ｋ次高調波のフーリエ係数が入る
yf = sp.fftpack.fft(y[:N])
plt.plot(np.abs(yf)) #含まれる周波数成分を描画
plt.show()

yif = sp.fftpack.ifft(yf[:N]).real #ifftで元の波形に戻してみる
plt.plot(t, yif) #戻ったか、確認。最初と同じグラフのはず
plt.show()

yf[10] = 0 #第10次高調波をカット
yf[N-10] = 0 #第10次高調波をカット
#t = np.arange(0, 1, 1/fs)
#1/fs秒刻みの時刻の配列。基本周期T = 1。
yif = sp.fftpack.ifft(yf[:N]).real
plt.plot(t, yif)
plt.show()

yf1 = sp.fftpack.fft(yif[:N]) #10次はカットしてある
plt.plot(np.abs(yf1)) #周波数成分を調べる
plt.show()
#==================

4－8　SciPyのFFT、IFFTを使ったわかりやすい例（4）

#==================
import numpy as np
import scipy as sp
from scipy.fftpack import fft
import matplotlib.pyplot as plt

pi = np.pi

#周期T = 3秒。基本周波数=1/3Hz。
fs = 40 #サンプリング周波数。データ1個で1/40秒。
N = 120 #データ数は3/(1/40)=120
#ｔは1/40秒刻みの時刻の配列。1/40秒ごとに3秒間、サンプリング
t = np.arange(0, 3, 1/40)
tt = np.arange(0, 120) #要素は0から1刻みで119まで


y = np.sin(2*pi*4*t)+2*np.sin(2*pi*10*t)    #周波数は4Hz，10Hz
plt.plot(t, y)     #元の波形を描画
plt.show()

#yf[k]には第ｋ次高調波のフーリエ係数が入る
yf = sp.fftpack.fft(y[:N])
#配列 tt には0，1，2，……, 119 が入っている
#yf[k]には第ｋ次高調波（周波数はｋ／3 Hz）の成分が入るので、 
#横軸の目盛りを= tt / 3.0 とすれば周波数が表示される
plt.plot(tt / 3.0, np.abs(yf)) #含まれる周波数成分を描画
plt.show()
#==================

5－1　ノイズ入りのサウンドからノイズを消す

1本目

#==================
#ノイズ入りのwavファイルを作る。左は440Hz、右は700Hzで10秒間。
#
import scipy as sp
import sys
import numpy as np
import wave
from scipy.fftpack import fft
import matplotlib.pyplot as plt
pi = np.pi
fs = 44100 #サンプリング周波数
s = 10 #10秒の長さのwavファイルを作る
#s秒分のサウンドデータを作るには、fs*sフレームが必要。
N = fs * s #フレーム数（左右ひと組で1フレーム）
t = np.arange(0, s, 1/fs) #横軸の配列を作る
#sin(2π*□＊t) の周波数は□の部分。16ビットで録音するので、
#－32768～32767がデータの範囲。だから－15000～15000の範囲の
#大きさの音を作ることにした。
#
#440Hzの音+ノイズ
ldata = 15000 * np.sin(2*pi*440*t) + 2500 *np.random.randn(t.size) 
#700Hzの音+ノイズ
rdata = 15000 * np.sin(2*pi*700*t) + 2500 * np.random.randn(t.size) 
#配列のサイズは左、右のサイズの和
sounddata = np.zeros(N * 2, dtype = "int16")
sounddata[0::2] = ldata
sounddata[1::2] = rdata
writewave = wave.Wave_write("keshitene.wav") #保存先ファイル
#setparamの引数は
#(nchannels, sampwidth, framerate, nframes, comptype, compname)
#先頭から順にステレオ、サンプルのサイズ2バイト、
#サンプリング周波数fs、フレーム数N。
writewave.setparams((2, 2, fs, N, 'NONE', 'NONE'))
writewave.writeframes(sounddata) #編集したサウンドデータを保存
writewave.close()
sys.exit()
#==================

2本目

#==================
#ノイズ入りのwavファイルを作る。左は440Hz、右は700Hzで10秒間。
#
import scipy as sp
import sys
import numpy as np
import wave
from scipy.fftpack import fft
import matplotlib.pyplot as plt
pi = np.pi
fs = 44100 #サンプリング周波数
s = 10 #10秒の長さのwavファイルを作る
#s秒分のサウンドデータを作るには、fs*sフレームが必要。
N = fs * s #フレーム数
t = np.arange(0, s, 1/fs) #横軸の配列を作る
#sin(2π*□＊t) の周波数は□の部分。16ビットで録音するので、
#－32768～32767がデータの範囲。だから－15000～15000の範囲の
#大きさの音を作ることにした。
#440Hzの音+ノイズ
ldata = 15000 * np.sin(2*pi*440*t) + 2500 *np.random.randn(t.size)
#700Hzの音+ノイズ
rdata = 15000 * np.sin(2*pi*700*t) + 2500 * np.random.randn(t.size)
#■■■
yf = sp.fftpack.fft(ldata[:N])
plt.plot(abs(yf))
plt.show()
#■■■
#配列のサイズは左、右のサイズの和
sounddata = np.zeros(N * 2, dtype = "int16")
sounddata[0::2] = ldata
sounddata[1::2] = rdata
writewave = wave.Wave_write("keshitene.wav") #保存先ファイル
#setparamの引数は
#(nchannels, sampwidth, framerate, nframes, comptype, compname)
#先頭から順にステレオ、サンプルのサイズ2バイト、
#サンプリング周波数fs、フレーム数N。
writewave.setparams((2, 2, fs, N, 'NONE', 'NONE'))
writewave.writeframes(sounddata) #編集したサウンドデータを保存
writewave.close()
sys.exit()
#==================

3本目

#==================
#ノイズ入りのwavファイルを作る。左は440Hz、右は700Hzで10秒間。
#
import scipy as sp
import sys
import numpy as np
import wave
from scipy.fftpack import fft
import matplotlib.pyplot as plt
pi = np.pi
fs = 44100 #サンプリング周波数
s = 10 #10秒の長さのwavファイルを作る
#s秒分のサウンドデータを作るには、fs*sフレームが必要。
N = fs * s #フレーム数
t = np.arange(0, s, 1/fs) #横軸の配列を作る
#sin(2π*□＊t) の周波数は□の部分。16ビットで録音するので、
#－32768～32767がデータの範囲。だから－15000～15000の範囲の
#大きさの音を作ることにした。
#440Hzの音+ノイズ
ldata = 15000 * np.sin(2*pi*440*t) + 2500 *np.random.randn(t.size)
#700Hzの音+ノイズ
rdata = 15000 * np.sin(2*pi*700*t) + 2500 * np.random.randn(t.size)
#■■■
yf = sp.fftpack.fft(ldata[:N])
yf[abs(yf) < 7500000] = 0 #yfの絶対値が750万未満のデータは0に
ldata = sp.fftpack.ifft(yf[:N]).real
yf = sp.fftpack.fft(rdata[:N])
yf[abs(yf) < 7500000] = 0 #yfの絶対値が750万未満のデータは0に
rdata = sp.fftpack.ifft(yf[:N]).real
#■■■
#配列のサイズは左、右のサイズの和
sounddata = np.zeros(N * 2, dtype = "int16")
sounddata[0::2] = ldata
sounddata[1::2] = rdata
writewave = wave.Wave_write("keshitene.wav") #保存先ファイル
#setparamの引数は
#(nchannels, sampwidth, framerate, nframes, comptype, compname)
#先頭から順にステレオ、サンプルのサイズ2バイト、
#サンプリング周波数fs、フレーム数N。
writewave.setparams((2, 2, fs, N, 'NONE', 'NONE'))
writewave.writeframes(sounddata) #編集したサウンドデータを保存
writewave.close()
sys.exit()
#==================

5－2　ノイズ入りの信号から周波数成分を見つける

#==================
import sys
import numpy as np
import scipy as sp
from scipy.fftpack import fft
import matplotlib.pyplot as plt
pi = np.pi
fs = 1000  #サンプリング周波数
N = fs * 4  #サンプル数
t = np.arange(0, 4, 1/fs) #1/fs秒刻みの時刻の配列(4秒間)
#100Hz,440Hz の信号にノイズを混ぜる
y=np.sin(2*pi*100*t)+1.5*np.sin(2*pi*440*t)+3* np.random.randn(t.size)
plt.plot(t, y) #元の波形を描画
plt.show()
#yf[k]には第ｋ次高調波のフーリエ係数が入る
yf = sp.fftpack.fft(y[:N])
#得られる最大の周波数成分は fs = 1000ヘルツ。
#つまりFFTの結果のグラフの右端は fs = 1000ヘルツの成分。
plt.plot(np.abs(yf)) #含まれる周波数成分を描画
plt.show()
sys.exit()
#==================

5－3　肉声を電話の声に変換する

#==================
import sys
import numpy as np
import scipy as sp
from scipy.fftpack import fft
import wave
N = 8192 #FFTするサイズ。2の13乗。
#
#--------------------------------
#サウンドデータの高域成分をカットする関数
#
#databuf: サウンドデータの入ったバッファ
#num1: 第num1次高調波から第num2高調波までを残して消す
def editsound(databuf, num1, num2):
    databuf1 = np.zeros(databuf.size)
    for i in range(int(databuf.size / N)):
        #databufの先頭から順にNデータずつ処理
        databuf2 = databuf[N * i: N * (i+1)] #Nデータを取り出す
        yf = sp.fftpack.fft(databuf2)
        for j in range(1, num1):
            yf[j] = 0 #低周波成分を0に
            yf[N - j] = 0
        for j in range(num2+1, int(N/2)):
            yf[j] = 0 #高周波成分を0に
            yf[N - j] = 0
        #ifftの結果は複素数。もとが実データなので結果の虚数部分は0。
        #コピー先は複素数の配列ではないので実数部分だけをコピーする。
        databuf1[N * i: N * (i+1)] = sp.fftpack.ifft(yf).real
    return databuf1
#--------------------------------
#メインルーチン
wavefile = wave.open("人の声.wav","rb")
#サウンドデータ部分の読み取り
buf = wavefile.readframes(wavefile.getnframes())
buf = np.frombuffer(buf, dtype= "int16") #16bit符号付き整数に変換
params = wavefile.getparams()
wavefile.close()
print(params)
if (params.sampwidth != 2) or (params.nchannels != 2):
    sys.exit() #16ビットで記録していないか、ステレオでなければ終了
sounddata = buf.copy() #サウンドデータ用バッファ
sounddata[0::2] = editsound(sounddata[0::2], 56, 631) #左の音を編集
sounddata[1::2] = editsound(sounddata[1::2], 56, 631) #右の音を編集
writewave = wave.Wave_write("電話の声.wav") #保存先ファイル名
writewave.setparams(wavefile.getparams()) #同じパラメータをセット
writewave.writeframes(sounddata) #編集したデータをファイルへ保存
writewave.close()
sys.exit()
#==================

5ー4　FFTは信用できるのか

#==================
import sys
import numpy as np
import scipy as sp
from scipy.fftpack import fft
import matplotlib.pyplot as plt
print("▼")
pi = np.pi
fs = 200 #サンプリング周波数
N = 600 #サンプル数
t = np.arange(0, N/fs, 1/fs) #1/fs秒刻みの時刻の配列
y = np.sin(2*pi*40*t)+2*np.sin(2*pi*90*t)

#yf[k]には第ｋ次高調波のフーリエ係数が入る
yf = sp.fftpack.fft(y[:N])
plt.plot(np.abs(yf)) #含まれる周波数成分を描画
plt.show()

print("▲")
sys.exit()
#==================

5－5　Pythonでリアルタイムに動くグラフを描く

#==================
#リアルタイムにグラフを描画
#
import sys
import numpy as np
import matplotlib.pyplot as plt
pi = np.pi
t = np.arange(0, 3, 1/40) #横軸
for i in range(1,1000):
    plt.ylim(-5, 5) #y軸の最小値、最大値をセット
    y = 3*np.sin(2*pi*4*t+0.3*i)
    plt.plot(t, y) #波形を描画
    plt.pause(0.1) #pyplot.show()では1回描画して止まってしまう
    plt.clf() #今描いたグラフを消す
sys.exit()
#==================

5－6　変化する周波数成分を表示する

#==================
#リアルタイムに周波数成分を描画
#
import sys
import numpy as np
import scipy as sp
from scipy.fftpack import fft
import matplotlib.pyplot as plt
pi = np.pi
N = 120    #データ数は3/(1/40)=120
# t は1/40秒刻み、0秒以上3秒未満の時刻の配列。
#周期T = 3秒。基本周波数=1/3Hz。サンプリング周波数は40Hz
t = np.arange(0, 3, 1/40)
for i in range(1,100):
    y = np.sin(2*pi*(1/3*30)*t) + np.sin(2*pi*(1/3*10 + 0.1*i)*t)
    yf = sp.fftpack.fft(y[:N])
#yf[k]には第ｋ次高調波のフーリエ係数が入る
    plt.plot(np.abs(yf))    #含まれる周波数成分を描画
    plt.pause(0.1)
    plt.clf()
sys.exit()
#==================

6－1　離散フーリエ係数の値を変える際の注意点

1本目

#==================
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]には第ｋ次高調波のフーリエ係数が入る
for i in range(8):
    print("■yf", i, yf[i])
yif = sp.fftpack.ifft(yf)
for i in range(8):
    print("■■yif", i, yif[i])
sys.exit()
#==================

2本目

#==================
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]には第ｋ次高調波のフーリエ係数が入る
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)
for i in range(8):
    print("■■yif", i, yif[i])
sys.exit()
#==================

3本目

#==================
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]には第ｋ次高調波のフーリエ係数が入る
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)
for i in range(8):
    print("■■yif", i, yif[i])
sys.exit()
#==================

6－2　最も基本的なイコライザを作る

#==================
import sys
import numpy as np
import scipy as sp
from scipy.fftpack import fft
import wave
import datetime
N = 22050 #FFTするサイズ
#--------------------------------
#サウンドデータの、指定した周波数成分をrate倍する
#
#databuf: サウンドデータの入ったバッファ
#num1～num2次高調波の成分をrate倍にする。maxはN/2－1
def editSound(databuf, num1, num2, rate):
    databuf1 = np.zeros(databuf.size)
    for i in range(int(databuf.size / N)):
        #Nデータずつ処理
        yf = sp.fftpack.fft(databuf[N * i: N * (i+1)])
        for j in range(num1, num2+1):
            yf[j] = yf[j] * rate #周波数成分をrate倍に
            yf[N - j] = yf[N - j] * rate
        #ifftの結果は複素数。もとが実データなので結果の虚数部分は0。
        #だがコピー先は実数の配列ではないので実数部分だけをコピーする。
        databuf1[N * i: N * (i+1)] = sp.fftpack.ifft(yf).real
    return databuf1
#--------------------------------
wavefile = wave.open("C:\\Users\\xxxx\\Desktop\\BGM.wav","rb")
#サウンドデータ部分の読み取り
buf = wavefile.readframes(wavefile.getnframes())
buf = np.frombuffer(buf, dtype= "int16") #16bit符号付き整数に変換
params = wavefile.getparams()
wavefile.close()
print(params)
if (params.sampwidth != 2) or (params.nchannels != 2):
    sys.exit() #16ビットで記録していないか、ステレオでなければ終了
sounddata = buf.copy()  #サウンドデータ用バッファ
#高音を編集
sounddata[0::2] = editSound(sounddata[0::2], 500, 5000, 0) #左を編集
sounddata[1::2] = editSound(sounddata[1::2], 500, 5000, 0) #右を編集
#低音を編集
sounddata[0::2] = editSound(sounddata[0::2], 1, 499, 1.2)  #左を編集
sounddata[1::2] = editSound(sounddata[1::2], 1, 499, 1.2)  #右を編集
#保存先
writewave = wave.Wave_write("C:\\Users\\xxxx\\Desktop\\keshitene.wav")  
#保存先ファイルに元と同じパラメータをセット
writewave.setparams(wavefile.getparams())
writewave.writeframes(sounddata)  #編集したサウンドデータを保存
writewave.close()

dt = datetime.datetime.now()
print('▲▲▲▲ {}:{}:{}'.format(dt.hour, dt.minute, dt.second))

sys.exit()
#==================