Równanie różnicowe filtru dolnoprzepustowego pierwszego stopnia

Funkcja przejścia filtru dolnoprzepustowego 1 stopnia dana ogólnie

$H(s) = \displaystyle\frac{1}{s R C + 1}$

Korzystając z informacji w http://en.wikipedia.org/wiki/Low-pass_filter równanie różnicowe takiego filtra możemy zapisać jako:

$y(n) = b_0x(n) + b_1x(n-1) + b_2x(n-2) - a_1y(n-1) - a_2y(n-2)$

, gdzie:

$b_0 = \displaystyle\frac{T}{RC+T}$
$b_1 = 0$
$b_2 = 0$
$a_1 = b_0 - 1$
$a_2 = 0$

Jeśli skorzystamy tutaj z metody takiej jak przy filtrze drugiego stopnia to nasze wynikowe blepy będą się nieznacznie różnić. Maksymalny błąd względny będzie rzędu 0.001. Podana tutaj wyżej metoda jest pozbawiona tego błędu.

Istnieje wiele metod konwersji filtrów ciągłych w czasie na dyskretne. W zależności od charakterystyki funkcji H(s) metody te mogą dawać mniej lub bardziej dokładną postać H(z).

Kolejna uwaga: porównując wersję lsim do wersji bq widzimy, że powyższa metoda nie generuje takich samych wyników jak lsim. Takie same wyniki jak lsim generuje metoda taka jak dla drugiego rzędu. Jeśli uznać, że lsim jest najdokładniejszy, zaś metody różnicowe są mniej dokładne to metoda z filtru drugiego stopnia jest dokładniejsza. Czemu więc w WinUAE zastosowano metodę mniej dokładną.

Porównując charakterystykę częstotliwością sygnału lsim z bq wygenerowanym przez metodę taką jak dla filtra drugiego rzędu widzimy, że choć różnice w czasie są znaczne (0.001), różnice w częstotliwości są znacznie mniejsze (0.0001).

Kolejna uwaga: oryginalny kod w pythonie wykorzystuje tak zwane frequency prewarping, brak tego w kodzie C# powoduje błąd rzędu $10^{-5}$. Przy czym nie wiem czy ten błąd jest spowodowany głównie brakiem FQ PW. Tak więc wydaje się, że stosowanie tej techniki mija się tutaj z celem.

Window function - Kaiser w Python

Tutaj implementacja w Python z wykorzystaniem odwołania do scipy.special.iv, czyli Modified Bessel function of first order. Mam nadzieję, że w następnym poście podam kod na scipy.special.iv. Wtedy wszystko będzie już gotowe do przeniesienia w C#. By tam generować BLEP-y programowo jakie się nam wymarzy.

def kaiser(samples, beta):
    result = [0] * samples
    M = samples - 1
    for i in range(samples):
        num = special.iv(1, beta*math.sqrt(1-(2*i/M - 1)**2))
        den = special.iv(1, beta)
        result[i] = num / den
    return result

Transformata Z i równanie różnicowe filtru dolnoprzepustowego drugiego stopnia

Filtr dolnoprzepustowy pierwszego stopnia zostanie omówiony w innym poście. Nie możemy użyć poniższych wzorów podstawiając po prostu A=0. Ale za to podstawienie A=0 i B=0 daje poprawne wyniki.

Filtru dolnoprzepustowego drugiego stopnia:

$H(s) = \displaystyle\frac{1}{s^2 R_1 R_2 C_1 C_2 + s(R_1 C_1 + R_2 C_2) + 1}$

Ogólnie:

$H(s) = \displaystyle\frac{1}{s^2 A + s B + 1}$

Filtr pierwszego stopnia to filtr dla którego A = 0.

Korzystając z wzoru na związek pomiędzy transformatą Laplaca a Z:

$s =\displaystyle\frac{2}{T} \frac{(z-1)}{(z+1)}$

Po przekształceniach otrzymujemy:

$H(z) = \displaystyle\frac{T^2 + 2T^2z^{-1} + T^2 z^{-2}}{4A+2TB+T^2 + (2T^2-8A)z^{-1}+(4A-2TB+T^2)z^{-2}}$

Co możemy zapisać jako:

$H(z)=\displaystyle\frac{b_0+b_1z^{-1}+b_2z^{-2}} {1+a_1z^{-1}+a_2z^{-2} }$

, gdzie:

$a_0=4A+2TB+T^2$
$a_1=(2T^2-8A)/{a_0}$
$a_2=(4A-2TB+T^2)/{a_0}$
$b_0={T^2}/{a_0}$
$b_1={2T^2}/{a_0}$
$b_2={T^2}/{a_0}$

T to okres co jaki chcemy otrzymywać kolejne próbki, w naszym przypadku jest to częstotliwość z jakiej wielokrotnością custom chip jest w stanie zmieniać wyjście.

Filtr zapisany w takiej postaci to Digital biquad filter.

Takie funkcji przejścia odpowiada następujące równianie różnicowe:

$y(n) = b_0x(n) + b_1x(n-1) + b_2x(n-2) - a_1y(n-1) - a_2y(n-2)$

Wszystko to pozwala nam na generowania przebiegów BLEP wprost z równania różnicowego, którego parametry są bezpośrednio określone parametrami elektrycznymi filtra i częstotliwością z jaką chcemy ciągły sygnał próbkować.

Rekonstrukcja dźwięku - notatki

Jak jest rozdzielczość DACa w Pauli. Wydaje się, że nawet jak na ówczesne czasy 14 bitów nie było poza ich zasięgiem, ale czy tak jest w rzeczywistości. Jak skalowane jest 8-bitów z sampla + 6-bitów z głośności.

Czy na A1200 której chipy są dwa razy szybsze Paula jest w stanie częściej zmieniać stan wyjścia DACa. To samo tyczy się A500. Istnieje tryb w którym CPU, nie DMA wypełnia rejestry dźwiękowe Pauli.

Na sztywno założyłem, że częstotliwość próbkowania hosta to 44100. Sam AudioReconstructor może pracować z dowolną. Tylko, że dla każdej takiej częstotliwości trzeba przygotować osobne BLEPy.

Generacja BLEP

Po ustaleniu, że nie filtrujemy naszego sygnału w dziedzinie częstotliwości ani czasu, tylko wykorzystujemy odpowiedzi na skok pozostaje nam ustalić jakie odpowiedzi na skok potrzebujemy.

Po pierwsze mamy nasze oryginalne BLEPy z WinUAE. Mamy BLEPy WinUAE generowane przez oryginalny skrypt. Są one takie same jak generowane przez skrypt w genblep. Wygenerowane przez nasz BLEPy różnią się od tych oryginalnych z WinUAE. Czemu nie wiem.

Orygialne BLEPy nie wykorzystały większości funkcji z pakietu scipy: lsim, tworzenie filtrów. Wszystko zostało napisane ręcznie.

Najważniejsze było dla mnie potwierdzenie, że generowane przeze mnie BLEPy są takie same w jak WinUAE. Ponieważ BLEPy WinUAE różnią zarówno od tych generowanych przez oryginalny skrypt jak i od generowanych przez mnie (wszystkie 3 grupy różnią się od siebie tak naprawdę). Napisałem więc prosty kod, które zadaniem było znalezienie takich stałych czasowych filtrów dla których moje wykresy wyglądałyby jak te w Pythonie. Stałe te są ujęte pod nazwami A500 WinUAE i A1200 WinUAE.

Dzięki temu metodą z wykorzystaniem lsim jestem w stanie wygenerować zarówno moje BLEPy jak te oryginalne z WinUAE.

W stosunku do WinUAE postanowiłem jeszcze wyprodukować wersję NTSC BLEPów, różnica jest niewielka, ale zawsze. Oczywiście można zapisać tylko jedną wersję i odpowiednio skalować i może nawet interpolować by pobrać prawidłową wartość dla NTSC. Na razie zapisuje di XMLa obie wersje, jak będzie później to się okaże.

Druga sprawa to dokładność BLEPa. W WinUAE jest to 16 bitów. U mnie na razie liczba zmiennoprzecinkowa.

Trzecie sprawa to dithering. Czy musimy go zastosować ?

Impuls służący do generacji naszego BLEPa jest potraktowany filtrem dolnoprzepustowym o częstotliwości 21KHz. Dzięki czemu próbkowanie wyjścia dźwięku z częstotliwością 44KHz powinno wywołać niewielki aliasing. Autor skryptu w Pythonie twierdzi, że powinien się też nadać do próbkowania z 48KHz. Teoretycznie jest to prawda. Na pewno nie nadaję się on do próbkowania z niższymi częstotliwościami, a takie wspierał WinUAE. Osobiście wydaje mi się, że współcześnie będę obsługiwał tylko jeden tryb 44KHz.

Częstotliwość 21KHz została wybrana zamiast 22KHz gdyż dla tej drugiej wyższe częstotliwości nie były zbyt dobrze tłumione.

Na większość tych pytań będę musiał odpowiedzieć kiedy wygeneruje WAVE z próbek wyjść Pauli ( kanał, wartość, cykl). I zrobię jakąś analizę spektralną. Dla porównania mam jedną melodię, do której posiadam także plik MOD, którą ktoś nagrał w 4 podstawowych konfiguracjach filtrów (A1200/A500, LED on/off).

No i ostatnia kwestia - wpływ starzenia się na elementy filtrów. Innymi słowy jaki był wiek Amig z których mam zgrane sygnały i czy ma on znaczenie.

Na razie mam wygenerowane BLEPy dla wszystkich parametrów filtrów zgodnie z posiadanymi schematami elektrycznymi. Oraz tak dobrane parametry by odtwarzać przebiegi z jakich korzysta WinUAE.

Kompilacja UADE pod Cygwin

Krok pierwszy to instalacja Cygwin. Wraz Z GCC4, Make, Install, Sed. Po pierwszej instalacji odpalamy instalację jeszcze raz i dorzucamy pkg-config.

Pobieramy źródła UADE, wypakowujemy je, uruchamiamy Cygwin. Montujemy katalog ze źródłami:

mount -f "D:\Programowanie\C++\Moje programy\winuae\docs\audio\uade 2.13" /src

Okazuje się, że nie mam dostępnego polecenia Make, uruchamiamy instalator i doinstalowujemy Make.

Uruchamiamy: ./configure

Dowiemy się, że nie zostaną zainstalowane jakieś pluginy dla Linuxowych playerów - nieważne. Co ważniejsze nie mamy GCC C. Znowu instalator i doinstalowujemy GCC4 C Core.

Uruchamiamy: ./configure. Tym razem brakuje libao. Doinstalowujemy.

Tak na marginesie wygląda na to, że instalator całkiem sprawnie sobie radzi przy włączonej konsoli Cygwin.

Uruchamiamy: ./configure. Bez błędów.

Teraz make. I make soundcheck.

Przy okazji warto się zapoznać z plikiem INSTALL.readme.

Teraz jeszcze dopieszczamy wszystko. Robimy by cała skompilowana zawartość lądowała w katalogu /src/bin/.

Wywołujemy:

make clean
./configure --prefix=/src/bin
make
make soundcheck
make install

Teraz czas tak dostosować proces by z poziomu Windows dało się uruchomić nasz UADE. W sumie musiałem skopiować trzy biblioteki dll z katalogu Cygwin.

Teraz musimy sobie poradzić z systemem plików. Odnajdujemy plik uadeconfig.h i modyfikujemy ścieżki na względne tak jak poustawiało nam pliki ostatnie polecenie make install.

Ostatecznie poradziłem sobie z tym modyfikując główny Makefile:

BINDIR = /src/bin
DATADIR = /src/bin/share/uade2
DOCDIR = {DOCDIR}
MANDIR = /src/bin/share/man/man1
LIBDIR = /src/bin/lib/uade2

Modyfikując src\frontends\uade123\Makefile:

BINDIR = /src/bin

I modyfikując plik uadeconfig.h:

#define UADE_CONFIG_BASE_DIR "share/uade2"
#define UADE_CONFIG_UADE_CORE "lib/uade2/uadecore.exe"

Podane zmiany są kasowane przez polecenie ./configure.

Rekonstrukcja sygnału wyjściowego na podstawie odpowiedzi na skok

Mamy funkcję odpowiedzi na skok. Każda zmian na na wyjściu ADC generuje na wyjściu filtra odpowiedź na skok. Za każdym razem naszą funkcję odpowiedzi na skok musimy przeskalować by wartości skoku były sobie równe.

Nasz filtr posiada następujące właściwości. Jest liniowy. Czyli odpowiedź na sumę dwóch sygnałów wyjściowych jest sumą odpowiedzi na każdy z sygnałów wejściowych z osobna. Jest time-invariant. Odpowiedź na sygnał wejściowy zawsze następuje po takim samym czasie.

Dzięki temu możemy powiedzieć, że sygnał wyjściowy z naszego filtru, będący odpowiedzią na pewien zbiór skoków jednostkowych na wejściu jest złożeniem odpowiedzi na skok tego filtru na każdy wejściowy skok z osobna.

Ponieważ oscylację naszego skoku z czasem stabilizują się potrzebujemy pamiętać na tyle długo na ile się zmienia. WinUAE wykorzystuje długość skoku 2048 cykli i takiej długości mniej więcej musi być bufor time-ahead.

Przykładowe odpowiedzi na skok tak jak w WinUAE:



Do przeczytania polecam ten pdf: Hard Sync Without Aliasing.

Metoda ta nada się także do symulacji dźwięku w starszych systemach, np. w C64.

Metoda ta jest dużo szybsza od stosowania splotu.

Odpowiedź na skok, a odpowiedź na impuls

Mając odpowiedź na impuls układu poszukujemy jego odpowiedzi na skok.

Splot impulsu z funkcją skoku po stronie czasu jest równoznaczny z mnożeniem po stronie częstotliwości, czyli w tym wypadku filtrowaniu.

$(f * g)(t) = \displaystyle\int_0^t f(\tau)g(t-\tau)\,d\tau$

Podstawiając za g(t) funkcję skoku otrzymujemy:

$\displaystyle\int_0^t f(\tau)\, d\tau$$

, gdzie f(t) to odpowiedź impulsowa filtra.


Przykład funkcji w matlabie do generacji odpowiedzi na skok, posiadając odpowiedź na impuls:

function [ blep ] = step_calc( blip )

step = ones(1, size(blip, 2));
blep = conv(step, blip);
blep = blep(1:1:size(blip, 2));
blep = blep / blep(end);

end

I funkcja odwrotna do powyższej:

function [ blip ] = destep_calc( blep )

blip = [blep, 0] - [0, blep];
blip = blip(1:1:end-1);
blip = blip / trapz(blip);

end

Jak we mgle

Mając podane rozwiązanie na patelni w postaci gotowego kodu w Pythonie, który to jest używany do generacji odpowiedzi na skok, które są używane do symulacji dźwięku w Amidze, Ja postanowiłem odkryć na nowo Amerykę. Główna przyczyna: nie rozumiałem tego kodu w Pythonie i wydawało mi się, że można to zrobić lepiej, bardziej elegancko.

Zacząłem od wygenerowania odpowiedzi impulsowej filtru. Z charakterystyki częstotliwościowej filtru otrzymanej z LTSpice, gdzie narysowałem schemat filtru i dokonałem symulacji pracy takiego układu. Bezpośrednio z równań funkcji przejścia, wyliczając ich transformaty odwrotne. Z funkcji MatLAB impulse. Z charakterystyki częstotliwościowej filtru wyliczonej z funkcji przejścia. Dla porównania miałem także wyniki skryptu w pythonie i tabel z WinUAE. No i jeszcze ostani sposób filtrowanie odbywa się w domenie częstotliwości i poprzez IFT generujemy impuls.

Wszystkie te rozwiązania sprowadzały się do stworzenia odpowiedzi impulsowej filtru. W oparciu o nią pokazanie jak zachowuje się filtr przepuszczając przez niego impuls, którego widmo zawiera częstotliwości do 21KHz (prostokąt w domenie częstotliwości). Wyższych nie ma co analizować raz ze względu na częstotliwości odcięcia filtrów no i tego, że mamy do czynienia z dźwiękiem.

Splot tych dwóch impulsów daje nam wynik - odpowiedź impulsową naszego filtru na sygnał zawierający częstotliwości do 21KHz.

O tym dlaczego 21KHz powiem prawdopodobnie więcej później.

Taka odpowiedź impulsowa jest podstawą do dalszej zabawy.

Zanim się połapałem o co chodzi tak naprawdę w kodzie zawartym w Pythonie, z jakiś postów których kopie mam w dokumentacji myślałem, że bierzemy taką odpowiedź impulsową i splatamy ją z naszym sygnałem wyjściowym z Pauli otrzymując przefiltrowany sygnał wyjściowy.

I to pewnie by działało, ale jest wolniejsze od rozwiązania zastosowanego w WinUAE.

Zamiast odpowiedzi impulsowej generujemy odpowiedź na skok. I teraz każda zmiana sygnału wyjściowego z Pauli generuje odpowiedź na skok - z reguły coś w rodzaju gasnącej oscylacji. Sygnał wyjściowy jest superpozycją odpowiedzi na wszystkie skoki z osobna. Z uwagi na to, że sygnał odpowiedzi na skok szybko zanika w WinUAE analizie poddaje się tylko skoki z 2048 ostatnich cykli.

Poza tym wygenerowane przeze mnie impulsy wynikowe generowane są poprzez splot, co samo w sobie powoduje zniekształcenia (oba sygnały są skończone w czasie).

Filtr drugiego rzędu

Filtr dolnoprzepustowy aktywny drugiego rzędu w konfiguracji Sallen-Key:



Filtr ten jest włączany programowo. Wraz z filtrem pierwszego rzędu powoduje jeszcze silniejsze tłumienie sygnału.

Funkcja przejścia takiego filtru:

$H(s) = \displaystyle\frac{1}{s^2 R_1 R_2 C_1 C_2 + s(R_1 C_1 + R_2 C_2) + 1}$

Częstotliwość odcięcia:

$f_{CUT} = \displaystyle\frac{1}{2\pi\sqrt{R_1 R_2 C_1 C_2}}$

Dla filtru z A500 wynosi ona 3.1kHZ.

Zgodnie z tym co czytałem tranzystor odcinający ten filtr przepuszcza w stanie odcięcia niewielki prąd. Z uwagi na to, że jest to wartość stała powinna ona zostać odfiltrowana na kondensatorach filtrujących.

Funkcja przejścia filtru pierwszego stopnia

Typowy schemat aktywnego filtru dolnoprzepustowego:



Filtr ten jest używany we wszystkich modelach Amigi. Nie da się go wyłączyć.

Funkcja przejścia ma postać:

$H(s) = \displaystyle\frac{1}{s R_2 C + 1}$

Częstotliwość odcięcia takiego filtru:

$f_{CUT} = \displaystyle\frac{1}{2 \pi R_2 C}$

Dla typowej A500 będzie to 4.4KHz.

Filtr audio

Typowy filtr audio w Amidze - po jednym na każdy kanał składa się w połączonych kaskadowo dwóch filtrów dolnpoprzepustowych: filtru pierwszego rzędu i filtru drugiego rzędu w topologii Sallen–Key. Oba filtry są operarte na wzmacniaczu LF347. Filtr drugiego stopnia jest sterowany programomo. W niektórych modelach o statnie filtru informuje dioda LED zasilania, jeśli jest przygaszona to filtr ten jest wyłączony. Stanem filtru sterujemy poprzez jedno z wyjść układu CIA.

Typowy schemat, który nie zmienia się od A1000 do A4000:

Filtry dla różnych modeli Amig różnią tylko wartością elementów RC.

Z posiadanych przeze mnie schematów wynika, że CD32 nie była wyposażona w filtr drugiego stopnia. Zaś niektóry modele A1000 i A2000 były wyposażone w zupełnie inna konfigurację elementów RC, ale pewnie dające podobną charakterystykę filtru.

Zestawione w tabele wartości z różnych schematów prezentują się następująco:

Model	C331	R331	R332	R333	C332	C333
"A500 R6a7"	100nF	360	10k	10k	6.8nF	3.9nF
"A500"	100nF	360	10k	10k	7.5nF	3.9nF
"A500+"	100nF	360	10k	10k	6.8nF	3.9nF
"A600"	100nF	360	10k	10k	6.8nF	3.9nF
"A1000"
"A1200"	6.8nF	680	10k	10k	6.8nF	3.9nF
"A1200 #2"	3.9nF	1500	10k	10k	6.8nF	3.9nF
"A2000 #1"	100nF	360	10k	10k	6.8nF	3.9nF
"A2000 #2"	47nF	750	10k	10k	6.8nF	3.9nF
"A2000 #3"
"A3000"	100nF	360	10k	10k	6.8nF	3.9nF
"CD32"	6.8nF	330	10k	10k	6.8nF	3.9nF
"CDTV"	47nF	750	10k	10k	6.8nF	3.9nF
"A4000"	47nF	750	10k	10k	6.8nF	3.9nF
"A500 winuae"	90nF	360	10k	10k	6.3nF	3.7nF
"A1200 winuae"	7.3nF	680	10k	10k	6.3nF	3.7nF

Dwa puste rzędy to wspomniane wyżej modele A1000 i A2000 z inną topologią RC.

Co do dwóch ostatnich filtrów. Tutaj trochę wybiegnę w przyszłość. W sieci znalazłem dokument który zawiera kod w pythonie z wykorzystaniem SciPy, który generuje przebiegi odpowiedzi na skok tego filtru. W kodzie tym nie ma podanych parametrów RC wprost, są za to podane parametry funkcji przejścia. Z nich wyliczyłem teoretyczne parametry C zakładając, że parametry R są takie same. Co więcej zostały one odpracowane tak by filtr programowy dawał podobne efekty jak ten w Amidze.

Tutaj do uwzględnienia jest kolejna kwestia. Starzenie się elektroniki i jej wpływ na pracę filtrów. Na razie nie potrafię oszacować, czy i jak powinna się zmienić charakterystyka 5-letniego filtru.

Generacja dźwięku - wstęp

Za generację dźwięku w Amidze, odpowiedzialny jest układ Paula. Układ ten posiada dwa analogowe wyjścia audio - kanał lewy i prawy. Rozdzielczość ADC to 8-bitów.

Dodatkowo istnieją jakieś tryby kiedy jeden kanał moduluje częstotliwość lub/i amplitudę drugiego. W takim trybie możliwe jest uzyskanie sygnału o 14-bitowej rozdzielczości. Na pewno jest to bardzo rzadko wykorzystywane, nie jestem pewny czy WinUAE to nawet emuluje. Gry z tego nie korzystają, ewentualnie jakieś programy muzyczne. Na razie tymi trybami się nie przejmuje.

Paula wewnętrznie posiada 4 kanały, po dwa na każdy kanał stereo.

Amiga to komputer, którego działanie silnie zależy od zegara taktującego. Cały proces generacji obrazu, dźwięku odbywa się w takt zegara. Ponieważ obraz generowany jest w PAL i NTSC dla obu systemów mamy inne zegary taktujące, odpowiednio $FREQ_{PAL} = 7093790$, $FREQ_{NTSC} = 7159090$. Stąd wynika maksymalna częstotliwość sygnału jaką Paula jest w stanie wygenerować: połowa powyższego.

Z uwagi na czas w jakim powstała Amiga zastosowanie DSP nie było możliwe. Mimo to Amiga radzi sobie z generacją zdigitalizowanego dźwięku, ale odbywa się to kosztem jakości.

Paula generuje sygnał o maksymalnej częstotliwości około 3.5 MHz, który jest filtrowany i wzmacniany.

Paula w każdej linii pobiera 4 16-bitowe słowa, czyli po dwie próbki na każdy z 4 kanałów. Ponieważ moment uzyskania dostępu do pamięci w linii ekranu nie jest ustalony Paula posiada bufor 2x4 bajty. Kiedy 4 próbki są odtwarzane, 4 kolejne mogą zostać załadowane.

Częstotliwości wyświetlania obrazu są następujące:

$FPS_{PAL} = 50$
$FPS_{NTSC} = 59.94$

Liczba linii ma klatkę w obu systemach:

$LINES_{PAL} = 625/ 2$
$LINES_{NTSC} = 525/2$

525 i 625 to liczba linii w przeplocie czyli na dwie klatki.

Stąd możemy policzyć ile maksymalnie próbek na sekundę w obu systemach jest w stanie pobrać i wygenerować Paula:

$SAMPLESMAX_{PAL} = 2 * LINES_{PAL} * FPS_{PAL} = 31250$
$SAMPLESMAX_{NTSC} = 2 * LINES_{NTSC} * FPS_{NTSC} = 31468.5$

2 to liczba próbek na kanał jakie maksymalnie może Paula pobrać z pamięci podczas jednej linii.

Tutaj trochę mi mieszają. Zgodnie z dokumentacją liczba próbek została ograniczona do $SAMPLESMAX = 28867$.

Tak więc teraz możemy policzyć co ile minimalnie cykli może się zmienić odtwarzana próbka by zachować to kryterium:

$CYCLESPERSAMPLE_{PAL} = FREQ_{PAL} / SAMPLESMAX / 2 = 123$
$CYCLESPERSAMPLE_{NTSC} = FREQ_{NTSC} / SAMPLESMAX / 2 = 124$

2 gdyż, częstotliwość procesora to jedno, a cykle to drugie. CPU posiada 16-bitową magistralę danych, zaś operuje na słowach 32-bitowych, czyli dostęp do każdej danej zajmuje mu dwa cykle zegarowe.

Paula posiada rejestr 16-bitowy definiujący przez ile cykli może być utrzymywana próbka na wyjściu. Jeśli wartość jest mniejsza od minimalnej nie ma gwarancji, że następna próbka zostanie uzyskana na czas. Odtwarzana jest wtedy próbka poprzednia.

Z pewnością jest tutaj jeszcze wiele problemów do rozwiązania. Skąd te 28867. Choć liczba linii w PAL to 625 maksymalna rozdzielczość to 400 linii. Reszta linii tworzy ramkę. Co się konkretnie dzieje jak nie mamy następnej próbki, a następną trzeba odtworzyć, co się dzieje jak właściwa próbka nadejdzie, czy następuje jakaś desynchronizacja.

Z drugiej strony możemy policzyć minimalną częstotliwość z jaką może być odtwarzany dźwięk. Zakładamy, że w rejestr jest wpisana wartość 63355.

$SAMPLESMIN_{PAL} = FREQ_{PAL} / 65535 / 2 = 54$
$SAMPLESMIN_{NTSC} = FREQ_{NTSC} / 65535 / 2 = 55$

Wartości te to ilość próbek generowanych na sekundę. Taka liczba próbek pozwala odtworzyć sinusoidę o częstotliwości 27 i 27.5. Z drugiej strony maksymalna częstotliwość to 14.4Khz. Człowiek jest w stanie usłyszeć dźwięki od około 20 do 20000Hz. 14.4KHz to ciągle jest dużo, gdyż 200000 to bardziej wartość teoretyczna.

Generowanie jednokanałowego dźwięku z maksymalną częstotliwością wymaga strumienia danych rzędu 28KB na kanał.

Przypuśćmy, że mamy zarejestrowaną próbkę jakiegoś instrumentu z tą częstotliwością. I chcemy ją odegrać na niższym tonie. Odtwarzamy więc ją wolniej wpisując do rejestru wartość większą niż 124. Konkretna wysokość tonu zależy od częstotliwości z jaką dźwięk został zarejestrowany.

Normalnie współcześnie melodia, efekty dźwiękowe są przygotowane w postaci strumienia WAV (MP3). Na Amidze z uwagi na wielkość pamięci maszyn z tamtych czasów nie było technicznie możliwe by odtwarzać muzykę jako WAV.

Próbki instrumentów były rejestrowane każdy z osobna, a następnie odtwarzane z wybraną częstotliwością na 4 kanałach (4 głosowa polifonia).

W tamtych czasach nie było praktycznie możliwe użycie DSP to konwersji częstotliwościowej próbki. Było to zbyt drogie. CPU był o wiele za wolny na tą operację. Taka konwersja, czyli resampling, to proces podczas którego zmieniamy częstotliwość sygnału i jednocześnie filtrujemy dolno-przepustowo sygnał by wyeliminować zjawisko aliasingu. Z braku takiego układu dźwięk generowany jest zniekształcony.

Częściowo można temu zaradzić przygotowując więcej niż jedną próbkę o kilku pośrednich częstotliwościach. Prawdopodobnie było to rzadko używane w grach, ale częściej przy odtwarzaniu modułu.

Rozwiązaniem problemu było filtrowanie sygnału. Każdy kanał z osobna. Filtr taki składał się z kaskadowo połączonych dwóch filtrów, jeden z nich włączony na stałe, drugi włączany programowo. Filtry te i ich parametry zostaną omówione później.

Z uwagi na to, że filtr włączony na stałe ma częstotliwość odcięcia 4.4KHz, generalnie nie jesteśmy w stanie usłyszeć niczego powyżej 7KHz. Z tego powoduje próbki powinny być rejestrowane z częstotliwością wyższą niż 14KHz. Czyli pomimo tego, że teoretycznie maksymalna częstotliwość jaką jesteśmy w stanie wygenerować to 14KHz, filtr ucina wszystko co jest powyżej 7KHz. Związku z tym (za dokumentacją) dla próbkowania większego niż 320 (4KHz) cykli zaczynamy gubić generowane wyższe częstotliwości.

Ich zadaniem jest eliminacja aliasingu, który charakteryzuje się w przypadku dźwięku z Amigi zakłóceniami na wysokich częstotliwościach.

Każdy z czterech kanałów jest wyposażony w 6 bitową regulację głośności.

Mimo wszystko było to prawdopodobnie wszystko co najlepsze co dało się uzyskać w tamtych czasach.

Resampling

To kombinacja upsamplingu i downsamplingu. Pozwala nam ona zmienić próbkowanie sygnału o niecałkowitą skale. Ponieważ jest to kombinacja upsamplingu i downsamplingu dla pewnych współczynników skalowania może się okazać, że współczynnik upsamplingu może być bardzo duży. Takie przypadki w moim przykładzie staram się odrzucać. No i dla współczynników skalowania całkowitych albo ich odwrotności musimy tylko skorzystać z upsamplingu lub downsamplingu.

Funkcja:

function [ tout, fout ] = resampling( fin, fsin, fsout )

m = lcm(fsin, fsout);
u = m / fsin;
d = m / fsout;

if (u > 1000)
    error('up scale factor too big');
end

if (u > 1)
    [tout, fout] = upsampling(fin, fsin, fsin*u);
    if (d > 1)
        [tout, fout] = downsampling(fout, fsin*u, fsout);
    end
elseif (d > 1)
    [tout, fout] = downsampling(fin, fsin, fsout);
else
    tout = (1/fsin) * linspace(0, size(fin, 2), size(fin, 2)); 
    fout = fin;
end

end

Ograniczenie na błąd powinniśmy sobie sami ustawić, w zależności jakiej długości sygnały przetwarzamy i ile mamy pamięci, ... czasu.

Upsampling

Upsampling czyli zwiększenie częstotliwości próbkowania sygnału o zadany całkowity współczynik przy zachowaniu jego widma.

Funkcja:

function [ tout, fout ] = upsampling( fin, fsin, fsout )

scale = fsout / fsin;

if (fix(scale) ~= scale)
    error('scale should be integer');
end

if (scale <= 0)
    error('scale should be greater than zero');
end

fin = fin * scale;
fout = upsample(fin, scale);
tout =  linspace(0, (1 / fsout) * (size(fout, 2) - 1), size(fout, 2)); 

[~, k] = sinc_gen(fsout, fsin/2);
if (mod(size(k, 2), 2) == 0)
    [~, k] = sinc_gen(fsout, fsin/2, fix(size(k, 2) / 2) * 2 + 1);
end

fout = conv(fout, k, 'same');

fout = fout(1:1:end-scale+1);
tout = tout(1:1:end-scale+1);

end

Funkcja upsample uzupełnia fin zerami. Zera w ilości scale-1 są wstawiane między próbki i na koniec. My te końcowe zera usuwamy. Dzięki czemu ostatnia próbka sygnału wyjściowego odpowiada ostatniej próbce sygnału wejściowego.

Okno sinc musi mieć nieparzystą długość by uniknącć przesunięcia wyniku.

Sam splot to tzw. idealna rekonstrukcja, tutaj o tyle prostsza i szybsza, że współczynnik skalowania jest całkowity i sinc rekonstruujący może być wyliczony raz i używamy dla każdego punktu z osobna.

W stosunku do wbudowanej w Matlab funkcji interp błąd przez nas popełniany jest na tyle niewielki, że możemy powiedzieć, że wynika on z innej metody interpolacji.

Przykład:

%% input data

fsin = 6000;
fsout = 30000;
scale = fsout / fsin;

func = @(ti) sin(2*pi*650*ti-pi/4) + sin(2*pi*1500*ti+pi/2);
tend = 1/500*2;

% build input func.
Tin = 1/fsin;
tin = (0:Tin:tend);
fin = func(tin);

% reconstruct
[fout, tout] = upsampling(fin, fsin, fsout);

% calculate ideal reconstruction
forg = func(tout);

fout1 = interp(fin, scale);
fout1 = fout1(1:1:end-scale+1);

%% compare them

hold on;
plot(tout, forg, 'r');
plot(tin, fin, 'g');
plot(tout, fout, 'b');
plot(tout, fout1, 'm');
legend('ideal','in', 'out', 'interp');

%% compare fq

FIN = fftshift(abs(fft(fin)));
FOUT = fftshift(abs(fft(fout)));
FOUT1 = fftshift(abs(fft(fout1)));
s = size(FIN, 2) / size(FOUT, 2);
FORG = s*fftshift(abs(fft(forg)));
FOUT = s*FOUT;
FOUT1 = s*FOUT1;

fqout = fsout/2*linspace(-1, 1, size(FORG, 2));
fqin = (1/Tin)/2*linspace(-1, 1, size(FIN, 2));

hold on;
plot(fqout, FORG, 'r');
plot(fqin, FIN, 'g');
plot(fqout, FOUT, 'b');
plot(fqout, FOUT1, 'm');
legend('ideal','in', 'out', 'interp');

%% error

error_calc(fout, forg);
error_calc(fout1, forg);

Downsampling

Downsampling czyli zmniejszenie częstotliwości próbkowania sygnału o zadany całkowity współczynnik przy zachowaniu jego widma (poniżej częstotliwości do której zmniejszamy). Ponieważ tutaj rekonstrukcja odbywa się z częstotliwością mniejszą o całkowitą wartość w przeciwieństwie do rekonstrukcji wystarczy nam tylko jeden impuls sinc. Cały kod jest znacznie szybszy.

Funkcja:

function [ fout, tout ] = downsampling( fin, fsin, fsout )

scale = fsin / fsout;

if (fix(scale) ~= scale)
    error('scale should be integer');
end

if (scale <= 0)
    error('scale should be greater than zero');
end

% works only with odd window.
k = sinc_gen(fsin, fsout/2);
if (mod(size(k, 2), 2) == 0)
    k = sinc_gen(fsin, fsout/2, size(k, 2) + 1);
end

tout = linspace(0, (1/fsin) * (size(fin, 2)-1), ... 
                ceil(size(fin, 2)/scale));
fout = zeros(1, size(tout, 2));

for m=1:1:size(fout, 2)
    fout(m) = conv_point(fin, k, (m - 1) * scale + 1);
end

end

Przykład:

%% input data

% resampling
fsin = 30000;
fsout = 6000;
scale = fsin / fsout;

func = @(ti) sin(2*pi*650*ti-pi/4) + sin(2*pi*1500*ti+pi/2);
tend = 1/500*2.4; % czym dluzszy tym lepsze widmo

% build input func.
Tin = 1/fsin;
tin = (0:Tin:tend);
fin = func(tin);

% downsampling
[fout, tout] = downsampling(fin, fsin, fsout);

fout1 = decimate(fin, scale);

% calculate ideal reconstruction
forg = func(tout);

%% compare them

fig = figure;
hold on;
plot(tout, forg, 'r');
plot(tin, fin, 'g');
plot(tout, fout, 'b');
plot(tout, fout1, 'm');
legend('ideal','in', 'out', 'decimate');

%% compare fq

FIN = fftshift(abs(fft(fin)));
FOUT = fftshift(abs(fft(fout)));
FOUT1 = fftshift(abs(fft(fout1)));
s = size(FIN, 2) / size(FOUT, 2);
FORG = s*fftshift(abs(fft(forg)));
FOUT = s*FOUT;
FOUT1 = s*FOUT1;

fqout = fsout/2*linspace(-1, 1, size(FORG, 2));
fqin = (1/Tin)/2*linspace(-1, 1, size(FIN, 2));

hold on;
plot(fqout, FORG, 'r');
plot(fqin, FIN, 'g');
plot(fqout, FOUT, 'b');
plot(fqout, FOUT1, 'm');
legend('ideal','in', 'out', 'decimate');

%% error;

error_calc(fout, forg);
error_calc(fout1, forg);

Błąd mojej wersji jest porównywalny do tej wbudowanej w Matlab. Ilość próbek wejściowych jest równa wielokrotnością skalowania. Jak widać obie wersje rozmijają się, Matlab na początku, ja na końcu. Dla innych liczb próbek obie funkcję potrafią zwrócić bardzo porównywalne wyniki.

Idealna rekonstrukcja sygnału

Mamy sygnał fin, który został spróbkowany z częstotliwością fsin. Chcemy wygenerować sygnał, który został spróbkowany z częstotliwością fsout. Jeśli fsout > fsin w sygnale wyjściowym nie powinno być żadnej częstotliwości powyżej fsin. Jeśli fsout < fsin w sygnale wyjściowym powinny się znaleźć tylko częstotliwości mniejsze od fsout z sygnału wejściowego.

Wszelkie interpolacje (spline, bicubic, itp) w dziedzinie czasu odpadają, gdyż wprowadzimy tym częstotliwości do sygnału, których w nim nie było oryginalnie. Jedne sprawują się gorzej, inne lepiej, ale generalnie interpolacje odtworzą sygnał wizualnie niby dobrze, częstotliwościowo nie będzie to to samo. Ma to szczególne znaczenie kiedy pracujemy z sygnałem dźwiękowym. No dobra za wyjątkiem podanej niżej interpolacji. Rekonstrukcję możemy też przeprowadzić z dziedzinie częstotliwości.

Aby odtworzyć oryginalny sygnał potrzeba nam minimum dwie próbki na na okres, w praktyce powinno być ich więcej. Możemy próbkować dokładnie w miejscach gdzie sinusoida przechodzi przez zero. A nawet jeśli jesteśmy minimalnie minięci podczas próbkowania z zerem, pobierane przez nas wartości są małe i błąd przez to jest duży.

Teoria mówi, że sygnał oryginalny może zostać dokładnie zrekonstruowany ze swych próbek, jeśli próbki były pobierane w odstępach czasu T, zaś sygnał próbkowany nie zawiera częstotliwości większych niż połowa częstotliwości próbkowania. Czyli jeśli próbkujemy z częstotliwością 44100Hz, to sygnał nie powinien zawierać częstotliwości większych niż 22050Hz. W przeciwnym razie częstotliwości wyższe od limitu pojawią się w sygnale zrekonstruowanym zniekształcając go (aliasign). Twierdzenie Whittaker–Shannon'a mówi:

$x(t) = \sum_{n=-\infty}^{\infty} x[n] \cdot {\rm sinc}\left(\frac{t - nT}{T}\right)\,$

Inaczej sygnał oryginalny to splot impulsów z próbek z funkcją sinc o odpowiednim charakterze.

Teraz mając odtworzony sygnał ciągły, próbkujemy go ponownie z nową częstotliwością. Jeśli z zmniejszą to najpierw pozbywamy się z niego częstotliwości większych od połowy nowej częstotliwości próbkowania.

Całość robi poniższy kod. Jeśli mamy zrekonstruować sygnał z niższą częstotliwością to filtrowane są próbki wejściowe. Taki sygnał podlega rekonstrukcji zgodnie z wyżej podaną formułą tylko w interesujących nas czasach.

Kod funkcji:

function [ fout, tout ] = reconstruction( fin, fsin, fsout )

Tout = 1 / fsout;
Tin = 1 / fsin;

if (fsin == fsout)
    fout = fin;
    tout = Tout* linspace(0, size(fout, 2), size(fout, 2)); 
    return;
end

scale = fsout / fsin;

if (scale < 1)
    % for compatibility with downsampling and decimate
    fout = zeros(1, ceil(size(fin, 2)* scale));
else
    fout = zeros(1, fix(size(fin, 2) * scale));
end

tout = linspace(0, (1/fsin) * (size(fin, 2) - 1), size(fout, 2));            

% when downsampling cut frequencies below fsout/2 from fin 
% function.
if (fsin > fsout)
    k = sinc_gen(fsin, fsout/2);
    fin = conv(fin, k, 'same');
end

W = fix(160 / (fsout*2 / fsin)); % experimental
W = min(W, size(fin, 2));

tin = (0:1:size(fin, 2)-1) * Tin;

if (W > size(fin, 2))
    W = size(fin, 2);
end

for m=1:size(fout, 2)
    [~, min_index] = min(abs(tout(m) - tin));
    min_index = fix(min_index - W/2);
    if (min_index < 1)
        min_index = 1;
    end
    max_index = min_index + W - 1;
    if (max_index > size(fin, 2))
        max_index = size(fin, 2);
    end
    sc = sinc((tout(m) - (min_index-1:1:max_index-1)*Tin) / Tin);
    fout(m) = sum(fin(min_index:1:max_index) .* sc);
end

% for compatibility with upsampling
if (scale > 1)
    fout = fout(1:1:end-scale+1);
    tout = tout(1:1:end-scale+1);
end

Dla idealnej rekonstrukcji okno W powinno być równe długości funkcji, w praktyce może być znacznie krótsze. Jeśli mamy do czynienia z downsamplingiem to najpierw sygnał trzeba przefiltrować filtrem dolnoprzepustowym. Inaczej częstotliwości powyżej fsout, które były w fin odbiją się i pojawią się jako zakłócenia w sygnale. Z tak zakłóconym sygnałem już niż nie możemy zrobić - nie da się stwierdzić co jest oryginałem, a co zakłóceniem.

Rekonstrukcja taka jest dokładna, ale wyjątkowo czasochłonna, dla każdej próbki wyjściowej trzeba policzyć osobny impuls sinc.

Normalnie okno sinc powinno być równe wektorowi wejściowemu, tutaj jest krótsze, o długości powiedzmy eksperymentalnej dla której nie widać różnicy.

Jest też naniesionych parę poprawek na kompatybilność z innymi funkcjami dotyczącymi długości wektora wyjściowego.

Przykład:

%% input data

% resampling
fsin = 30000;
fsout = 7300;
Tin = 1/fsin;
scale = fsin / fsout;

func = @(ti) sin(2*pi*650*ti-pi/4) + sin(2*pi*1500*ti+pi/2);
tend = 1/500*2;

%build input func.
tin = (0:Tin:tend);
fin = func(tin);

% reconstruct
[fout, tout] = reconstruction(fin, fsin, fsout);

%resample
[fout1, tout1] = resampling(fin, fsin, fsout);

% calculate ideal reconstruction
forg = func(tout);

%% compare them

fig = figure;
hold on;
plot(tout, forg, 'r');
plot(0:Tin/20:tend, func(0:Tin/20:tend), 'g');
plot(tout, fout, 'b');
plot(tout1, fout1, 'black');
legend('ideal','in', 'out reconstr', 'out resampl');

%% compare fq

FORG = fftshift(abs(fft(forg)));
FIN = fftshift(abs(fft(func(0:Tin:tend))));
FOUT = fftshift(abs(fft(fout)));
FOUT1 = fftshift(abs(fft(fout1)));

fqout = fsout/2*linspace(-1, 1, size(FORG, 2));
fqout1 = fsout/2*linspace(-1, 1, size(FOUT1, 2));
fqin = (1/Tin)/2*linspace(-1, 1, size(FIN, 2));

s = size(FIN, 2) / size(FOUT, 2);

hold on;
plot(fqout, s*FORG, 'r');
plot(fqin, FIN, 'g');
plot(fqout, s*FOUT, 'b');
plot(fqout1, s*FOUT1, 'black');
legend('ideal','in', 'out reconstr', 'out resampl');

%% error

error_calc(fout, forg);
error_calc(fout1, forg);

Filtr FIR o minimalnym przesunięciu fazowym

Weźmy wszystkie filtry o takiej samej charakterystyce częstotliwościowej G, ale różniących się charakterystyką przesunięcia fazowe. Filtr o minimalnym przesunięciu fazowym to taki filtr, że w dowolnym okresie czasu $<0, t>$ energia sygnału wyjściowego jest większa lub równa dowolnemu innemu filtrowi z G. Czyli sygnał w takim filtrze najszybciej narasta i najszybciej maleje. Sygnał jest maksymalnie przesunięty ku zeru.

Dla filtrów FIR Matlab posiada funkcję do znajdowania MPF (minimum-phase filter) rceps. Niestety nie działa ona za dobrze. Rezultat nie wygląda tak jak powinien, a na dodatek dla niektórych długości wektora policzenie wyniku jest niemożliwe.

Moja wersja:


function [ minblip ] = fir_min_phase( impulse )

impulse1 = [impulse, zeros(1, size(impulse, 2)*7)];

y = ifft(log_safe(abs(fft(impulse1))));

%y(1) = y(1) * 1;
y(2:fix(end/2)) = y(2:fix(end/2)) * 2;
%y(fix(end/2+1) = y(fix(end/2+1) * 1;
y(fix(end/2)+2:end) = 0;

minblip = real(ifft(exp(fft(y))));

minblip = minblip(1:1:size(impulse, 2));

minblip = minblip / trapz(minblip);

end

function [ r ] = log_safe( v )

r = log(v);

for k=1:1:size(r, 2)
if (r(k) == -inf)
r(k) = -1e-100;
end
end

end


Kompletna magia dla mnie - matematycznie nie wiem jak to działa. Zmniejszając początkową rozdzielczość okna, albo źle dobierając okno +/-1 możemy otrzymać wersję Matlabową.

Porównanie impulsu sinc i jego wersji MPF.


fs = 7093790/2;
d1 = firwin(fs, 21000, 2047);
d2 = fir_min_phase(d1);
[~, d3] = rceps(d1);

hold on;
plot(d1, 'r');
plot(d2, 'g');
plot(d3, 'b');
legend('org', 'my', 'matlab');

D1 = abs(fftshift(fft(d1)));
D2 = abs(fftshift(fft(d2)));
D3 = abs(fftshift(fft(d3)));
fx = linspace(-fs/2, fs/2, 2047);

hold on;
plot(fx, D1, '*');
plot(fx, D2, 'g');
plot(fx, D3, 'o');
legend('org', 'my', 'matlab');


Widzimy, że we wszystkich 3 przypadkach charakterystyka częstotliwościowa filtrów jest taka sama.

Idealny filtr dolnoprzepustowy

Idealny filtr dolnoprzepustowy powinien po stronie częstotliwości być linią prostą o wartości dodatniej od f=0 do częstotliwości odcięcia i powyżej nie przyjmować wartość zero. Ponieważ w transformacie Fouriera odwrotnej chcąc po stronie czasu uzyskać wartości rzeczywiste prostokąt ten powinien być symetryczny względem f=0.

Najpierw zdefiniujmy sobie tzw. funkcję prostokątną:

$\text{rect}(t) = \begin{cases}0 & \text{if } \; |t| \;> \;^1/_2 \\ ^1/_2 & \text{if }\; |t| \;= \;^1/_2 \\1 & \text{if }\; |t| \;<\; ^1/_2 \\\end{cases}$

W naszym wypadku przyjmując za częstotliwość odcięcia B filtr po stronie częstotliwości powinien mieć postać $\displaystyle rect(\frac{f}{2B})$.

Jego transformatą po stronie czasu jest:

$\displaystyle f(t)=\mathcal{F}^{-1}\{ rect(\frac{f}{2B})\} = 2Bsinc(2Bt)$

Funkcja generująca taki filtr:

function [ k, t ] = sinc_gen( sampling_rate, cutfq, samples )

  T = 1 / sampling_rate; 

  if (nargin == 2)
    samples = 160 / (cutfq*2 / sampling_rate); % experimental
  end

  tout = T * (1:samples); 
  tout = tout - T*(samples+1)/2;
  B = 2*cutfq;
  k = sinc(B*tout);
  k = k / trapz(k);

  if (nargout == 2)
      t = tout;
  end

end

Jeśli nie podamy ilości próbek jest ona ustalana na najlepszą (według mnie) wartość, która przy filtrowaniu nie powoduje błędów.

Zmienna sampling_rate to częstotliwość próbkowania sygnału. sampling_rate/2 to największa częstotliwość w sygnale. Czyli np: jeśli sampling_rate=44100, to $f_{max}=22050$. Na samym końcu impuls jest normalizowany, tak by jego pole powierzchni było równe jeden, dzięki czemu splot z sygnałem zachowuje jego energię.

Zauważmy, że dla takiej samej wartości samples jeśli stosunek sampling_rate do cutoff_fq FIR będzie taki sam to odpowiedź impulsowa będzie taka sama (w sensie wartości w tablicy).

Przykładowy filtr sinc:

[~, w3] = sinc_gen(44100, 4000, 200);
tx = linspace(0, 200/44100, 200);
plot(tx, w3);
title('sampling rate=44100Hz, cutoff fq=4kHz, samples=200');

Filtr jest symetryczny i zmierza do zera jak czaszmierza do nieskończoności. Dla parzystej liczby sampling_rate dwie środkowe próbki mają takie same wartości. Dla nieparzystej największą wartość ma środkowa próbka (ten jest bardziej pożądany, gdyż nie przesuwa sygnału podczas splotu).

Wraz ze zwiększaniem sampling rate i utrzymywaniem wysokiej cutoff_fq kształt filtru będzie coraz bardziej przypominać impuls.

W dziedzinie częstotliwości filtr taki wygląda tak:

f3 = abs(fftshift(fft(w3)));
fx = linspace(-44100/2, 44100/2, 200);
plot(fx, f3);
title('sampling rate=44100Hz, cutoff fq=4kHz, samples=200');

Oscylacje na brzegach są spowodowane tym, że nas filtr jest gwałtownie się urywa na brzegach. Aby tego uniknąć na filtry nakłada się okno, jest ich wiele rodzajów, i nie zamierzam wnikać, który się do czego się nadaje.

w4 = w3 .* blackman(200)';
f3 = abs(fftshift(fft(w4)));
fx = linspace(-44100/2, 44100/2, 200);
plot(fx, f3);
title('sampling rate=44100Hz, cutoff fq=4kHz, samples=200');

Moja wyspecjalizowana funkcja generująca sinca:

function [ k, t ] = firwin( fs, fcut, len, prevent_shift )

k = sinc_gen(fs, fcut);

if (sum(size(len) - [1 1]) == 0)
    L = len;
else
    L = size(len, 2);
end

if (size(k, 2) > L)
    k = sinc_gen(fs, fcut, L);
end;

p = true;
if (nargin == 4)
    p = prevent_shift;
end

% prevent time shifting
if (p)
    if (mod(size(k, 2), 2) == 0)
        k = sinc_gen(fs, fcut, size(k, 2) - 1);
    end
end

% apply window
w = kaiser(size(k, 2), 9);
k = k .* w';

% normalize
k = k / trapz(k);

if (nargout == 2)
    t = linspace(0, (size(k, 2)-1)/fs, size(k, 2));
end

end

Korzysta ona z ustawionego na stałe okna (Kaiser, 9). Można także zapobiegać generacji sinca o parzystej liczbie próbek, splot z takim sincem w czasie przesuwa nam rezultat o jeden.

Porównanie do funkcji wbudowanej w Matlaba:

spectrum_cut = 50000;
Fs = 350000;
len = 100;

k1 = fir1(len-1, 2 * spectrum_cut / Fs, kaiser(len, 9));
k1 = k1 / trapz(k1);

k2 = firwin(Fs, spectrum_cut, len, false);

hold on;
plot(k1, '*');
plot(k2, 'r');

err = error_calc(k1, k2);

Funkcja pomocnicza:

function [ r ] = error_calc( v1, v2 )

res = sum(abs((v1 - v2) / v2) / size(v1, 2));

if (nargout == 0)
    display(sprintf('error: %f', res));
else
    r = res;
end

end

Splot

Poniżej podane są funkcje w Matlabie realizujące splot, działające jak ich odpowiedniki w Matlabie:

function [ fout ] = conv_full( fin, k )

  if (size(fin, 1) ~= 1)
    error('bad fin size');
  end
  
  if (size(k, 1) ~= 1)
    error('bad k size');
  end
  
  s1 = size(fin, 2);
  s2 = size(k, 2);
  
  fout = zeros(1,s1+s2-1);
  
  for m=1:size(fout, 2)
    fin1 = max(m-s2+1, 1);
    fin2 = min(m, s1);
    k1 = max(m-s1+1, 1);
    k2 = min(m, s2);
    fout(m) = sum(fin(fin1:fin2) .* k(k2:-1:k1));
  end

end

function [ fout ] = conv_same( fin, k )

  if (size(fin, 1) ~= 1)
    error('bad fin size');
  end
  
  if (size(k, 1) ~= 1)
    error('bad k size');
  end
  
  s1 = size(fin, 2);
  s2 = size(k, 2);
  
  L = fix(size(k,2)/2)+1;
  H = L + size(fin, 2) - 1;
  
  fout = zeros(1, H-L+1);
  
  for m=L:H
    fin1 = max(m-s2+1, 1);
    fin2 = min(m, s1);
    k1 = max(m-s1+1, 1);
    k2 = min(m, s2);
    fout(m-L+1) = sum(fin(fin1:fin2) .* k(k2:-1:k1));
  end

end

Funkcja conv_same przydaje się w przetwarzaniu sygnałów, wtedy z reguły u traktujemy jako sygnał zaś v jako odwrotną transformatę Fouriera jakiegoś filtru, którego chcemy użyć. Splotowi u i v w czasie odpowiada mnożenie po stronie częstotliwości. Wynik conv_same będzie takiego samego rozmiaru jak u.

Funkcja v może być także odpowiedzią impulsową funkcji przejścia H(s) dowolnego układu. Znając H(s) liczymy jaka jest odpowiedź układu na impuls (1 po stronie s), czyli jest to transformata odwrotna 1*H(s). Mówimy wtedy o odpowiedzi impulsowej filtru. Po stronie czasu mamy układ o funkcji przejścia H(s), podajemy na jego wejście impuls (po stronie częstotliwości widmo impulsu jest ciągłe od minus nieskończoności do plus nieskończoności i o stałej amplitudzie). To co otrzymujemy na wyjściu to odpowiedź w czasie przefiltrowanych wszystkich możliwych częstotliwości. Transformatą Fouriera tej odpowiedzi wyznacza nam charakterystykę częstotliwością filtru.

I jeszcze przykład funkcji realizującej splot dla danego punktu:

function [ r ] = conv_point( fin, k, point )

  if (size(fin, 1) ~= 1)
    error('bad fin size');
  end
  
  if (size(k, 1) ~= 1)
    error('bad k size');
  end
  
  if (point < 1)
    error('bad point index');
  end
  
  if (point > size(fin, 2))
    error('bad point index');
  end
  
  s1 = size(fin, 2);
  s2 = size(k, 2);
  
  L = fix(size(k,2)/2)+1;
  H = L + size(fin, 2) - 1;
  
  m = point + L - 1;
  
  fin1 = max(m-s2+1, 1);
  fin2 = min(m, s1);
  k1 = max(m-s1+1, 1);
  k2 = min(m, s2);
  r = sum(fin(fin1:fin2) .* k(k2:-1:k1));

end

Porównanie wersji same i full:

fsin1 = 500; % czestotliwosc pierwsze sinusoidy
fsin1_phase = pi/2; %-pi/4.14; % przesuniecie fazowe pierwszej sinusoidy
fsin2 = 2000; % czestotliwosc drugiej sinusoidy
fsin2_phase = -pi/4; %pi/2.11; % przesuniecie fazowe drugiej sinusoidy
fs = 24000; % czestotliwosc probkowania

% funkcja wejsciowa.
func_in = @(ti) sin(2*pi*fsin1*ti+fsin1_phase) + ...
sin(2*pi*fsin2*ti+fsin2_phase);

fcut = 1000; % czestotliwosc odciecia
tend = 3/fsin1; % dlugosc funkcji
sinc_samples = 200;

% przygotuj wektor czasu.
ts = 1/fs; % odstep czasu pomiedzy probkowaniami
tx = (0:ts:tend);

if (mod(size(tx, 2), 2) == 1)
tx = tx(1:1:end-1);
end

% wejscie
fin = func_in(tx)/15;

% filtr
w1 = sinc_gen(fs, fcut, sinc_samples);

fout1 = conv(fin, w1, 'same');
fout2 = conv(fin, w1, 'full');

figure(1);
hold on;
plot(fin, 'r');
plot(w1, 'g');
plot(fout1, 'b');
plot(fout2, 'c');
legend('wejscie', 'filtr', 'conv full', 'conv same');

Przykład:

fs = 30000; % czestotliwosc probkowania
fcut = 1000; % czestotliwosc odciecia
fsin1 = 500; % czestotliwosc pierwsze sinusoidy
fsin1_phase = pi/2; %-pi/4.14; % przesuniecie fazowe pierwszej sinusoidy
fsin2 = 2000; % czestotliwosc drugiej sinusoidy
fsin2_phase = -pi/4; %pi/2.11; % przesuniecie fazowe drugiej sinusoidy
tend = 1/fsin1*3.35; % dlugosc funkcji
sinc_samples = 1000;

% funkcja wejsciowa.
func_in = @(ti) sin(2*pi*fsin1*ti+fsin1_phase) + ...
sin(2*pi*fsin2*ti+fsin2_phase);

% funkcja wyjsciowa - idealne odciecie.
func_out = @(ti) sin(2*pi*fsin1*ti+fsin1_phase);

% przygotuj wektor czasu.
ts = 1/fs; % odstep czasu pomiedzy probkowaniami
tx = (0:ts:tend);

% oblicz wartosci obu funkcji w czasie.
fin = func_in(tx);
fout = func_out(tx);

%[~, fin] = sinc_gen(fs, (fsin1+fsin2)*1.12, size(tx, 2));
%[~, fout] = sinc_gen(fs, fsin1, size(tx, 2));

% przygotuj sinca
[~, w1] = sinc_gen(fs, fcut, sinc_samples);

% oba okna w czestotliwosci.
fw1 = fft(w1);
fwx = linspace(-fs/2, fs/2, size(fw1, 2));

% funkcja wejsciowa w czestotliwosci.
ffin = fftshift(abs(fft(fin)));
ffinx = linspace(-fs/2, fs/2, size(ffin, 2));

% przyciecie poprzez splot
fout2 = conv(fin, w1, 'same');

figure(1);
hold on;
plot(fwx, fftshift(abs(fw1)), 'r');
legend('sinc');
title('filtr dolnoprzepustowy w przestrzeni czestotliwosci');

figure(2);
hold on;
plot(fin, 'c');
plot(fout);
plot(fout2, 'r');
title('porowanie w czasie');
legend('wejscie', 'idealne wyjscie', 'przyciecie przez splot');

figure(3);
hold on;
plot(ffinx, fftshift(abs(fft(fin))), 'c');
plot(ffinx, fftshift(abs(fft(fout))));
plot(ffinx, fftshift(abs(fft(fout2))), 'r');
title('porowanie w czestotliwosci');
legend('wejscie', 'idealne wyjscie', 'przyciecie przez splot');

Po odkomentowaniu linii zastępujących nasze funkcje sincami, wyniki przedstawiają się tak:

Test naszych funkcji:

% przetestujmy poprawnosc dzialania naszej funkcji conv_cut.

a = [1,2,3,4];
b = [1,2,3,4,5];
c = [1,2,3,4,5,6];
d = [1,2,3,4,5,6,7];

x11 = conv(a, a, 'same');
x12 = conv(a, b, 'same');
x13 = conv(a, c, 'same');
x14 = conv(a, d, 'same');

x21 = conv(b, a, 'same');
x22 = conv(b, b, 'same');
x23 = conv(b, c, 'same');
x24 = conv(b, d, 'same');

x31 = conv(c, a, 'same');
x32 = conv(c, b, 'same');
x33 = conv(c, c, 'same');
x34 = conv(c, d, 'same');

x41 = conv(d, a, 'same');
x42 = conv(d, b, 'same');
x43 = conv(d, c, 'same');
x44 = conv(d, d, 'same');

y11 = conv_same(a, a);
y12 = conv_same(a, b);
y13 = conv_same(a, c);
y14 = conv_same(a, d);

y21 = conv_same(b, a);
y22 = conv_same(b, b);
y23 = conv_same(b, c);
y24 = conv_same(b, d);

y31 = conv_same(c, a);
y32 = conv_same(c, b);
y33 = conv_same(c, c);
y34 = conv_same(c, d);

y41 = conv_same(d, a);
y42 = conv_same(d, b);
y43 = conv_same(d, c);
y44 = conv_same(d, d);

test = @(a,b) assert((size(a,2) == size(b,2)) && (sum(abs(a - b)) == 0));

test(x11, y11);
test(x12, y12);
test(x13, y13);
test(x14, y14);

test(x21, y21);
test(x22, y22);
test(x23, y23);
test(x24, y24);

test(x31, y31);
test(x32, y32);
test(x33, y33);
test(x34, y34);

test(x41, y41);
test(x42, y42);
test(x43, y43);
test(x44, y44);

% przetestujmy poprawnosc dzialania naszej funkcji conv_full.

a = [1,2,3,4];
b = [1,2,3,4,5];
c = [1,2,3,4,5,6];
d = [1,2,3,4,5,6,7];

x11 = conv(a, a, 'full');
x12 = conv(a, b, 'full');
x13 = conv(a, c, 'full');
x14 = conv(a, d, 'full');

x21 = conv(b, a, 'full');
x22 = conv(b, b, 'full');
x23 = conv(b, c, 'full');
x24 = conv(b, d, 'full');

x31 = conv(c, a, 'full');
x32 = conv(c, b, 'full');
x33 = conv(c, c, 'full');
x34 = conv(c, d, 'full');

x41 = conv(d, a, 'full');
x42 = conv(d, b, 'full');
x43 = conv(d, c, 'full');
x44 = conv(d, d, 'full');

y11 = conv_full(a, a);
y12 = conv_full(a, b);
y13 = conv_full(a, c);
y14 = conv_full(a, d);

y21 = conv_full(b, a);
y22 = conv_full(b, b);
y23 = conv_full(b, c);
y24 = conv_full(b, d);

y31 = conv_full(c, a);
y32 = conv_full(c, b);
y33 = conv_full(c, c);
y34 = conv_full(c, d);

y41 = conv_full(d, a);
y42 = conv_full(d, b);
y43 = conv_full(d, c);
y44 = conv_full(d, d);

test = @(a,b) assert((size(a,2) == size(b,2)) && (sum(abs(a - b)) == 0));

test(x11, y11);
test(x12, y12);
test(x13, y13);
test(x14, y14);

test(x21, y21);
test(x22, y22);
test(x23, y23);
test(x24, y24);

test(x31, y31);
test(x32, y32);
test(x33, y33);
test(x34, y34);

test(x41, y41);
test(x42, y42);
test(x43, y43);
test(x44, y44);

%% przetestujmy poprawnosc dzialania naszej funkcji conv_point.

a = [1,2,3,4];
b = [1,2,3,4,5];
c = [1,2,3,4,5,6];
d = [1,2,3,4,5,6,7];

x11 = conv(a, a, 'same');
x12 = conv(a, b, 'same');
x13 = conv(a, c, 'same');
x14 = conv(a, d, 'same');

x21 = conv(b, a, 'same');
x22 = conv(b, b, 'same');
x23 = conv(b, c, 'same');
x24 = conv(b, d, 'same');

x31 = conv(c, a, 'same');
x32 = conv(c, b, 'same');
x33 = conv(c, c, 'same');
x34 = conv(c, d, 'same');

x41 = conv(d, a, 'same');
x42 = conv(d, b, 'same');
x43 = conv(d, c, 'same');
x44 = conv(d, d, 'same');

y11 = zeros(1, size(a, 2));
y12 = zeros(1, size(a, 2));
y13 = zeros(1, size(a, 2));
y14 = zeros(1, size(a, 2));
for i=1:1:size(a, 2)
y11(i) = conv_point(a, a, i);
y12(i) = conv_point(a, b, i);
y13(i) = conv_point(a, c, i);
y14(i) = conv_point(a, d, i);
end

y21 = zeros(1, size(b, 2));
y22 = zeros(1, size(b, 2));
y23 = zeros(1, size(b, 2));
y24 = zeros(1, size(b, 2));
for i=1:1:size(b, 2)
y21(i) = conv_point(b, a, i);
y22(i) = conv_point(b, b, i);
y23(i) = conv_point(b, c, i);
y24(i) = conv_point(b, d, i);
end

y31 = zeros(1, size(c, 2));
y32 = zeros(1, size(c, 2));
y33 = zeros(1, size(c, 2));
y34 = zeros(1, size(c, 2));
for i=1:1:size(c, 2)
y31(i) = conv_point(c, a, i);
y32(i) = conv_point(c, b, i);
y33(i) = conv_point(c, c, i);
y34(i) = conv_point(c, d, i);
end

y41 = zeros(1, size(d, 2));
y42 = zeros(1, size(d, 2));
y43 = zeros(1, size(d, 2));
y44 = zeros(1, size(d, 2));
for i=1:1:size(d, 2)
y41(i) = conv_point(d, a, i);
y42(i) = conv_point(d, b, i);
y43(i) = conv_point(d, c, i);
y44(i) = conv_point(d, d, i);
end

test = @(a,b) assert((size(a,2) == size(b,2)) && (sum(abs(a - b)) == 0));

test(x11, y11);
test(x12, y12);
test(x13, y13);
test(x14, y14);

test(x21, y21);
test(x22, y22);
test(x23, y23);
test(x24, y24);

test(x31, y31);
test(x32, y32);
test(x33, y33);
test(x34, y34);

test(x41, y41);
test(x42, y42);
test(x43, y43);
test(x44, y44);

Funkcja pomocnicza:

function [ k, t ] = sinc_gen( sampling_rate, cutfq, samples )

T = 1 / sampling_rate;

if (nargin == 2)
samples = 160 / (cutfq*2 / sampling_rate); % experimental
end

tout = T * (1:samples);
tout = tout - T*(samples+1)/2;
B = 2*cutfq;
k = sinc(B*tout);
k = k / trapz(k);

if (nargout == 2)
t = tout;
end

end