Down/up-sampling

Jako dane weźniemy przebieg ze spotkania 7. W pierwszej części pomiary są co godzinę, a później w nieregularnych odstępach:

t = czytaj_historie("temperatury.csv");
plot(t.last_changed,t.state)

size(t)

ans = 1x2
   514     2

Czas pomiaru

delta_t = hours(t.last_changed(end) - t.last_changed(1))

delta_t = 305.1252

Popatrzmy jak często robione są pomiary. Można użyć do tego funkcję diff(), która wylicza różnice pomiędzy kolejnymi warrtościami tablicy

delta = diff(t.last_changed);
min(delta)

ans = duration
   00:00:15

max(delta)

ans = duration
   02:39:32

mean(delta)

ans = duration
   00:35:41

duration oznacza tu odstęp czasu.

Można skorzystać z histogramu, ale, nie będzie to miało zbyt wielkiego sensu

% histogram(delta)
delta(1)

ans = duration
   01:00:00

Na początku co godzinę, a póxniej jakoś dziwnie

delta(end)

ans = duration
   00:18:45

Interesuje nas sytuacja gdy będziemy mieli pomiary z N punktów. Alternatywnie możemy zażądać, żeby pomiary były co godzinę…

Wyznaczony przedział czsowy to

N = 64;
delta_t

delta_t = 305.1252

Zatem chcemy aby pomiary były robione co

dt = delta_t/(N - 1)

dt = 4.8433

Dt = hours(dt)

Dt = duration
   4.8433 hr

czas = t.last_changed(1):Dt:t.last_changed(end);
size(czas)

ans = 1x2
     1    64

Możemy teraz użyć interpolacji:

Dane wejściowe: t.last_changed i t_state; czas to „nowe punkty pomiarowe”.

temperatury = interp1(t.last_changed, t.state, czas);

I wyrysować dwa wykresy (ograniczam nieco zakres, żeby lepiej było widać)

plot(t.last_changed, t.state,'-o', czas, temperatury, '+')
xlim([datetime(2024,11,7,18,40,00)...
      datetime(2024,11,8,12,00,00)])
ylim([7 10])
legend("przebieg oryginalny","interpolacja")

Teraz możemy zrobić FFT

Okres pomiarów to:

dt

dt = 4.8433

zatem częstotliwość

Fs = 24/dt % razy na dobę

Fs = 4.9553

Liczba pomiarów

L = size(temperatury, 2)

L = 64

Y = fft(temperatury);
plot(abs(fftshift(Y)))

Zachowuję nazwy wszystkich zmiennych, żeby nie tłumaczyć procedury skalowania…

P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = Fs*(0:(L/2))/L;
stem(f(f<2),P1(f<2),'r')
title("Współczynniki Fouriera")
ylabel("Amplituda")
xlabel("razy na dobę")

Wyraźnie widać, że występuje jakaś skłądowa harmoniczna o częstości zbliżónej do „raz na dobę”, ale uzyskana rozdzielczość nie bardzo pozwala dokładnie określić jaka to częstość.

Żeby zwiększyć rozdzielczość trzeba przedłużyć czas pomiarów… Ale takich danych nie mamy