# How can I generate a sine wave with time varying frequency that is continuous, in C?

How can I generate a sine wave with time varying frequency that is continuous? How can I resolve the following problem? I want a continuous graph.

I'm generating this simply like this:

for(int i = 0; i < pcm_buffer_size - 1; i += 2) {
float sample = gain * sin(((float) t * M_PI * 2 * 500) / (float) (sample_rate));

printf("%f\n", sample);

t++;
}

for(int i = 0; i < pcm_buffer_size - 1; i += 2) {
float sample = gain * sin(((float) t * M_PI * 2 * 1000) / (float) (sample_rate));

printf("%f\n", sample);

t++;
}


Where 500 and 1000 is the frequency.

• @Peter K I think closing this for being a coding question might be a bit premature. The OP's attempted solution is not a varying frequency, but a changed one. All that needs to be done for that is to make sure the signal is continuous at the junction by picking the correct phase offset. Smoothness can be achieved by a small overlap and a sliding average. For a varying frequency the answer is simpler. Just use a non-linear function in your trig argument. A quadratic is the next easiest one. E.g. $A \cos( C_0 + C_1 \theta + C_2 \theta^2 )$. – Cedron Dawg Feb 11 at 18:36
• @CedronDawg OK! Reopening.... – Peter K. Feb 11 at 20:17
• The smoothly varying version is known as a "chirp," if you want to look up various popular chirps like a quadratic chirp or a cubic chirp. Each has interesting useful properties. – Cort Ammon Feb 13 at 1:18
• This seems a bit related to the question on Math StackExchange: Terminology re: continuity of discrete a*sin(t). – Mark S. Feb 13 at 22:14

Well lets go lol

@AnonSubmitter85 give to you a nice answer, but let me show my way to do it in matlab, and this maybe can be very easy to port to C:

First I'm creating 256 samples in 500hz sampled at 44100hz take a look how I accumulate the phase and in the end of the first loop I put the phase between the interval 0 and 2pi...

Nice now lets go to the second loop to create more 256 samples of 1000Hz and to do it continuous I'm using the last phase :-), here is the code:

m_phase=0;

signal=[];

f=500;
fs=44100

phaseInc = 2*pi*f/fs;

for i=1:256
signal(i) = sin(m_phase);
m_phase = m_phase + phaseInc;
end

%place the phaser between the 0 and 2pi range
m_phase = mod(m_phase, 2*pi);

f=1000
phaseInc = 2*pi*f/fs;

for i=257:256*2
signal(i) = sin(m_phase);
m_phase = m_phase + phaseInc;
end

%place the phaser between the 0 and 2pi range
m_phase = mod(m_phase, 2*pi);


here is the plot of the above code::

It's a little weird to see so many answers but none that presents an actual answer in C or explains how and why to do it.

The general idea is to maintain a phase that is incremented by a step size that is calculated from the frequency and sample rate. This way you will never get a phase discontinuity.

When doing this, one has to be very careful with accumulated round-off errors in floating point variables because, while the relative round-off error stays mostly the same no matter how large the numbers are, the absolute round-off error increases depending on the magnitude of the number. (There is an in-depth and heavy explanation in the article What Every Computer Scientist Should Know About Floating-Point Arithmetic.)

If you simply keep adding your step size to the phase, it will soon reach magnitudes that are too large, so it has to be kept in check. Since we are dealing with a single waveform we can limit the phase to between 0.0 and 1.0 exclusive by wrapping, or modulo. The easiest and likely most efficient way to do this is to use an unsigned integer for the phase, and let the C compiler take care of the wraps. C is annoying because a lot of arithmetic operations are not defined, but unsigned integer wrap-around is defined the way we want it to, so we can take advantage of that.

The program below outputs a waveform on standard output using the technique described above. The phase is kept in an unsigned integer, incremented with another integer calculated based on the requested frequency and sample rate. Whenever the phase goes above whatever limit set by the integer size, it will automatically wrap to zero.

The integer phase is then scaled to the wanted index into the sine function. This could easily be changed to use a wave-table or an interpolating lookup.

Some parts are unusually complicated because of portability and compatibility. In a fixed setting, a lot can be simplified for readability. Likewise there are places where you might want to add better rounding.

#include <math.h>
#include <stdio.h>

#define M_TWOPI 6.283185307179586476925286766559

/*
* The phase must be an unsigned integer.
* 'maxphase' is for example 65536 if phase_t is 16 bits.
*/
typedef unsigned long phase_t;
double maxphase = (double)((phase_t)0-(phase_t)1)+1.0;

double fs = 44100;

phase_t hz_to_delta( double hz )
{
return maxphase*hz/fs+0.5;
}

float sample_phase( phase_t phase )
{
return sin( phase/maxphase*M_TWOPI );
}

int main( void )
{
long i;
phase_t delta, iphase = 0;

delta = hz_to_delta( 500.0 );
for( i=0; i<fs; ++i )
printf( "%e\n", sample_phase( iphase += delta ) );

delta = hz_to_delta( 1000.0 );
for( i=0; i<fs; ++i )
printf( "%e\n", sample_phase( iphase += delta ) );

return 0;
}

• FYI, this post was originally closed because asking for striclty coding to be done is against SE policy. The OP's code can be fixed with the addition of one line of code between the two loops and a modification to one line in the second loop. I'm not going to code it, but simply set a variable for the phase offset at the junction and add it into the argument in the second loop. I just chimed in to make a distinction between "varying" (from the title) and "change in" (from the code example). – Cedron Dawg Feb 12 at 18:16
• @CedronDawg I see, though I still think that the "integer phase" trick is fundamental enough to practical DSP implementations to at least warrant a mention. – pipe Feb 12 at 18:30
• It should always be considered, but isn't always determinative. The extra divides and type conversions you pay may not be worth it, nor a bunch of extra function calls. On a side note, your code would benefit by replacing your divides with multiplications of their stored reciprocals. – Cedron Dawg Feb 12 at 18:37

The easiest way to do this is to note that frequency is by definition the derivative of the phase. Thus, you can define what the frequency of each sample is and then integrate that to get the phase. This keeps the phase continuous. For example, in matlab this would look something like so:

% The two frequencies of interest.
f1 = 500;
f2 = 1000;

% Sample frequency and locations.
fs = 50 * max( f1, f2 );
T = 0.01;
t = 0 : 1/fs : T;

% Make the first half of the samples have frequency f1
% and the rest have frequency f2.
f = zeros( size( t ) );
f( 1 : floor(length(t)/2) ) = f1;
f( floor(length(t)/2) + 1 : end ) = f2;

% Now integrate the above to get the phase.
phi = cumsum( 2*pi * f ) / fs;

x = cos( phi );
plot( t, x );


• “frequency is by definition the derivative of the phase” – I would disagree with that. Frequency is by definition the number of repeats in a signal per time unit. The derivative of phase is – by definition – the phase velocity. It so happens that frequency is equal to ω/2π, though arguably that's specific to sinusoidals. – leftaroundabout Feb 13 at 15:11
• I can't say that I have ever heard the term phase velocity before. Perhaps instantaneous frequency would be a more accurate term for how I am using it. Every field has its own argot I suppose. – AnonSubmitter85 Feb 13 at 15:22
• Actually I meant angular frequency (the German term is Winkelgeschwindigkeit, lit. “angle velocity”). So, ok, it is called a “frequency” in English, but I would argue frequency and angular frequency should still be kept apart. If nothing else, then at least because of the 2·π factor. – leftaroundabout Feb 13 at 16:03
• @leftaroundabout: fun fact: phase velocity is a term in physics, for waves propagating through space or a medium (e.g. an electrical signal in a transmission line). Phase velocity = angular velocity / wavenumber, and has dimensions of distance / time, not angle / time. en.wikipedia.org/wiki/Phase_velocity = how fast the apparent "crest" of a wave moves. There are some related terms like "group velocity", see wikipedia for an animation. – Peter Cordes Feb 14 at 3:16

In addition to incrementing phase (instead of incrementing time and multiplying by frequency, potentially causing jaggies), also note that the input to your trig function might need to be range limited (or wrapped) to prevent range reduction loss of precision in the trig function implementation.

After incrementing by delta phase, I usually limit phase (by circular wrapping) to between -2pi and 2pi, or -pi to +pi, etc.

• Thanks for your answer. What do you mean by range reduction loss of precision? – deltafft Feb 12 at 11:30
• If the argument gets large, relative to the number of bits in a floating point mantissa, the fractional bits of the phase (small phase angle deltas) get increasingly quantized or even disappear. – hotpaw2 Feb 12 at 15:06

In general, a pure steady real sinusoidal signal can be modeled as:

$$x(t) = A \cos( \theta t + \phi )$$

The argument $$\theta t + \phi$$ is usually called the "phase" of the signal which can be confusing as the same term is also used for the $$\phi$$, so let's call it the phase function.

For a steady frequency, the phase function is linear, and the instantaneous frequency is the same as the constant frequency. In the stipulated situation, you have to make sure the value of the phase function matches at your junction and your signal will be "continuous" (actually a meaningless term in a discrete context.)

For a varying signal, you want a phase function that is not linear. The simplest, as mentioned in my comment, is to use a quadratic. Generalizing the coeffiecient names, it can be formulated as:

$$x(t) = A \cos( C_0 + C_1 t + C_2 t^2 )$$

This cause the frequency to vary at a linear rate. Other functions are also possible. Matching this to your situation and coding it is your responsibility.

For embedded devices the high level answers do not suffice, this would be my approach, a lookup table with phase counting. Interpolation on a different timing loop can be used to smooth out values between lookups. By using fixed lookups a step change in frequency is accomodated smoothly for free

You only need 1/4 of the wave lookup table to generate a reasonable sine wave, but this example has a memory inneficient full wave lookup in order to simplify code

The frequency is determined by the next tick timeout

#define F1 500 // in hz
#define F2 1000
#define SAMPLES 255
#define CUTOVER 10*SAMPLES // in phase cycles

#define PHASE_STEP(f) 1.0/(f*SAMPLES)

uint32_t sine[SAMPLES] = {0, ... }; //fixed point sine samples
uint32_t output; // used to drive DAC/DSP output value

void isr()
{
static uint8_t phase = 0;
output = sine[phase];
phase++;
}

void main()
{
// can be a timer or otherwise
for(unsigned int cycle = 0; cycle < CUTOVER; cycle++) {
isr()
sleep(PHASE_STEP(F1));
}

while (1) {
isr()
sleep(PHASE_STEP(F2));
}
}


• A lookup table does work, sure, but it inevitably has some discontinuities at the nodes which can be a problem. CORDIC has discontinuities too, and isn't ideally suited to floating-point. Approximation to a polynomial can give very good results for surprisingly little calculation though - OlliW's "fast functions" appears to be the first good reference for this, but I've derived higher-order versions with the same basic principles and got very good results. – Graham Feb 13 at 14:03
• @Grahan good points Certainly, lookup based example is extremely simplified and in practice many optimizations and optimized data structures can be done, but yes fundamentally if the function can be polynomialized (e.g. only sine output) this is inexpensive in fixed point as well. – crasic Feb 14 at 1:22

Maybe you'd be interested in some ASCII graphics as well.

See my complete C solution to generate plots in C using printfs here https://stackoverflow.com/questions/53898625/plot-of-the-exponential-function/53899222#53899222

The result is something like this:

Change the following as per your need:

double equation(const double x)
{
return exp(-0.1 * x) * sin(0.5 * x);
}


EDIT:

To get a sinusoid of growing frequency, one could use the following equation as a starting point:

double equation(const double x)
{
const double xx = x - LIMIT_X_MIN;

return sin(0.05 * xx * xx);
}


• I got curious, but I must be doing something wrong, because the program just hangs for me with a maxed CPU core. I literally copy-pasted your code, tried less width, height, points, same thing. I use Archlinux and compile it with gcc test.c -lm -o test, no errors or messages. It's the same with normal 100x40 terminal window, maximized, in console without X. – a concerned citizen Feb 12 at 21:39
• The code you've shown changes the amplitude of the signal over time, not the requency, as was asked. – ambitiose_sed_ineptum Feb 12 at 21:59
• @ambitiose_sed_ineptum This code is meant as an example from one of my previous answers. One needs to change the equation to get the correct results. In this case that would be removing the exponential. – John Feb 13 at 10:48
• @aconcernedcitizen Good catch. I had a bug in there. Please see the revised code. The only change is instead of get_next_x(x) you need x = get_next_x(x). – John Feb 13 at 10:59
• It works now. As for the downvote, whoever did it probably considered that it does not directly address the OP, instead of this answer being an addition to it. Just a guess. – a concerned citizen Feb 13 at 11:27

Little bit of automation

clear all;
close all;

m_phase         = 0;
frequency       = 500;
fs              = 44100

signal          = [];

samples         = 2^8;

initial         = 1;

for k=1:3

frequency     = k*frequency
phaseInc      = 2*pi*frequency/fs;

if(k == 1)
initial     = 1;
else
initial     = ((k-1)*samples) + 1;
endif

for i=initial:k*samples

signal(i)  = sin(m_phase);
m_phase    = m_phase + phaseInc;

end

m_phase       = mod(m_phase, 2*pi);

endfor

subplot(211);
plot(signal);

N               = 512;

X               = fftshift(fft(signal,N));

df              = fs / N;
sampleIndex     = -N/2:N/2-1;
relative_f      = sampleIndex * df

subplot(212);
plot(relative_f, abs(X));
axis([min(relative_f) max(relative_f)])



• But I wonder what can we infer from FFT? – jomegaA Feb 13 at 20:19
• I think there are some things missing to show the frequencies in your FFT – ederwander Feb 13 at 22:06