Wykład nr 10
Wstęp do transformaty falkowej Haara
Przygotowała Klaudia Wrzask

Oknowa transformata Fouriera i główna motywacja dla transformaty falkowej

Jakie informacje otrzymujemy wykonując transformatę Fouriera?

Clear[recwindow,sygnal,sygnalwin1,spectogram1,sygnalwin2,spectogram2]

Rozważmy sygnał, w którym częstotliwość jest zmienna w czasie np. Let's consider a signal in which frequency changes in time.

sygnal[t_]:=Sin[2 Pi t] UnitStep[t (3-t)]+Sin[4 Pi t] UnitStep[(3-t)(t-6)]+Sin[6 Pi t]UnitStep[(6-t)(t-9)]
plotsyg=Plot[sygnal[t], {t,0,10},PlotRange->All]

Transformata Fouriera daje nam odpowiedź na pytanie jakie częstości wchodzą w skład sygnału, nie wiadomo jednak w jakim czasie występuje dana częstotliwość:
Fourier Transform gives us answers about the frequencies in a signal,but doesn't give us answer about their time.

plotFTsyg=Plot[Evaluate[Abs[FourierTransform[sygnal[t],t,s]]], {s,0,30},PlotRange->All]

A teraz weźmy pod uwagę ten sam sygnał "puszczony od tyłu":

sygnal2[t_]:=Sin[6 Pi t] UnitStep[t (3-t)]+Sin[4 Pi t] UnitStep[(3-t)(t-6)]+Sin[2 Pi t]UnitStep[(6-t)(t-9)]
plotsyg2=Plot[sygnal2[t], {t,0,10}]

plotFTsyg2=Plot[Evaluate[Abs[FourierTransform[sygnal2[t],t,s]]], {s,0,30}]

plotFTArgsyg=Plot[Evaluate[Arg[FourierTransform[sygnal[t],t,s]]], {s,0,30},PlotRange->All];
plotFTArgsyg2=Plot[Evaluate[Arg[FourierTransform[sygnal2[t],t,s]]], {s,0,30}];

Jak widać oba sygnały mają jednakowe widmo wartości bezwzglednej. Mając zatem tylko transformatę Fouriera tych sygnałów nie byli byśmy w stanie stwierdzić jak wygląda rozkład częstotliwości w czasie. Nie wiedzielibyśmy np."czy melodia grana jest do przodu czy od tyłu".
As You can see both signals have the same Fourier Transform. Knowing only this transform we wouldn't know how the frequencies change in time. It's like not knowing if the melody is played forward or backward.

GraphicsGrid[{{plotsyg,plotsyg2},{plotFTsyg,plotFTsyg2},{plotFTArgsyg,plotFTArgsyg2}}]

Oknowa Transformata Fouriera (WFT):
Krok 1
Definiujemy okno: w tym przypadku będzie to okno prostokątne o szerokości L i początku a.

recwindow[t_,L_,a_]:= UnitStep[(t-a) (L-(t-a))]
Plot[{recwindow[t,1,2],recwindow[t,2,4]},{t,0,10}]

Krok 2
Patrzymy na sygnał przez okno tzn. wycinamy kawałek sygnału. W tym przypadku nasze okno ma szerokość 3:

sygnalwin1[t_,a_]:=sygnal[t] recwindow[t,3,a]
Plot[sygnalwin1[t,0],{t,0,5}]

Krok 3
Wykonujemy transformatę Fouriera takiego wycietego sygnału. W kolejnych etapach będziemy przesuwać okno i wykonywać kolejne transformaty. W ten sposób WFT będzie zależna od parametru a, który jest po prostu początkiem naszego przesuwającego się okna.

spectogram1[s_,a_]:=Abs[Integrate[sygnalwin1[t,a] E^(I s t),{t, a, a+3},Assumptions->Im[a]==0]]

Activate[spectogram1[4,4]]

Krok 4
Wykonujemy spektogram,czyli wykres zależności częstotliwości sygnału od czasu:w tym przypadku czas określa początek okna,czyli parametr a.

Date[]
DensityPlot[spectogram1[s,a],{a,0,12},{s,0,24},PlotPoints->12]
Date[]

Zasada nieoznaczoności

Sprawdźmy jaki otrzymamy spectogram przyjmując węższe  okno np. L=1: Lets see how a spectogram will look like for a more narrow window L=1.

sygnalwin2[t_,a_]:=sygnal[t] recwindow[t,1,a]
Plot[sygnalwin2[t,0],{t,0,5}]

spectogram2[s_,a_]:=Abs[Integrate[sygnalwin2[t,a] E^(I s t),{t, a, a+1},Assumptions->Im[a]==0]]
DensityPlot[spectogram2[s,a],{a,0,12},{s,0,24},PlotPoints->12]

Widać, że przy zwężeniu okna, czyli dokładniejszym określeniu czasu tracimy pewną infirmację o częstościach (większe rozmycie na osi rzędnych). Jest to oczywiście zgodne z teorią nieoznaczoności.
As we can see, by taking a more narrow window, i.e. defining time more precise, we loose some information about the frequency. This of course can be explained by the uncertainty theory.

Podstawowe definicje:

Filtr górnoprzepustowy - Highpass filter (H0 - oznaczenie na diagramie)

yHP[x_List,n_]:=x[[n-1]]/2-x[[n]]/2

Clear[x];
x=Insert[Table[0,{n,1,19}],1,10]

Table[yHP[x,n],{n,2,20}]

Filtr dolnoprzepustowy - Lowpass filter (G0 - oznaczenie na diagramie)

yLP[x_List,n_]:=x[[n-1]]/2+x[[n]]/2

Clear[x];
x=Insert[Table[0,{n,1,19}],1,10]

Table[yLP[x,n],{n,2,20}]

Dwukrotne zmiejszenie liczby próbek - Downsampling

DownSample[x_List]:=Module[{y},y=Take[x,{2,Length[x],2}];y];

DownSample[{1,2,4,7,9,3,5,20}]

Procedura transformaty falkowej Haara Krok po kroku
Haar wavelet transform Step by Step

Krok 1 - Dany jest sygnał wejściowy x; wstawiany 0 na początek by móc skorzystać z procedury yLP i yHP
Step 1 - The input signal is x; we put a 0 at the begining of x so we can use the yLP and yHP.

Clear[x];
x={3,5,9,2};
x=PrependTo[x,0]

Krok 2 - Sygnał wejściowy przepuszczamy przez filtry górno i dolno przpustowe
Step 2 - We put the input signal through the lowpass and highpass filter

d1=Table[yHP[x,n],{n,2,5}]

a1=Table[yLP[x,n],{n,2,5}]

Krok 3 - Dwukrotnie zmiejszamy liczbę próbek tak otrzymanych sygnałów wyjsciowych
Step 3 - Downsampling

d12n=DownSample[d1]

a12n=DownSample[a1]

Krok 4 - Teraz bierzemy pod uwagę sygnał wejściowy  przepuszczony przez filtr dolnoprzepustowy o dwukrotnie zmiejszonej liczbie próbek i powtarzamy całą procedurę :
Step 4 - Now we consider only downsampled lowpass filtered signal:

a12n=PrependTo[a12n,0]

d2=Table[yHP[a12n,n],{n,2,Length[a12n]}]

a2=Table[yLP[a12n,n],{n,2,3}]

d22n=DownSample[d2]

a22n=DownSample[a2]

Krok 5 - Podsumowując współczynniki rozwiniecia falkowego to {a2(2n) , d2(2n) , d1(2n)}
Krok 5 - These are the wavelet transform coefficients {a2(2n) , d2(2n) , d1(2n)}

Wspoczynniki={19/4,-3/4,-1,(7/2)}

Wspoczynniki=Join[a22n,d22n,d12n]

Procedura odwrotnej transformaty falkowej Haara Krok po kroku
Inverse Haar wavelet transform Step by Step

Krok 1 - Tym razem mamy dane współczynniki transformaty Haara:
Step 1 - This time our data are the wavelet transform coefficients:

Wspoczynniki={19/4,-3/4,-1,(7/2)}

Highpass filter (H1-oznaczenie na diagramie)

xHP[y_List,n_]:=y[[n]]-y[[n-1]]

Lowpass filter (G1-oznaczenie na diagramie)

xLP[y_List,n_]:=y[[n-1]]+y[[n]]

Dopełnianie zerami:
Upsampling

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]

UpSample[{1,2,4,20}]

Krok 2 - Dopełniamy zerami dwa pierwsze współczynniki transformaty
Step 2 - We upsample two first coefficients

a2n=UpSample[a22n];
Dimensions[a2n]
PrependTo[a2n,0]
Dimensions[a2n]

d2n=UpSample[d22n];
PrependTo[d2n,0]

Krok 3 - Otrzymane wartości odfiltrowujemy wg. procedur xHP i xPL. Otrzymane winiki dodajemy do siebie.
Step 3 - The results we filter with xHP i xPL. The results we add together.

lp=Table[xLP[a2n,n],{n,2,3}]

hp=Table[xHP[d2n,n],{n,2,3}]

a12n=lp+hp

Krok 4 - Otrzymany w ten sposób sygnał znowu dopełniamy zerami i odfiltrowujemy procedurą xLP.
Step 4 - The result we upsample and filter with xLP

a1n=UpSample[a12n];
PrependTo[a1n,0]

lp=Table[xLP[a1n,n],{n,2,5}]

Krok 5 - Współczynniki Transformaty trzeci i czwatry (d12n) dopełniamy zerami i odfiltrowujemy prz pomocy procedury xHP
Step 5 - We upsample the third and fourth coefficient and filter them with xHP

d1n=UpSample[d12n];
PrependTo[d1n,0]

hp=Table[xHP[d1n,n],{n,2,5}]

Krok 6
Aby otrzymać sygnał wejsciowy nalezy już tylko dodac do siebie dwa ostatnie wyniki z kroku 4 i 5.
Step 6
Now we add the results form step 4 and 5. These are the coefficients of the inverse Haare wavelet transform.

x=lp+hp

Filtry dolno i górnoprzepustowe interpretacja graficzna

Przykład 1
W jaki sposób przy pomocy filtrów górno i dolno przepustowych wkonać analizę funkcji:
We are going to analyse a function using lowpass and highpass filters:

Do tego ćwiczenia użyję funkcji f(t)=(e^(t/5)) sin((t^2)) spróbkowanej 256 razy i zaszumionej.
For this exercise we are going to use  f(t)=(e^(t/5)) sin((t^2)) sampled 256 times with a noise.

f[t_]:=E^(t/5) Sin[t^2]

Plot[f[t],{t,0,4}, PlotRange-> All]

Clear[x]
x=Table[N[f[n/2^6]]+Random[], {n,1,2^8}];

sygnal=ListPlot[x,Joined->True]

Sygnał wejściowy przepuszczamy przez filtry górno i dolno przpustowe
We lowpass and highpass filter our input data

x=PrependTo[x,0];
d1=Table[yHP[x,n],{n,1,(2^8)+1}];
a1=Table[yLP[x,n],{n,1,(2^8)+1}];

Dwukrotnie zmiejszamy liczbę próbek tak otrzymanych sygnałów wyjsciowych. Rysujemy wykresy otrzymanych próbek. Jak widać filtr dolnoprzepustowy powoduje wygładzenie funkcji i odpowiada za "wygląd ogólny" funkcji natomiast filtr górnoprzepustowy za szczegóły.
We downsample the results and plot them. As You can see the lowpass filter smoothes the function and shows the general form and the hightpass filter shows the detales.

d1sygnal=ListPlot[DownSample[d1] , Joined->True]
a1sygnal=ListPlot[DownSample[a1], Joined->True]

Analizując dolny wykres można dojść do wnoisku, że funkcja została skompresowana, a pewne informacje zatracone. Zatem do pełnego jej odtworzenia potrzebne są dodatkowe informacje i są to właśnie współczynniki DownSample[d1].
Analizę naszej funkcji możemy także przeprowadzić na kolejnym poziomie próbkowania:

Analysing the upper plot we can come to a conclusion that we lost some information during such a compresion.  This means that for the reconstruction we need some extra information which are the DownSample[y1] coefficients.
We can also analyse our function on other levels of resolution:

a12n=DownSample[a1];
PrependTo[a12n,0];

d2=Table[yHP[a12n,n],{n,1,2^7+1}];
a2=Table[yLP[a12n,n],{n,1,2^7+1}];

d2sygnal=ListPlot[DownSample[d2], Joined->True]
a2sygnal=ListPlot[DownSample[a2], Joined->True]

I jeszcze na kolejnym...
And one more time...

a22n=DownSample[a2];
PrependTo[a22n,0];
a3=Table[yLP[a22n,n],{n,1,2^6+1}];
d3=Table[yHP[a22n,n],{n,1,2^6+1}];
d3sygnal=ListPlot[DownSample[d3], PlotJoined->True]
a3sygnal=ListPlot[DownSample[a3], PlotJoined->True]

GraphicsGrid[{{sygnal},{a1sygnal,d1sygnal},{a2sygnal,d2sygnal},{a3sygnal,d3sygnal}}]

Przykład 2

f[t_]:=(0.2 Sin[2 Pi 8 t]+Sin[2 Pi 2 t]) UnitStep[-t(t-1)]+(0.2 Sin[2 Pi 16 t]+Sin[2 Pi 2 t]) UnitStep[-(t-2)(t-1)]+(0.2 Sin[2 Pi 32 t]+Sin[2 Pi 2 t]) UnitStep[-(t-2)(t-3)]+(0.2 Sin[2 Pi 64 t]+Sin[2 Pi 2 t]) UnitStep[-(t-3)(t-4)];
Plot[f[t],{t,0,4}, PlotRange-> All]

Clear[x]
x=Table[N[f[n/2^8]], {n,1,2^10}];
sygnal=ListPlot[x,PlotJoined->True];
x=PrependTo[x,0];
d1=Table[yHP[x,n],{n,1,(2^10)+1}];
a1=Table[yLP[x,n],{n,1,(2^10)+1}];
d1sygnal=ListPlot[DownSample[d1] , Joined->True];
a1sygnal=ListPlot[DownSample[a1], Joined->True];
a12n=DownSample[a1];
PrependTo[a12n,0];
d2=Table[yHP[a12n,n],{n,1,2^9+1}];
a2=Table[yLP[a12n,n],{n,1,2^9+1}];
d2sygnal=ListPlot[DownSample[d2], Joined->True];
a2sygnal=ListPlot[DownSample[a2], Joined->True];
a22n=DownSample[a2];
PrependTo[a22n,0];
d3=Table[yHP[a22n,n],{n,1,2^8+1}];
a3=Table[yLP[a22n,n],{n,1,2^8+1}];
d3sygnal=ListPlot[DownSample[d3], Joined->True];
a3sygnal=ListPlot[DownSample[a3], Joined->True];
a32n=DownSample[a3];
PrependTo[a32n,0];
d4=Table[yHP[a32n,n],{n,1,2^7+1}];
a4=Table[yLP[a32n,n],{n,1,2^7+1}];
d4sygnal=ListPlot[DownSample[d4], Joined->True];
a4sygnal=ListPlot[DownSample[a4], Joined->True];

GraphicsGrid[{{sygnal},{a1sygnal,d1sygnal},{a2sygnal,d2sygnal},{a3sygnal,d3sygnal},{a4sygnal,d4sygnal}}]

Przykład 3

f[t_]:=(0.2 Sin[2 Pi 32 t]+Sin[2 Pi 2 t]) UnitStep[-t(t-1)]+(0.4 Sin[2 Pi 32 t]+Sin[2 Pi 2 t]) UnitStep[-(t-2)(t-1)]+(0.6 Sin[2 Pi 32 t]+Sin[2 Pi 2 t]) UnitStep[-(t-2)(t-3)]+(0.8 Sin[2 Pi 32 t]+Sin[2 Pi 2 t]) UnitStep[-(t-3)(t-4)];
Plot[f[t],{t,0,4}, PlotRange-> All]

Clear[x]
x=Table[N[f[n/2^8]], {n,1,2^10}];
sygnal=ListPlot[x,Joined->True];
x=PrependTo[x,0];
d1=Table[yHP[x,n],{n,1,(2^10)+1}];
a1=Table[yLP[x,n],{n,1,(2^10)+1}];
d1sygnal=ListPlot[DownSample[d1] , Joined->True];
a1sygnal=ListPlot[DownSample[a1], Joined->True];
a12n=DownSample[a1];
PrependTo[a12n,0];
d2=Table[yHP[a12n,n],{n,1,2^9+1}];
a2=Table[yLP[a12n,n],{n,1,2^9+1}];
d2sygnal=ListPlot[DownSample[d2], Joined->True];
a2sygnal=ListPlot[DownSample[a2], Joined->True];
a22n=DownSample[a2];
PrependTo[a22n,0];
d3=Table[yHP[a22n,n],{n,1,2^8+1}];
a3=Table[yLP[a22n,n],{n,1,2^8+1}];
d3sygnal=ListPlot[DownSample[d3], Joined->True];
a3sygnal=ListPlot[DownSample[a3], Joined->True];
a32n=DownSample[a3];
PrependTo[a32n,0];
d4=Table[yHP[a32n,n],{n,1,2^7+1}];
a4=Table[yLP[a32n,n],{n,1,2^7+1}];
d4sygnal=ListPlot[DownSample[d4], Joined->True];
a4sygnal=ListPlot[DownSample[a4], Joined->True];

GraphicsGrid[{{sygnal},{a1sygnal,d1sygnal},{a2sygnal,d2sygnal},{a3sygnal,d3sygnal},{a4sygnal,d4sygnal}}]

Pytania na kolokwium

1. Przeprowadzić procedurę transrormaty Haara dla sygnału np. {16,16,16,16,8,8,4,4} - p.p  100% - w pytaniu pojawi się inny sygnał
2. Przeprowadzić procedurę odwrotnej transrormaty Haara dla wyniku z poprzedniego pytania - p.p 100%
3. Jak przebiega procedura transformaty Haara - można narysowac drzewo (diagram) i opisać filtry i oznaczenia - p.r
4. Jak przebiega procedura odwrotnej transformaty Haara - można narysowac drzewo (diagram) i opisać filtry i oznaczenia - p.r
5. Jaka jest interpretacja graficzna filtrów górno i dolnoprzepustowych - p.r

Zadanie na laboratorium

Zadanie nr 12

Należy:
1. Wybrać sobie dowolną funkcję taką, w której występują różne częstotliwości w różnych przedziałach czasowych
3. Spróbkować funkcje - tak, aby było widać wszystkie częstotliwości wchodzące w skład sygnału
4. Dokonać analizy funkcji na czterech poziomach rozdzielczości używając filtrów górno i dolno przepustowych Haara - może być tak jak to zostało przedstawione w trzech przykładach na końcu wykładu.
5. Plus krótki opis nowo poznanej funkcji Mathematici + wlasne przykłady