clc
clear
format(8)
close()

//Wczytanie plików

Meteo = fscanfMat('J:\ćwiczenie 2\xd\meteo_dane_ok.txt')
Nivel = fscanfMat('J:\ćwiczenie 2\xd\Nivel_poz3_ok.txt')
Pkt407 = fscanfMat('J:\ćwiczenie 2\xd\Punkt_407_ok.txt')

//Zamiana daty i godziny na czas GPS i część dnia

data_m=datenum(Meteo(:,3),Meteo(:,2),Meteo(:,1),Meteo(:,4),Meteo(:,5),Meteo(:,6))-datenum(Meteo(:,3),1,1,00,00,00)-1
data_n=datenum(Nivel(:,3),Nivel(:,2),Nivel(:,1),Nivel(:,4),Nivel(:,5),Nivel(:,6))-datenum(Nivel(:,3),1,1,00,00,00)-1
data_p=datenum(Pkt407(:,8),Pkt407(:,7),Pkt407(:,6),Pkt407(:,9),Pkt407(:,10),Pkt407(:,11))-datenum(Pkt407(:,8),1,1,00,00,00)-1

Meteo(:,1:6)=[]
Nivel(:,1:6)=[]
Pkt407(:,6:11)=[]

M =[data_m Meteo]
N =[data_n Nivel]
P =[data_p Pkt407(:,6:$)]


//metoda regresji liniowej dla danych z tachimetru - wielkość Hz

OUTLIER = [];
[wp, kp] = size(P);
mV = 0.0003;
for i = 1: wp-30
    OKNO = P(i:i+30, :);
    BADANA = P(i+15, :);
    A = [OKNO(:, 1), ones(31,1)];
    L = [OKNO(:, 2)];
    x = pinv(A'*A) * A'*L;
    Vteor = x(1,1) *OKNO(:, 1)+x(2,1);

    V = Vteor - OKNO(:,2)
    VV = V.^2;
    m0 = sqrt(sum(VV)/(length(OKNO)-1));

    if m0 > 2*mV
        //rownianie regresji liniowe nie pasuje - nie ma trendu liniowego
        //znalezc inna metode na szukanie outlierow

    else //wystepuje trend liniowy
    [maxval, maxrow]=max(VV);
        if sqrt(maxval) > 5*mV
        //OUTLIER (odstająca)
            OUTLIER=[OUTLIER; P(i+maxrow-1,:)];
            unique_outliers = unique(OUTLIER, 'r')
        else
        //dane ok, trend w oknie bez outlierow
        end
        //wyswietlenie outlierów na wykresie

    end
end

//wyrzucanie outlierów

Removed = P //kopia macierzy P potrzebna do wykresu

//tworzenie macierzy indeksów wierszy z kolumny P do wyrzucenia
for i = 1:size(unique_outliers,1)
    DEL(i) = vectorfind(P,unique_outliers(i,:),'r')
end

P(DEL,:) = [] //wyrzucenie outlierów

figure()
//wyswietlenie outlierów na wykresie
subplot(2,1,1)
plot(Removed(:, 1), Removed(:, 2), unique_outliers(:,1), unique_outliers(:,2), '*r')
//wyswietlenie wykresu bez outlierow
subplot(2,1,2)
plot(P(:, 1), P(:, 2))