Jak wdrożyć oscylator cyfrowy?

Mam zmiennoprzecinkowy system cyfrowego przetwarzania sygnałów, który działa ze stałą częstotliwością próbkowania próbek na sekundę, zaimplementowany przy użyciu procesora x86-64. Zakładając, że system DSP jest synchronicznie zablokowany na czymkolwiek, co jest najlepszym sposobem na wdrożenie oscylatora cyfrowego przy jakiejś częstotliwości ? $f_s = 32768$ $f$

W szczególności chcę wygenerować sygnał: gdzie dla próbki o numerze .

y (t) = \sin (2 π f t)

$y(t) = \sin(2\pi f t)$

t = n / f_{s}

$t=n/f_s$

n

$n$

Jednym z pomysłów jest śledzenie wektora który obracamy o kąt w każdym cyklu zegara. $(x,y)$ $\Delta\phi = 2\pi f/f_s$

Jako implementacja pseudokodu Matlab (rzeczywista implementacja znajduje się w C):

%% Initialization code

f_s = 32768;             % sample rate [Hz]
f = 19.875;              % some constant frequency [Hz]

v = [1 0];               % initial condition     
d_phi = 2*pi * f / f_s;  % change in angle per clock cycle

% initialize the rotation matrix (only once):
R = [cos(d_phi), -sin(d_phi) ; ...
     sin(d_phi),  cos(d_phi)]

Następnie, w każdym cyklu zegara, obracamy nieco wektor wokół:

%% in-loop code

while (forever),
  v = R*v;        % rotate the vector by d_phi
  y = v(1);       % this is the sine wave we're generating
  output(y);
end

Pozwala to na obliczenie oscylatora tylko przy 4 multiplikacjach na cykl. Martwię się jednak o błąd fazy i stabilność amplitudy. (W prostych testach zdziwiłem się, że amplituda nie umarła ani nie wybuchła natychmiast - być może sincosinstrukcja gwarantuje ?). $\sin^2+\cos^2=1$

Jak to zrobić?

software-implementation oscillator

— nibot
źródło

Odpowiedzi:

Masz rację, że podejście ściśle rekurencyjne jest podatne na gromadzenie się błędów w miarę wzrostu liczby iteracji. Jednym z bardziej niezawodnych sposobów, w jaki jest to zwykle wykonywane, jest użycie sterowanego numerycznie oscylatora (NCO) . Zasadniczo masz akumulator, który śledzi chwilową fazę oscylatora, zaktualizowany w następujący sposób:

δ = \frac{2 π f}{f_{s}}

$\delta = \frac{2\pi f}{f_s}$

ϕ [n] = (ϕ [n - 1] + δ) \mod 2 π

$\phi[n] = (\phi[n-1] + \delta) \mod 2\pi$

W każdej chwili pozostaje Ci zatem zamiana zakumulowanej fazy w NCO na pożądane sinusoidalne wyjścia. To, jak to zrobisz, zależy od wymagań dotyczących złożoności obliczeniowej, dokładności itp. Jednym z oczywistych sposobów jest po prostu obliczenie wyników jako

x_{c} [n] = \cos (ϕ [n])

$x_c[n] = \cos(\phi[n])$

x_{s} [n] = \sin (ϕ [n])

$x_s[n] = \sin(\phi[n])$

używając dowolnej implementacji sinus / cosinus, jaką masz dostępną. W systemach o dużej przepustowości i / lub wbudowanych mapowanie z wartości fazowych na wartości sinusoidalne / cosinusowe jest często wykonywane za pomocą tabeli przeglądowej. Rozmiar tabeli odnośników (tj. Wielkość kwantyzacji, którą wykonujesz w argumencie fazowym dla sinusa i cosinusa) może być wykorzystany jako kompromis między zużyciem pamięci a błędem aproksymacji. Zaletą jest to, że ilość wymaganych obliczeń jest zwykle niezależna od wielkości tabeli. Ponadto, w razie potrzeby, możesz ograniczyć swój rozmiar LUT, korzystając z symetrii właściwej funkcji cosinus i sinus; naprawdę musisz przechowywać tylko jedną czwartą okresu próbkowanej sinusoidy.

Jeśli potrzebujesz większej dokładności niż LUT o rozsądnej wielkości, zawsze możesz spojrzeć na interpolację między próbkami tabeli (na przykład interpolacja liniowa lub sześcienna).

Inną zaletą tego podejścia jest to, że włączenie modulacji częstotliwości lub fazy do tej struktury jest banalne. Częstotliwość wyjściową można modulować odpowiednio zmieniając , a modulację fazową można wdrożyć, po prostu dodając bezpośrednio do . $\delta$ $\phi[n]$

— Jason R.
źródło

Dziękuję za odpowiedź. Jak czas wykonania sincosróżni się od garstki mnożeń? Czy są jakieś pułapki, na które należy uważać podczas modoperacji?

— nibot

To atrakcyjne, że dla wszystkich oscylatorów w systemie można zastosować tę samą LUT między fazami i amplitudą.

— nibot

Jaki jest cel modu 2pi? Widziałem także implementacje, które wykonują mod 1.0. Czy możesz rozwinąć, do czego służy operacja modulo?

— BigBrownBear00

@ BigBrownBear00: Operacja modulo utrzymuje w możliwym do zarządzania zakresie. W praktyce, jeśli nie miałbyś modulo, z czasem stałby się on bardzo dużą liczbą dodatnią lub ujemną (łączna ilość skumulowanej fazy). Może to być złe z kilku powodów, w tym z ewentualnego przepełnienia lub utraty precyzji arytmetycznej oraz zmniejszonej wydajności ocen funkcji cosinus i sinus. Typowe implementacje są szybsze, jeśli nie muszą najpierw przeprowadzać redukcji argumentów w zakresie .

ϕ [n]

$\phi[n]$

[0, 2 π)

$[0, 2\pi)$

— Jason R

Współczynnik porównaniu do 1,0 to szczegół implementacji. Zależy to od dziedziny funkcji trygonometrycznych platformy. Jeśli oczekują wartości z zakresu (tj. Kąt jest mierzony w cyklach), wówczas równanie dla zostanie dostosowane tak, aby odzwierciedlało tę różną jednostkę. Wyjaśnienie w powyższej odpowiedzi zakłada, że użyto typowej jednostki kątowej radianów.

2 π

$2\pi$

[0, 1.0)

$[0, 1.0)$

ϕ [n]

$\phi[n]$

— Jason R

To, co masz, to bardzo dobry i wydajny oscylator. Potencjalny problem dryfu numerycznego można faktycznie rozwiązać. Twoja zmienna stanu v składa się z dwóch części, z których jedna jest częścią rzeczywistą, a druga częścią urojoną. Zadzwońmy wtedy ri ja. Wiemy, że r ^ 2 + i ^ 2 = 1. Z czasem może to dryfować w górę i w dół, jednak można to łatwo skorygować przez pomnożenie przez współczynnik korekcji wzmocnienia taki jak

g = \frac{1}{\sqrt{r^{2} + i^{2}}}

$g = \frac{1}{\sqrt{r^{2}+i^{2}}}$

Oczywiście jest to bardzo kosztowne, jednak wiemy, że korekcja wzmocnienia jest bardzo bliska jedności i możemy to przybliżyć za pomocą zwykłego rozszerzenia Taylora do

g = \frac{1}{\sqrt{r^{2} + i^{2}}} \approx \frac{1}{2} \cdot (3 - (r^{2} + i^{2}))

$g = \frac{1}{\sqrt{r^{2}+i^{2}}}\approx \frac{1}{2}\cdot \left ( 3-(r^{2}+i^{2}) \right )$

Co więcej, nie musimy tego robić na każdej pojedynczej próbce, ale raz na każde 100 lub 1000 próbek wystarczy, aby utrzymać tę stabilność. Jest to szczególnie przydatne, jeśli wykonujesz przetwarzanie oparte na ramkach. Aktualizacja raz na klatkę jest w porządku. Oto szybki Matlab oblicza 10 000 000 próbek.

%% seed the oscillator
% set parameters
f0 = single(100); % say 100 Hz
fs = single(44100); % sample rate = 44100;
nf = 1024; % frame size

% initialize phasor and state
ph =  single(exp(-j*2*pi*f0/fs));
state = single(1 + 0i); % real part 1, imaginary part 0

% try it
x = zeros(nf,1,'single');
testRuns = 10000;
for k = 1:testRuns
  % overall frames
  % sample: loop
  for i= 1:nf
    % phasor multiply
    state = state *ph;
    % take real part for cosine, or imaginary for sine
    x(i) = real(state);
  end
  % amplitude corrections through a taylor exansion aroud
  % abs(state) very close to 1
  g = single(.5)*(single(3)-real(state)*real(state)-imag(state)*imag(state) );
  state = state*g;
end
fprintf('Deviation from unity amplitude = %f\n',g-1);

— Hilmar
źródło

Hilmar wyjaśnia tę odpowiedź w innym pytaniu: dsp.stackexchange.com/a/1087/34576

— sircolinton

Możesz uniknąć niestabilnego dryfu wielkości, jeśli nie sprawisz, że rekursywnie zaktualizujesz wektor v. Zamiast tego obróć prototypowy wektor v do bieżącej fazy wyjściowej. Nadal wymaga to niektórych funkcji wyzwalania, ale tylko raz na bufor.

Brak dryfu wielkości i arbitralnej częstotliwości

pseudo kod:

init(freq)
  precompute Nphasor samples in phasor
  phase=0

gen(Nsamps)
    done=0
    while done < Nsamps:
       ndo = min(Nsamps -done, Nphasor)
       append to output : multiply buf[done:done+ndo) by cexp( j*phase )
       phase = rem( phase + ndo * 2*pi*freq/fs,2*pi)
       done = done+ndo

Możesz zrezygnować z mnożenia, funkcji wyzwalania wymaganych przez cexp, a moduł pozostanie powyżej 2pi, jeśli możesz tolerować kwantowane tłumaczenie częstotliwości. np. fs / 1024 dla 1024 próbnego bufora fazorów.

— Mark Borgerding
źródło