function wave = oscillator(varargin)
%OSCILLATOR generates a standard waveform, click train, or noise
% OSCILLATOR(wavetype,duration,frequency)
%
% Input arguments:
%
% wavetype (string):
% 'Sinusoid'
% 'Triangle'
% 'Square'
% 'Sawtooth'
% 'Reverse Sawtooth'
% 'Linear Sweep'
% 'Log Sweep'
% 'Click Train'
% 'White Noise'
% 'Pink Noise'
% 'Brown Noise'
% 'Grey Noise'
% 'Speech Noise'
%
% duration (in seconds)
% frequency (in Hz)
% NOTE: linear and log sweeps require [start stop] frequency vector
%
% Optional input arguments:
% OSCILLATOR(wavetype,duration,frequency,gate,phase,sample_freq)
% gate (in seconds): duration of a raised cosine on/off ramp
% phase (in radians): starting phase of the waveform.
% sample_freq (in samples): 44100 is default, custom rates are possible
%
% Examples:
%
% wave = OSCILLATOR('Sinusoid',1,1000); % simple pure tone at 1000 Hz.
% wave = OSCILLATOR('Sawtooth',2,440); % 2 second sawtooth at 440 Hz.
% wave = OSCILLATOR('Pink Noise',1); % 1 second of pink (1/F) noise
% wave = OSCILLATOR('Linear Sweep',2,[440 880]); % linear sweep from 440 to 880 Hz.
% wave = OSCILLATOR('Log Sweep',2,[20 20000],.01); % ramped on/off log sweep.
% wave = OSCILLATOR('White Noise',1,[],0.1); %ramped on and off noise
% wave = OSCILLATOR('Sinusoid',1,220,0,pi/2,48000); %pure tone with a
% starting phase of 90 degrees and sample rate set to 48000.
%
% Omitting 'wavetype' sets it to sinusoid, omitting 'duration' sets it to
% one second, and omitting 'frequency sets it to 440 Hz. Gate is set to 0,
% phase to 0 and sample rate to 44100; All output waves are scaled from -1
% to 1.
%
% Note that the noise signals (pink, brown, grey, etc) are only
% approximations arrived at through filtering, rather than through the more
% accurate iterative processes found elsewhere. Suggestions for better (but
% just as fast) methods or improved coefficients are heartily welcomed.
% (c) W. Owen Brimijoin - MRC Institute of Hearing Research
% Tested on Matlab R2011b and R14
% Version 1.0 18/05/12 - original
% Version 1.1 29/06/12 - added pink and speech-shaped noise options
% Version 1.2 26/06/13 - added swept sine and brown and grey noise options
%input handling:
switch nargin,
case 0, wavetype='Sinusoid';duration=1;frequency=440;gate=0;phase=0;sample_rate=44100;
case 1, wavetype=varargin{1};duration=1;frequency=440;gate=0;phase=0;sample_rate=44100;
case 2, wavetype=varargin{1};duration=varargin{2};frequency=440;gate=0;phase=0;sample_rate=44100;
case 3, wavetype=varargin{1};duration=varargin{2};frequency=varargin{3};gate=0;phase=0;
sample_rate=44100;
case 4, wavetype=varargin{1};duration=varargin{2};frequency=varargin{3};gate=varargin{4};
phase=0;sample_rate = 44100;
case 5, wavetype=varargin{1};duration=varargin{2};frequency=varargin{3};gate=varargin{4};
phase=varargin{5};sample_rate = 44100;
case 6, wavetype=varargin{1};duration=varargin{2};frequency=varargin{3};gate=varargin{4};
phase=varargin{5};sample_rate = varargin{6};
otherwise, error('incorrect number of input arguments');
end
%for noise, frequency is likely omitted, set to default to avoid error:
if isempty(frequency),frequency=440;end
%start input checking:
if ~isstr(wavetype),
error('1st argument ''wavetype'' must be a string.')
end
if ~isnumeric(duration) || numel(duration) ~= 1 || duration < 0 || ~isreal(duration),
error('2nd argument ''duration'' must be a single positive number.')
end
if ~isnumeric(frequency) || ~ismember(numel(frequency),[1 2]) || sum(frequency <= 0) || ~isreal(frequency) || sum(frequency>=sample_rate/2),
error('3rd argument ''frequency'' must be 1 or 2 (for sweeps) positive numbers less than or equal to the Nyquist.')
end
if ~isnumeric(gate) || numel(gate) ~= 1 || gate < 0 || ~isreal(gate) || 2*gate>duration,
error('optional 4th argument ''gate'' must be a positive number less than or equal to half the duration.')
end
if ~isnumeric(phase) || numel(phase) ~= 1 || ~isreal(phase),
error('optional 5th argument ''phase'' should be a single number between -2pi and 2pi.')
end
if ~isnumeric(sample_rate) || numel(sample_rate) ~= 1 || sample_rate < 0 || ~isreal(sample_rate),
error('optional 6th argument ''sample_rate'' must be a single positive number.')
end
if numel(frequency)>1,
if ~any(strcmpi(wavetype,{'log sweep','linear sweep'})),
error('For signals apart from sweeps, ''frequency'' must be 1 positive number.')
end
else
if any(strcmpi(wavetype,{'log sweep','linear sweep'})),
error('For swept sine signals, ''frequency'' must be 2 positive numbers.')
end
end
%end input checking
num_samples = floor(sample_rate*duration);%duration in samples
phase = mod(phase,2*pi)/(2*pi); %phase rescaled from 2*pi to 1
switch lower(wavetype), %generate the chosen waveform:
case 'sinusoid',
% use in-built sine function, offset by the phase argument:
wave = sin((2*pi*phase)+(2*pi*frequency*linspace(0,duration,num_samples)))';
case 'triangle',
% use modulo to fold a line from 0 to frequency, sign-reverse positive values and rescale:
wave = 2*mod(phase+.25+linspace(0,duration*frequency,num_samples)',1)-1;
wave(wave>0) = -wave(wave>0);wave = 2*(wave+0.5);
case 'square',
% same as sinusoid...
wave = sin((2*pi*phase)+(2*pi*frequency*linspace(0,duration,num_samples)))';
% all positive values set to 1 and all negative values to -1:
wave(wave>=0) = 1;wave(wave<0)=-1;
case 'sawtooth',
% simple modulo with phase added linearly:
wave = 2*mod(phase+.25+linspace(0,duration*frequency,num_samples)',1)-1;
case 'reverse sawtooth',
% same as sawtooth...
wave = 2*mod(phase+.25+linspace(0,duration*frequency,num_samples)',1)-1;
% but sign-reversed:
wave = -wave;
case 'linear sweep',
% based on a linearly increasing frequency vector, generate swept sine:
t = linspace(0,duration,num_samples);
wave = sin(2*pi*((frequency(2)-frequency(1)).*(duration.^-1)./2.*(t.^2) + frequency(1).*t + phase))';
case 'log sweep',
% based on the natural logarithm of linearly increasing time and frequency vectors, generate swept sine:
t = linspace(0,duration,num_samples);
freq = duration/log(frequency(2)/frequency(1))*(frequency(1)*(frequency(2)/frequency(1)).^(t/duration)-frequency(1));
wave = sin(2*pi * (freq + phase))';
case {'click train','click'},
% change phase argument to a fraction of the signal period:
phase_offset = mod(round(-phase/frequency*sample_rate),round(1/frequency*sample_rate));
% use the index 'phase offset' to specify location of 1 values:
wave(phase_offset+round([1:1/frequency*sample_rate:num_samples])) = 1;
% ensure that the signal is the correct duration:
wave(1+num_samples)=0; wave = wave(1:num_samples)';
case {'white noise','white','noise'},
% use rand to create noise, wave is then normalized to a max of 1:
wave = 2*(rand(num_samples,1)-0.5);wave = wave./max(abs(wave));
case {'pink noise','pink'},
% 1/f filter coefficients computed using prony's method:
B = [1 -0.69935907 -1.3088024 0.89838254 0.31198493 -0.19603952 0.0017025764];
A = [1 -1.0701074 -1.1525308 1.3802168 0.12319231 -0.31148431 0.031216025];
% filter white noise and normalize:
wave = filter(B,A,2*(rand(num_samples,1)-0.5));wave = wave./max(abs(wave));
case {'brown noise','brown'},
% 1/f^2 filter coefficients computed using prony's method:
B = [1 -0.63757782 -1.2342674 0.75352979 0.25330774 -0.1247317 0.00063048409];
A = [1 -1.4817211 -0.75304159 1.7951497 -0.3000657 -0.33985127 0.079550214];
% filter white noise and normalize:
wave = filter(B,A,2*(rand(num_samples,1)-0.5));wave = wave./max(abs(wave));
case {'grey noise','grey'},
% generic grey noise filter coefficients based on the ISO 60-phon
% Equal-loudness contour and computed using prony's method:
B = [1 -0.8248 0.0617 0.3488 -0.5826 0.5134 -0.2323];
A = [1 -0.8546 0.1445 0.2604 -0.7555 0.4847 -0.2747];
% filter white noise and normalize:
wave = filter(B,A,2*(rand(num_samples,1)-0.5));wave = wave./max(abs(wave));
case {'speech noise','speech'},
% Bernard Seeber's (MRC Institute of Hearing Research) filter
% coefficients. Note that these coefficients were generated based
% on the CCITT standard G.227 telephone signal. This is similar to
% the long term spectrum of speech.
B = [0.00396790391508 0.00032556793042 -0.00314367152058 -0.00104604251859 -0.0000887591994];
A = [1 -3.39268359295324 4.31295903323020 -2.43473845585969 0.51493759484342];
% filter white noise and normalize:
wave = filter(B,A,2*(rand(num_samples,1)-0.5));wave = wave./max(abs(wave));
otherwise,
display('Unrecognized waveform type');wave = [];return
end
if gate>0,
% create the raised cosine ramp:
t = linspace(0,pi/2,floor(sample_rate*gate))';
gate = (cos(t)).^2;
% ramp the beginning of the signal on:
wave(1:length(gate)) = wave(1:length(gate)).*flipud(gate);
% ramp the end of the signal off:
wave(end-length(gate)+1:end) = wave(end-length(gate)+1:end).*gate;
end
