ブログ

scipyのinterpolateは今まで使ったことがなかったが，事前にx, yの２つの数値列を与えておくと補間するための関数を返してくれるという仕様になっている．この関数として返すという仕様がfminなどで使うときの便利さにつながる．
さらに，interpolateはデータの範囲が超えたときにエラーにしない仕様がある．例えば以下のような例である

# -*- coding: utf-8 -*-
import numpy as np
import scipy as scipy
from scipy.optimize import fmin
from scipy import interpolate
import matplotlib
import matplotlib.pyplot as plt

x = np.linspace(0, 100, num=1001, endpoint=True);
y = np.zeros(x.shape);

for i in range(y.shape[0]):
    if i < 500 and i > 50:
        y[i] = 1;
f4 = scipy.interpolate.interp1d(x, y, 'cubic', fill_value='extrapolate')
xi = np.linspace(-100, 200, num=121)

plt.plot(x, y, linewidth='3');
plt.plot(xi, f4(xi))
plt.grid(True)

plt.savefig('20190727test.png', format='png', dpi=100)
plt.show();

さらに数値最小化の追加の例までついかすると例えば以下のようになる．

# -*- coding: utf-8 -*-
import numpy as np
import scipy as scipy
from scipy.optimize import fmin
from scipy import interpolate
import matplotlib
import matplotlib.pyplot as plt

x = np.linspace(0, 100, num=1001, endpoint=True);
y = (x - 40) ** 2

f2 = scipy.interpolate.interp1d(x, y, 'cubic', fill_value='extrapolate')

minimum = scipy.optimize.fmin(f2, 1)

xi = np.linspace(-100, 200, num=121)

plt.plot(x, y, linewidth='3');
plt.plot(xi, f2(xi))
plt.plot(minimum, f2(minimum), '+')
plt.grid(True)

plt.savefig('20190727test2.png', format='png', dpi=100)
plt.show();

導関数を使わない数値最小化ならこれらでもよいが，微分係数を使って勾配を下に下らせるように最適化を目指すアルゴリズムを使用したい場合には，計算が増えてしまう．そこで，UnivariateSplineというものもある．まず上の置き換えだけなら，

from scipy.interpolate import UnivariateSplie
f2 = scipy.interpolate.UnivariateSpline(x, y, s=0, ext=0)

にするだけである．さらにf2.derivatives()なるものが使えるようになる
f2.derivatives(minimum)は[ほぼ0, ほぼ0, 2, ほぼ0]になるはず．順にf2の補間結果，1次微分係数，2次微分係数,...

音を入力出力する例としては，sd.playrec(z, samplerate=fs, channels=1)に変えるだけである．
下の例ではz2が録音されたデータの列になっている．

注意すべき点は，秒単位ではタイミングがずれることである．タイミング関係は使用している環境や負荷によるかもしれない．

import numpy as np
import sounddevice as sd
import matplotlib.pyplot as plt

fs = 44100
fbase = 4000;
t = np.linspace(0, 44099, num=44100, endpoint=False);
z = 0.2 * np.sin(t * 2 * np.pi * (fbase)) + 0.2 * np.sin(t * 2 * np.pi * (fbase + 1)) + 0.2 * np.sin(t * 2 * np.pi * (fbase + 2)) + 0.2 * np.sin(t * 2 * np.pi * (fbase + 3));

z2 = sd.playrec(4 * z, fs, channels=1)#sd.rec(int(duration * fs), samplerate=fs, channels=1)
sd.wait()

sd.play(z2, fs);
# 表示
Fz = np.fft.fft(z2) / z2.shape[0] * 2; # 折り返すのでパワーが２分の１になっている．
Fz[0] = Fz[0] / 2; # 平均成分は折り返さない．
Fz_abs = np.abs(Fz);

freq = np.fft.fftfreq(len(z), d=1.0/fs);
tim = np.arange(0, z2.shape[0]) * 1.0 / fs;
plt.figure(1)
plt.subplot(211);
plt.plot(tim, z2);

plt.subplot(212);
plt.plot(freq, Fz_abs * 1e6);
plt.xlim([0, 20000])

plt.show();

参考： sounddeviceでPythonを使って録音する

以下のようになる．実は，折り返しが始まった周波数より上はnumpy.fft.fftfreqで作った方が正しい．というのはlinspaceで周波数リストを作ると折り返しは再現されないのだ．
fftfreqで作成すると絶対値だけは周波数を再現する．負符号で折り返しを表現する．

import numpy as np
import matplotlib.pyplot as plt

# 例の信号を作成
t = np.linspace(0.001, 4, 4000);
z = 0.1 + 0.2 * np.sin(t * 10 * 2 * np.pi) + 0.2 * np.sin(t * 33 * 2 * np.pi);

# サンプリング周波数
fsmp = 1 / (t[1] - t[0]);
# 解析時間
t_len = t.max() - t.min();

Fz = np.fft.fft(z) / z.shape[0] * 2; # 折り返すのでパワーが２分の１になっている．
Fz[0] = Fz[0] / 2; # 平均成分は折り返さない．
Fz_abs = np.abs(Fz);

freq = np.linspace(0, fsmp - 1/(t_len + t[1] - t[0]), num=len(t));
freq2 = np.fft.fftfreq(len(t), d=1.0/fsmp);

plt.figure(1)
plt.plot(freq, freq2 - freq);

plt.show();

音を出力する例としては以下のようになる．
zがnumpy.arrayでこの例では正弦波の重ね合わせになっている．

import numpy as np
import sounddevice as sd

fs = 44100

t = np.linspace(0, 44099, num=44100, endpoint=False);
z = 0.2 * np.sin(t * 2 * np.pi * 11000) + 0.2 * np.sin(t * 2 * np.pi * 11001) + 0.2 * np.sin(t * 2 * np.pi * 11002) + 0.2 * np.sin(t * 2 * np.pi * 11003) + 0.2 * np.sin(t * 2 * np.pi * 11004) + 0.2 * np.sin(t * 2 * np.pi * 11005);
sd.play(z, fs)

この例は再生だけなので後で，録音も以下のように行える．

import numpy as np
import sounddevice as sd

fs = 44100

duration = 5  # seconds
z2 = sd.rec(int(duration * fs), samplerate=fs, channels=1)
sd.wait()

sd.play(z2, fs); #　録音したものを再生

参考： sounddeviceでPythonを使って録音する

Visual Studio 2019ではExpressは提供されない．
Communityへの切り替えを促されている．SharpDevelopの開発が終了した今Expressが無償で手に入る.NETアプリケーションの商用開発の最新ツール的な位置づけだったが，また状況が変化した．

学習用途では無料で問題ないが，ライセンスが若干異なるので，仕事で使う分にはライセンスの比較が必要だ．

必要に応じてProfessionalのライセンスが必要になる案件が生じるということかもしれない．
最近知ったが，Visual Studio MarketplaceにてVisual Studio Professionalの月額サブスクリプションを取り扱っているらしい．本当に必要な時にはこちらも調査したい．現在の為替だと月５千円強となる．
通常Visual Studio Professionalを購入する場合に比べると一度に支払う金額が抑えられるか...