Wykład nr 11
Zastosowanie transformaty Haara do filtrowania i kompresji sygnału 
Przygotowała Klaudia Wrzask

Procedura Transformaty Haara i konstrukcja sygnału przy jej użyciu
Haar Transform procedure

Procedura transformaty Haara

Clear[output]
DWT[input_List]:=Module[{ output={}},
x=input;
yLP[x_List,n_]:=0.5 (x[[n-1]]+x[[n]]);
yHP[x_List,n_]:=0.5 (x[[n-1]]-x[[n]]);
While[1≤Length[x],
x=PrependTo[x,0];
y0=Table[yLP[x,n],{n,2,Length[x]}];
y0=Take[y0,{2,Length[y0],2}];
y1=Table[yHP[x,n],{n,2,Length[x]}];
y1=Take[y1,{2,Length[y1],2}];
x=y0;
If[Length[y0]==1,y1=Join[y0,y1]];
output=Flatten[AppendTo[y1,output]];
];
output];

DWT[{8,4,2,8,6,6,2,4}]

Odwrotna transformata Haara

UpSample[x_List]:=Module[{i=1,y=x},While[i≤2 Length[x],If[Mod[i,2]==0,y=Insert[y,0,i];i=++i,i=++i]];y];

IDWT[input_List]:=Module[{j=1, output={}},
temp=input;
xLP[y_List,n_]:=y[[n-1]]+y[[n]];
xHP[y_List,n_]:=y[[n]]-y[[n-1]];
While[2^j≤ Length[input],
y0=temp[[Range[1,2^(j-1)]]];
y1=temp[[Range[2^(j-1)+1,2^j]]];
y=UpSample[y0];
y=PrependTo[y,0];
y0=Table[xLP[y,n],{n,2,1+2^j}];
y=UpSample[y1];
y=PrependTo[y,0];
y1=Table[xHP[y,n],{n,2,1+2^j}];
output=temp;
output=Take[output,-(Length[output]-2^j)];
output=Join[y0+y1,output];
temp=output;
j=++j];
output];

IDWT[{5.`,0.5`,0.5`,1.5`,2.`,-3.`,0.`,-1.`}]

Krok 0 (wykonujemy tylko po to aby przekonać się, że dokompozycja Haara prowadzi do tej samej funkcji)
Przybliżanie funkcji za pomocą funkcji schodkowej (skalującej):
Step 0
Aproximating by step functions:

Poniżej przedstawioną funkcję schodkową będziemy nazywali podstawową funkcją skalującą Haara

ϕ[x_]:=UnitStep[x(1-x)]
Plot[ϕ[x],{x,-0.5,1.5}]

Funkcja ϕ[x,j,k]=ϕ[(2^j)x-k] to funkcja skalująca o skali j i translacji (przesunięciu) k. Podstawowa funkcja skalująca jest zatem postaci ϕ[x, 0, 0]. Funkcję ϕ[x,j,k] nazywamy także matką falką.
Poniżej przykładowa rekonstrukcja funkcji sinus przy pomocy funkcji skalujących o skali 5.

ϕ[x_,j_,k_]:=ϕ[2^j x-k];
xp=Table[x,{x,0,1,1/(2^5-1)}];
yp=Table[Sin[2 Pi xp]];
Dimensions[xp]
Dimensions[yp]
ys=Table[ϕ[t,5,k],{k,0,2^5-1}];
Plot[yp.ys,{t,0,1}]

Krok 1
Konstrukcja funkcji Haara przy pomocy funkcji skalującej:
Step 1
Defining Haar wavelet by scaling function:

Funkcję ψ[x]=ϕ[2x]-ϕ[2x-1] nazywamy podstawową falką Haara. Została ona skonstruowana przy pomocy funkcji skalującej.

ψ[x_]:=ϕ[2 x]-ϕ[2 x-1]
Plot[ψ[x],{x,-1,2}]

Funkcję ψ[x,j,k]=ψ[(2^j)x-k] nazywamy falką Haara o skali j i translacji k. Podstawowa falka Haara jest zatem postaci ψ[x,0,0].

ψ[x_,j_,k_]:=ψ[(2^j) x-k]
Plot[ψ[x,1,1],{x,0,1.5}]
Plot[ψ[x,2,2],{x,0,1.5}]

Krok 2
Wykonujemy transformatę Haara spróbkowanej funkcji
Step 2
Haar wavelet traansform of the sampled function:

THaar=DWT[yp];

Krok 3
Wykonyjemy dekompozycję naszej funkcji przy pomocy współczynników rozwinięcia transformaty Haara i funkcji Haara. (przemnażamy odpowiednie współczynniki Haara przez odpowiednie funkcje Haara)
Step 3
Now we can decompose our function with the Haar transform coeffcients and the Haar wavelets

funkcje=Flatten[Table[{ϕ[t],Table[ψ[t,j,k],{j,0,4},{k,0,2^j-1}]}]];

Plot[THaar.funkcje,{t,0,1}]

Własności funkcji Haara - trochę matematyki :)

Niech V0 będzie przestrzenią wektorową funkcji stałych na przedziale półotwartym [0; 1). Przestrzeń (Subscript[V, j]) - jest to przestrzeń zawierająca wszystkie funkcje chwilo stałe, o odcinkach długości 1/(2^j) w przedziale [0; 1).
Przykłady funkcji żyjących na przestrzeni (Subscript[V, 2]):

Plot[2 UnitStep[-x(x-0.25)]+4 UnitStep[-(x-0.25)(x-0.5)]+3 UnitStep[-(x-0.5)(x-0.75)]+1 UnitStep[-(x-0.75)(x-1)],{x,-0.5,1.5}]
Plot[1 UnitStep[-x(x-0.25)]+-4 UnitStep[-(x-0.25)(x-0.5)]-2 UnitStep[-(x-0.5)(x-0.75)]+2 UnitStep[-(x-0.75)(x-1)],{x,-0.5,1.5}]

Potocznie można powiedzieć, że przestrzeń wektorowa (Subscript[V, 0]) posiada jeden przedział [0; 1).
Przestrzeń wektorowa (Subscript[V, 1]) posiada dwa przedziały [0; 1/2) i  [1/2; 1).
Przestrzeń wektorowa (Subscript[V, 2]) posiada cztery przedziały [0; 1/4) i  [1/4; 1/2) i  [1/2; 3/4) i  [3/4; 1), przestrzeń wektorowa (Subscript[V, 3]) posiada osiem przedziałów i.t.d...
Prosta baza dla przestrzeni (Subscript[V, j]) tworzona jest przez zbiór przeskalowanych i przesuniętych funkcji prostokątnych ϕ[x,j,k]. W takim przypadku bazę nazywamy bazą Haara.
Poniżej przykład wektorów bazowych (tutaj wektory bazowe są funkcjami) przestrzeni (Subscript[V, 2]):

phi20=Plot[ϕ[x,2,0],{x,-0.5,1.5}];
phi21=Plot[ϕ[x,2,1],{x,-0.5,1.5}];
phi22=Plot[ϕ[x,2,2],{x,-0.5,1.5}];
phi23=Plot[ϕ[x,2,3],{x,-0.5,1.5}];

GraphicsRow[{phi20,phi21,phi22,phi23}]

(*) Każdy wektor w przestrzeni (Subscript[V, j]) znajdzie swoją reprezentacje w przestrzeni (Subscript[V, j+1]). Twierdzenie odwrotne nie jest prawdziwe.

Definiujemy iloczyn skalarny
⟨f|g⟩=(∫_0^1) f(x) g(x) dx

Teraz można zdefiniować nową przestrzeń wektorową Wj, która będzie ortogonalnym dopełnieniem przestrzeni (Subscript[V, j]) w przestrzeni (Subscript[V, j+1]). Innymi słowy stworzymy przestrzeń wektorową Wj, która będzie przestrzenią wszystkich funkcji w przestrzeni (Subscript[V, j+1]), takich, które będą ortogonalne (iloczyn skalarny równy zero) do wszystkich funkcji w przestrzeni (Subscript[V, j]). Nieformalnie możemy myśleć o (Subscript[W, j]) jako o przestrzeni zawierającej informacje o szczegółach w (Subscript[V, j+1]), które nie mogą być przedstawione w przestrzeni (Subscript[V, j]). Liniowo niezależne funkcje ψ[x,j,k] (wektory) rozpinające przestrzeń (Subscript[W, j]) nazywane są falkami.

Te funkcje bazowe mają dwie ważne właściwości:
(**)   Funkcje bazowe przestrzeni Wj razem z funkcjami bazowymi przestrzeni (Subscript[V, j]) tworzą bazę przestrzeni (Subscript[V, j+1]).
(***) Każda funkcja bazowa przestrzeni Wj jest ortogonalna do każdej funkcji bazowej przestrzeni (Subscript[V, j]).
(** i ***)  Każdy wektor w przestrzeni (Subscript[V, j+1]) znajdzie swoją reprezentacje w przestrzeni (Subscript[V, j]) ⋃ (Subscript[W, j]).

Przykład:
(Subscript[V, 3]) -> (Subscript[V, 2]) ⋃ (Subscript[W, 2])  ->  (Subscript[V, 1])⋃(Subscript[W, 1]) ⋃ (Subscript[W, 2])  ->  (Subscript[V, 0])⋃(Subscript[W, 0])⋃(Subscript[W, 1]) ⋃ (Subscript[W, 2]) 

Przypuśćmy, że mamy sygnał w przedziale [0; 1) spróbkowany (2^3)=8 razy. 
Możemy go przedstawić za pomocą kombinacji liniowej funkcji skalujących ϕ[x,3,k], k=0,(2^3)-1 z odpowiednimi współczynnikami.
Możemy go przedstawić za pomocą kombinacji liniowej funkcji skalujących ϕ[x,2,k] i funkcji Haara ψ[x,2,k], k=0,(2^2)-1 z odpowiednimi współczynnikami (oczywiście nie tymi samymi co powyżej).
Możemy go przedstawić za pomocą kombinacji liniowej funkcji skalujących ϕ[x,1,k] i funkcji Haara ψ[x,1,k], k=0,1 i funkcji Haara ψ[x,2,k], k=0,(2^2)-1 z odpowiednimi współczynnikami.
Możemy go przedstawić za pomocą kombinacji liniowej funkcji skalującej ϕ[x,0,0] i funkcji Haara ψ[x,0,0] i funkcji Haara ψ[x,1,k], k=0,1 i funkcji Haara ψ[x,2,k], k=0,(2^2)-1 z odpowiednimi współczynnikami. 
Zwróćmy uwagę na to, że w każdym z czterech powyższych przypadków mamy osiem funkcji bazowych.
Poniżej możliwe funkcje bazowe dla przestrzeni (Subscript[V, 3])

GraphicsRow[{Plot[ϕ[x,3,0],{x,-0.5,1.5}],Plot[ϕ[x,3,1],{x,-0.5,1.5}],Plot[ϕ[x,3,2],{x,-0.5,1.5}],Plot[ϕ[x,3,3],{x,-0.5,1.5}], Plot[ϕ[x,3,4],{x,-0.5,1.5}],Plot[ϕ[x,3,5],{x,-0.5,1.5}],Plot[ϕ[x,3,6],{x,-0.5,1.5}],Plot[ϕ[x,3,7],{x,-0.5,1.5}]}]

GraphicsRow[{Plot[ϕ[x,2,0],{x,-0.5,1.5}],Plot[ϕ[x,2,1],{x,-0.5,1.5}],Plot[ϕ[x,2,2],{x,-0.5,1.5}],Plot[ϕ[x,2,3],{x,-0.5,1.5}], Plot[ψ[x,2,0],{x,-0.5,1.5}],Plot[ψ[x,2,1],{x,-0.5,1.5}],Plot[ψ[x,2,2],{x,-0.5,1.5}],Plot[ψ[x,2,3],{x,-0.5,1.5}]}]

GraphicsRow[{Plot[ϕ[x,1,0],{x,-0.5,1.5}],Plot[ϕ[x,1,1],{x,-0.5,1.5}],Plot[ψ[x,1,0],{x,-0.5,1.5}],Plot[ψ[x,1,1],{x,-0.5,1.5}], Plot[ψ[x,2,0],{x,-0.5,1.5}],Plot[ψ[x,2,1],{x,-0.5,1.5}],Plot[ψ[x,2,2],{x,-0.5,1.5}],Plot[ψ[x,2,3],{x,-0.5,1.5}]}]

GraphicsRow[{Plot[ϕ[x,0,0],{x,-0.5,1.5}],Plot[ψ[x,0,0],{x,-0.5,1.5}],Plot[ψ[x,1,0],{x,-0.5,1.5}],Plot[ψ[x,1,1],{x,-0.5,1.5}], Plot[ψ[x,2,0],{x,-0.5,1.5}],Plot[ψ[x,2,1],{x,-0.5,1.5}],Plot[ψ[x,2,2],{x,-0.5,1.5}],Plot[ψ[x,2,3],{x,-0.5,1.5}]}]

Odfiltrowanie sygnału przy pomocy transformaty Haara krok po kroku

Krok 0
Przybliżanie funkcji za pomoca funkcji schodkowej (skalującej) dla funkcji spróbkowanej 2^10 razy:
Step 0
Aproximation by step functions of a 2^10 times sampled function:

f[t_]=E^(t/5) Sin[t^2] ;
fsample=Table[N[f[n/2^8]]+0.3 Random[], {n,1,2^10}];
ListPlot[fsample,PlotJoined->True]

Krok 1
Wykonujemy transformatę Haara spróbkowanej funkcji
Step 1
Haar wavelet traansform of the sampled function:

Timing[haarf=DWT[fsample];]

Krok 2
Wykonyjemy dekompozycję naszej funkcji przy pomocy współczynników rozwinięcia transformaty Haara i funkcji Haara.
Step 2
Now we can decompose our function with the Haar transform coeffcients and the Haar wavelets

funkcje=Flatten[Table[{ϕ[t],Table[ψ[t,j,k],{j,0,9},{k,0,2^j-1}]}]];
Plot[haarf.funkcje,{t,0,1}]

Dimensions[funkcje]

Krok 3
Filtrowanie: czyli wyzerowanie współczynników transformaty Haara dla j powyżej pewnej liczby, czyli pominięcie ostatnich współczynników transformaty. 
Technicznie na wykresach poniżej, współczynniki te nie zostały wyzerowane tylko po prostu pominięte. 
Step 3
Filtering: setting Haar transform parameters 0 for j above certain value.

syg=Take[haarf,2^8];
funkcje1=Flatten[Table[{ϕ[t],Table[ψ[t,j,k],{j,0,7},{k,0,2^j-1}]}]];
Plot[syg.funkcje1,{t,0,1}]

syg=Take[haarf,2^7];
funkcje1=Flatten[Table[{ϕ[t],Table[ψ[t,j,k],{j,0,6},{k,0,2^j-1}]}]];
Plot[syg.funkcje1,{t,0,1}]

Kompresja sygnału krok po kroku
Step by step signal compression

Reprezentacja sygnału przy pomocy jego transformaty falkowej, w naszym przypadku transformaty Haara, ma wiele zalet. Jedną z nich jest to, że duża liczba współczynników transformaty okazuje się mało istotna z punktu widzenia amplitudy sygnału. Opuszczenie lub wyzerowanie mało istotnych współczynników wprowadza jedynie małe błędy w sygnale rekonstruowanym, dając tym samym kompresję danych. Taka metoda kompresji jest kompresją stratną.

Krok 0
Przybliżanie funkcji za pomoca funkcji schodkowej (skalującej) dla funkcji spróbkowanej 2^10 razy:
Step 0
Aproximation by step functions of a 2^10 times sampled function:

Clear[f,fsample]
f[t_]=E^(t/5) Sin[t^2]+0.1 Sin[1000 n] ;
fsample=Table[N[f[n/2^8]], {n,1,2^10}];
ListPlot[fsample,PlotJoined->True]

Krok 1
Wykonujemy transformatę Haara sprobkowanej funkcji
Step 1
Haar wavelet traansform of the sampled function:

haarf=DWT[fsample];

Krok 2
Wykonyjemy dekompozycję naszej funkcji przy pomocy współczynników rozwinięcia transformaty Haara i funkcji Haara.
Step 2
Now we can decompose our function with the Haar transform coeffcients and the Haar wavelets

funkcje=Flatten[Table[{ϕ[t],Table[ψ[t,j,k],{j,0,9},{k,0,2^j-1}]}]];
Plot[haarf.funkcje,{t,0,1}]

Krok 3
Kompresja: Wyzerowanie wszystkich współczynników transformaty Haara poniżej pewnej wartości
Step 3
Compression: setting Haar transform parameters 0 below a certain threshold

Zacznijmy od wyzerowania wszystkich współczynników poniżej wartości 0.01. Takich współczynników jest 120, zatem jest to kompresja ok. 12 %. Jak widać taka kompresja nie zmienia widocznie postaci sygnału.
Lets start from setting 0 all the parameters below 0.01. There are 120 such parameters, so our compression is about 12%. As we can see this commpression does not visibly change the signal.

Dimensions[Select[haarf,Abs[#]<0.01&]]
N[Dimensions[Select[haarf,Abs[#]<0.01&]]/Dimensions[haarf]]
sygcomp=Chop[haarf,0.01];
Plot[sygcomp.funkcje,{t,0,1}]

Tak wygląda kompresja ok. 54%:

Dimensions[Select[haarf,Abs[#]<0.04&]]
N[Dimensions[Select[haarf,Abs[#]<0.04&]]/Dimensions[haarf]]
sygcomp=Chop[haarf,0.05];
Plot[sygcomp.funkcje,{t,0,1}]

A tak ok. 88%

Dimensions[Select[haarf,Abs[#]<0.08&]]
N[Dimensions[Select[haarf,Abs[#]<0.08&]]/Dimensions[haarf]]
sygcomp=Chop[haarf,0.08];
Plot[sygcomp.funkcje,{t,0,1}]

Pytania na kolokwium

1. Funkcje Haara i funkcje skalujące - podać definicje i narysować przykładowe funkcje -  p.p 
2. Własności matematyczne funkcji Haara i skalujących funkcji Haara - p.r. 
3. Na czym polega rekonstrukcja sygnału przy pomocy transformaty Haara? - p.p 
4. Jak przebiega filtrowanie (odszumianie) sygnału przy pomocy transformaty Haara? -  p.r
5. Jak przebiega kompresja sygnału przy pomocy transformaty Haara? - p.r

Zadanie na laboratorium

Zadanie nr 13 - 4p.

Należy
1. Wymyślić dowolną funkcje f(t) (pewien sygnał - każdy powinien mieć inny) i ją zaszumić (dowolną metodą)
2. Narysować wykres f(t) i zaszumioną f(t)
3. Spróbkować funkcję (sygnał), tak aby otrzymać przynajmniej 2^10 próbek
4. Odfiltrować sygnał na pierwszych trzech poziomach rozdzielczości korzystając z transformaty Haara
5. Narysować wykresy otrzymanych funkcji
6. Proszę opisać nowo poznaną funkcję Mathematici plus własne przykłady