%% ########################## MDF obrada CWT-om ###########################
% Skripta (Funkcija) kojoj vršimo glavnu CWT 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

% Morletov wavelet
%Fc=812.5;

Ts=1/1000;
Fs=1/Ts; % frekvencija otipkavanja

% uèitanje EMG signala
load('gk002.mat');
signal= gk002(:,1);
kut= gk002(:,4);

% definiramo ime WAVELETA, npr. 'morl', 'cmor1-1' 
waveime= 'morl'; %cmor"Fb"-"Fc", Fb=širina, Fc=centralna freq
%waveime= 'cmor1-1';
%waveime= 'db4';
%waveime= 'sym4';

% za brojPodjelaskale= 25 dobivamo gotovo iste rezultate kao za
% brojPodjelaskale= 490, a vrijeme izvoðenja je neusporedivo kraæe
brojPodjelaSkale= 25;
Fd= 10; % donja frekvencija
Fg= 250; % gornja frekvencija

% za F=Fg, skala A je najmanja 
A1 = scal2frq(Fg, waveime, Ts);

% za F=Fd, skala A je najveæa
A2 = scal2frq(Fd, waveime, Ts);

logKoef= log(A2)/log(A1);

raspon=linspace(1,logKoef,brojPodjelaSkale);
skale= A1.^raspon;
frekvencije= scal2frq(skale, waveime, Ts);
frekvencije= fliplr(frekvencije);

% izraèunavanje ili uèitavanje wavelet koeficijenata
%load('cwt_morl.mat');
%load('cwt_cmor1-1.mat');
X = cwt(signal, skale, waveime);

% 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 CWTpMDF_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


CWTpMDF_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 CWTpMDF su poredani kaskadno u retcima
    CWTpMDF_L{end+1}= CWTobradaMDF(ctrlPrikaza,signal,kut,X,Fs,kutLMA,k,R,LMA,frekvencije,waveime);
    
end 


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

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

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



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

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

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



%% ################ Obrada MDF sa L kao parametrom, uz R=L ################
% izlazna varijabla CWTpMDF_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

CWTpMDF_RL= [];
% 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
    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 CWTpMDF su poredani kaskadno u retcima
    CWTpMDF_RL{end+1}= CWTobradaMDF(ctrlPrikaza,signal,kut,X,Fs,kutLMA,k,R,LMA,frekvencije,waveime);
    
end 


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

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

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



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

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

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