%% Digital Signal Processing and Digital Audio
%
% Table of Contents
% 1 Introduction
%     Digital Signal Processing (DSP)
%     Digital Audio Applications
% 2. What is Sound?
% 3. Digitising Sound: Digtial Sampling
%     Computer Manipulation of Sound
% 4.Digital Sampling
%     Sample Rates and Bit Size
%     Digital Sampling
%     Nyquist's Sampling Theorem
%     If you get Nyquist Sampling Wrong?
%     What does aliasing sound like?
%         Sine Wave
%         Piano
%     Implications of Sample Rate and Bit Size
%     Practical Implications of Nyquist Sampling Theory
% 5. Basic Digital Audio Signal Processing
%     Simple Waveforms
%         The Sine Wave and Sound
%         MATLAB Sine Wave Frequency
%         MATLAB Sine Wave Plot of  cycles
%         MATLAB Square and Sawtooth Waveforms
%         Recap: Relationship Between Amplitude, Frequency and Phase
%         Phase of a Sine Wave
%     Digital Audio Effects Example
% 6. Basic DSP Concepts and Definitions
%         The Decibel (dB)
%         Why Use Decibel Scales?
%         Decibel and Chillies!
%         Decibel and acoustics
%         Examples of dB measurement in Sound
%         Signal to Noise Ratio
% Supporting functions
%
%% 1 Introduction
% Issues to be covered (Over next few lectures):
% Digital Signal Processing and Digital Audio
% Sampling Theorem
% Digital Audio Signal Processing
% Digital Audio Effects
%% Digital Signal Processing (DSP)
% DSP includes many different topics, such as:
% Digital filters
% Analysis of signals and systems (especially in terms of frequency)
% Synthesis of signals
% Detection of signals and estimation of signal and system parameters
% Data compression
% and on and on ...
% Related Topics DSP is the intersection of a number of different areas of study:
% Mathematics
% Electrical engineering
% Signals and systems
% Analog circuit theory
% Computer architecture, (and more)
% Probability and statistics
% Computer programming
% Strong Link to Image and Video Processing — we will study some introductory aspects of theses topics in a few weeks time.
%% Digital Audio Applications
% Application of Digital Audio — A few selected examples:
% Music Production
% Hard Disk Recording
% Sound Synthesis
% Samplers
% Effects Processing
% Video
% Audio Important Element: Music and Effects
% Web
% Many uses on Web:
% Streaming Audio
% Spotify
% Listen to Web Radio
% Element of a Web Page
%% 2. What is Sound?
%
% Sound is transmitted via air pressure waves: some form of oscillation needs to move air to generate sound.
% How is sound generated:
% Sound Sources— A few ways to generate sound:
% Air Pressure changes
% Electrical — Loud Speaker
% Acoustic — Direct Pressure Variations
%
% Sound Reception: Air pressure waves revieved and need to translate/interpret these waves as sound.
% Sound Destination— Receives Sound. Examples:
% Electrical — Microphone produces electric signal
% Ears — Responds to pressure hear sound (MPEG Audio — exploits this fact)
%% 3. Digitising Sound: Digtial Sampling
%
%
% Microphone:
% Receives sound
% Converts to  analog signal.
% Computer like discrete entities
% Need to convert Analog-to-Digital — Dedicated Hardware (e.g. Soundcard)
% Also known as Digital Sampling
%% Computer Manipulation of Sound
% Digital Audio Examples Digital Signal Processing routines range from being trivial to highly complex :
% Volume
% Fading in/out, Cross-fading
% Reverse
% Repeating/Looping
% Echo/Reverb/Delay
% Filtering
% Signal Analysis
% Change Pitch and Tempo
% We'll meet a few examples in due course (and lot more in CM3106 Multimedia, if studied next year).
%% 4.Digital Sampling
%% Sample Rates and Bit Size
% -------------------------------------------------------
% Bit Size —  How do we store each sample value (Quantisation)?
% 8 Bit Value(0-255)
% 16 Bit Value(Integer) (0-65535)
% -------------------------------------------------------
% -------------------------------------------------------
% Sample Rate — How many Samples to take?
% 11.025 KHz— Speech (Telephone 8 KHz)
% 22.05 KHz— Low Grade Audio (WWW Audio, AM Radio)
% 44.1 KHz— CD Quality
% -------------------------------------------------------
%% Digital Sampling
% Sampling process basically involves:
% Measuring the analog signal at regular discrete intervals
% Recording the value at these points
%
%
%
%
%% Nyquist's Sampling Theorem
%
% The Sampling Frequency is critical to the accurate reproduction of a digital version of an analog waveform
% -------------------------------------------------------
% Nyquist's Sampling Theorem:
% The Sampling frequency for a signal must be at least twice the highest frequency component in the signal.
% -------------------------------------------------------
% The following graphs are taken and adapted from Nikhil Prathapani (2020). Sampling Theorem Illustration (https://www.mathworks.com/matlabcentral/fileexchange/45109-sampling-theorem-illustration), MATLAB Central File Exchange. Retrieved October 1, 2020)

% "Continuous" time interval
t=0:0.01:20;
T=4;fm=1/T;

% Continuous sinewave
x=cos(2*pi*fm*t);

plot(t,x,'LineStyle',"--",'linewidth',2,"Color",[0.4940, 0.1840, 0.5560]);
xlabel('time');ylabel('amplitude');
grid;
title('Original sinewave');
% Sampling at Signal Frequency
fs1=fm;

% Discrete time interval
n1=0:(1/fs1):20;

% Sampled sinewave
x1=cos(2*pi*fm*n1);

stem(n1,x1,'linewidth',3);
xlabel('number of samples');ylabel('amplitude');
hold on;
plot(n1,x1,'linewidth',3);
plot(t,x,'LineStyle',"--","Color",[0.4940, 0.1840, 0.5560],'linewidth',2);
xlabel('time');ylabel('amplitude');
grid;
title('Sampling at Signal Frequency');
legend('samples','reconstructed wave','original wave')
hold off

%

% Sampling at Twice Nyquist Frequency
fs2=2*fm;

% Discrete time interval
n2=0:(1/fs2):20;

% Sampled sinewave
x2=cos(2*pi*fm*n2);

stem(n2,x2,'linewidth',3);
xlabel('number of samples');ylabel('amplitude');
hold on;
plot(n2,x2,'linewidth',3);
plot(t,x,'LineStyle',"--","Color",[0.4940, 0.1840, 0.5560],'linewidth',2);
xlabel('time');ylabel('amplitude');
grid;
title('Sampling at Twice Nyquist Frequency');
legend('samples','reconstructed wave','original wave')
hold off

%

% Sampling at above Nyquist Frequency
fs3=20*fm;

% Discrete time interval
n3=0:(1/fs3):20;

% Sampled sinewave
x3=cos(2*pi*fm*n3);

stem(n3,x3,'linewidth',3);
hold on;
plot(n3,x3,'linewidth',3);
plot(t,x,'LineStyle',"--","Color",[0.4940, 0.1840, 0.5560],'linewidth',2);
xlabel('number of samples');ylabel('amplitude');
xlabel('time');ylabel('amplitude');
grid;
title('Sampling at above Nyquist Frequency');
legend('samples','reconstructed wave','original wave')
hold off

%
%% If you get Nyquist Sampling Wrong?

% Frequencies below Nyquist
fs4=1.2*fm;
fs5=1.6*fm;
fs6=1.8*fm;

% Discrete time intervals
n4=0:(1/fs4):20;
n5=0:(1/fs5):20;
n6=0:(1/fs6):20;

% Sampled sinewave
x4=interp1(n4,cos(2*pi*fm*n4),t,'pchip');
x5=interp1(n5,cos(2*pi*fm*n5),t,'pchip');
x6=interp1(n6,cos(2*pi*fm*n6),t,'pchip');

figure
xlabel('number of samples');ylabel('amplitude');
hold on;
plot(t,x,'LineStyle',"--","Color",[0.4940, 0.1840, 0.5560],'linewidth',2);
plot(t,x4,'linewidth',3,"Color",[0.6350, 0.0780, 0.1840]);
plot(t,x5,'linewidth',3,"Color",[0.8500, 0.3250, 0.0980]);
plot(t,x6,'linewidth',3,"Color",[0.9290, 0.6940, 0.1250]	);
xlabel('time');ylabel('amplitude');
grid;
hold off;

%
%
% Digital Sampling Artefacts Arise — Effect known as Aliasing
% Aliasing Affects Audio, Imagery and Video
%
%% What does aliasing sound like?
%% Sine Wave

% Downloading the videos, if not alredy in the folder
folderContent=dir;
if ~any(strcmp({folderContent.name},'sin_alias.mov'))
    websave('sin_alias.mov','http://users.cs.cf.ac.uk/Dave.Marshall/Multimedia/PDF/Lecture_Media/sin_alias.mov');
end
if ~any(strcmp({folderContent.name},'sin_no_alias.mov'))
    websave('sin_no_alias.mov','http://users.cs.cf.ac.uk/Dave.Marshall/Multimedia/PDF/Lecture_Media/sin_no_alias.mov');
end

playWithoutAliasing=1;

if playWithoutAliasing
    [y,Fs]=audioread('sin_no_alias.mov');
else
    [y,Fs]=audioread('sin_alias.mov');
end
    
player=audioplayer(y,Fs);
play(player)
% Execute stop(player) to stop the reproduction

%
%% Piano

if ~any(strcmp({folderContent.name},'aliasing_piano.mov'))
    websave('aliasing_piano.mov','http://users.cs.cf.ac.uk/Dave.Marshall/Multimedia/PDF/Lecture_Media/aliasing_piano.mov');
end
[y,Fs]=audioread('aliasing_piano.mov');

player=audioplayer(y,Fs);
play(player)
% Execute stop(player) to stop the reproduction

% MATLAB Code for Sine Demos above: Plot Version , Audio Version
% More on image and video sampling artefacts later.
%
%% Implications of Sample Rate and Bit Size
% Affects Quality of Audio
% Ears do not respond to sound in a linear fashion
% Decibel (dB) a logarithmic measurement of sound
% 16-Bit has a signal-to-noise ratio of 98 dB — virtually inaudible
% 8-bit has a signal-to-noise ratio of 50 dB
% Therefore, 8-bit is roughly 8 times as noisy
% 6 dB increment is twice as loud
% Audio Sample Rate and Bit Size Examples

freq= '44';
if freq=='11'
    bitSample='8';
else
    bitSample= '8';
end
websave([ freq 'k-' bitSample 'bit.aiff'],['http://users.cs.cf.ac.uk/Dave.Marshall/Multimedia/audio/' freq 'k-' bitSample 'bit.aiff']);
[y,Fs]=audioread([freq 'k-' bitSample 'bit.aiff']);
player=audioplayer(y,Fs);
play(player)
% Execute stop(player) to stop the reproduction

% Web Link: Click Here to Hear Sound Examples
% Affects Size of Data: Memory Required for 1 Minute of Digital Audio
%
%
%% Practical Implications of Nyquist Sampling Theory
% Filtering of Signal
% Must (low pass) filter signal before sampling:
%
% Otherwise strange artefacts from high frequency (above Nyquist Limit) signals would appear in the sampled signal.
% -------------------------------------------------------
% Example: Why are CD Sample Rates 44.1 KHz?
% Upper range of human hearing is around 20-22 KHz — Apply Nyquist Theorem
% -------------------------------------------------------
%% 5. Basic Digital Audio Signal Processing
%
% In this section, we look at some basic aspects of Digital Audio Signal Processing:
% Some basic definitions and principles
% Filtering
% Basic Digital Audio Effects
%% Simple Waveforms
% Simple Waveforms
%
% Frequency is the number of cycles per second and is measured in Hertz (Hz)
% Wavelength is  inversely proportional to frequencyi.e. Wavelength varies as \frac{1}{frequency}
%% The Sine Wave and Sound
%
% The general form of the sine wave we shall use (quite a lot of) is as follows: y = A. sin( 2\pi.n.F_w/F_s)
% where:
% Ais the amplitude of the wave,
% F_w is the frequency of the wave,
% F_s is the sample frequency,
% n is the sample index.
% MATLAB function: sin() used — works in radians
% For more information look at:

% MATLAB Doc and Help

showDoc = 0;
if showDoc
    doc sin;
end

showHelp = 1;

if showHelp
    help sin;
end

%
% Example – Simple Basic sin()function and plot, run the following code:
%

amp = 0.5;
freq = 2;
Fs = 256;
n = 0:floor(Fs/freq);

y = amp.*sin(2*pi*n*freq/Fs);

plot(y); % Figure above used stem() instead of plot();
title ('1 Cycle Duration Sine Wave');

% Note some convenience functions are supplied (see code at bottom, also in zip supplied) to plot sin() and other waveforms
%
% As we will often call sin() with a need for specific paremeters (amplitude,freqeuncy,phase,duration and sample_freqeincy) we will also often use a convenient fuction, mysin():

% mysin() example call

amp = 0.5;
freq = 2;
phase = 0; % See below for example of phase
samp_freq = 256;
dur = samp_freq/freq; % Number of Samples so on cycle again here.

y = mysin(amp,freq,phase,dur, samp_freq);

stem(y);

%
% We sill also create other simple waveforms which conveniently create (synthesise) with synth():

% synth() example call
amp = 0.5;
freq = 2;
samp_freq = 256;
dur = 1/freq; % Note: Duration here is in Seconds so 1 cycle is 1 wavelength

y = synth(freq,dur,amp,samp_freq,'sin');

plot(y)

%
%% MATLAB Sine Wave Frequency
%
% To Set the  Frequency of Sine Wave:

F_s = 11025; % Sample Frequency
F_w = 440; %  Frequency of Sine Wave
amp = 0.5;

nsec = 2; 
dur= nsec*F_s;

n = 0:dur;

y = amp*sin(2*pi*n*F_w/F_s);

figure(3) 
plot(y(1:500)); title('Sine Wave of Freqency F_w');

%
%% MATLAB Sine Wave Plot of n cycles
% Plotting of n cycles of a Sine Wave:

F_s = 11025; % Sample Frequency
F_w = 440; %  Frequency of Sine Wave
amp = 0.5;

ncyc = 2;
n=0:floor(ncyc*F_s/F_w);
y = amp*sin(2*pi*n*F_w/F_s);
figure(4) 
plot(y); 
title('N Cycle Duration Sine Wave');

%
%% MATLAB Square and Sawtooth Waveforms
% MATLAB Square and Sawtooth Waveforms MATLAB Square and Sawtooth Waveforms:

F_s = 11025; % Sample Frequency
F_w = 440; %  Frequency of Sine Wave

amp=1;
ncyc = 2;
n=0:floor(ncyc*F_s/F_w);

ysq = amp*square(2*pi*n*F_w/F_s); 

ysaw = amp*sawtooth(2*pi*n*F_w/F_s);

figure(6); 

hold on;

plot(ysq,'b'); 
plot(ysaw,'r'); 
title('Square (Blue)/Sawtooth (Red) Waveform Plots'); 
hold off;

%
%% Recap: Relationship Between Amplitude, Frequency and Phase
%
%% Phase of a Sine Wave
%
% Recall the following code example from CM2208: sinphasedemo.m

% Simple Sin Phase Demo
samp_freq = 400;
dur = 800; % 2 seconds
amp = 1;
phase = 0;
freq = 1;
s1 = mysin(amp,freq,phase,dur,samp_freq);

axisx = (0:dur-1)*360/samp_freq; % x axis in degrees

plot(axisx,s1,'b');
set(gca,'XTick',[0:90:axisx(end)]);


% change phase

phase = 90;

s2 = mysin(amp,freq,phase,dur,samp_freq);
hold on
plot(axisx,s2,'r');
set(gca,'XTick',[0:90:axisx(end)]);
hold off;

%%
%
%% Digital Audio Effects Example
% A very simple example, with MATLAB Code: A basic ampliier:
% height 0.4in mic.pdf
%
%

% Read in some waveform: Could use audioread for a local file.

[u1,FS]=webread('http://users.cs.cf.ac.uk/Dave.Marshall/CM2208/LECTURES/TON2.WAV');
y=interp(u1(1:30),4);

% "Analog" Waveform: y is an emulation/visualisation of a the analog wave
% We use interp to resample the data at a higher rate using lowpass interpolation

subplot(2,4,1);
plot((0:96)*5,y(1:97));
title('x(t)');
axis([0 500 -0.05 0.05]);

% ADC (Analog to digital): Obviously u1 is already a digirised waveform, now we plot u1
xlabel('t in \musec \rightarrow');
subplot(2,4,2);
stem((0:24),u1(1:25),'.');axis([0 25 -0.05 0.05]);
xlabel('n \rightarrow');
title('x(n)');

% Apply some DAFX: In this case, a trivial case of amplification by a
% factor of 0.5 – This is simply multipliction, of course.
subplot(2,4,3);
stem((0:24),0.5*u1(1:25),'.');axis([0 25 -0.05 0.05]);
xlabel('n \rightarrow');
title('y(n)');

%DAC: Visualisation of DAC – we could upsample x(t) but for ease of illustration we
% just scale (amplify) y(t)

subplot(2,4,4);
plot((0:96)*5,0.5*y(1:97));axis([0 500 -0.05 0.05]);
xlabel('t in \musec \rightarrow');
title('y(t)');

%
%
% The code above explained in more detail:
% An analog signal, x(t) with signal amplitude continuous over time, t.
% Following ADC the signal is converted into a a discrete-time and quantised amplitude signal, x(n) — a stream of samples over discrete time index, n
% The time distance between two consecutive samples, the sample interval, T (or sampling period).
% The the sampling frequency is  f_s = \frac{1}{T} — the number of samples per second measured in Hertz (Hz).
% Next we apply some simple DAFX — E.g here we multiply the signal by a factor of 0.5 to produce  y(n) = 0.5.x(n).
% The signal y(n) is then forwarded to the DAC which reconstruct an analog signal y(t)
%
%
%% 6. Basic DSP Concepts and Definitions
%% The Decibel (dB)
% When referring to measurements of power or intensity, we express these in decibels (dB):
% X_{dB} = 10 \log_{10} \left( \frac{X}{X_0} \right)
%
% where:
% X is the actual value of the quantity being measured,
% X_0 is a specified or implied reference level,
% X_{dB} is the quantity expressed in units of decibels, relative to X_0.
% X and X_0 must have the same dimensions — they must measure the same type of quantity in the the same units.
% The reference level itself is always at 0 dB — as shown by setting X = X_0 (note: \log_{10}(1) = 0).
% Matlab functions:
%
% db(x) converts the elements of x to decibels (dB), see:

openExample('signal/DecibelsFromVoltageAndPowerExample')

%
% mag2db(y) expresses in decibels (dB) the magnitude measurements specified in y, see:

openExample('signal/MagnitudeResponseOfAHighpassFilterExample')

%
%% Why Use Decibel Scales?
%
% When there is a large range in frequency or magnitude, logarithm units often used.
% If X is greater than X_0 then X_{dB} is positive (Power Increase)
% If X is less than X_0 then X_{dB} is negative (Power decrease).
% Power Magnitude = |X(i)|^2so (with respect to reference level)
%                 X_{dB} =  10 \log_{10} ( |X(i)|^2 )
%
%  =  20  \log_{10} ( |X(i)| )
%                which is an expression of dB we often come across.
%
% pow2db(y) expresses in decibels (dB) the power measurements specified in y.

% pow2db(y) expresses in decibels (dB) the power measurements specified in y. The relationship between power and decibels is ydb = 10 log10(y).

% from openExample('signal/PowerSpectrumOfANoisySinusoidExample')

% Generate 1024 samples of a noisy sinusoid having a normalized frequency of  2\pi/3 rad/sample. Estimate the power spectrum of the signal using pwelch. Express the estimate in decibels and plot it.

n = 0:1024-1;
x = cos(2*pi*n/3) + randn(size(n));

[pxx,w] = pwelch(x,'power');

dB = pow2db(pxx);

plot(w/pi,dB)
xlabel('\omega / \pi')
ylabel('Power (dB)')

%
%% Decibel and Chillies!
% Decibels are used to express wide dynamic ranges in a many applications (not just in signal processing!):
%
%% Decibel and acoustics
% dB is commonly used to quantify sound levels relative to some 0 dB reference.
% The reference level is typically set at the  threshold of human perception
% Human ear is capable of detecting a very large range of sound pressures --- hence the need for a logarithmic decibel sclae.
%% Examples of dB measurement in Sound
%
% Threshold of Pain:
% The ratio of sound pressure that causes permanent damage from short exposure to the limit that (undamaged) ears can hear is above a million:
% The ratio of the maximum power to the minimum power is above one (short scale) trillion ( 10^{12}).
% The log of a trillion is 12, so this ratio represents a difference of 120 dB.
% 120 dB is the quoted Threshold of Pain for Humans.
% Speech Sensitivity Human ear is not equally sensitive to all the frequencies of sound within the entire spectrum:
% Maximum human sensitivity at noise levels at between 2 and 4 kHz (Speech)
% These are factored more heavily into sound descriptions using a process called frequency weighting.
% Filter (Partition) into frequency bands concentrated in this range.
% Used for Speech Analysis
% Mathematical Modelling of Human Hearing
% Audio Compression (E.g. MPEG Audio)
%          More on this Later in the module
% Digital Noise increases by 6dB per bit:
% In digital audio sample representation (linear pulse-code modulation (PCM)),
% The first bit (least significant bit, or LSB) produces residual quantization noise (bearing little resemblance to the source signal)
% Each subsequent bit offered by the system doubles the  resolution, corresponding to a 6  (= 10*\log_{10}(4)) dB.
% So a 16-bit (linear) audio format offers 15 bits beyond the first, for a dynamic range (between quantization noise and clipping) of (15 x 6) = 90 dB, meaning that the maximum signal is 90 dB above the theoretical peak(s) of quantisation noise.
% 8-bit linear PCM similarly gives (7 x 6) = 42 dB.
% 48 dB difference between 8- and 16-bit which is (48/6 (dB))  --- 8 times as noisy.
%            More on this Later too.
%
%% Signal to Noise Ratio
%
% Signal-to-noise ratio is a term for the power ratio between a signal (meaningful information) and the background noise:
%     SNR = \frac{P_{signal}}{P_{noise}} = \left(\frac{A_{signal}}{A_{noise}}\right)^2
% where P is average power and A is RMS amplitude.
% Both signal and noise power (or amplitude) must be measured at the same or equivalent points in a system, and within the same system bandwidth.
% Because many signals have a very wide dynamic range, SNRs are usually expressed in terms of the logarithmic decibel scale:
%     SNR_{dB} = 10 \log_{10} \left( \frac{P_{signal}}{P_{noise}} \right)= 20 \log_{10} \left(\frac{A_{signal}}{A_{noise}}\right)

% MATLAB has a built snr() function - see doc snr for more examples

showDoc = 0;
if showDoc
    doc snr;
end

%
% Generate a 20-millisecond rectangular pulse sampled for 2 seconds at 10 kHz:

Tpulse = 20e-3;
Fs = 10e3;
t = -1:1/Fs:1;
x = rectpuls(t,Tpulse);

% Embed the pulse in white Gaussian noise such that the signal-to-noise ratio (SNR) is 53 dB. Reset the random number generator for reproducible results
%

rng default

SNR = 53;
y = randn(size(x))*std(x)/db2mag(SNR);

plot(x,'b');

figure

plot(y,'r')
s = x+y;


figure

plot(s,'g')

% Use the snr() function to compute the SNR of the noisy signal.

pulseSNR = snr(s,s-x)

%
% See also psnr():

%Read image and create a copy with added noise. The original image is the reference image.

ref = imread('pout.tif');

imshow(ref);

noiseref = imnoise(ref,'salt & pepper', 0.02);

imshow(noiseref);

% Calculate the PSNR.

[peaksnr, snr] = psnr(noiseref, ref);
  
fprintf('\n The Peak-SNR value is %0.4f', peaksnr);
 % The Peak-SNR value is 22.6437
fprintf('\n The SNR value is %0.4f \n', snr);

%
%% Supporting functions
% mysin():

function s = mysin(amp,freq,phase,dur, samp_freq)
% example function to so show how amplitude,frequency and phase
% are changed in a sin function
% Inputs: amp - amplitude of the wave
%         freq - frequency of the wave
%         phase - phase of the wave in degree
%         dur - duration in number of samples
%         samp_freq - sample frequncy


x = 0:dur-1;
phase = phase*pi/180;

s = amp*sin(2*pi*x*freq/samp_freq + phase);

end

% synth():

function y=synth(freq,dur,amp,Fs,type)
% y=synth(freq,dur,amp,Fs,type)
%
% Synthesize a single note
%
% Inputs:
%  freq - frequency in Hz
%  dur - duration in seconds
%  amp - Amplitude in range [0,1]
%  Fs -  sampling frequency in Hz
%  type - string to select synthesis type
%         current options: 'fm', 'sine', 'sin' or 'saw'
%


    if nargin<5
      error('Five arguments required for synth()');
    end
    
    N = floor(dur*Fs);
    n=0:N-1;
    if (strcmp(type,'sine') || strcmp(type,'sin'))
      y = amp.*sin(2*pi*n*freq/Fs);
    
    elseif (strcmp(type,'saw'))
    
      T = (1/freq)*Fs;     % period in fractional samples
      ramp = (0:(N-1))/T;
      y = ramp-fix(ramp);
      y = amp.*y;
      y = y - mean(y);
    
    elseif (strcmp(type,'fm'))
    
      t = 0:(1/Fs):dur;
      envel = interp1([0 dur/6 dur/3 dur/5 dur], [0 1 .75 .6 0], 0:(1/Fs):dur);
      I_env = 5.*envel;
      y = envel.*sin(2.*pi.*freq.*t + I_env.*sin(2.*pi.*freq.*t));
      
    else
      error('Unknown synthesis type');
    end
    
    % smooth edges w/ 10ms ramp
    if (dur > .02)
      L = 2*fix(.01*Fs)+1;  % L odd
      ramp = bartlett(L)';  % odd length
      L = ceil(L/2);
      y(1:L) = y(1:L) .* ramp(1:L);
      y(end-L+1:end) = y(end-L+1:end) .* ramp(end-L+1:end);
    end
end