sine wave that exponentialy changes between frequencies f1 and f2 at given time/amount of samples

Question

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 code:

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)

How to change this linear accumulative phase/delta approach to expotential frequency sweep and be smooth to human ear.

Abduce

I think this question is more about math than programming, see: math.stackexchange.com – Disease

This question has the math, and a python implementation: possible duplicate of sine wave that slowly ramps up frequency from f1 to f2 for a given time – Roquelaure

This platform is for asking programming-related question. If you ask your question like as you have asked above, it can be just answered by somebody who knows the particular branch of engineering. Please try to make it as much generic as possible. – Ineluctable

if the linear question was acceptable, why isn't the log one? i'll answer this when i have some free time. i don't think it can be that hard. #11200009 – Islam

@andrewcooke: I don't see any point in saying this question do not fit here but I strongly suggest that OPs should be encouraged to ask questions that reflects the question in programming terms so that anybody willing to answer will less have to deal with the obscure physics of the problem and focus more and to the point to the programming aspect. – Ineluctable

One difference, andrew cooke, would be that in the case of the linear question, the OP added his/her own code to do it, and was actually asking about a problem he/she had with it. No code is provided here. In any case, I find it a nice problem to think a little ;). – Retirement

Answer 1

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

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.

Islam

I did all the math and now I understand it, but I think I'm missing some nuance. When I call sweep(16000.0, 16500.0, 256.0/48000.0, 256) I get strange beginning of sweep: tinypic.com/r/5bpob9/5 – Abduce

[please ignore my earlier replies if you saw them] i am pretty sure this is just an artefact of the sampling. the 16kHz signal is quite close to the nyquist frequency of 24kHz so you're getting poorly sampled sine waves which are giving that strange envelope. if you reduce the frequency range to 8000 to 8500 then it looks better, but you still see something. so it's just "normal" sampling artefacts and not a problem. – Islam

see update - if you oversample the same scan it looks much better. – Islam

Thank You for such comprehensive answer! Now I am facing something different. I want to start my signal at different time because I want to sweep it back and forth f1 to f2 to f1. When I start from 0 to t the k-parameter (in f_inst(t) = f1 * exp(k * t)) is chosen such that f_inst(t2) = f2 and automatically for f_inst(0) = f1 because of exp(k*0) == 1. Now when t > 0 there is no option to start exactly from initial frequency. My second question is if I want to continue with end frequency do I have to add the phase calculated from sweep ? – Abduce

i don't understand the question. but if you want the curves to "join up" then yes, just add the phases as required (see comment just before Update) – Islam

I tested it thoroughly and it seems like you made one crucial mistake. The phase is an integral of frequency, and you put the (exponentially changing) frequency straight into sine. This way it doesn't even start at y=0 for x=0 on your plot. Which @Bas did, both in theory and implementation. – An

It is really weird, even with integral it does not behave as it should. The frequency is just wrong, see edge case: desmos.com/calculator/womya2hq4k – An

Answer 2

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 frequency f_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.

Impervious

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

Endogen