Generate sine wave/tone Matlab

Download
genTone.m


How-to
Generating tones or sine waves of specific frequencies in Matlab is conceptually simple, but generating the correct frequency over the correct time period at the correct sampling rate can be a bit confusing at first. This article deals mainly to generating sine wave tones, and genTone.m is designed to simplify this process, but the timing aspects also apply to any signal output from a digital system.

In order to get the output frequency of the tone correct, the sampling rate of the (digital) output device must be known; 44100 Hz is a typical value for a basic sound card (for clarity of plotting the examples here use lower sampling rates). Remember that the minimum sampling rate (Hz, samples/second) must be at least double the highest frequency in the output tone (Hz, or  wave cycles/second, see: Nyquist).

For timing, we need to know 3 things:

• Desired duration of the tone (in seconds)
• Desired frequency of the tone (in Hz)
• Sampling rate of the output device (in Hz)

From this information we can calculate how many points we need to generate a tone of a certain frequency and make it last for the correct amount of time. For example, for a 2 second tone at 1000 Hz sampling rate needs 2*1000 = 2000 "points" or "samples" to last for 2 seconds when output at 1000 samples per second.

The first thing to do is generate a time vector with the correct number of points. This can be done in a couple of ways, for example

Fs = 1000 % Sampling rate, Hz

dur = 1 % Duration, S

freq = 10 % Frequency, Hz

amp = 1 % Amplitude, for example, Volts
phase = 0 % Start point of wave, rads

timeVec = 0 : 1/Fs : dur - (1/Fs)

This command generates a vector from 0 (:) in steps of 1/Fs (:) up to duration (minus 1 step, as it's starting from 0). Here step is the size of the time step so that we end up with dur*Fs points in the vector. 

An alternative command would be

timeVec = linspace(0, dur-1/Fs, dur*Fs)

Linspace creates a vector from input 1 to input 2, with the number of steps specified in input 3. Again starting from 0, so end point is the duration minus one time step.

Next step is generating the sine wave. Note that Matlab's sine function expects radians as an input. This means that:

sin(-pi) = 0

sin(-pi/2) = -1

sin(0) = 0

sin(pi/2) = +1

sin(pi) =  0

We want a wave that runs from 0 to -1 to 0 to 1 to 0, so the range we need to use is -pi to pi - that's a total range of 2*pi. This pattern continues for values above pi, so it doesn't actually matter where we start, just so long as whatever timeVec we have runs through one full cycle every second.

tone1HzRads = timeVec * 2 * pi

Additionally, the pi range used can be shifted to shift the phase of the wave (in radians) simply by adding it:

tone1HzRads = (timeVec * 2 * pi) + phase

"Phase" is in radians and refers to the point in the cycle the wave is at at any given time, adjusting it effectively shifts the start point.

tone1HzRads now contains the values we need to take the sine of the generate one cycle every second with the sampling rate and duration specified (and starting with the phase specified). ie: 

tone1Hz = sin(tone1HzRads)

To manipulate amplitude, multiply outside the sine, to manipulate frequency multiply within the sine.

tone = amp * sin(freq * tone1HzRads)

Or in other "words"

tone = amp * sin(2 * pi * freq .* timeVec + phase)

If the output is intended to be played from a speaker, there's another technical consideration. It's generally advisable to add a brief rise and fall envelope (20 ms usually works well) so the amplitude at the start and end doesn't change so suddenly. This stops the speaker from clicking at the onset and offset.

An easy way to do this is to create an envelope for the wave. This will be a vector of 0 to 1 of the same length as the wave. Multiplying the wave by the envelope will attenuate the wave's amplitude anywhere where the envelope is <1. To add a rise and fall, we'll create an envelope that starts at 0, rises to 1, remains at 1 for the majority of the wave, then falls to 0 near the end.

For the rise and fall regions, a quarter of a cosine cycle will be added to the envelope. For the rise we want a rising quarter. For the fall we want a falling quarter. The rest of the points we don't want to change. 

First we need to calculate how many points we need - this depends on the length of the rise/fall and the sampling rate:

riseDur = 20/1000; % Rise duration (20 ms) in S
nPtsRise = round(riseDur * Fs); % Corresponding number of points

Then generate a cosine wave over the range -pi/2 to 0, in steps of the length of this range/the number of points (minus 1 again):

tVecRise = -pi/2 : abs(-pi/2-0)/(nPtsRise-1) : 0; 

Or, in linspace terms, the equivalent would be:

tVecRise = linspace(-pi/2, 0, nPtsRise)

The rise can then be generated in a similar way as the sine wave previously, except this time we don't need to adjust the phase, amplitude, or frequency. The fall is simply the reverse of the rise, the end:-1:1 index in to rise returns rise starting with the last point and ending at the first, reversing it.

rise = cos(tVecRise);
fall = rise(end:-1:1);

The envelope itself will be mostly 1s, and will be the same length as the tone

nPts = length(tone);
env = ones(1, nPts);

The it's simply a case of dropping the rise and fall in to the correct parts of the envelope:

env(1:nPtsRise) = rise;
env(end-(nPtsRise-1):end) = fall;

The end-(nPtsRise-1):end index on the second line selects the end of the envelope. For example, if the rise was 20 points long (nPtsRise = 20) and the envelope was 100 points long, then the end index in to env would refer to point 100. ie. 100 - (20-1): 100, or 81:100. Note the minus 1, this is because 80:100 would be 21 points

Finally the envelope is applied to the tone by doing an element-wise multiplication.

tone = env .* tone;

In the examples below the tones generated will be plotted visually. The way to output a signal through speakers depends on the hardware and software setup, but the simplest way is to use Matlab's audio functions - see wavplay in older versions of Matlab and audioplayer in more recent versions. If reliable or short latency is required, external toolboxes such as Psychtoolbox are more useful.

Try this:

p = audioplayer(tone, Fs);
play(p)

genTone.m examples

genTone is designed to automate the steps above. It takes 6 inputs, amp, freq, duration, rise duration, phase and sampling rate. It can return up to 3 outputs, the generated tone, the time vector used to generate the tone and the envelope used to modify the tone.

%% Simple sine wave

amp = 2;
freq = 10;
dur = 1;
riseDur = 0;
Fs = 200;
phase = 0.5;
[tone, tVec] = genTone(amp, freq, dur, riseDur, phase, Fs);

clf
plot(tVec, tone)
ylabel('Amp, V')
xlabel('Time, S')
hold on
scatter(tVec, tone)

%% Sine wave with rise/fall
amp = 2;
freq = 8;
dur = 1;
riseDur = 0.5;
Fs = 300;
phase = 0
[tone, tVec, env] = genTone(amp, freq, dur, riseDur, phase, Fs);

clf
plot(tVec, tone)
hold on
scatter(tVec, tone)
plot(tVec, amp*env)
ylabel('Amp, V')
xlabel('Time, S')

%% Tone at realistic sampling rate
amp = 1;
freq = 800;
dur = 1;
riseDur = 0.02;
Fs = 44100;
phase = 0;
[tone, tVec, env] = genTone(amp, freq, dur, riseDur, phase, Fs);

clf
plot(tVec, tone)
hold on
plot(tVec, amp*env)
ylabel('Amp, V')
xlabel('Time, S')

%% Multiple waves
amp = [1, 3, 3];
freq = [4, 7, 15];
dur = 0.5;
riseDur = [0.02, 0.03, 0.04];
Fs = 500;
phase = [-1, 0.2, 1.6];

nTones = 3;
tones = NaN(3, round(dur*Fs));
for n = 1:nTones
    [tones(n,:), tVec] = genTone(amp(n), freq(n), dur, riseDur(n), phase(n), Fs);
end
comboTone = sum(tones)./nTones;

clf
plot(tVec, tones')
hold on
scatter(tVec(1:2:end), comboTone(1:2:end))
ylabel('Amp, V')
xlabel('Time, S')

%% Breaking Nyquist
figure
for f = 300:-1:2
    amp = 2;
    freq = 8;
    dur = 1;
    riseDur = 0.25;
    Fs = f; % Decrease sampling rate
    phase = 0;
    [tone, tVec, env] = genTone(amp, freq, dur, riseDur, phase, Fs);
    
    clf
    plot(tVec, tone)
    title(['Sampling rate = ', num2str(f), 'Hz'])
    ylabel('Amp, V')
    xlabel('Time, S')
    hold on
    scatter(tVec, tone)
    drawnow
    
end

genTone.m Code

function [tone, tVec, env] = genTone(amp, freq, dur, riseDur, phase, Fs)

% amp = amplitude, arb. (eg. Volts)

% freq = frequency, Hz

% dur = duration, S

% riseDur = rise/fall time, in S

% Fs = Sampling frequency, Hz

% Calculate number of points required

dt = 1/Fs;

% Length in pts

% npts (pts) = dur (S) * Fs (pts/S)

nPts = round(dur * Fs);

% Generate time vector

tVec = 0:dt:dur-dt;

% Generate tone

tone = amp * sin(pi*2*freq.*tVec + phase);

% Generate rise/fall

nPtsRise = round(riseDur * Fs);

tVecRise = -pi/2 : abs(-pi/2-0)/(nPtsRise-1) : 0;

% Equiv to linspace(-pi/2, 0, nPtsRise)

% Generate rise

rise = cos(tVecRise);

fall = rise(end:-1:1);

% Add to envelope

env = ones(1, nPts);

env(1:nPtsRise) = rise;

env(end-(nPtsRise-1):end) = fall;

% Apply envelope to tone

tone = env .* tone;