% furulya_elemzes.m
% Frekvenciamérés, spektrogram számítás

% Felvétel betöltése
clear;
[y, fs] = audioread('furulya_004.wav');
% Csak a mirkofoncsatornát tartjuk meg
y = y(:, 2);
% Idõtengely
n = size(y, 1);
t = (0 : n-1).' / fs;

% Felüláteresztõ szûrés
[b, a] = butter(3, 50 / (fs/2), 'high');
y = filter(b, a, y);

% Ábrázolás
figure;
plot(t, y);
xlabel('Idõ [s]');
ylabel('Jel [-]');
%% Jelszakasz kiválasztása
t_lim = 5.65 + [1.6 1.9];
sel = t > t_lim(1) & t < t_lim(2);
y_sel = y(sel);
n_sel = size(y_sel, 1);
t_sel = (0 : n_sel-1).' / fs;

figure;
plot(t_sel, y_sel);
xlabel('Idõ [s]');
ylabel('Jel [-]');

%% Frekvencia meghatározása modulációval
% Tippelt frekvencia (pl. FFT-bõl)
f1_hat = 900;
% Egyoldali moduláció
y_m = exp(-2i*pi*f1_hat*t_sel);
y_mod = y_sel .* y_m;

% Aluláteresztõ szûrés
[b, a] = butter(3, 100 / (fs/2));
y_filt = filter(b, a, y_mod);

% Szûrt eredmény
% A szûrt jel komplex értékû, lassan változó jel lesz
figure;
plot(t_sel, [real(y_filt), imag(y_filt), abs(y_filt)]);
xlabel('Idõ [s]');
ylabel('Modulált jel [-]');
legend({'Valós rész', 'Képzetés rész', 'Abszolút érték'});
%% Paraméterillesztés
% A szûrt jel bármely két szomszédos mintájára teljesül (ideális esetben):
%    gamma * y_filt[k] = y_filt[k+1]
% Mivel sok mintánk van, sok egyenletet írhatunk fel az egyetlen gamma
% ismeretlenre. Ezt azután a legkisebb négyzetek módszerével oldjuk meg.
% (Utóbbit a matlabban a "\" mûvelet végzi el automatikusan.)

N0 = 5000;                      % Eltolás, a szûrõ beállási ideje miatt
B = y_filt(N0+1 : end-0);       % Jobb oldali vektor
A = y_filt(N0+0 : end-1);       % Baloldali vektor

% Egyenletrendszer megoldás
gamma = A \ B;
% delta_f a tényleges és a tippelt frekvencia eltérése, mely a gamma
% komplex értékû paraméter fázisából és a mintavételi frekibõl adódik
delta_f = angle(gamma) * fs / (2*pi);

fprintf('Tippelt frekvencia: %.1f Hz\n', f1_hat);
fprintf('Tényleges frekvecia: %.1f Hz\n', f1_hat + delta_f);

%% Spektrogram számítása
% FFT pontszám (a frekvencia és idõfelbontást határozza meg)
n_fft = 4096;               
% Hann ablakfüggvény a spektrális szivárgás csökkentéséhez
win = hann(n_fft);
% Átlapolás - így látszólag jobb idõfelbontást kapunk
n_ol = n_fft * 7 / 8;

% Spektrogram számítása a matlab beépített függvényével
% kimenet: s - spektrogram, w - frekvencia, ts - spektrogram idpvektor
[s, w, ts] = spectrogram(y, win, n_ol, n_fft, fs);
%% Spektrogram ábrázolása
figure;
pcolor(ts, w, 20*log10(abs(s)));
shading interp;
ylim([0 8000]);
ylabel('Frekvencia [Hz]');
xlabel('Idõ [s]');
cb = colorbar;
ylabel(cb, 'Magnitúdó [dB]');
