% Piano Sustain Effect Script
% Method 2: IR

% Parameters
MIDI = 60;
METHOD = 2;
start = 4437 + 10000;    % Start after attack
window_size = 3000;      % Analysis window size
NFFT = 4096;            % FFT size
resonant_bins = 1000;     % Number of frequency bins to analyze
scaling_factor = 0.2;    % Reconstruction scaling

% Load audio files
noPedal_path = sprintf('../data/no_pedal_%d.wav', MIDI);
pedal_path = sprintf('../data/long_pedal_%d.wav', MIDI);
[noPedal, Fs1] = audioread(noPedal_path);
[pedal, Fs2] = audioread(pedal_path);

% Convert to mono if stereo
if size(noPedal, 2) > 1
    noPedal = mean(noPedal, 2);
end
if size(pedal, 2) > 1
    pedal = mean(pedal, 2);
end

% Get analysis chunks
pedal_chunk = pedal(start:start+window_size);
noPedal_chunk = noPedal(start:start+window_size);

% Calculate spectra
pedal_spec = abs(fft(pedal_chunk, NFFT));
noPedal_spec = abs(fft(noPedal_chunk, NFFT));

% Create frequency difference profile
freq_diff = pedal_spec(1:resonant_bins) - noPedal_spec(1:resonant_bins);
freq_diff = freq_diff(:);
freq_diff = max(freq_diff, 0);  % Keep only positive differences

% Build impulse response
ir_freq = zeros(NFFT, 1);
ir_freq(1:resonant_bins) = freq_diff;
ir_freq(NFFT:-1:NFFT-resonant_bins+2) = conj(freq_diff(2:end));
ir = real(ifft(ir_freq));

% Apply decay window
window = exp(-(0:length(ir)-1)'/1000);
ir = ir .* window;

% Test reconstruction on short segment
test_len = 1000;
noPedal_test = noPedal(start:start+test_len-1);
pedal_test = pedal(start:start+test_len-1);

filtered = conv(noPedal_test, ir, 'same');
reconstructed = noPedal_test + scaling_factor * filtered;

% Plot test reconstruction
figure;
plot(pedal_test, 'b');
hold on;
plot(reconstructed, 'r');
title('Original vs Reconstructed Waveform from Deconvolution)');
legend('Original', 'Reconstructed');
xlabel('Sample');
ylabel('Amplitude');
grid on;

% Plot impulse response
figure;
plot(ir);
title('Impulse Response from Deconvolution');
xlabel('Sample');
ylabel('Amplitude');
grid on;

% Reconstruction
filtered = conv(noPedal, ir);
reconstructed_full = noPedal + scaling_factor * filtered(1:length(noPedal));

% Normalize to prevent clipping
reconstructed_full = reconstructed_full / max(abs(reconstructed_full));

% Save and play
sound(reconstructed_full, Fs1);

% Save result
% synth_path = sprintf('../result/%d_%d.wav', METHOD, MIDI);
% audiowrite(synth_path, reconstructed_full, Fs1);
% disp('Processed audio saved');