%% ######################### MDF obrada STFT-om ###########################
% Skripta (Funkcija) kojoj vršimo glavnu STFT analizu MDF parametara sa 
% promjenjivim vrijednostima L i R
% L= [51 101 151 201 251 301 351 401 451 501]
% R=20 i R=L
%
%Ulazne varijable:
%   - Ts,Fs - frekvencija otipkavanja EMG signala
%   - signal - EMG signal
%   - kut - signal kuta pomoæu kojega odreðujemo vremenske odsjeèke
%   kontrakcije
%Izlazne varijable:
%   - STFTpMDF_L - vektor koji sadrži poèetnu frekvenciju i nagib
%   regresijskog pravca uz fiksni parametar R=20
%   - STFTpMDF_RL - vektor koji sadrži poèetnu frekvenciju i nagib
%   regresijskog pravca uz parametar R=L
% Prikaz:
%   - ctrlPrikaza - kontrolni bit koji odreðuje da li za svaku vrijednost 
%   parametara L i R crtamo prikaz dobivenih vrijednosti MDF 
% Pozivi funkcija:
%   - Tsrednje
%   - STFTobradaMDF
%
% Vedran Srhoj-Egekher, 2008-2009 
%

close all;
clear all;

%%
% uèitavanje signala i njihovih parametara

Ts=1/1000;
Fs=1/Ts;

load('gk002.mat');
signal= gk002(:,1);
kut= gk002(:,4);

% kontrolni bit koji odreðuje da li se prikazuje prikaz za svaku vrijednost
% parametra L, 0=nema prikaza, 1=ima prikaza.
ctrlPrikaza= 0;

L= [51 101 151 201 251 301 351 401 451 501];

%% ############## Obrada MDF sa L kao parametrom, uz R=20ms ###############
% izlazna varijabla STFTpMDF_L

kutLMA= 5; % broj èlanova u pomiènom usrednjavanju pri demodulaciji signala kuta
Tsr= Tsrednje(kut,kutLMA);
%L=floor(0.1*Tsr) + ~mod(floor(0.1*Tsr),2); % širina vremenskog prozora
R=20; % korak pomaka prozora u [ms], tj. R=20ms, i skaliranje u formuli za LMA
LMA= floor(0.1*Tsr/R) + ~mod(floor(0.1*Tsr/R),2); % broj èlanova u pomiènom usrednjavanju


STFTpMDF_L= [];
% prozor L uvijek simetrièan sa neparnim brojem uzoraka
%L= [63 127 255 511 767 1023 1279 1535 1791 2047];

for k=L
    
    % prozor uvijek simetrièan sa neparnim brojem uzoraka
    k % samo da bih mogao pratiti tijek izvoðenja
    
    % elementi u STFTpMDF su poredani u retku
    STFTpMDF_L{end+1}= STFTobradaMDF(ctrlPrikaza,signal,kut,Fs,kutLMA,k,R,LMA);
    
end 


% prikaz promjene nagiba REGRESIJSKOG pravca ovisno o širini prozora L
pom= cell2mat(STFTpMDF_L);

% u STFTpMDF_L1 su spremljeni nagibi REGRESIJSKOG pravca
STFTpMDF_L1= pom(1:2:length(pom))*1000*60;

figure, plot(L,STFTpMDF_L1);
%xlim([0 2100])
xlabel('Širina prozora u [ms]'); ylabel('Nagib REGRESIJSKOG pravca [Hz/min]');
title('STFT:Prikaz promjene nagiba REGRESIJSKOG pravca ovisno o širini prozora L, uz R=20ms');
hold on; 
plot(L,STFTpMDF_L1,'rd');
hold off;



% prikaz promjene poèetne frekvencije ovisno o širini prozora L
pom= cell2mat(STFTpMDF_L);

% u STFTpMDF_L2 su spremljene poèetne frekvencije
STFTpMDF_L2= pom(2:2:length(pom));

figure, plot(L,STFTpMDF_L2);
%xlim([0 2100])
xlabel('Širina prozora u [ms]'); ylabel('Pocetna frekvencija [Hz]');
title('STFT:Prikaz promjene poèetne frekvencije ovisno o širini prozora L, uz R=20ms');
hold on; 
plot(L,STFTpMDF_L2,'rd');
hold off;



%% ################ Obrada MDF sa L kao parametrom, uz R=L ################
% izlazna varijabla STFTpMDF_RL 

kutLMA= 5; % broj èlanova u pomiènom usrednjavanju pri demodulaciji signala kuta
Tsr= Tsrednje(kut,kutLMA);
%R= 20; % R=20ms samo da imamo argument za izraèunavanje LMA u sljedeæem retku
%LMA= floor(0.1*Tsr/R) + ~mod(floor(0.1*Tsr/R),2); % broj èlanova u pomiènom usrednjavanju


STFTpMDF_RL= [];
% prozor L uvijek simetrièan sa neparnim brojem uzoraka
%L= [63 127 255 511 767 1023 1279 1535 1791 2047];

for k=L
    
    % R=L tako da se spekrogrami ne preklapaju
    % zapravo uz ovaj uvjet R=L imamo i zadovoljen (granièni) uvjet 
    % rekonstrukcije GABOROVE EKSPANZIJE T*Omega <= 2*pi
    R=k
    
    % #################################################################
    % pokušaj da se prilagodi broj uzoraka usrednjavanja sa brojem
    % elemenata u vektoru sa frekvencijama medijana MDF
    LMA= floor(0.1*Tsr/R) + ~mod(floor(0.1*Tsr/R),2) %isto kao i LMA=1;
    %
    % #################################################################
    
    % elementi u STFTpMDF su poredani u retku
    STFTpMDF_RL{end+1}= STFTobradaMDF(ctrlPrikaza,signal,kut,Fs,kutLMA,k,R,LMA);
    
end 


% prikaz promjene nagiba REGRESIJSKOG pravca ovisno o širini prozora L i
% pomaku R=L
pom= cell2mat(STFTpMDF_RL);

% u STFTpMDF_RL1 su spremljeni nagibi REGRESIJSKOG pravca
STFTpMDF_RL1= pom(1:2:length(pom))*1000*60;

figure, plot(L,STFTpMDF_RL1);
%xlim([0 2100])
xlabel('Širina prozora u [ms]'); ylabel('Nagib REGRESIJSKOG pravca [Hz/min]');
title('STFT:Prikaz promjene nagiba REGRESIJSKOG pravca ovisno o širini prozora L, uz R=L');
hold on; 
plot(L,STFTpMDF_RL1,'rd');
hold off;



%prikaz promjene poèetne frekvencije ovisno o širini prozora L i pomaku R=L
pom= cell2mat(STFTpMDF_RL);

% u STFTpMDF_RL2 su spremljene poèetne frekvencije
STFTpMDF_RL2= pom(2:2:length(pom));

figure, plot(L,STFTpMDF_RL2);
%xlim([0 2100])
xlabel('Širina prozora u [ms]'); ylabel('Pocetna frekvencija [Hz]');
title('STFT:Prikaz promjene poèetne frekvencije ovisno o širini prozora L, uz R=L'); 
hold on; 
plot(L,STFTpMDF_RL2,'rd');
hold off;


