% DFT and CQT computer exercise

clear                 % clear all variables
close all             % close all figures

% intialize parameters and signals
f_1=2^-2*440;         % frequency of first tone
f_2=2^-(5/3)*440;     % frequency of second tone
fs=44.1E3;            % sampling frequency

N_1=fs/f_1;           % number of samples for one period of signal 1
k=0:1*fs;             % generate sampling index vector for 1 second   
y_1=sin(k*2*pi/N_1);  % generate audio signal 1

N_2=fs/f_2;           % number of samples for one period of signal 2
y_2=sin(k*2*pi/N_2);  % generate audio signal 2

y = y_1;              % input signal to analyse

% [y,fs] = audioread('audiofile_ex4.4.wav'), y = y.';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute DFT spectrum of a short segment
M=1024;                     % DFT length = number of frequency bins
w = hanning(M)';            % generate window function
%!!!!!! X_DFT =             % compute DFT/FFT

% generate vector of analysis frequencies
freqvec = 1:M;        %!!!!!

figure(1)                                 % open figure (1)
subplot(1,2,1)                            % create sub-plot
plot(freqvec, 20*log10(abs(X_DFT)));      % plot DFT spectrum
xlabel('Frequency (Hz)');
ylabel('dB')
title('DFT Spectrum');

return                      % remove this command to continue into CQT

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute CQT spectrum of a short segment
b = 1;              % number of CQT bins per semi-tone
f_min = 55;         % minimum analysis frequency for CQT
f_max = fs/2;       % compute Nyquist frequency

% compute number of frequency bins, M is the highest mu
M = floor(12*b*log2(f_max/f_min))+1;

% compute Q factor
Q = 1./(2^(1/(12*b))-1);

% compute center frequencies of CQT
f_mu=f_min.*2.^((0:M-1)/(12*b));

% compute lengths of analysis windows
N_mu=Q*fs./f_mu; N_mu=round(N_mu);

% compute shift of segments = sub-sampling factor
R=floor(min(N_mu)/2);

% number of signal segements to process
numframes = 1;                                      % just one frame
% numframes = ceil((length(y) - max(N_mu))./R);     % all frames

X_CQT = zeros(numframes,M);

for mu = 1:M
  % for all frequency bins mu
  % compute window with approriate length
  % !!!! w = 
  % generate index vector
  % !!!! K =  
  % precompute segment-independent parts
  W = w.* exp(-1i*(2*pi*Q*(K.'-1))/N_mu(mu));
  % compute CQT sum
  for lambda = 1:numframes 
    X_CQT(lambda,mu) = sum(y(K+(lambda-1)*R).*W.')/N_mu(mu);
  end
end

figure(1)                 % open figure
subplot(1,2,2)            % create sub-plot
freqvec_log = 1:M;        % !!!! create frequency axis
plot(freqvec_log,20 .* log10(abs(X_CQT)));   % plot CQT spectrum
% axis([0 M-1 -120 0])
title('CQT Spectrum');    
xlabel('Frequency (Hz)');
ylabel('dB')


figure(10)
x_axis=0:(length(y)-max(N_mu))/R;
x_axis=x_axis*R/fs;
mu=0:M-1;
imagesc(x_axis, freqvec_log, db(X_CQT'),[-100 0])
xlabel ('t / s')




