% elemzes_furulya.m
% Analyze the sound produced by a recorder

clear;

%% Load the input
[x_in, fs] = audioread('recorder_C5_110Pa.wav');

%% Plot the loaded data
figure;
n_samp = size(x_in, 1);
t = (0 : n_samp-1).' / fs;
plot(t, x_in);
xlabel('Time [s]');
ylabel('Signal [-]');

%% Clean the signal using a low-pass filter
f_cut = 50;     % Cut-off frequency of the filter
[b, a] = butter(3, f_cut/(fs/2), 'high');
x_filt = filter(b, a, x_in);

figure;
plot(t, x_in, t, x_filt);

%% Cut the signal
t_lims = [5 10];         % One second for stationary analysis
sel = t >= t_lims(1) & t < t_lims(2);

t_sel = t(sel);
t_sel = t_sel - t_sel(1);
x_sel = x_filt(sel);

figure;
plot(t_sel, x_sel);

%% Listen
soundsc(x_sel, fs);

%% Calculate the spectrum
X = fft(x_sel);
n_sel = length(x_sel);
f = (0 : n_sel-1) / n_sel * fs;

plot(f, 20*log10(abs(X)));
xlim([0 10000]);
xlabel('Frequency [Hz]');
ylabel('Amplitude [dB]');

%% Find the fundamental frequency
[~, i_max] = max(abs(X));

f_est = 530;       % Estimated frequency
fprintf('Estimated frequency: %.2f Hz\n', f_est);

% Complex modulation signal
x_mod = exp(-2i*pi*f_est*t_sel);
% Transform signal by modulation
x_trans = x_sel .* x_mod;

% Create a low-pass filter
[b, a] = butter(3, 20/(fs/2), 'low');
x_lpf = filter(b, a, x_trans);

% Plot the result of modulation
figure;
plot(t_sel, [real(x_lpf), imag(x_lpf), abs(x_lpf)]);
xlabel('Time [s]');
ylabel('Modulated and filtered signal [-]');
legend({'Real', 'Imag', 'Abs'});

% Take the second half of the signal
n_half = ceil(n_sel/2);

lhs = x_lpf(n_half : end-1);
rhs = x_lpf(n_half+1 : end);

% Find gamma = exp(1i*delta_f*Ts)
gamma = lhs \ rhs;
% Compute delta_f from gamma
delta_f = angle(gamma) / (2*pi) * fs;
% Find the correctred frequency 
f1 = f_est + delta_f;
fprintf('Corrected estimation: %.2f Hz\n', f1);

%% Resample the original signal
ratio = 64;                     % Ratio for resampling
fs_re = f1 * ratio;             % The new sampling frequency

[p, q] = rat(fs_re / fs);
x_res = resample(x_sel, p, q);
% Take a window length such that coherent sampling is ensured
n_res = floor(length(x_res) / ratio)*ratio;
x_res = x_res(1 : n_res);

w = hann(n_res);
X_res = fft(x_res .* w);
f_res = (0 : n_res-1).' / n_res * fs_re;

% Plot the spectrum of the resampled signal
figure;
plot(f_res, 20*log10(abs(X_res)));
n_harm = 10;
harm_idx = (1 : n_harm)*(n_res/ratio)+1;
hold on;
plot(f_res(harm_idx), 20*log10(abs(X_res(harm_idx))));
xlim([0 10000]);
xlabel('Frequency [Hz]');
ylabel('Amplitude [dB]');
legend({'Spectrum', 'Spectral envelope'});
%% Synthesize
T_sim = 5;
ampls = abs(X_res(harm_idx));
t = (0 : T_sim*fs-1).'/fs;
freqs = f1*(1 : n_harm);

x_synth = sin(2*pi*t*freqs) * ampls;

t_ADSR = [0 .1 .2 T_sim-.1 T_sim];
A_ADSR = [0 1.1 1.0 1.0 0];
a_ADSR = interp1(t_ADSR, A_ADSR, t, 'linear');

x_adsr = x_synth .* a_ADSR;

%% Compare
t_cmp = 5+[0 T_sim];
t = (0 : n_samp-1).' / fs;
sel = t >= t_cmp(1) & t < t_cmp(2);
x_cmp = x_in(sel);
