Sine Wave That Exponentialy Changes Between Frequencies F1 And F2 At Given Time/amount Of Samples

I'm trying to implement Python method that generates sine wave, which ramps up between two freq exponentially. Linear change was solved in [this question] with following Python cod

Solution 1:

The trick in this sort of problems is to understand the relation between frequency modulation and phase modulation, these two are closely related. A sine with a constant frequency f and amplitude A can be described as (formulas, not python code):

x(t) = A sin(2 * pi * f * t)

but a different way to write this is by first defining a phase phi as a function of time:

phi(t) = 2 * pi * f * t
x(t) = A sin(phi(t))

The thing to note here is that frequency f is the derivative of the phase, divided by 2*pi: f = d/dt(phi(t)) / (2*pi).

For a signal which has a frequency that is varying in time, you can similarly define an instantaneous frequencyf_inst:

x(t) = A sin(phi(t))
f_inst = d/dt(phi(t)) / (2*pi)

What you want to do is the opposite of this, you have a given instantaneous frequency (your logarithmic sweep), which you need to convert into a phase. Since the opposite of derivation is integration, you can calculate the appropriate phase like this (still formulas):

phi(t) = 2 * pi * Integral_0_to_t {f_inst(t) dt}
x(t) = A sin(phi(t))

What you are doing here is some sort of phase modulation of a signal (with zero frequency) to obtain the required instantaneous frequency. This is pretty easy to do in 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()

Resulting image:

But (as also noted in the dupe mentioned by Dave), you probably don't want a logarithmic sweep, but an exponential one. Your ear has a logarithmic perception of frequency, so a smooth/linear musical scale (think the keys on a piano) are therefore spaced exponentially. This can be achieved by simply redefining your instantaneous frequency f_inst(t) = f1 * exp(k * t), where k is chosen such that f_inst(t2) = f2.

If you want to use amplitude modulation at the same time, you can simply change A to a time dependent A(t) in the formulas.

Solution 2:

Bas's answer is great, but doesn't actually give an analytic solution, so here's that part...

As far as I can tell, you want something like sin(Aexp(Bt)) where A and B are constants. I'll assume time starts at 0 and continues to C (if it starts at some other time, subtract that from both).

Then, as Bas said, I think, if we have sin(g(t)) frequency f is such that 2 * pi * f = dg / dt. And we want that to be f0 at time 0 and fC at time C.

If you go through the maths, which is easy (it really is - last year of school level), you get:

B = 1/C * log(fC/f0)
A = 2 * pi * f0 / B

and here's some code that goes from 1 to 10Hz in 5 seconds using 1000 samples:

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)

which gives:

(and you can add in a constant - sin(g_t + k) - to get the starting phase wherever you want).

Update

To show that the issue you are seeing is an artefact of sampling, here's a version that does oversampling (if you set it as an argument):

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 inrange(n_steps):for oversample inrange(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

If you check the code you'll see that all that setting n_oversample to 4 does (the second call) is add a higher resolution to the timesteps. In particular, the code when oversample = 0 (ie fractional_step = 0) is identical to before, so the second plot includes the points in the first plot, plus extra ones that "fill in" the missing data and make everything look much less surprising.

Here's a close-up of the original and the oversampled curve near the start, showing what is happening in detail:

Finally, this kind of thing is completely normal and does not indicate any kind of error. When an analogue signal is generated from the digital waveform you'll get "the right" result (assuming the hardware is working right). This excellent video will explain things if they are not clear.

Solution 3:

This old question caught my eye when it showed up in the right margin of another question in the list of related questions. The answers by Bas Swinckels and andrew cooke are great; I'm adding this answer to point out that SciPy has the functionscipy.signal.chirp for generating such a signal. For the exponential variation of the frequency, use method="logarithmic". (The argument method can also be "linear", "quadratic" or "hyperbolic". For a chirp using an arbitrary polynomial, you can use scipy.signal.sweep_poly.)

Here's how you can use chirp to generate the sweep from 1 Hz to 10 Hz in 1000 samples over 5 seconds:

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()