%Compute the response of a vowel; ECE-537 HW05, 
%Oct 23, 2018: Fix formula for f0 (t/Duration) not t)/Duration
%Sept 16, 2006
%Sept 12, 2004
%copyright: Jont B Allen

clear all;close all; clf; hold off
Fm=22050; %cd player maximum frequency
Fmax=10*Fm; %Upsample for computations
Fs=2*Fmax;
T=1/Fs;
Duration=1.5; %maximum duration in seconds
NFT=Duration/T; %Length of data

t=T*(0:NFT-1); %time vector

%Define pitch: linear change from 150 to 350;, 1 Hz variation; small 20 Hz verbrato
f0=150+100*t/Duration + 25*sin(2*pi*t/Duration) + 0.2*sin(2*pi*20*t);
%f0=150+100*t/Duration + 25*sin(2*pi*t/Duration);
%f0=150+100*t/Duration + .2*sin(2*pi*20*t);

%Generate a glottal pulse
ug=max(0,sin(2*pi*f0.*t));
ug=[diff(ug),0];
ip=find(ug>0);
ug(ip)=ug(ip)-ug(ip); %Remove these points
ug=abs(ug).^3.0; %Take power to modify the spectrum
ug=.99*ug/max(abs(ug)); %Renormalize to .99

%ug=max(0,sin(2*pi*f0.*t)) +.5*max(0,sin(250+2*pi*f0.*t)).^4;
%Design a filter to be used as a glottal pulse
%[b,a]=butter(2,500/Fm);
  % ug=sin(2*pi*1000*t); %debug with a 1 kHz tone

d0=0.008/T; %6 ms of waveform
t0=1:10:d0; % 10 ms, downsampled, time index for graphics
ud=ug(t0); td=t(t0); %Decimated wavefrom for graphics

figure(1);
subplot(211); plot(1000*t,f0);
ylabel('f_0(t)'); xlabel('time [ms]')

Ud=fft(ud); %compute the spectrum of the glottal pulse
Ud=Ud/Ud(1);

%Graphics of Glottal spectrum pulse
figure(2);
subplot(211); plot(1000*td,ud);
ylabel('u_g(t)'); xlabel('time [ms]')

%Plot power spectrum
subplot(212);
M=length(Ud);
dF=2*Fm/(M-1); %Freq. increment
f=0:dF:(2*Fm);
m=2:((M-1)/2);

loglog(f(m),abs(Ud(m)));
xlabel('Freq. [kHz]'); ylabel('Spectrum |Ud(f)|');
xlim([100,5000]); ylim([1e-2,1])
hold on;

%Plot average spectrum of speech, from Dunn and White (1939)
HFS=29.067; %HFS from 8.75 dB/oct
IF1=max(find(f<500));     %f(IF1) = 500 Hz
SS=ones(size(f)); %Spectral level 1 volt for f<500 Hz
SS(IF1:end) =(f(IF1)./f(IF1:end)).^(HFS/20);
plot(f,SS,'r--');

sprintf('Length of array %g [points], Duration %g[s]',NFT,max(t))

%END of example code

%write it out:
%Down sample entire waveform
wavwrite(ug(1:10:end),2*Fm,'ug.wav');
%!play -v 5 ug.wav %linux
%soundsc('ug.wav',2*Fm); %windows
system('play -v 5 ug.wav')
