Wykład nr 6
Filtrowanie przy pomocy TF

Filtrowanie sygnału krok po kroku:
Filtering step by step

Krok 1
Dany jest sygnał: s(t)=Exp[-t^2/10](Sin[2 t]+2 Cos[4 t]+0.4 Sin[t] Cos[50 t]) in<0,2π>. Zanalizujmy tą funkcję w przedziale <0,2π>.
Step 1
We are ginen a signal: s(t)= Exp[-t^2/10](Sin[2 t]+2 Cos[4 t]+0.4 Sin[t] Cos[50 t]) in<0,2π>.

sygnal[t_]:=Exp[-t^2/10](Sin[2 t]+2 Cos[4 t]+0.4 Sin[t] Cos[50 t]);
Plot[sygnal[t],{t,0,7}]

Krok 2
Zastępujemy funkcję ciągłą w czasie jej wartościami zmierzonymi w określonych odstępach czasu, czyli wykonujemy operację próbkowania. Generujemy w ten sposób listę 256 jednakowo odległych punktów sygnału. Wyniki zapisujemy w tablicy.
Step 2
We change the analog function with it's discrete values in certain steps of time, meaning we sample the signal.
We generate this way a list o 256 equally distanced points. Our result we write in a table.

sygnalprob[n_]:=sygnal[2 Pi n/256]
tablica=Table[N[sygnal[2 Pi n/256]],{n,256}];
ListPlot[tablica,PlotStyle->PointSize[0.01]]

Krok 3
Wykonujemy transformatę Fouriera spróbkowanej funkcji. wyniki zapisujemy w tablicy:
Step 3
We take a Fourier Transform of our sampled function. Results we write down in a table:

Fourier[tablica];

Krok 4
Wykonujemy wykresy (części rzeczywistej, urojonej i modułu) otrzymanych współczynników. Analizując wykresy możemy wywnioskować, że zaszumienie ma częstotliwości odpowiadające n≥6.
Step 4
We plot our table (it's imaginary, real part and absolute value). Analyzing the plot we find that the noise has  frequencies bigger than for n=6.

{ListPlot[Re[Fourier[tablica]],Joined->True,PlotRange->All],
ListPlot[Im[Fourier[tablica]],Joined->True,PlotRange->All],
ListPlot[Abs[Fourier[tablica]],Joined->True,PlotRange->All]}

Krok 5
Konstrujemy taki filtr, który przepuści częstotliwości dla n≤6 i n≥251, natomiast wyzeruje pozostałe. Przypomnij sobie zad 4 i własności transformat Fouriera  (OverHat[f])[n-k]=(OverBar[(OverHat[f])])[k].
Step 5
We consider a filter which will pass all the frequencies for n≤6 i n≥ 251 and band all others. Remind Yourself exercise 4 and the symmetry properties: (OverHat[f])[n-k]=(OverBar[(OverHat[f])])[k].

filtr=Table[UnitStep[(n-6)(n-251)],{n,1,256}]

Krok 6
Filtrujemy transformatę Fouriera i wyniki zapisujemy w tablicy.
Step 6
We filter the Fourier Transform. Results we write down in a table:

odfiltrowanatablica=filtr Fourier[tablica]

Krok 7
Dla sprawdzenia filtru rysujemy wykresy otrzymanych współczynników.
Step 7
Just to find out if we chose the right filter we plot the (real part, imaginary part and the absolute value) of the coefficients

GraphicsRow[{ListPlot[Re[odfiltrowanatablica],Joined->True,PlotRange->All],
ListPlot[Im[odfiltrowanatablica],Joined->True,PlotRange->All],
ListPlot[Abs[odfiltrowanatablica],Joined->True,PlotRange->All]}]

Krok 8
Wykonujemy transformatę odwrotną odfiltrowanych danych.
Step 8
We take the inverse transform of the filtered data.

InverseFourier[odfiltrowanatablica];

Krok 9
Wykres odfiltrowanej (wygładzonej) funkcji:
Step 9
This is the plot of the filtered out (smooth) function.

ListPlot[Re[InverseFourier[odfiltrowanatablica]],Joined->True,PlotRange->All]

Clear[tablica, nowatablica, sygnal, sygnalprob, odfiltrowanatablica, filtr]

Filtrowanie sygnału krok po kroku II (przy użyciu splotu):
Filtering step by step II (with use of convolution)

Krok 1
Dany jest sygnał: s(t)=Exp[-t^2/10](Sin[2 t]+2 Cos[4 t]+0.4 Sin[t] Cos[50 t]). t∈<0,2π>.
Step 1
We are ginen a signal: s(t)=Exp[-t^2/10](Sin[2 t]+2 Cos[4 t]+0.4 Sin[t] Cos[50 t]) in<0,2π>.

sygnal[t_]:=Exp[-t^2/10](Sin[2 t]+2 Cos[4 t]+0.4 Sin[t] Cos[50 t])

Plot[sygnal[t],{t,0,7}]

Krok 2
Zastępujemy funkcję ciągłą w czasie jej wartościami zmierzonymi w określonych odstępach czasu, czyli wykonujemy operację próbkowania. 
Generujemy w ten sposób listę  256 jednakowo odległych  punktów sygnału. Wyniki zapisujemy w tablicy.
Step 2
We change the analog function with it's discrete values in certain steps of time, meaning we sample the signal.
We generate this way a list o 256 equally distanced points. Our result we write in a table.

sygnalprob[n_]:=sygnal[2 Pi n/256]
tablica=Table[N[sygnal[2 Pi n/256]],{n,1,256}];
ListPlot[tablica,PlotStyle->PointSize[0.01]]

Krok 3
Wykonujemy transformatę Fouriera spróbkowanej funkcji. Wyniki zapisujemy w tablicy:
Step 3
We take a Fourier Transform of our sampled function. Results we write down in a table:

Fourier[tablica];

Krok 4
Wykonujemy wykresy (części rzeczywistej, nierzeczywistej i modułu) otrzymanych współczynników. Analizujac wykresy możemy wywnioskować, że zaszumienie ma częstotliwości dla n≥6.
Step 4
We plot our table (it’s imaginary, real part and absolute value). Analyzing the plot we find that the noise has  frequencies bigger than for n=6.

GraphicsRow[{ListPlot[Re[Fourier[tablica]],Joined->True,PlotRange->All],
ListPlot[Im[Fourier[tablica]],Joined->True,PlotRange->All],
ListPlot[Abs[Fourier[tablica]],Joined->True,PlotRange->All]}]

Krok 5
Konstrujemy taki filtr, który przepuści częstotliwości dla n≤6 i n≥ 251, natomiast wyzeruje pozostałe. Przypomnij sobie zadanie 4 i własność symetrii transformat Fouriera: (OverHat[f])[n-k]=(OverBar[(OverHat[f])])[k]. 
Step 5
We consider a filter which will pass all the frequencies for n≤6 i n≥ 251 and band all others. Remind Yourself exercise 4 and the symmetry properties: (OverHat[f])[n-k]=(OverBar[(OverHat[f])])[k].

filtr=Table[UnitStep[(n-6)(n-251)],{n,1,256}]

Krok 6
Wykonujemy transformatę odwrotną filtru. 
Step 6
We take the inverse transform of the filter.

InverseFourier[filtr];

Krok 7
Obliczamy splot spróbkowanej funkcji i odwrotnej transformaty Fouriera filtra. W tym celu użyłam ListConvolve[] 
With (Subscript[kernelK, r])and (Subscript[lista, s]),(RefLink[ListConvolve, paclet:ref/ListConvolve])[ker,list]computes(∑_r)(Subscript[K, r])(Subscript[a, s-r]),where the limits of the sum are such that the kernel never overhangs either end of the list.. 
Step 7
We take the convolution of our function and the inverse transform of the filter.

splot=ListConvolve[tablica,InverseFourier[filtr],{2,2}];
Dimensions[splot]
Dimensions[tablica]

Krok 8
Wykonyjemy wykres splotu czyli naszej odfiltrowanej (wygładzonej) funkcji: 
Step 8
We plot the convolution that is our filtered (smooth) function

ListPlot[Re[splot],Joined->True,PlotRange->All]

Odtwarzanie dźwięku przy pomocy programu Mathematica:

Tony dobrze nastrojonego instrumentu mają ściśle określone wysokości, czyli ściśle określone częstotliwości ich tonów podstawowych. Tak np. ton a1 ma częstotliwość podstawową 440 Hz, c1 częstotliwość 261 Hz.

a1[t_]:=Sin[2 Pi 440 t]
Plot[a1[t],{t,0,1}]

Play[a1[t],{t,0,3}]

c[t_]:=Sin[2 Pi 261 t]
Plot[c[t],{t,0,1}]

Play[c[t],{t,0,3}]

Zastępumemy funkcję a1[t]  ciągłą w czasie jej wartościami zmierzonymi w określonych odstępach czasu, czyli wykonujemy operację próbkowania. 
Próbkujemy z cz ęstością 4000. Wyniki zapisujemy w tablicy i odtwarzamy dźwięk.

tablicaa1=Table[N[a1[n/4000]],{n,1,12000}];
ListPlot[tablicaa1,PlotStyle->PointSize[0.01]]

ListPlay[tablicaa1]

Coś się nie zgadza. Wnikliwe ucho stwierdzi zapewne,że otrzymany sygnał jest o oktawę wyższy, czyli jest to druga harmoniczna tonu a o częstotliwości 880 Hz.

a2[t_]:=Sin[2 Pi 880 t]
Plot[a2[t],{t,0,1}]

Play[a2[t],{t,0,3}]

Options[Play]

Skąd ten błąd? 
W Play domyślnie ustawione jest SampleRate->8000.

tablica2a1=Table[N[a1[n/8000]],{n,1,24000}];
ListPlot[tablica2a1,PlotStyle->PointSize[0.01]]

ListPlay[tablica2a1]

Teraz udało się!
Musimy także uważać na zjawisko aliasingu.  Tw. Nyquista mówi o tym, że sygnał możemy odtworzyć z wektora jego wartości dysktretnych nδt jeżeli nie było w nim częstości wyższych niż (1/(2δt)). W przypadku dźwięku a1 ta częstotliwość to 440, zatem powinniśmy próbkować przynajmniej z częstotliwością 880.

k=889;
tablicaa1=Table[N[a1[n/k]],{n,1,2400}];
ListPlot[tablicaa1,PlotStyle->PointSize[0.01]];

ListPlay[tablicaa1,SampleRate->k]

ListPlot[Abs[Fourier[tablicaa1]],Joined->True,PlotRange->All]

Kompresja sygnału

W zależności od rodzaju transmisji stosuje się próbkowanie:
Częstotliwość próbkowania       Zastosowania
8,00 kHz                               Telefonia
16 kHz                                  Komunikacja multimedialna
22,05 kHz                              Komputery osobiste
32,00 kHz                              Cyfrowe radio i telewizja
44,10 kHz                              Płyty CD – Audio
48,00 kHz                              Magnetofony DAT,HDTV
96,00 kHz                              Płyty DVD

Zmiejszenie liczby próbek
Kwantyzacja jest nieodwracalnym odwzorowaniem zmniejszającym dokładność danych przez ograniczenie ich zbioru wartości. 
Procedura:
1. Wejściowy zbiór wartości sygnału dzielony jest na rozłączne przedziały. 
2. Każda wartość wejściowa wypadająca w określonym przedziale jest w wyniku kwantyzacji odwzorowana na jedną wartość wyjściową przypisaną temu przedziałowi, czyli tak zwany poziom reprezentacji. Poziomem reprezentacji może być górna bądź dolna granica przedziału, jednak najczęściej jest nią wartość ze środka przedziału.

k=8000;
tablicaa1=Table[N[a1[n/k]],{n,1,24000}];
n=2
tablicaa1[[1]]
tabl=Table[Round[(2^n) tablicaa1[[m]]]/(2^n),{m,1,24000}];
ListPlot[tabl,PlotStyle->PointSize[0.01]];

ListPlay[tabl,SampleRate->k]

Pytania na kolokwium

1. W jakich krokach przebiega filtrowanie sygnału przy użyciu Transformaty Fouriera - p.p. 
2. W jakich krokach przebiega filtrowanie sygnału przy użyciu splotu - p.r.

Ćwiczenia na laboratorium

Ćwiczenie 9 - 3p.
Należy
1. Wymyślić sobie dowolny sygnał w przedziale <0,2π>, który składa się z wysokich i niskich częstotliwości (można to zrobić dość analogicznie jak w wykładzie).
2. Proszę przeprowadzić procedurę odfiltrowania dla trzech zaprojektowanych przez siebie filtrów (tzn. pobawić się filtrami, możecie np. usunąć wysokie lub tylko niskie częstotliwości lub tylko trochę wysokich...)
3. Proszę krótko opisać wnioski z ćwiczenia
4. Opisać jedną funkcję Mathematici jak w poprzednich ćwiczeniach.