6 WAVEファイルの内容を FFT して加工 (高い周波数のカット) してみる

一つ一つコピペして実行するのが面倒という人用に、 Mathematica 9 で作ったノートブック

http://nalab.mind.meiji.ac.jp/~mk/lecture/fourier-2017/piano-cutoff.nb

も提供します (保存してから Mathematica で開いて下さい)。

(* ピアノでドミソの和音を弾いてものを録音したデータを取り込む *)
snd = Import["http://nalab.mind.meiji.ac.jp/~mk/fourier/piano-do-mi-so.wav","Sound"]

{left, right} = snd[[1, 1]];
sr = snd[[1, 2]]

(* データを切り出すための関数. t1〜t2, t1〜(t1+1) *)
take[tbl_, t1_, t2_] := Take[tbl, {Floor[t1*sr], Floor[t2*sr]}]
take1[tbl_, t_] := Take[tbl, {Floor[t*sr], Floor[t*sr] + sr - 1}]

(* 1秒後から1秒間のデータを tbl に代入して、プロット *)
g = ListPlot[tbl = take1[left, 1.0], PlotRange -> All]

(* サンプルレートを変えて、再生 *)
ListPlay[tbl]
ListPlay[tbl, SampleRate -> sr]
ListPlay[tbl, SampleRate -> sr/2]
ListPlay[tbl, SampleRate -> 2*sr]
ListPlay[tbl, SampleRate -> Floor[sr*1.5]]

(* 離散Fourier変換してみる *)
c = Fourier[tbl, FourierParameters -> {-1, -1}];
g = ListPlot[Abs[c], Joined -> True, PlotRange -> All]

(* 周波数 f Hz 以上をカットオフするための関数 *)
cutoff[f_] := 
  Join[Table[1, {n, f + 1}], Table[0, {n, sr - 2*f - 1}],  Table[1, {n, f}]];

(* ド、ミ、ソの基本周波数を表示 *)
440.0*2^(-{9, 5, 2}/12)
(* 500ヘルツ以上をカットオフして鳴らしてみる *)
c2 = c*cutoff[500];
g3 = ListPlot[Abs[c2], Joined -> True, PlotRange -> All]
tbl2 = Re[InverseFourier[c2, FourierParameters -> {-1, -1}]];
ListPlay[tbl2, SampleRate -> sr]

(* カットオフする周波数を色々変えてみる *)
do = Re[InverseFourier[c*cutoff[300], FourierParameters -> {-1, -1}]];
domi = Re[InverseFourier[c*cutoff[360], FourierParameters -> {-1, -1}]];
domiso = Re[InverseFourier[c*cutoff[450], FourierParameters -> {-1, -1}]];
domiso2 = Re[InverseFourier[c*cutoff[900], FourierParameters -> {-1, -1}]];
original = Re[InverseFourier[c, FourierParameters -> {-1, -1}]];

(* 鳴らしてみる *)
ListPlay[do, SampleRate -> sr]
ListPlay[domi, SampleRate -> sr]
ListPlay[domiso, SampleRate -> sr]
ListPlay[domiso2, SampleRate -> sr]
ListPlay[original, SampleRate -> sr]