2つの周波数の間で指数関数的に増加する正弦波を生成するPythonメソッドを実装しようとしています。次のPythonコードを使用して、 [この質問]で線形変化が解決されました。
from math import pi, sin
def sweep(f_start, f_end, interval, n_steps):
for i in range(n_steps):
delta = i / float(n_steps)
t = interval * delta
phase = 2 * pi * t * (f_start + (f_end - f_start) * delta / 2)
print t, phase * 180 / pi, 3 * sin(phase)
sweep(1, 10, 5, 1000)
この線形累積位相/デルタ アプローチを指数周波数スイープに変更し、人間の耳に滑らかにする方法。
この種の問題の秘訣は、周波数変調と位相変調の関係を理解することです。これら 2 つは密接に関連しています。f一定の周波数と振幅を持つ正弦波は、次のAように記述できます (Python コードではなく数式)。
x(t) = A sin(2 * pi * f * t)
しかし、これを書く別の方法は、最初にフェーズphiを時間の関数として定義することです。
phi(t) = 2 * pi * f * t
x(t) = A sin(phi(t))
ここで注意すべきことは、周波数fは位相の微分を 2*pi: で割ったものであるということですf = d/dt(phi(t)) / (2*pi)。
時間とともに変化する周波数を持つ信号の場合、同様に瞬時周波数 f_instを定義できます。
x(t) = A sin(phi(t))
f_inst = d/dt(phi(t)) / (2*pi)
あなたがしたいことはこれの反対です。あなたは位相に変換する必要がある特定の瞬間周波数（対数スイープ）を持っています。導出の反対は積分であるため、次のように適切なフェーズを計算できます (式はそのままです)。
phi(t) = 2 * pi * Integral_0_to_t {f_inst(t) dt}
x(t) = A sin(phi(t))
ここで行っているのは、必要な瞬時周波数を取得するための信号 (ゼロ周波数) のある種の位相変調です。これは numpy で非常に簡単に行うことができます:
from pylab import *
n = 1000 # number of points
f1, f2 = 10, 30 # frequency sweep range in Hertz
t = linspace(0,1,1000)
dt = t[1] - t[0] # needed for integration
# define desired logarithmic frequency sweep
f_inst = logspace(log10(f1), log10(f2), n)
phi = 2 * pi * cumsum(f_inst) * dt # integrate to get phase
# make plot
plot(t, sin(phi))
xlabel('Time (s)')
ylim([-1.2, 1.2])
grid()
show()
結果の画像:
しかし (Dave が言及しただましにも記されているように)、おそらく対数スイープではなく、指数スイープが必要です。あなたの耳は周波数を対数的に知覚するため、滑らかで直線的な音階 (ピアノの鍵盤を考えてください) は指数関数的に配置されます。これは、瞬時周波数f_inst(t) = f1 * exp(k * t)を再定義するだけで実現できます。kf_inst(t2) = f2
同時に振幅変調を使用したい場合は、式をA時間依存に変更するだけです。A(t)
Basの答えは素晴らしいですが、実際には分析的な解決策を提供していないため、その部分は次のとおりです...
私が知る限り、sin(Aexp(Bt))whereAやBare 定数のようなものが必要です。0時間が から始まり、続くと仮定しますC(別の時間に始まる場合は、両方からそれを差し引きます)。
それから、バスが言ったように、もし私たちがそのようなsin(g(t))頻度を持っているなら、私は思う. そして、私たちはそれが時と時にあることを望んでいます。f2 * pi * f = dg / dtf00fCC
簡単な計算を行うと (実際には - 学校の最後の年レベルです)、次のようになります。
B = 1/C * log(fC/f0)
A = 2 * pi * f0 / B
1000 サンプルを使用して 5 秒で 1 Hz から 10 Hz に変化するコードを次に示します。
from math import pi, sin, log, exp
def sweep(f_start, f_end, interval, n_steps):
b = log(f_end/f_start) / interval
a = 2 * pi * f_start / b
for i in range(n_steps):
delta = i / float(n_steps)
t = interval * delta
g_t = a * exp(b * t)
print t, 3 * sin(g_t)
sweep(1, 10, 5, 1000)
与える：
(そして、定数を追加することができます - sin(g_t + k)- 開始フェーズを好きな場所に取得できます)。
アップデート
表示されている問題がサンプリングのアーティファクトであることを示すために、オーバーサンプリングを行うバージョンを次に示します (引数として設定した場合)。
from math import pi, sin, log, exp
def sweep(f_start, f_end, interval, n_steps, n_oversample=1):
b = log(f_end/f_start) / interval
a = 2 * pi * f_start / b
for i in range(n_steps):
for oversample in range(n_oversample):
fractional_step = oversample / float(n_oversample)
delta = (i + fractional_step) / float(n_steps)
t = interval * delta
g_t = a * exp(b * t)
print t, 3 * sin(g_t)
sweep(16000.0, 16500.0, 256.0/48000.0, 256) # looks strange
sweep(16000.0, 16500.0, 256.0/48000.0, 256, 4) # looks fine with better resolution
n_oversampleコードを確認すると、4 への設定 (2 番目の呼び出し) がタイムステップに高い解像度を追加するだけであることがわかります。特に、 コード when oversample = 0(ie fractional_step = 0) は以前と同じであるため、2 番目のプロットには、最初のプロットの点に加えて、欠落しているデータを「埋め」、すべてがそれほど驚くべきものにならないようにする追加の点が含まれます。
以下は、元のカーブとオーバーサンプリングされたカーブの開始付近のクローズアップで、何が起こっているかを詳細に示しています。
最後に、この種のことは完全に正常であり、エラーを示すものではありません。デジタル波形からアナログ信号が生成されると、「正しい」結果が得られます (ハードウェアが正しく機能していると仮定します)。 この優れたビデオは、明確でない場合に説明します。
この古い質問は、関連する質問のリスト内の別の質問の右マージンに表示されたときに私の目を引きました。Bas Swinckels と andrew cooke による回答は素晴らしいです。scipy.signal.chirpこの回答を追加して、SciPy にはそのような信号を生成する機能があることを指摘します。周波数の指数関数的変化には、 を使用しますmethod="logarithmic"。(引数は、またはの場合methodもあります。任意の多項式を使用するチャープの場合は、 を使用できます。)"linear""quadratic""hyperbolic"scipy.signal.sweep_poly
chirp以下は、5 秒間で 1000 サンプルで 1 Hz から 10 Hz までのスイープを生成するために使用できる方法です。
import numpy as np
from scipy.signal import chirp
import matplotlib.pyplot as plt
T = 5
n = 1000
t = np.linspace(0, T, n, endpoint=False)
f0 = 1
f1 = 10
y = chirp(t, f0, T, f1, method='logarithmic')
plt.plot(t, y)
plt.grid(alpha=0.25)
plt.xlabel('t (seconds)')
plt.show()
