How can I generate a sine wave with consistent "vibrato"

I am trying to create a .wav file which contains a 440Hz sine wave tone, with 10Hz vibrato that varies the pitch between 430Hz and 450Hz. Something must be wrong with my approach, because when I listen to the generated .wav file, it sounds like the "amplitude" of the vibrato (e.g. the highest/lowest pitch reached by the peaks and troughs of the waveform of the vibrato) just progressively increases over time, instead of staying between 430-450Hz. What is wrong with my approach here? Here is some minimal python code which illustrates the issue:

import math
import wave
import struct

SAMPLE_RATE = 44100

NOTE_PITCH_HZ = 440.0        # Note pitch, Hz
VIBRATO_HZ = 10.0             # Vibrato frequency, Hz
VIBRATO_VARIANCE_HZ = 10.0    # Vibrato +/- variance from note pitch, Hz

NOTE_LENGTH_SECS = 2.0      # Length of .wav file to generate, in seconds

NUM_SAMPLES = int(SAMPLE_RATE * NOTE_LENGTH_SECS)

# Generates a single point on a sine wave
def _sine_sample(freq: float, sine_index: int):
    return math.sin(2.0 * math.pi * float(freq) * (float(sine_index) / SAMPLE_RATE))

samples = []
for i in range(NUM_SAMPLES):
    # Generate sine point for vibrato, map to range -VIBRATO_VARIANCE_HZ:VIBRATO_VARIANCE_HZ
    vibrato_level = _sine_sample(VIBRATO_HZ, i)
    vibrato_change = vibrato_level * VIBRATO_VARIANCE_HZ

    # Mofidy note pitch based on vibrato state
    note_pitch = NOTE_PITCH_HZ + vibrato_change
    sample = _sine_sample(note_pitch, i) * 32767.0

    # Turn amplitude down to 80%
    samples.append(int(sample * 0.8))

# Create mono .wav file with a 2 second 440Hz tone, with 10Hz vibrato that varies the
# pitch by +/- 10Hz (between 430Hz and 450Hz)
with wave.open("vibrato.wav", "w") as wavfile:
    wavfile.setparams((1, 2, SAMPLE_RATE, NUM_SAMPLES, "NONE", "not compressed"))

    for sample in samples:
        wavfile.writeframes(struct.pack('h', sample))

A more straight forward approach that will accomplish what you want is to use a phasor (linear ramp that goes from 0 to 1 then shoots back down to 0) to look up the sin of that value. Then, you can control the amount the phasor increments (the frequency of vibrato).

Here is the code. I lowered the sampling rate to make it easier to look at:

import math
import matplotlib.pyplot as plt

SAMPLE_RATE = 10000

NOTE_PITCH_HZ = 100.0        # Note pitch, Hz
VIBRATO_HZ = 20.0             # Vibrato frequency, Hz
VIBRATO_VARIANCE_HZ = 20.0    # Vibrato +/- variance from note pitch, Hz

NOTE_LENGTH_SECS = 2.0      # Length of .wav file to generate, in seconds

NUM_SAMPLES = int(SAMPLE_RATE * NOTE_LENGTH_SECS)

# Generates a single point on a sine wave
def _sine_sample(freq: float, sine_index: int):
    return math.sin(2.0 * math.pi * float(freq) * (float(sine_index) / SAMPLE_RATE))

phasor_state = 0
phasored_samples = []
samples = []
unmodulated_samples = []
for i in range(NUM_SAMPLES):

    # Generate sine point for vibrato, map to range -VIBRATO_VARIANCE_HZ:VIBRATO_VARIANCE_HZ
    vibrato_level = _sine_sample(VIBRATO_HZ, i)
    vibrato_change = vibrato_level * VIBRATO_VARIANCE_HZ

    # Mofidy note pitch based on vibrato state
    note_pitch = NOTE_PITCH_HZ + vibrato_change
    samples.append(_sine_sample(note_pitch, i)+5)
    unmodulated_samples.append(_sine_sample(NOTE_PITCH_HZ, i))
    phasored_samples.append(math.sin(2*math.pi*phasor_state)+10)
    phasor_inc = note_pitch/SAMPLE_RATE
    phasor_state += phasor_inc
    if phasor_state>=1:
        phasor_state -=1
plt.plot(unmodulated_samples, label='unmodulated')
plt.plot(samples, label='not working')
plt.plot(phasored_samples, label='using phasor')
plt.legend()
plt.show()

A zoom in on the output shows you the difference between these approaches:

Keep in mind though, that this still isn't quite right. A violinist or vocalist will vibrate up and down in a more or less linear trajectory, not a sinusoidal one. To be more 'correct' (if that is what you are going for, that is) would be to compute the change in phase increment as a triangle wave, not a sinusoidal one.

Your mistake was not integrating.

You can't just take instantaneous frequency and multiply it by elapsed time, nor can you just take the sine of a frequency without any time component in it. When frequency changes, simply multiplying would ignore the history of the oscillation (how much it has spun already). Integrating takes care of the "distance traveled" (spun).

Here's a script using numpy. I used smaller values so the plotting looks readable.

Approach:

• generate per-sample array of time. generally useful.
• generate per-sample array of frequency. that's going to wobble.
• integrate the frequency over time, to get phase.
• scale to get radians (it used to be cycles)
• slap a sin() on it to get amplitude

Some parameters:

NOTE_LENGTH_SECS = 0.2      # Length of .wav file to generate, in seconds
SAMPLE_RATE = 10000

NOTE_PITCH_HZ = 100.0        # Note pitch, Hz
VIBRATO_HZ = 20.0             # Vibrato frequency, Hz
VIBRATO_VARIANCE_HZ = 20.0    # Vibrato +/- variance from note pitch, Hz

Generating the oscillation:

t = np.arange(0, NOTE_LENGTH_SECS, 1 / SAMPLE_RATE)

freq = NOTE_PITCH_HZ + VIBRATO_VARIANCE_HZ * np.sin(2 * np.pi * VIBRATO_HZ * t)

# integrate. `∫ f dt` where dt = 1/fs.
phase = np.cumsum(freq) / SAMPLE_RATE

# convert phase from cycles to radians
phase *= 2 * np.pi

# generate waveform
signal = np.sin(phase)

Plotting:

plt.figure(figsize=(12, 4))
plt.plot(t, signal * 20, label='Signal') # amplification for visibility
plt.plot(t, freq, label='Frequency')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.show()

Sanity check:

• I asked for 100 Hz with some wobble, 0.2 seconds, so that should be 20 cycles. Looks like it.

• The wobble wobbles at 20 Hz. Over 0.2s, that is four wobbles. Looks like it.

• Wobble spans 120 to 80 Hz, so the space between peaks should be something like 3:2. Looks like it.