Analiza sygnałów

Transformata Fouriera na komputerze działa na skończonej, spróbkowanej porcji danych — i każde z tych trzech ograniczeń zostawia ślad. Skończoność daje przeciek widmowy, próbkowanie daje aliasing, a dyskretność zamienia całkę w sumę, którą FFT liczy setki razy szybciej niż definicja. Ten moduł pokazuje wszystkie trzy zjawiska na policzonych danych, a nie na rysunku poglądowym.

Moduł łączy zakres „analizy sygnałów” (P2) i „analizy harmonicznej” (P3) ze specyfikacji — obie pozycje pokrywały się na sekcji o aliasingu i twierdzeniu o próbkowaniu, więc zostały scalone.

O(N log N)

DFT i FFT

Dyskretna transformata z definicji i szybka — ten sam wynik, setki razy szybciej

X_k = Σ xₙ e^(−2πikn/N)
f_s > 2f_max

Próbkowanie i aliasing

Dwie różne częstotliwości dające identyczne próbki — twierdzenie Nyquista-Shannona

f_alias = |f − k·f_s|
|H(ω)|

Filtracja cyfrowa

Filtr FIR i IIR, charakterystyka częstotliwościowa liczona ze wzoru i zmierzona

y = h ∗ x ↔ Y = H·X
PRZECIEK

Przeciek widmowy i okna

Dlaczego jedna czysta sinusoida rozmazuje się na całe widmo — i jak temu zaradzić

okno Hanna, Hamminga
K = P/(P+R)

Filtr Kalmana

Estymacja trajektorii z zaszumionych pomiarów, z porównaniem błędów

predykcja → korekta
ZADANIA

Przykłady krok po kroku

Rozdzielczość DFT, częstotliwość pozorna, wzmocnienie Kalmana

3 rozwiązane zadania

Dyskretna transformata Fouriera

X_k = Σ_{n=0}^{N−1} xₙ · e^(−2πikn/N)
Wprost z definicji to N² mnożeń. FFT wykorzystuje symetrię e^(−2πi/N) i dzieli zadanie na pół, uzyskując N·log₂N — dla N = 1024 to różnica rzędu 200 razy.
rozdzielczość widma Δf = f_s/N–
zakres widma (Nyquist)–
max |FFT − DFT z definicji|–
czas: DFT / FFT–
Parseval: Σ|xₙ|² kontra (1/N)Σ|X_k|²–
sygnał w dziedzinie czasu
widmo amplitudowe |X_k|

Twierdzenie o próbkowaniu i aliasing

Nyquist-Shannon: sygnał o paśmie ograniczonym przez f_max da się dokładnie odtworzyć z próbek, o ile f_s > 2·f_max. Powyżej tej granicy wysokie częstotliwości podszywają się pod niskie i nic już ich nie rozróżni.
f_alias = |f − round(f/f_s)·f_s|
granica Nyquista f_s/2–
częstotliwość pozorna f_alias–
czy występuje aliasing–
max różnica próbek f kontra f_alias–
szczyt w widmie DFT–
sygnał prawdziwy, próbki i sygnał pozorny
widmo — gdzie naprawdę ląduje energia

Filtracja cyfrowa

Filtr o odpowiedzi impulsowej h działa przez splot: y = h ∗ x. W dziedzinie częstotliwości to zwykłe mnożenie Y = H·X — dlatego charakterystyka |H(ω)| mówi wszystko o tym, co filtr robi z którą częstotliwością.
częstotliwość odcięcia (−3 dB)–
max ||H| zmierzone − wzór|–
RMS szumu przed / po–
RMSE wobec sygnału czystego–
opóźnienie grupowe–
sygnał zaszumiony i po filtracji
charakterystyka |H(ω)| — wzór i pomiar

Przeciek widmowy i funkcje okna

DFT zakłada, że sygnał jest okresowy z okresem równym długości okna. Jeśli nie mieści się w nim całkowita liczba okresów, powstaje skok na granicy — a skok ma szerokie widmo. Stąd przeciek: jedna czysta sinusoida rozmazuje się na wszystkie prążki.
czy k jest całkowite–
energia poza listkiem głównym–
szerokość listka głównego–
poziom pierwszego listka bocznego–
sygnał i obwiednia okna
widmo w decybelach — widać przeciek

Filtr Kalmana

Filtr na przemian przewiduje stan z modelu ruchu i koryguje go pomiarem. Wzmocnienie K decyduje, komu bardziej ufać:
K = P⁻/(P⁻ + R)
x = x⁻ + K·(z − x⁻)
P = (1 − K)·P⁻
Duże R (zaszumiony czujnik) → małe K → filtr trzyma się modelu. Duże Q (nieprzewidywalny ruch) → duże K → filtr ufa pomiarom.
RMSE surowych pomiarów–
RMSE filtru Kalmana–
RMSE średniej ruchomej (odniesienie)–
poprawa względem pomiarów–
wzmocnienie K w stanie ustalonym–
wariancja P w stanie ustalonym–
trajektoria prawdziwa, pomiary i estymata
wariancja P i wzmocnienie K w czasie

Przykład 1 — rozdzielczość widma i czas pomiaru

Rejestrujemy sygnał z częstotliwością próbkowania f_s = 1000 Hz przez N = 500 próbek. Jakie częstotliwości rozróżnimy?

1. Czas trwania rejestracji: T = N/f_s = 500/1000 = 0,5 s.
2. Rozdzielczość widma: Δf = f_s/N = 1000/500 = 2 Hz. Prążek k odpowiada częstotliwości k·Δf.
3. Równoważnie Δf = 1/T — i to jest sedno: rozdzielczość zależy wyłącznie od czasu obserwacji, nie od częstotliwości próbkowania. Chcąc rozróżnić 440 Hz od 441 Hz, trzeba mierzyć co najmniej 1 sekundę, choćby próbkować milion razy na sekundę.
4. Zakres widma: od 0 do f_s/2 = 500 Hz (granica Nyquista). Prążki k = 0…250 niosą informację; wyższe są lustrzanym odbiciem (dla sygnału rzeczywistego X_{N−k} = konjugata X_k).
5. Sygnał 300 Hz trafi w prążek k = 300/2 = 150. Sygnał 301 Hz nie trafi w żaden prążek dokładnie — jego energia rozleje się na sąsiednie (przeciek widmowy, zakładka „Przeciek widmowy i okna”).
6. Sygnał 700 Hz przekracza Nyquista i pojawi się jako 1000 − 700 = 300 Hz, czyli w tym samym prążku co prawdziwe 300 Hz. Po spróbkowaniu nie da się ich rozróżnić — dlatego przed przetwornikiem A/C zawsze stoi analogowy filtr antyaliasingowy.
7. Koszt obliczeń: DFT z definicji to N² = 250 000 mnożeń zespolonych, FFT to N·log₂N ≈ 500 · 9 = 4 500. Zakładka „DFT i FFT” mierzy oba czasy i porównuje wyniki — różnią się o ok. 10⁻¹², czyli tylko o błąd zaokrągleń.

Przykład 2 — częstotliwość pozorna

Sygnał sin(2π·18·t) próbkujemy z f_s = 20 Hz. Co zobaczymy?

1. Granica Nyquista to f_s/2 = 10 Hz, a 18 > 10 — będzie aliasing.
2. Próbki: x_n = sin(2π·18·n/20) = sin(2π·0,9·n).
3. Korzystamy z okresowości sinusa: 2π·0,9·n = 2π·n − 2π·0,1·n, a sin(2πn − θ) = −sin(θ) dla całkowitego n. Zatem
   x_n = −sin(2π·0,1·n) = −sin(2π·2·n/20).
4. Próbki są identyczne jak dla sygnału o częstotliwości 2 Hz (z odwróconą fazą). Ogólnie f_alias = |f − round(f/f_s)·f_s| = |18 − 20| = 2.
5. Konsekwencja: mając tylko próbki, nie ma żadnego sposobu odróżnić 18 Hz od 2 Hz. To nie jest niedoskonałość algorytmu — to utrata informacji przy próbkowaniu. Aplikacja mierzy różnicę między próbkami obu sygnałów: wychodzi rzędu 10⁻¹³, czyli zero.
6. Znane skutki: koła wozu kręcące się „do tyłu” w filmie (f_s = 24 klatki/s), pasy mory na zdjęciu drobnej kratki, fałszywy niski ton w źle spróbkowanym nagraniu.
7. Przypadek graniczny f = f_s/2 = 10 Hz jest szczególny: próbki wypadają wtedy w tych samych fazach co pół okresu, więc amplituda zależy od fazy początkowej i może wyjść nawet zero. Dlatego twierdzenie wymaga f_s > 2f_max ostro, a nie ≥.

Przykład 3 — wzmocnienie Kalmana w stanie ustalonym

Rozważmy najprostszy przypadek: stała wielkość x, model bez dynamiki (x_{k+1} = x_k), szum procesu Q i szum pomiaru R. Jak zachowuje się wzmocnienie K?

1. Predykcja: P⁻_k = P_{k−1} + Q (niepewność rośnie o szum procesu).
2. Korekta: K_k = P⁻_k/(P⁻_k + R), a nowa niepewność P_k = (1 − K_k)·P⁻_k.
3. W stanie ustalonym P_k = P_{k−1} = P. Podstawiając:
   P = (1 − K)(P + Q), gdzie K = (P+Q)/(P+Q+R).
4. Oznaczmy A = P + Q. Wtedy 1 − K = R/(A+R), więc P = A·R/(A+R), czyli A − Q = A·R/(A+R). Mnożąc przez (A+R): A² + AR − QA − QR = AR, stąd A² − QA − QR = 0.
5. A = [Q + √(Q² + 4QR)]/2 (bierzemy pierwiastek dodatni), a K = A/(A+R).
6. Sprawdźmy granice:
   • Q → 0 (idealny model): A → 0, więc K → 0 — filtr przestaje słuchać pomiarów i uśrednia coraz dłużej.
   • R → 0 (idealny czujnik): A → Q, więc K = Q/(Q+R) → 1 — filtr bierze pomiar bez poprawki.
7. Dla Q = R = 1: A = (1 + √5)/2 = φ ≈ 1,618, a K = φ/(φ+1) = φ/φ² = 1/φ ≈ 0,618. Złota proporcja pojawia się tu z tego samego równania kwadratowego, co w ciągu Fibonacciego — miły przypadek.
8. Sens praktyczny: K to waga pomiaru w średniej ważonej między przewidywaniem a obserwacją. Filtr Kalmana jest optymalny w sensie najmniejszego błędu średniokwadratowego spośród wszystkich estymatorów liniowych — zakładka mierzy jego RMSE i porównuje z surowymi pomiarami oraz ze średnią ruchomą.