Laboratorium Cyfrowego Przetwarzania Sygnałów
Ćwiczenie 5 – Wieloczęstotliwościowe przetwarzanie sygnałów – interpolacja i decymacja
- 1 -
Wydział Elektryczny
Zakład Automatyki
LABORATORIUM CYFROWEGO PRZETWARZANIA SYGNAŁÓW
Ćwiczenie 5
Wieloczęstotliwościowe przetwarzanie sygnałów – interpolacja i decymacja
1. Cel ćwiczenia
• Zapoznanie się z metodami i efektami cyfrowej zmiany częstotliwości próbkowania sygnałów
dyskretnych.
• Przeprowadzenie dekompozycji i rekonstrukcji podpasmowej sygnału za pomocą
zaprojektowanego banku filtrów pasmowoprzepustowych.
2. Podstawy teoretyczne
Przetwarzanie wieloczęstotliwościowe (multirate processing) jest ważnym obszarem cyfrowego
przetwarzania sygnałów dyskretnych, który nie ma odpowiednika w dziedzinie analogowej. Podstawowe
idee polegają na takiej konwersji sygnału, aby mógł on być przetwarzany z możliwie najmniejszą
częstotliwością, dopasowanie częstotliwości do określonego standardu, pojemności kanału
transmisyjnego lub mocy obliczeniowej procesora DSP. Operacja zmniejszania częstotliwości
próbkowania sygnału (o czynnik całkowity) nazywana jest decymacją. Jeżeli potrzebne jest zwiększenie
częstotliwości przetwarzania, np. w celu dokładniejszej reprezentacji sygnału, stosuje się operację
odwrotną, czyli interpolację. Kombinacja interpolacji i decymacji w połączeniu z odpowiednią filtracją
dolnoprzepustową umożliwia zmianę częstotliwości próbkowania o czynnik ułamkowy.
Oprócz podstawowych operacji interpolacji i decymacji przetwarzanie wieloczęstotliwościowe
obejmuje bardziej zaawansowane tematy, takie jak synteza banków filtrów pasmowoprzepustowych do
rozkładu i rekonstrukcji sygnałów czy rozkład sygnału za pomocą transformacji falkowych (wavelets).
2.1. Interpolacja sygnału
Proces interpolacji o czynnik U polega w zasadzie na estymacji lub rekonstrukcji U-1 wartości
sygnału x(n), próbkowanego z częstotliwością f
sx
, pomiędzy istniejącymi próbkami. Pierwszym
elementem interpolatora (rys. 1) jest ekspander, który pomiędzy każde sąsiednie próbki sygnału
wejściowego wstawia U-1 próbek zerowych (upsampling). Sygnał wyjściowy z ekspandera
/
( ) |
( / ) dla
0,
, 2 ,
( )
0
dla pozostalych
n m U
e
x n
x m U
m
U
U
x m
m
=
=
= ± ±
⎧
= ⎨
⎩
…
(5.1)
ma większą częstotliwość próbkowania f
sy
=U
⋅
f
sx
. W dziedzinie częstotliwości zespolonej X
e
(z)=X(z
U
).
Drugim elementem interpolatora jest filtr dolnoprzepustowy (LP) H
I
(z), który nadaje dodanym
próbkom wartości niezerowe przez interpolację pomiędzy próbkami oryginalnymi generując na wyjściu
sygnał interpolowany y(m). Idealny filtr LP interpolatora ma charakterystykę widmową
dla 0 |
|
(
)
0
dla pozostalych
j
I
U
H e
U
π
Ω
⎧
≤ Ω ≤
⎪
= ⎨
⎪
Ω
⎩
.
(5.2)
Laboratorium Cyfrowego Przetwarzania Sygnałów
Ćwiczenie 5 – Wieloczęstotliwościowe przetwarzanie sygnałów – interpolacja i decymacja
- 2 -
i zerowe opóźnienie grupowe
τ
(
Ω). Wzmocnienie równe U zapewnia takie same amplitudy obwiedni
sygnałów x(n) i y(m) (tylko co U-ta próbka sygnału wejściowego filtra jest niezerowa). Skala
unormowanej pulsacji kątowej
Ω
=
Ω
y
=
Ω
x
/U, ponieważ
2
2
,
y
x
x
y
sx
sy
f
f
f
f
π
π
Ω =
Ω =
.
(5.3)
Charakterystyka H
I
(e
j
Ω
) odpowiada nieskończonej i nieprzyczynowej (tzn. o wartościach niezerowych
dla czasu m<0) charakterystyce impulsowej:
sin(
/ )
dla
0
( )
/
1
dla
0
I
m
U
m
h m
m
U
m
π
π
⎧
≠
⎪
= ⎨
⎪
=
⎩
.
(5.4)
gdzie sin(
πx)/(πx)=sinc(x). Podane warunki mogą być w dobrym przybliżeniu spełnione tylko przez
nieprzyczynowy filtr odpowiednio wysokiego rzędu. Filtr interpolatora projektuje się zwykle zadając
pulsację końca -3 dB pasma przepustowego
Ω
pass
=
απ/U, gdzie 0<α≤1 jest wypełnieniem pasma Nyquista
przez ograniczone widmo sygnału wejściowego.
W dziedzinie czasu wyjście interpolatora jest splotem:
( )
( )
( )
(
) ( )
(
) ( / )
I
e
I
e
I
k
k
y m
h m
x m
h m k x k
h m k x k U
∞
∞
=−∞
=−∞
=
∗
=
−
=
−
∑
∑
,
(5.5)
a w dziedzinie transformaty Z:
( )
( )
( )
( )
( )
U
I
e
I
Y z
H z X z
H z X z
=
⋅
=
⋅
. (5.6)
Rys. 1. Schemat interpolacji o czynnik U
Rys. 2. Przebiegi sygnałów w procesie interpolacji o czynnik U=3: a) sygnał wejściowy,
b) sygnał wyjściowy z ekspandera (upsampling), c) sygnał wyjściowy po filtracji wygładzającej LP
Ekspander
↑ U
Filtr LP
H
I
(z)
x(n)
y(m)
x
e
(m)
f
sy
f
sx
f
sy
=U
⋅
f
sx
x(n)
x
e
(m)
y(m)
m
m
n
0
c)
b)
a)
Laboratorium Cyfrowego Przetwarzania Sygnałów
Ćwiczenie 5 – Wieloczęstotliwościowe przetwarzanie sygnałów – interpolacja i decymacja
- 3 -
h
I
(0-n)
x
e
(n)
h
I
(1-n)
h
I
(2-n)
h
I
(3-n)
h
I
(4-n)
Rys. 3. Widma sygnałów w procesie interpolacji o czynnik U=3: a) sygnał wejściowy, b) sygnał z
ekspandera (upsampling) i charakterystyka |H
I
(f)| filtra LP, c) sygnał wyjściowy interpolatora po filtracji
Rys. 4. Interpolacja jako proces filtracji liniowej
– splotu xe(n) z charakterystyką impulsową filtra
hI(n). Odwrócona charakterystyka hI(-n) jest
przesuwana na tle sygnału xe(n). Dla każdego
przesunięcia y(n) to suma iloczynów
niezerowych próbek i współczynników filtra.
Zwrócić uwagę, kiedy zera xe(n) nakładają się na
hI(k-n) i kiedy zera hI(k-n) nakładają sę na xe(n).
0
3f
sx
|X(f)|
f
2f
sx
f
sx
-f
sx
-2f
sx
-3f
sx
0
f
sy
=
3f
sx
|X
e
(f)|
f
-f
sy
=-3f
sx
0
f
sy
|Y(f)|
f
-f
sy
|H
I
(f)|
½f
sx
-½f
sx
a)
b)
c)
Laboratorium Cyfrowego Przetwarzania Sygnałów
Ćwiczenie 5 – Wieloczęstotliwościowe przetwarzanie sygnałów – interpolacja i decymacja
- 4 -
W wielu przypadkach wystarczająco dobre wyniki daje zastosowanie prostej interpolacji liniowej lub
ekstrapolacji zerowego rzędu ZOH (zero order hold; stosowana w przetwornikach C/A).
W przypadku interpolacji ZOH każda próbka sygnału x(n) jest po prostu powtarzana U razy:
(0)
0,1, ,
1
( )
,
1, , 2
1
( )
(2 )
2 , 2
1, ,3
1
e
e
e
x
m
U
x U
m U U
U
y m
x
U
m
U U
U
=
−
⎧
⎪
=
+
−
⎪
= ⎨
=
+
−
⎪
⎪⎩
…
…
…
,
(5.7)
co daje schodkowy sygnał wyjściowy i odpowiada filtrowi o przyczynowej charakterystyce impulsowej:
( )
( )
(
1)
(
(
1))
ZOH
h
m
m
m
m
U
δ
δ
δ
=
+
− + +
−
−
…
. (5.8)
Interpolacja liniowa jest realizowana jako operacja nieprzyczynowa, której może odpowiadać
charakterystyka impulsowa:
1 | | /
dla | |
1
( )
0
dla pozostalych
lin
m U
m U
h m
m
−
≤ −
⎧
= ⎨
⎩
(5.9)
(interpolacja tylko pomiędzy dwiema próbkami sygnału oryginalnego).
2.2. Decymacja sygnału
Decymacja jest operacją odwrotną do interpolacji. Decymacja o czynnik D polega na D-krotnym
zmniejszeniu częstotliwości próbkowania sygnału. Pierwszym elementem decymatora jest filtr
dolnoprzepustowy H
D
(z) (rys. 5), którego zadaniem jest ograniczenie pasma sygnału wejściowego x(n) do
zakresu D razy mniejszego. Charakterystyka widmowa idealnego filtra LP decymatora
1 dla 0 |
|
(
)
0
dla pozostalych
j
D
H e
D
π
Ω
⎧
≤ Ω ≤
⎪
= ⎨
⎪
Ω
⎩
.
(5.10)
gdzie skala
Ω
=
Ω
x
=
Ω
y
⋅
D.
W praktyce filtr decymatora projektuje się zadając pulsację początku pasma
zaporowego
Ω
stop
=
π/D.
Filtracja LP zapobiega aliasingowi po następującej dalej operacji redukcji częstotliwości próbkowania
polegającej na pozostawieniu co D-tej próbki sygnału przefiltrowanego x
f
(n): f
sy
=f
sx
/D. Zmniejszenie
szybkości próbkowania ogranicza na wyjściu reduktora D razy pasmo Nyquista, czyli zakres
częstotliwości możliwych do reprezentowania w sygnale wyjściowym y(m). Sygnał wyjściowy:
( )
( ) |
( )
( ) |
( ) (
)
f
n mD
D
n mD
D
k
y m
x n
h n
x n
h k x mD k
∞
=
=
=−∞
=
=
∗
=
⋅
−
∑
,
(5.11)
gdzie h
D
(n) jest charakterystyką impulsową filtra H
D
(z). W przypadku idealnej filtracji LP w dziedzinie
transformaty Z:
1/
1
( )
(
)
D
Y z
X z
D
=
.
(5.12)
Filtracja LP sygnału w decymatorze może być realizowana w sposób nieprzyczynowy (off-line),
podobnie jak w interpolatorze, lub przyczynowy (on-line) tylko na bazie próbek dostępnych do momentu
aktualnego, co wprowadza przesunięcie fazowe. W przypadku filtra SOI rzędu N=2L o symetrycznej
odpowiedzi h
D
(n), charakterystyki przyczynowa i nieprzyczynowa są wzajemnie przesunięte:
h
nieprzycz
(n) = h
przycz
(n+L),
H
nieprzycz
(z) = z
L
H
przycz
(z).
(5.13)
Ze względu na oszczędności obliczeniowe i łagodniejsze wymagania co do filtrów LP decymacja dla
dużych czynników D realizowana jest jako kombinacja dwóch lub kilku decymacji z mniejszymi
wartościami D.
Laboratorium Cyfrowego Przetwarzania Sygnałów
Ćwiczenie 5 – Wieloczęstotliwościowe przetwarzanie sygnałów – interpolacja i decymacja
- 5 -
Rys. 5. Schemat decymacji o czynnik D
Rys. 6. Przebiegi sygnałów w procesie decymacji o czynnik D=3: a) sygnał wejściowy, b) sygnał po
filtracji LP, c) sygnał wyjściowy po redukcji częstotliwości próbkowania (downsampling)
Rys. 7. Widma sygnałów w procesie decymacji o czynnik D=3: a) sygnał wejściowy i charakterystyka
|H
D
(f)| filtra LP ograniczającego pasmo, b) sygnał wyjściowy decymatora po filtracji i redukcji
częstotliwości próbkowania
2.3. Połączenie interpolacji i decymacji
Połączenie interpolacji i decymacji umożliwia zmianę szybkości próbkowania o czynnik wymierny
U/D (rys. 8). Oba filtry LP połączone szeregowo pracują na częstotliwości wyjściowej ekspandera, co
umożliwia zastąpienie ich pojedynczym filtrem o transmitancji H
I
(z)H
D
(z).
y(m)
x
f
(n)
m
c)
n
a)
n
b)
x(n)
0
0
f
sx
|
X(f)|
f
-f
sx
0
3f
sy
=f
sx
|Y(f)|
f
½f
sx
-½f
sx
2f
sy
f
sy
-f
sy
-2f
sy
-3f
sy
=-f
sx
|H
D
(f)|
½f
sy
-½f
sy
B
a)
b)
Filtr LP
H
D
(z)
x(n)
y(m)
x
f
(n)
f
sx
f
sx
f
sy
=f
sx
/D
Reduktor
↓ D
Laboratorium Cyfrowego Przetwarzania Sygnałów
Ćwiczenie 5 – Wieloczęstotliwościowe przetwarzanie sygnałów – interpolacja i decymacja
- 6 -
Rys. 8. Przetworzenie częstotliwości próbkowania o czynnik U/D=4/3: a) kombinacja interpolacji i
decymacji, b) zastąpienie dwóch filtrów LP jednym o transmitancji H
I
(z)H
D
(z), c) widmo sygnału
wejściowego, d) widmo sygnału wyjściowego filtra po interpolacji, e) ostateczne widmo sygnału
wyjściowego po decymacji. Oznaczenia:
f
sx
=f
s_old
f
sy
=f
s_new
2.4. Próbkowanie sygnału pasmowego
Standardowe próbkowanie dolnopasmowe sygnału o pasmowym widmie o szerokości B wokół
częstotliwości nośnej (środkowej) f
c
(rys. 9a) teoretycznie wymaga dużej częstotliwości próbkowania
f
s
≥f
c
+B/2. Zastosowanie redukcji szybkości próbkowania (downsamplingu, decymacji bez filtracji)
umożliwia uzyskanie powielenia oryginalnego widma wokół częstotliwości zerowej wskutek aliasingu,
bez straty informacji (rys. 9b). Jest to znacznie oszczędniejsza metoda reprezentacji sygnałów
pasmowych stosowana w dekompozycji sygnału za pomocą banków filtrów pasmowoprzepustowych.
W celu uniknięcia nakładania się powieleń widma sygnału zredukowana częstotliwość próbkowania
musi spełniać warunek
1
2
/ 2
/ 2
c
c
s
f
B
f
B
f
k
k
+
−
≤
≤
+
, f
s
≥2B
(5.14)
gdzie k jest liczbą naturalną.
Odtworzenie sygnału oryginalnego jest możliwe po przeprowadzeniu interpolacji pasmowo-
przepustowej: powrocie do większej (pierwotnej) częstotliwości próbkowania i filtracji BP wydzielającej
tylko oryginalne pasmo sygnału.
Laboratorium Cyfrowego Przetwarzania Sygnałów
Ćwiczenie 5 – Wieloczęstotliwościowe przetwarzanie sygnałów – interpolacja i decymacja
- 7 -
Rys. 9. Próbkowanie sygnału pasmowego: a) widmo oryginalnego sygnału ciągłego, b) powielenia
widma sygnału próbkowanego, B – szerokość pasma sygnału, f
c
– częstotliwość nośna (centralna),
f
s
– częstotliwość próbkowania
2.5. Podpasmowa dekompozycja i rekonstrukcja sygnału za pomocą banku filtrów
Przetwarzanie sygnału w podpasmach częstotliwości jest stosowane w wielu praktycznych
aplikacjach, m.in. algorytmach kompresji audio (MPEG-1,2 z poziomem mp3 oraz MPEG-4) i obrazów
(JPEG2000). Jedną z metod dekompozycji sygnału na składowe pasmowe jest zastosowanie banku
(zespołu) filtrów BP pokrywających zakres widma sygnału, a następnie decymacja składowych sygnału w
podpasmach, co umożliwia np. kompresję sygnału (w zależności od wariancji składowych sygnału w
podpasmach reprezentuje się je w formie zawierającej mniej lub więcej bitów).
M-kanałowy bank filtrów do podpasmowej dekompozycji sygnału x(n) jest przedstawiony na rys. 10a,
a charakterystyki amplitudowo-częstotliwościowe jego filtrów BP otrzymane w wyniku rzeczywistej
modulacji kosinusowej jednego prototypu - na rys. 11. Przy jednakowych szerokościach pasm wszystkie
współczynniki reduktorów i ekspanderów są jednakowe N
i
=N, i=0,1,...,M-1. Kiedy liczba pasm jest
równa współczynnikowi redukcji, N=M, mamy do czynienia z próbkowaniem krytycznym, kiedy N>M
występuje podpróbkowanie, a kiedy N<M występuje nadpróbkowanie. Odpowiednie umiejscowienie
charakterystyk H
k
(e
j
Ω
) filtrów na osi częstotliwości uzyskuje się w tym przypadku jako efekt modulacji
kosinusowej (w dziedzinie czasu) charakterystyki impulsowej p(n), 0≤n≤L-1, filtra prototypowego LP o
pulsacji granicznej
Ω
p
=
π/(2M), który należy zaprojektować. Charakterystyki impulsowe poszczególnych
filtrów BP banku analizy:
1
1
( ) 2 ( ) cos
( 1)
,
0,1, ,
1
2
2
4
k
k
L
h n
p n
k
n
k
M
M
⎡ π
−
π⎤
⎛
⎞⎛
⎞
=
⋅
+
−
+ −
=
−
⎜
⎟⎜
⎟
⎢
⎥
⎝
⎠⎝
⎠
⎣
⎦
…
(5.15)
Modulacja kosinusowa charakterystyki impulsowej prototypu powoduje przesuwanie charakterystyki
widmowej do odpowiedniego położenia na osi pulsacji
Ω, ponieważ:
0
0
0
0
(
)
(
)
0
1
( ) cos
( )
(
)
(
)
2
2
j
n
j
n
j
j
e
e
h n
n h n
H e
H e
− Ω
Ω
Ω−Ω
Ω+Ω
+
⎡
⎤
⋅
Ω =
↔
+
⎣
⎦
(5.16)
Odtworzenie ˆ
( )
x n sygnału oryginalnego uzyskuje się po interpolacji, czyli przepuszczeniu
zdecydowanych składowych
y
k
(m) przez ekspandery i filtry BP rekonstrukcji (syntezy, rys. 10b).
Jeżeli prototyp
p(n) jest symetryczny, to odpowiedzi impulsowe filtrów syntezy g
k
(n) są odwróconymi w
czasie odpowiedziami h
k
(n).
|X
a
(f)|
|X(f)|
f
f
B
f
c
f
c
-f
c
-f
c
k
⋅
f
s
k
⋅
f
s
½B
-½B
0
a)
b)
f
s
f
s
Laboratorium Cyfrowego Przetwarzania Sygnałów
Ćwiczenie 5 – Wieloczęstotliwościowe przetwarzanie sygnałów – interpolacja i decymacja
- 8 -
a)
b)
Rys. 10. M-kanałowy bank (zespół) filtrów podpasmowej analizy (a) i syntezy (b) sygnału
(subband decomposition and reconstruction)
Rys. 11. Charakterystyki amplitudowo-częstotliwościowe M-pasmowego banku filtrów otrzymane w
wyniku rzeczywistej modulacji kosinusowej charakterystyki impulsowej p(n) prototypu LP
Zależności pomiędzy transformatami Z sygnałów z rys.10 są następujące:
- po filtracji :
( )
( )
( ),
0,1, ,
1
k
k
U z
H z X z
k
M
=
⋅
=
−
…
(5.17)
- po redukcji:
1
1/
2 /
0
1
( )
(
),
N
l
N
j
N
k
k
N
N
l
Y z
U w z
w
e
N
−
−
π
=
=
=
∑
(5.18)
Y
k
(z) jest sumą powielonych poprzesuwanych (zmodulowanych przez
w
N
l
) widm U
k
(z
1/N
).
- wyjścia ekspanderów:
( )
(
)
N
k
k
V z
Y z
=
(5.19)
- zrekonstruowane składowe:
ˆ ( )
( )
( )
k
k
k
X z
G z V z
=
⋅
(5.20)
Sygnał zrekonstruowany:
1
1
1
1
0
0
0
0
1
1
1
0
0
0
1
ˆ
ˆ
( )
( )
( ) (
)
( )
(
) (
)
1
(
)
(
)
( )
(
) ( )
M
M
M
N
N
l
l
k
k
k
k
k
N
N
k
k
k
l
N
M
N
l
l
l
N
k
N
k
N
l
l
k
l
X z
X z
G z Y z
G z
H w z X w z
N
X w z
H w z G z
X w z A z
N
−
−
−
−
=
=
=
=
−
−
−
=
=
=
⎛
⎞
=
=
=
=
⎜
⎟
⎝
⎠
⎛
⎞
=
=
⎜
⎟
⎝
⎠
∑
∑
∑
∑
∑
∑
∑
(5.21)
gdzie:
1
0
1
( )
(
)
( )
M
l
l
k
N
k
k
A z
H w z G z
N
−
=
=
∑
(5.22)
0
Ω
π/M
|H(e
j
Ω
)|
H
0
(e
j
Ω
)
H
1
(e
j
Ω
)
H
M-1
(e
j
Ω
)
−π/M
2π/M
π
−π
−2π/M
H
0
(e
j
Ω
)
H
1
(e
j
Ω
)
H
M-1
(e
j
Ω
)
...
...
x(n)
H
0
(z)
y
0
(m)
u
0
(n)
↓ N
1
H
1
(z)
y
1
(m)
u
1
(n)
↓ N
2
H
M-1
(z)
y
M-1
(m)
u
M-1
(n)
↓
N
M-1
filtry analizy
reduktory
filtry syntezy
G
0
(z)
0
ˆ ( )
x n
v
0
(n)
↑ N
1
y
0
(m)
G
1
(z)
1
ˆ ( )
x n
v
1
(n)
↑ N
2
y
1
(m)
G
M-1
(z)
1
ˆ
( )
M
x
n
−
v
M-1
(n)
↑N
M-1
y
M-1
(m)
ˆ( )
x n
ekspandery
Laboratorium Cyfrowego Przetwarzania Sygnałów
Ćwiczenie 5 – Wieloczęstotliwościowe przetwarzanie sygnałów – interpolacja i decymacja
- 9 -
Jeżeli spełnione będą warunki:
0
1
0
0
1
( )
( )
( )
( ) 0,
1, 2, ,
1
M
n
k
k
k
l
A z
H z G z
z
N
A z
l
N
−
−
=
=
= α ⋅
=
=
−
∑
…
,
(5.23)
to odtworzenie
0
ˆ( )
(
)
x n
x n n
= α ⋅
−
jest przeskalowanym i opóźnionym, ale niezniekształconym sygnałem
oryginalnym, a składowe zmodulowane są wytłumione. Dlatego (5.23) są warunkami dokładnej
rekonstrukcji (perfect reconstruction – PR), które powinny być spełnione, przynajmniej w przybliżeniu,
przez zaprojektowany bank filtrów. Niezerowe wartości A
l
(z) dla l≠0, spowodowane nieidealnymi
charakterystykami filtrów, powodują tzw. przecieki widmowe składowych pomiędzy podpasmami.
W algorytmach podpasmowego kodowania (kompresji) sygnałów, składowe podpasmowe
y
k
(m)
po
redukcji są jeszcze kwantowane co do wartości. Liczba poziomów kwantowania R
k
(czyli bitów
reprezentacji) w podpaśmie jest tym większa, im większa jest wariancja var(y
k
) składowej. Ponieważ
całkowita liczba bitów na zakodowaną próbkę R=
Σ
R
k
jest ograniczona i zależy np. od założonego bitrate
I=R
⋅
f
s
(bitów/s, f
s
– częstotliwość próbkowania), podpasma o najmniejszej wariancji mogą w ogóle nie
być reprezentowane w zakodowanym sygnale.
Literatura
1. Lyons R.G.: „Wprowadzenie do cyfrowego przetwarzania sygnałów”, WKŁ, 1999.
2. Zieliński T.P.: „Cyfrowe przetwarzanie sygnałów. Od teorii do zastosowań”, WKŁ, 2005.
3. Crochiere R.E., Rabiner L.R.: „Mltirate Signal Processing”, Prentice Hall, Englewood Cliffs, NJ, 1983.
Laboratorium Cyfrowego Przetwarzania Sygnałów
Ćwiczenie 5 – Wieloczęstotliwościowe przetwarzanie sygnałów – interpolacja i decymacja
- 10 -
3. Obliczenia komputerowe - zadania do wykonania
Do obliczeń wykorzystywane będą funkcje biblioteki Matlaba o nazwie Signal Processing Toolbox.
3.1. Interpolacja sygnału
A. Interpolacja nieprzyczynowa LP o ograniczonym widmie.
• Wygenerować sygnał wejściowy zawierający dwie składowe harmoniczne:
1
1
2
2
1
1
1
2
2
2
( )
cos
cos
1.0,
,
0.4,
x n
A
n A
n
A
f
A
f
=
Ω +
Ω
=
Ω = π
=
Ω = π
(częstotliwości f unormowane do częstotliwości Nyquista, tzn. f=1 odpowiada f
Nyq
=f
s
/2) za pomocą
funkcji:
>> x = c5_signal(512,[1 0.4],[0.1 0.9]);
% N=512 punktów
gdzie: [A1 A2]=[1 0.4], [f1 f2]=[0.1 0.9].
• Przeprowadzić interpolację sygnału (w skrypcie korzysta się z funkcji interp Matlaba):
>> nrys=[0 100];
% zakres osi OX przebiegów czasowych, domyślnie [0,100]
>> U=3;
% współczynnik rozszerzenia
>> L=4; alpha=0.4;
% rząd filtra interpolacyjnego N=2*L*U
Parametr
α jest zakładaną częstotliwością odcięcia sygnału wejściowego unormowaną do zakresu 0< α
≤1. Im większa jest jego wartość, tym szersze jest pasmo przepustowe i bardziej strome pasmo
przejściowe projektowanego w funkcji interp filtra interpolatora.
>> [xup,y,nup,b] = c5_interp(x, U, nrys, L, alpha);
% obliczenia bez wykresów
>> c5_interp(x, U, nrys, L, alpha);
% rysowanie wykresów
gdzie: xup – sygnał z ekspandera, y – sygnał wyjściowy interpolatora, nup – indeksy po rozszerzeniu,
b=[b
0
, b
1
, ..., b
N
]
– współczynniki (odpowiedź impulsowa) filtra SOI H(z)=b(z).
• Zaobserwować widma DFT sygnałów w procesie interpolacji:
>> c5_interpwid(x, xup, y ,U);
• Wyznaczyć i zaobserwować charakterystykę impulsową oraz amplitudowo-częstotliwościową filtra
interpolatora i porównać je z odpowiednimi charakterystykami filtra tego samego rzędu opisanego
wzorem (5.4)
>> c5_interpfilt(b,U);
• Powtórzyć powyższe obliczenia dla U=3 i 4 kombinacji parametrów: L=4, 8; alpha=0.4,
0.95
.
&
Zaobserwować różnice charakterystyk porównywanych filtrów. W których punktach charakterystyki
impulsowe pokrywają się i dlaczego?
&
Dla jakiego
α filtr interpolatora praktycznie pokrywa się z filtrem opisanym wzorem (5.4)? Zwrócić
uwagę na tłumienie składowej f
2
=0.9 w zależności od wartości
α.
B. Interpolacja nieprzyczynowa liniowa pomiędzy sąsiednimi próbkami – charakterystyka impulsowa
filtra interpolatora jest opisana wzorem (5.9).
• Sygnał wejściowy pozostaje taki sam jak w pkt.A (powinien być zapisany w zmiennej x w pamięci).
Przeprowadzić obliczenia dla U=5 i U=8:
>> [xup,y,nup,b] = c5_interplin(x, U);
% obliczenia, bez wykresów
>> c5_interplin(x, U);
% rysowanie wykresów
>> c5_interpwid(x, xup, y ,U);
>> c5_interpfilt(b,U);
&
Porównać wyniki z uzyskanymi w pkt.A.
C. Interpolacja sygnału dźwiękowego (*** sygnał na głośniki ***).
• Wczytać plik dźwiękowy w formacie .wav o podanej nazwie, np.:
>> [x,fs,nbits] = c5_sound('fire');
% tylko konwersja do zmiennej x
>>
c5_sound('fire');
% dźwięk i wykres
Laboratorium Cyfrowego Przetwarzania Sygnałów
Ćwiczenie 5 – Wieloczęstotliwościowe przetwarzanie sygnałów – interpolacja i decymacja
- 11 -
Liczba próbek ns, częstotliwość próbkowania fs i liczba bitów próbki nbits są podawane w oknie
komend Matlaba.
• Przeprowadzić interpolację sygnału:
>> nrys=[2000 2200];
% zakres osi OX przebiegów czasowych
>> U=2;
% współczynnik rozszerzenia
>> [xup,y,nup,b] = c5_interp(x, U, nrys);
% obliczenia, bez wykresów
>> c5_interp(x, U, nrys);
% rysowanie wykresów
>> c5_interpwid(x, xup, y ,U);
% widma częstotliwościowe
>> sound(y,U*fs)
% odegranie sygnału po interpolacji
Analogiczna funkcja wavplay(s,fs) używa urządzenia audio PC.
&
Czy słyszalne są różnice dźwięku przed i po interpolacji?
3.2. Decymacja sygnału
A. Decymacja bez przesunięcia fazowego (SOI ze zwierciadlanym uzupełnianiem próbek wejściowych w
celu uniknięcia przebiegów przejściowych lub nieprzyczynowa NOI).
• Wygenerować sygnał wejściowy zawierający dwie składowe harmoniczne:
>> x = c5_signal(512,[1 0.4],[0.04 0.4]);
% N=512 punktów
gdzie: [A1 A2]=[1 0.4], [f1 f2]=[0.04 0.4].
• Przeprowadzić decymację sygnału (w skrypcie korzysta się z funkcji decimate Matlaba, domyślnie
filtracja SOI rzędu N=30):
>> nrys=[0 100];
% zakres osi OX przebiegów czasowych
>> D=2;
% współczynnik redukcji
>> [xdown,y,ndown]=c5_dec(x,D,nrys);
% obliczenia, bez wykresów
>> c5_dec(x,D,nrys);
% rysowanie wykresów
gdzie: xdown – sygnał z reduktora bez filtracji LP, y – sygnał decymowany (z filtracją LP przed
reduktorem), ndown – indeksy po redukcji.
&
Zwrócić uwagę na brak przesunięcia fazowego sygnału po decymacji.
• Zaobserwować widma DFT sygnałów w procesie decymacji:
>> c5_decwid(x, xdown, y ,D);
&
Zwrócić uwagę, gdzie występuje skala częstotliwości
Ω
x
przed i
Ω
y
po redukcji.
• Powtórzyć obliczenia dla decymacji o czynnik 3:
>> D=3;
>> [xdown,y,ndown]=c5_dec(x,D,nrys);
>> c5_dec(x,D,nrys);
% rysowanie wykresów
>> c5_decwid(x, xdown, y ,D);
% widma częstotliwościowe
&
Wyjaśnić różnice przebiegów i widm sygnału zredukowanego (downsampled) bez filtracji LP w
porównaniu z sygnałem decymowanym z filtracją, biorąc pod uwagę D-krotne zawężenie pasma Nyquista
w wyniku redukcji. Co dzieje się ze składową sygnału wejściowego o większej częstotliwości? Na
wykresach widm wskazać efekt aliasingu.
B. Decymacja przyczynowa z przesunięciem fazowym.
• Przeprowadzić obliczenia dla decymacji o czynnik D=3 za pomocą funkcji:
>> c5_decprzycz(x, D, nrys, L);
dla L=8 i L=15, gdzie: N=2L – rząd przyczynowego filtra SOI decymatora o charakterystyce
impulsowej
(
)
( )
sinc
FIR
n L
h
n
D
D
α
α −
⎛
⎞
=
⎜
⎟
⎝
⎠
, (przyjęto
α=0.9).
Zarejestrować charakterystyki filtra.
&
Zwrócić uwagę na stan przejściowy i zależność jego czasu trwania (warunki początkowe są zerowe)
oraz przesunięcia fazowego w stanie ustalonym od L i wyjaśnić ten efekt.
C. Decymacja sygnału dźwiękowego (*** sygnał na głośniki ***).
• Wczytać plik dźwiękowy w formacie .wav o podanej nazwie, np.:
Laboratorium Cyfrowego Przetwarzania Sygnałów
Ćwiczenie 5 – Wieloczęstotliwościowe przetwarzanie sygnałów – interpolacja i decymacja
- 12 -
>> [x,fs,nbits] = c5_sound('fire');
% tylko konwersja do zmiennej x
>>
c5_sound('fire');
% dźwięk i wykres
• Przeprowadzić decymację sygnału:
>> nrys=[2000 2200];
% zakres osi OX przebiegów czasowych
>> D=2;
% współczynnik redukcji
>> [xdown,y,ndown]=c5_dec(x,D,nrys);
>> c5_dec(x,D,nrys);
% rysowanie wykresów
>> c5_decwid(x, xdown, y ,D);
>> sound(y,fs/D)
% odegranie sygnału po decymacji
• Powtórzyć eksperyment dla D=5.
&
Czy słyszalne są różnice dźwięku przed i po decymacji? Porównać zawartość widmową sygnału
wejściowego i wyjściowego.
• Przeprowadzić obliczenia dla decymacji z filtrem nieprzyczynowym i D=5.
>> [y,ndown]=c5_decprzycz(x, D, nrys);
% obliczenia, domyślnie L=15
>> c5_decprzycz(x, D, nrys);
% rysowanie wykresów
>> c5_decwid(x, xdown, y ,D);
D. Decymacja ze zwiększeniem rozdzielczości kwantowania.
• Wygenerować sygnał harmoniczny (jedna składowa) skwantowany do poziomów +1/-1 (1 bit):
zawierający dwie składowe harmoniczne:
>> x = sign(c5_signal(512,[1 0],[0.02 0]));
% funkcja signum
• Przeprowadzić decymację z przyczynowym filtrem uśredniającym o liniowej odpowiedzi impulsowej:
>> nrys=[0 200];
% zakres osi OX przebiegów czasowych
>> D=5;
% współczynnik redukcji
>> c5_decusred(x, D, nrys);
&
Jak sygnał wyjściowy odzwierciedla oryginalny sygnał harmoniczny? Wyjaśnić przebieg po
decymacji pamiętając, że filtr LP decymatora ma właściwości uśredniające.
3.3. Dekompozycja sygnału za pomocą banku filtrów pasmowoprzepustowych BP
A. Projektowanie filtra prototypowego LP SOI.
• Wprowadzić dane projektowe i zaprojektować prototyp:
>> M = 4;
% liczba kanałów projektowanego baku filtrów (decyduje o szerokości pasma
przepustowego
prototypu)
>> L = 201;
% dlugosc odpowiedzi impulsowej p(n) prototypu (= rząd filtra +1)
>> rs =96;
% min. tłumienie w paśmie zaporowym w dB (stopband rippling)
>> p = c5_prototyp(M,L,rs);
% tylko obliczenia
>> c5_prototyp(M,L,rs);
% wykresy
&
Zwrócić uwagę na zależność charakterystyki amplitudowej od M i rs (można powtórzyć obliczenia
dla innych wartości tych parametrów, ale pamiętać, że zmienia to p).
B. Projektowanie banku filtrów analizy (dekompozycji) i syntezy (rekonstrukcji) sygnału.
• Wyznaczenie macierzy h (wymiar MxL) odpowiedzi impulsowych filtrów analizy i g - filtrów
syntezy (k-ty wiersze macierzy jest charakterystyką filtra h
k
(n) (g
k
(n)), k=0,1,...,M-1) metodą
modulacji kosinusowej:
>> [h,g] = c5_bank1(M,p);
% tylko obliczenia
>> c5_bank1(M,p);
% wykresy
&
Przeanalizować skrypt (plik) c5_bank1.m. Zwrócić uwagę i wyjaśnić wstępne skalowanie
odpowiedzi impulsowej p(n) prototypu i różnice (na wykresach) zmodulowanych odpowiedzi
impulsowych filtrów analizy i syntezy. Czy są to filtry o liniowej fazie?
C. Przeprowadzić analizę (rozkład na M składowych pasmowych) sygnału dźwiękowego za pomocą
zaprojektowanych filtrów analizy.
• Podać nazwę pliku .wav z sygnałem do przetwarzania i wczytać sygnał do zmiennej:
>> wavfile='szum';
% sygnał z pliku szum.wav
Laboratorium Cyfrowego Przetwarzania Sygnałów
Ćwiczenie 5 – Wieloczęstotliwościowe przetwarzanie sygnałów – interpolacja i decymacja
- 13 -
>> [x,fs,nbits] = c5_sound(wavfile,nrys);
>> c5_sound(wavfile,nrys);
% dźwięk i wykres (***sygnał na głośniki***)
&
Zanotować parametry sygnału podane w oknie komend Matlaba.
• Określić lupę czasową rysowanych dalej przebiegów i obliczeń DFT:
>> nrys=[2000, 2500];
% wykresy i DFT w zakresie od próbki 2000 do 2500
• Przeprowadzić dekompozycję sygnału za pomocą filtrów BP h. Składowe y
k
(n), k=0,1,...,M-1, są
zapamiętywane w wierszach macierzy y:
>> [y,vary,varx] = c5_anal(h,x,fs,nrys);
% obliczenia
Odsłuchać sygnały składowych (***sygnał na głośniki***) naciskając dowolny klawisz:
>> c5_anal(h,x,fs,nrys);
% wykreślane są dwie pierwsze i ostatnia składowa: k=0,1,M-1
&
Podawane wariancje vary określają moc (a więc i poziom) składowych w podpasmach, co jest
wykorzystywane do alokacji bitów na podpasma przy kwantowaniu.
D. Kwantowanie poziomów składowych w podpasmach.
• Określić liczbę bitów R na próbkę sygnału po skwantowaniu (bitrate I=R
⋅
f
s
bitów/s) i wyznaczyć
alokację bitów na podpasma (zależną od wariancji vary):
>> R=8;
% do rozdziału jest R bitów na M pasm
>> Rk = c5_bitalloc(R,vary);
% wyliczone wartości Rk są podawane w oknie komend
• Przeprowadzić redukcję (downsampling) częstotliwości próbkowania składowych pasmowych o
czynnik M (powoduje to powielenie widm pasmowych):
>> yd=c5_analdown(y,M);
% yd – macierz składowych po downsamplingu
• Przeprowadzić kwantowanie zdecymowanych składowych pasmowych. Liczba poziomów składowej
po skwantowaniu jest przez liczbę bitów przypadającą na dane pasmo:
>> yq=c5_analquant(yd,fs,M,Rk,nrys);
% obliczenia, yq – składowe skwantowane
>> c5_analquant(yd,fs,M,Rk,nrys);
% wykresy
&
Sprawdzić, czy liczba poziomów na wykresach zgadza się z liczbą bitów przydzieloną na pasmo.
Porównać przebiegi czasowe i widma zarejestrowane przed i po redukcji częstotliwości próbkowania.
W jakim stopniu kwantowanie wpływa na widma składowych? Przeanalizować plik skryptowy
c5_analquant.m
pod kątem realizacji operacji kwantowania poziomów składowych.
3.4. Rekonstrukcja sygnału za pomocą banku filtrów
Uwaga: W tym punkcie wykorzystywane są zmienne wyliczone w pkt.3.3.
A. Rekonstrukcja sygnału ze składowych za pomocą zaprojektowanych w pkt.3.3 filtrów syntezy g:
>> xr=c5_synt(g,yq,fs,varx,nrys);
% obliczenia
>> c5_synt(g,yq,fs,varx,nrys);
% wykresy
Wariancja varx sygnału oryginalnego jest wykorzystywana do uzyskania sygnału zrekonstruowanego
ˆ( )
x n
(wektor xr) o takim samym poziomie. Funkcja dokonuje zwiększenia częstotliwości próbkowania
składowych do wartości pierwotnej (upsampling razy M), filtracji BP za pomocą filtrów g i zsumowania
wydzielonych składowych zrekonstruowanych. Przeanalizować plik skryptowy c5_synt.m.
&
Porównać ("na oko") przebiegi i widma sygnału odtworzonego
ˆ( )
x n i jego składowych ˆ ( )
k
x n
z
sygnałem wejściowym x(n) i składowymi analizy y
k
(n).
• Odsłuchać sygnał odtworzony (***sygnał na głośniki***):
>> sound(xr,fs);
B. Porównanie przebiegów sygnału oryginalnego i odtworzonego:
>> c5_compare(x,xr);
&
Zwrócić uwagę i wyjaśnić przesunięcie sygnału odtworzonego względem wejściowego (patrz:
charakterystyka fazowa w następnym podpunkcie).
• Zaobserwować sumaryczne charakterystyki impulsowe oraz widmowe filtrów analizy i syntezy
łącznie oraz tzw. przeciek widmowy związany z przenikaniem części składowej z jednego podpasma
do sąsiednich, a wynikający ze skończonego tłumienia filtrów BP w pasmach zaporowych:
>> c5_bank2(h,g);
Laboratorium Cyfrowego Przetwarzania Sygnałów
Ćwiczenie 5 – Wieloczęstotliwościowe przetwarzanie sygnałów – interpolacja i decymacja
- 14 -
3.5. Powtórzenie dekompozycji i rekonstrukcji dla innych danych wejściowych
A. Powtórzyć obliczenia od pkt.3.3.C dla sygnału nie będącego szumem, np.
>> wavfile='fire8';
B. (*)
Powtórzyć obliczenia od pkt.3.3.A dla innych wartości M, L (np. L=200 - parzyste, L=401),
rs
i zaobserwować wpływ zmian na wyniki.
4. Opracowanie sprawozdania
W sprawozdaniu należy zawrzeć zarejestrowane wyniki eksperymentów numerycznych z odpowiednimi
opisami oraz wyjaśnieniami wskazanych problemów.