Wave form function generator in Matlab, it also plays the wave sound.

oscillator.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Example:
%%%
%%% Create a 3000Hz triangel wave data for 3 seconds, and
%%% play the wave data at 44100Hz (default) sample rate, also
%%% save to .wav file.
%%% 
%%%
%%%      wave = OSCILLATOR('Triangle', 3, 3000);
%%%      soundsc(wave, 44100);
%%%      wavwrite(wave, 44100, 'noise.wav');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



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'
%    'Blue Noise'
%    'Violet 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_rate (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 to a
% peak absolute value of 1.0
%
% Note that the grey noise is *generic* grey and is based the ISO 66-phon
% Equal-loudness contour.

% (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
% Version 1.3 05/07/13 - added blue and violet and changed noise generation
%                        to a slower but more accurate ifft method. 

%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 or 0, set to default to avoid error:
if isempty(frequency) || any(frequency==0), frequency=1;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+.5+linspace(0,duration*frequency,num_samples)',1)-1;
        
    case 'reverse sawtooth',
        % same as sawtooth...
        wave =  2*mod(phase+.5+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_vals = mod(round(-phase/frequency*sample_rate),round(1/frequency*sample_rate));
        % use the index 'phase offset' to specify location of 1 values:
        wave(phase_vals+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'},
        b = -1; %set beta value
        %create frequency vector folded at center:
        freqs = [linspace(0,.5,floor(num_samples/2)),fliplr(linspace(0,.5,ceil(num_samples/2)))]';
        freqs = (freqs.^2).^(b/2);
        freqs(freqs==inf) = 0; %to prevent inf errors
        phase_vals = rand(num_samples,1);
        %now apply an inverse fft:
        wave = real(ifft(sqrt(freqs).*(cos(2*pi*phase_vals)+i*sin(2*pi*phase_vals))));
        wave = wave./max(abs(wave)); %normalize to +/- 1.
        
    case {'brown noise','brown'},
        b = -2; %set beta value
        %create frequency vector folded at center:
        freqs = [linspace(0,.5,floor(num_samples/2)),fliplr(linspace(0,.5,ceil(num_samples/2)))]';
        freqs = (freqs.^2).^(b/2);
        freqs(freqs==inf) = 0; %to prevent inf errors
        phase_vals = rand(num_samples,1);
        %now apply an inverse fft:
        wave = real(ifft(sqrt(freqs).*(cos(2*pi*phase_vals)+i*sin(2*pi*phase_vals))));
        wave = wave./max(abs(wave)); %normalize to +/- 1.
        
  case {'blue noise','blue'},
        b = 1; %set beta value
        %create frequency vector folded at center:
        freqs = [linspace(0,.5,floor(num_samples/2)),fliplr(linspace(0,.5,ceil(num_samples/2)))]';
        freqs = (freqs.^2).^(b/2);
        freqs(freqs==inf) = 0; %to prevent inf errors
        phase_vals = rand(num_samples,1);
        %now apply an inverse fft:
        wave = real(ifft(sqrt(freqs).*(cos(2*pi*phase_vals)+i*sin(2*pi*phase_vals))));
        wave = wave./max(abs(wave)); %normalize to +/- 1.
        
    case {'violet noise','violet'},
        b = 2; %set beta value
        %create frequency vector folded at center:
        freqs = [linspace(0,.5,floor(num_samples/2)),fliplr(linspace(0,.5,ceil(num_samples/2)))]';
        freqs = (freqs.^2).^(b/2);
        freqs(freqs==inf) = 0; %to prevent inf errors
        phase_vals = rand(num_samples,1);
        %now apply an inverse fft:
        wave = real(ifft(sqrt(freqs).*(cos(2*pi*phase_vals)+i*sin(2*pi*phase_vals))));
        wave = wave./max(abs(wave)); %normalize to +/- 1.
        
    case {'grey noise','grey'},
        %values from the ISO 66-phon Equal-loudness contour (adjusted for
        %optimal spline interpolation):
        freqs = [1 5 15 36 75 138 235 376 572 835 1181 1500 2183 2874 3718 ...
            4800 5946 7377 9051 10996 13239 15808 18735 22050].*sample_rate/44100;
        dB_vals = [61 61 56 40 25 17 11 7 5 4 6 8 3 1 1 4 9 14 17 16 10 5 2 1];
        
        %create level vector for use in inverse Fourier transform:
        freq = linspace(0,sample_rate/2*duration/2,floor(num_samples/2));
        spl = spline(freqs,dB_vals,freq); %upsample 
        levels = [spl,fliplr(spl)]; %reflect vector
        levels = 10.^(levels'./10); %change to power
        levels(levels==inf) = 0; %remove infinite values
        phase_vals = rand(length(levels),1); %generate random phase vals
        %now apply an inverse fft:
        wave = real(ifft(sqrt(levels).*(cos(2*pi*phase_vals)+i*sin(2*pi*phase_vals))));
        wave = wave./max(abs(wave));
        
    case {'speech noise','speech'},
        %values from Byrne et al (1994) JASA manuscript (adjusted for
        %optimal spline interpolation):
        freqs = [1 10 50 80 100 125 160 200 250 315 400 500 630 800 1000 1250 ...
            1600 2000 2500 3150 4000 5000 6300 8000 10000 12000 18000 22050].*sample_rate/44100;
        dB_vals = [35 36 38.6 43.5 54.4 57.7 56.8 60.2 60.3 59.0 62.1 62.1 60.5 ...
            56.8 53.7 53.0 52.0 48.7 48.1 46.8 45.6 44.5 44.3 43.7 43.4  43  41 38];
        
        %create level vector for use in inverse Fourier transform:
        freq = linspace(0,sample_rate*duration/2,floor(num_samples/2));
        spl = spline(freqs,dB_vals,freq);
        levels = [spl,fliplr(spl)]; %reflect vector
        levels = 10.^(levels'./10); %change to power
        levels(levels==inf) = 0; %remove infinite values
        phase_vals = rand(length(levels),1); %generate random phase vals
        %now apply an inverse fft:
        wave = real(ifft(sqrt(levels).*(cos(2*pi*phase_vals)+i*sin(2*pi*phase_vals))));
        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