Sine Wave That Slowly Ramps Up Frequency from F1 to F2 for a Given Time

sine wave that slowly ramps up frequency from f1 to f2 for a given time

if you want angular frequency (w=2 pi f) to vary linearly with time then dw/dt = a and w = w0 + (wn-w0)*t/tn (where t goes from 0 to tn, w goes from w0 to wn). phase is the integral of that, so phase = w0 t + (wn-w0)*t^2/(2tn) (as oli says):

void sweep(double f_start, double f_end, double interval, int n_steps) {
    for (int i = 0; i < n_steps; ++i) {
        double delta = i / (float)n_steps;
        double t = interval * delta;
        double phase = 2 * PI * t * (f_start + (f_end - f_start) * delta / 2);
        while (phase > 2 * PI) phase -= 2 * PI; // optional
        printf("%f %f %f", t, phase * 180 / PI, 3 * sin(phase));
    }
}

(where interval is tn and delta is t/tn).

here's the output for the equivalent python code (1-10Hz over 5 seconds):

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)

ps incidentally, if you're listening to this (or looking at it - anything that involves human perception) i suspect you don't want a linear increase, but an exponential one. but that's a different question...

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

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.

Generating a smooth sinusoidal wave

You can use the trigonometric theorems to get an iteration for the sequence of sine values.

sin(A+B) + sin(A-B) = 2*sin(A)*cos(B)

Thus if you want to generate the sequence of values sin(w*k*dt) then you only have to compute

s[0] = 0, s[1] = sin(w*dt), cc = 2*cos(w*dt)

and then iterate

s[k+1] = cc*s[k] - s[k-1]

---

This linear recursion has eigenvalues on the unit circle and thus accumulates floating point errors, which may lead to phase shift and changes in amplitude over very long time spans. However, locally it will always look like a sine wave.

The second effect can be avoided by iterating the cosine sequence c[k]=cos(w*k*dt) at the same time,

s[k+1] = c[1]*s[k] + s[1]*c[k]
c[k+1] = c[1]*c[k] - s[1]*s[k]

and periodically rescaling the pair c[k],s[k] to have euclidean length 1.

Recursively create a sine wave given a single sine wave value and the period

Lets start with some trigonometric identities:

sin(x)^2 + cos(x)^2 == 1
sin(x+y) == sin(x)*cos(y) + sin(y)*cos(x)
cos(x+y) == cos(x)*cos(y) - sin(x)*sin(y)

Given the sine and cosine at a point x, we can exactly calculate the values after a step of size d, after precalculating sd = sin(d) and cd = cos(d):

sin(x+d) = sin(x)*cd + cos(x)*sd
cos(x+d) = cos(x)*cd - sin(x)*sd

Given the initial sine value, you can calculate the initial cosine value:

cos(x) = sqrt(1 - sin(x)^2)

Note that there are two possible solutions, corresponding to the two possible square-root values. Also note that all the angles in these identities are in radians, and d needs to be negative if you're going back through the wave.

Plot exponential pulses over time

you need to use proper math equation y=f(t) for your signal instead. In that way you will obtain also the amplitude. I do not know the properties of your signal but it does not look like logarithmic. I would construct it like this:

y=( |sin(c0*t)| + sin(c0*t) )^c1 * sin(c2*t) * c3

where:

• t is time
• c0=frequency of gaps / 2*Pi
• c1 is sharpness of gap/signal transition
• c2=frequency of pulses / 2*Pi
• c3 is amplitude of signal

The first part just modulate the base signal with gaps and shaped amplitude shape, the last therm is your base signal.

The average pulses per second is not your signal frequency !!! instead it is around half of it (depends on the signal to gap ratio). Now just increment t in some loop and compute y(t) for your data. To avoid floating point rounding problems limit the t to common period between the 2 sin waves.

Weird glitch when creating sine audio signal

Digging out my rusty math I think it may be because:

Going in L steps from frequency F1 to F2 you have a frequency of

Fi = F1 + i * ( F2 - F1 ) / L

or with

( F2 - F1 ) / L = S

Fi = F1 + i * S

Now to find out how far we have progressed we need the integral over i, which would be

I(Fi) = i * F1 + 1/2 * i ^ 2 * S

Which give or take resembles the term inside your formula for the sine.

Note that you can gain efficiency by moving the constant part (S/2) out of the loop..

concatenate sinusoids of varying frequency

EDIT

I'm also a little confused by your choice of sine function, it doesn't look like the coefficients or duration are tied to real-world values, I would expect to see something like this:

p = np.concatenate((p, np.sin((2 * np.pi * f / sampling_rate) * np.arange(total_tone_time * sampling_rate) + phase)))
phase += 2 * np.pi * f * total_tone_time
phase %= 2 * np.pi  # strip off full cycles to avoid overflow

Where

• total_tone_time is the duration of the current tone being played
• sampling_rate is measured in hz (samples per second)

The purpose of adding in the phase between each term is to keep the resultant composite curve continuous, which might help avoid any pops in the waveform. If you can't hear anything, be sure the frequencies being output are in the audible range: http://en.wikipedia.org/wiki/Audio_frequency

I would try the above suggestion first before trying what I talk about below, which go much deeper down the rabbit hole.

/EDIT

Not an expert, but I think what you say about the discontinuity causing clicking could be correct. If it is, you could do a fast fade out/fade in of amplitude at the junction to avoid the click.

Here are some other possibilities I can think of.

Model and Fit

How are you sampling the frequencies? Are they sampled at specific times, so you can get frequency as a function of time? If so, you could try to fit a curve to the points (which is a non-trivial problem), then take the sine of the integral of the calculated curve:

Sin[2 * PI * Integral[freq(time), time, 0, current_time]]

Numerical Integration

As an alternative to fitting a curve, if your frequency sampling rate is fast enough to approximate a smooth curve, you can use the values directly in a numeric integral. For the below example I'm assuming your data is in the format of [[freq0, time0], [freq1, time1], ...] and that the frequency samples are evenly spaced over time at the same rate you'd like to sample your output waveform.

PI = 3.14159
waveform = [0] * len(data)
phase = 0
time_delta = data[1][1] - data[0][1]
for i, (f, t) in enumerate(data):
    if i != 0:
        phase += 0.5 * (f + last_f) * time_delta
        waveform[i] = sin(2 * PI * phase)
    last_f = f
    phase %= 2 * PI

Note that for the above I'm using the trapezoidal method, described here: http://www.mathworks.com/help/matlab/ref/trapz.html#bua4lsr

The one thing that concerns me about this method is that you are taking the absolute value of frequencies, which suggests to me that your frequencies may not be sampled from a data source yielding a well-behaved, continuous function.

Add a Ramp

Finally, if you are concatenating frequencies from a random sample with no direct time dependency, you could add a frequency ramp to bridge the gap.

You could try linear:

sin(PI * [(f2 - f1) / time_frame * t ** 2 + 2 * f1 * t])

Or exponential:

sin(2 * PI * time_frame * f1 * exp(log(f2 / f1) * t / time_frame) / log(f2 / f1))

Where

• time_frame is the duration you want the ramp to last
• t is chosen to start at zero and end at time_frame
• f1 is the frequency being left
• f2 is the frequency being entered

calculate multiple sine sweeps one after another

As @jaket pointed out in a comment, you must make the phase vary continuously from segment to segment (I'm paraphrasing a bit). Here's a variation of your code that shows one way you could do this. I don't have all your other code, so instead of self, the first argument of LinearSineSweep is a file to which the samples are written as text. (I've also tweaked the code to compensate for the fact that the requested interval will not, in general, be an exact multiple of the sampling period.) numpy are matplotlib are used to create the plot.

from __future__ import print_function, division

import math

def LinearSineSweep(f, fStart, fEnd, samplingTime, samplesPerSecond,
                    t0=0, phi0=0):
    nValues = int(samplesPerSecond * samplingTime)
    actualSamplingTime = nValues / samplesPerSecond
    for i in range(0, nValues):
        delta = float(i) / nValues
        t = actualSamplingTime * delta
        phase = 2 * math.pi * t * (fStart + (fEnd - fStart) * delta / 2)
        value = math.sin(phase + phi0)
        # Write the time and sample value to the output...
        print(t0 + t, value, file=f)
    phase = 2 * math.pi * actualSamplingTime * (fStart + (fEnd - fStart) / 2)
    return t0 + actualSamplingTime, phi0 + phase

if __name__ == "__main__":
    with open('out.csv', 'w') as f:
        t, phi = LinearSineSweep(f, 0, 1700, 0.001, 44100)
        t, phi = LinearSineSweep(f, 1700, 1700, 0.001, 44100, t, phi)
        t, phi = LinearSineSweep(f, 1700, 100, 0.001, 44100, t, phi)

    import numpy as np
    import matplotlib.pyplot as plt

    tvals, v = np.loadtxt('out.csv', unpack=True)
    plt.figure(figsize=(10, 4))
    plt.plot(tvals, v)
    plt.grid()
    plt.show()

Here's the plot:

---

Related Topics