% Demo script from Math 320 Tanya Leise
% Freq for notes from http://www.phy.mtu.edu/~suits/notefreqs.html
%
% Wavelet transform requires special toolbox (jlab)

clear

Fs=8192; % samples per second
toneFreq = [440 493.88 554.37 587.33 659.26];  % Tone frequency, in Hertz
Ntones=numel(toneFreq);
nSeconds = .5;   % Duration of the sound
y = [];
Nsteps=round(nSeconds*Fs);
t=(0:Nsteps-1)'/Fs;
for k=1:Ntones
    y=[y; sin(2*pi*toneFreq(k)*t)]; % audio signal with 5 notes of scale
end
sound(y,Fs)

L=Fs; 
n=numel(y);
t=(0:n-1)'/Fs; 
if mod(n,2)==0      % freq in Hertz
    k=(1/(nSeconds*length(toneFreq)))*[0:(n-2)/2 n/2 -(n-2)/2:-1];
else
    k=(1/(nSeconds*length(toneFreq)))*[0:(n-1)/2 -(n-1)/2:-1];
end
ks=fftshift(k);
S=y;
St=fft(S);  % Fourier transform of audio signal
figure('Position',[42 37 1140 669]);
plot(k,abs(St));xlim([0 1000]);
set(gca,'FontSize',24);xlabel('Frequency in Hertz (cycles per second) ');
ylabel('Magnitude of Fourier transform ')
title('Fourier transform of audio signal')

pause

figure('Position',[42 37 1140 669]); % movie
a=500;  % reducing window size blurs freq 
        % (increased time resolution -> decreased freq res)
        % but wider window may contain multiple freq's, so that
        % could also cause blurring of freq content
Sgt_spec=[];
tslide=t(1:100:end);
for j=1:length(tslide)
    g=exp(-a*(t-tslide(j)).^2); % Gabor
%     g=1.0*(abs(t-tslide(j))<a); % step fn causes lobes due to discont
    Sg=g.*S;
    Sgt=fft(Sg);
    Sgt_spec=[Sgt_spec abs(fftshift(Sgt))];
    subplot(3,1,1), plot(t,S,'k',t,g,'r');xlim([tslide(j)-.1 tslide(j)+.1]);
    title('Fourier transform of windowed signal')
    subplot(3,1,2), plot(t,Sg,'k');ylim([-1 1]);xlim([tslide(j)-.1 tslide(j)+.1]);
    subplot(3,1,3), plot(ks,abs(fftshift(Sgt))/max(abs(Sgt)));xlim([0 1000]);ylim([0 1]);
    drawnow 
    pause(0.02)
end

pause

figure('Position',[42 37 1140 669])  % spectrogram
pcolor(tslide,ks,Sgt_spec)
shading interp 
colormap(hot)
ylim([0 1000])
set(gca,'FontSize',24)
xlabel('Time in seconds ');ylabel('Frequency in Hertz ')

title('Spectrogram of audio signal')

pause

nvoice=64;gamma=3;beta=8;lowfreq=200;highfreq=800;
[fs,tau,qscaleArray] = CalcScaleForCWT(1/highfreq,1/lowfreq,nSeconds*Ntones,Fs,nvoice);
cwt=wavetrans(S,{1,gamma,beta,fs,'bandpass'},'mirror');

figure('Position',[42 37 1140 669]); % wavelet heat map
imagesc(abs(cwt'))
freqs=highfreq:-100:lowfreq+100;
for k=1:numel(freqs)
    ind(k)=find(1./tau<=freqs(k),1,'first');
end
set(gca,'FontSize',24,'ytick',ind,'yticklabel',freqs);colorbar;
xlabel('Time in seconds ');ylabel('Freq in Hertz ')
title('Wavelet transform of audio signal')

% wavelet ridge plot
[ir,jr,xr,fr]=ridgewalk(cwt,fs);
ind=find(~isnan(ir) & abs(xr)>.999);
figure('Position',[42 37 1140 669]);
plot(ir(ind)/Fs,1./tau(jr(ind)),'k.');
set(gca,'FontSize',24)
xlabel('Time in seconds ');ylabel('Freq in Hertz ')
title('Frequency estimate from wavelet transform')
text(.1,670,'(sharper time resolution of','Fontsize',24)
text(.2,640,'jumps than spectrogram)','Fontsize',24)