Implementing a sine wave with linearly changing frequency

Hi everyone

I was trying to implement a sine wave that changes its frequency linearly between f_start and f_end in a specific time t. So let's say:

f_start = 50 Hz

f_end = 100 Hz

t = 10 s

That means, i want the frequency of my sine to change linearly from 50 to 100 Hz in 10 seconds. In my mind, that seemed to be easy, so I started with the following approach:

clear;
t = 0:0.00001:10;
f_in_start = 50;
f_in_end = 100;
f_in = linspace(f_in_start, f_in_end, length(t));
y = sin(2*pi * f_in .* t);
plot(t,y)
f_out = 1./(2*diff(find_y_zeros(t,y)));
f_out_start = f_out(1);
f_out_end = f_out(end);

find_y_zeros is a function that I wrote to calculate the times where my y goes 0 (via interpolation). I tested this function seperately and everything seemed to work fine.

However, when I calculate the differences between the zeros, multiply that by 2 (to get the period length) and divide 1 by this, the frequencies at the beginning and the end of the oscillation are 50 and 150 Hz which was definitely not what I wanted above.

Then, by 'try-and-error', I found out that I could get what I wanted by changing line 4 of my code above to this:

f_in_end = f_in_start + 50/2;

Generally speaking, it seems that I get what I want when I divide the difference between my aimed f_start and f_end by 2. So, in this case, to get from 50 to 100 Hz, the difference is 50, and that is divided by 2, so I get f_in_start = 50 and f_in_end = 75.

From the distances of the zeros I can see that now, the frequency of my sine goes from 50 to 100 Hz and that seems to work on any other example (other frequencies and other times) too.

So... it finally works. But I want to know why. Could anybody tell my why the first naive approach did not work? I would be grateful for any hint on this.

Thanks in advance.

Answer 1

[Alfonso Nieto-Castanon]

the formula:

 sin( 2*pi*f(t)*t )

does not result in the desired sine wave with time-varying frequency. Rather the appropriate formula would use, instead of f(t)*t, the integral between 0 and t of f(t):

 sin( 2*pi* ∫f(t)dt )

Only when f(t) is a constant f value, its integral is f*t, and the sine wave is the familiar sin(2*pi*f*t).

The rational for this is that the sinusoid frequency is the rate-of-change / derivative of the sinusoid phase. If you are given a desired frequency profile f(t) you need to integrate that to compute the desired sinusoid phase (up to a constant additive term), from which you can then compute your desired waveform.

One possible way to produce your desired signal would be:

 Fs = 100000;
 t = 0:1/Fs:10;
 f_in_start = 50;
 f_in_end = 100;
 f_in = linspace(f_in_start, f_in_end, length(t));
 phase_in = cumsum(f_in/Fs);
 y = sin(2*pi*phase_in);
 plot(t,y)

[Alfonso Nieto-Castanon]

and to clarify your specific question, in your code, the instantaneous frequency signal f_in is:

 f(t) = f1+(f2-f1)*t/10

(a linear transition between f1 at t=0 and f2 at t=10). The phase of your sinusoid, as constructed by your code, is:

 f(t)*t = ( f1+(f2-f1)*t/10 )*t = f1*t + (f2-f1)*t^2/10

and its derivative at t=10 (the end frequency) is:

 f1 + (f2-f1)*t/5 = f1 + (f2-f1)*2 = 2*f2-f1

(i.e. 150Hz in your example with f1=50 and f2=100)

while the proper phase should look like:

 ∫f(t)dt = ∫f1+(f2-f1)*t/10 dt = f1*t + (f2-f1)*t^2/20

and its derivative at t=10 would be:

f1 + (f2-f1)*t/10 = f1 + (f2-f1) = f2

(100Hz in your example, with f1=50 and f2=100), which explains the strange "divide by 2" effect you were observing...

[Mauricio Galván García Luna]

What is making the phase in here

 Fs = 100000;
 t = 0:1/Fs:10;
 f_in_start = 50;
 f_in_end = 100;
 f_in = linspace(f_in_start, f_in_end, length(t));
 phase_in = cumsum(f_in/Fs);
 y = sin(2*pi*phase_in);
 plot(t,y)

Answer 2

[Alfonso Nieto-Castanon]

someone deleted his own answer including my own comments to it so, just in case these are interesting to anybody else, I am repeating a hopefully illustrative example below:

Using sin(2*pi*f(t)*t) instead of sin(2*pi*∫f(t)dt) is not really an uncommon mistake at all when building frequency-modulated signals. This is aggravated by the fact that, when f(t) is a polynomial function, the former leads to an "almost correct-looking" solution. Trying instead a non-polynomial function (e.g. f(t)=2+sin(t), a frequency modulated tone oscillating between 1Hz and 3Hz) one is able to see the differences between the two methods quite more clearly.

In the example below, we are trying to build a frequency modulated signal with frequency oscillating between 1Hz and 3Hz. The top graph shows the signal built using sin(2*pi*f(t)*t) which does not show the expected 1Hz to 3Hz frequency oscillation at all (things get even worse it you continue plotting beyond t=10s, as the signal frequency starts making wilder and wilder oscillations, then it starts aliasing, etc.), the middle graph shows the integral approximation using the cumsum trick sin(2*pi*cumsum(f(t))/Fs), and the bottom graph shows the exact-integral solution sin(2*pi*∫f(t)dt))

 t = 0:1/1000:10;
 f = 2+sin(t);
 y_wrong  = sin(2*pi*f.*t);
 y_approx = sin(2*pi*cumsum(f)/1000);
 y_exact  = sin(2*pi*(2*t-cos(t)));
 subplot(311); plot(t,y_wrong,'.-')
 subplot(312); plot(t,y_approx,'.-');
 subplot(313); plot(t,y_exact,'.-');

[André Costa]

Thanks a lot! I was geting almost the same problem you pointed!

[Farah Akhtar]

Thanks a lot for your clear description. Yet, I find that my code is still not working properly and I hoped you'd be able to see what I was doing wrong.

I'm trying to generate a cosine with varying frequency and plotting it against frequency. The physical significance of the problem is as follows:

I have a periodic function in frequency domain for which I'm trying to calculate its Fourier series. Consequently, I wish to plot the function in frequency domain, the calculate its Fourier Series and reconstruct it using the Fourier coefficeints. The function is as follows:

|V(f,ΔR)|=sqrt(1+Λ^2+2Λcos(2πfΔR/c))

I'm trying to plot this function from frequency 10MHz to 20MHz. Hence, it should be a sinosoid of increasing frequency as we increase f. For the first step in plotting the entire function, I wanted to render a varying frequency cosine of linearly varying frequencies but that was not successful.

As your comments say, the varying frequency cosine function used in this case would be a ramp. Therefore, the exact integral calculation would require an integral of the ramp which is a parabolic function. I used the following code to calculate and plot the function and get a results that I'm almost certain is incorrect. For one, the amplitude of the cosin is varying which, I believe, it should not and secondly, the result is also not quite smooth as your case. Following is my matlab code:

f_range = 10e6:0.1e6:20e6;
cos_exact = cos((2*pi*delta_R/c).*f_range(1) +(2*pi*delta_R/c).*(f_range(end) - f_range(1)).*f_range.^2./2);

figure(1);
f_range_label = f_range./1e6;
plot(f_range_label,cos_exact,'m');
title('VSWR vs Frequency');
xlabel('Frequency f (MHz)');
ylabel('VSWR (-)');

The resulting plot is posted here:

[Talmon Alexandri]

Thanks a lot for the nice explanation and good example. Helped me so much!!!!!

[pluton schmidt]

Very nice explanation! I am wondering whether there is a programming trick (or any other trick) able to accelerate the computation of the integral? It looks to me that in real-time applications, this integral is an issue and there might be a trick somewhere that would be inaudible to the human ear.