% elemzes_lemez.m
% Analysis of a plate hit by an impulse hammer

%% Load the data

clear;
load('lemez_meres.mat');

%% Plot the data
n_samp = size(data,1);              % Number of samples
t = (0 : n_samp-1).' / fs;          % Time axis [s]
figure;
% plot the data 
plot(t, data);
xlabel('Time [s]');
ylabel('Signals');
legend(channels);

%% Select a time limit to analyze
[t_lims, ~] = ginput;
%% Make selection
t_sel = t > t_lims(1) & t < t_lims(2);
t_hit = t(t_sel);
x_hit = data(t_sel, :);

%% Compute spectrum
n_hit = size(x_hit, 1);
f_hit = (0 : n_hit-1).' / n_hit * fs;
X_hit = fft(x_hit);
%% Create spectrum figure
figure;
ax1 = subplot(2,1,1);
plot(f_hit, 20*log10(abs(X_hit(:,2))));
ax2 = subplot(2,1,2);
plot(f_hit, 20*log10(abs(X_hit(:,3))));
linkaxes([ax1 ax2], 'x');
xlim([0 10000]);
xlabel(ax2, 'Frequency [Hz]');
ylabel(ax1, 'Acceleration');
ylabel(ax2, 'Sound pressure');
%% Create a high-pass filter
[b, a] = butter(3, 50/(fs/2), 'high');
x_filt = filter(b, a, x_hit);

%% Resample to allow narrowband filtering
fs_re = 10000;
[p, q] = rat(fs_re/fs);
x_re = resample(x_filt, p, q);      % Perform resampling
n_re = size(x_re, 1);               % Number of samples

%% Synthesize
f_sel = [270 486 873 1268];
n_sel = length(f_sel);
filt_order = 5;
filt_width = 20;

x_synth = zeros(size(x_re,1),1);
for i = 1 : n_sel
    f_cent = f_sel(i);
    [b, a] = butter(filt_order, (f_cent+filt_width*[-1 1])/(fs_re/2));
    x_re_filt = filter(b, a, x_re(:,3));
    x_synth = x_synth + x_re_filt;
end
%% Create a figure
figure;
t_re = (0 : n_re-1).'/fs_re;
plot(t_re, [x_re(:,3), x_synth]);
xlabel('Time [s]');
ylabel('Sound pressure signals [Pa]');
