Как реализовать цифровой генератор?

У меня есть система цифровой обработки сигналов с плавающей запятой, которая работает с фиксированной частотой дискретизации отсчетов в секунду, реализованная с использованием процессора x86-64. Предполагая, что система DSP синхронно связана с любыми вопросами, каков наилучший способ реализовать цифровой генератор на некоторой частоте ? $f_s = 32768$ $f$

В частности, я хочу сгенерировать сигнал: где для номера выборки .

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

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

t = n / f_{s}

$t=n/f_s$

n

$n$

Одна идея состоит в том, чтобы отслеживать вектор который мы поворачиваем на угол на каждом тактовом цикле. $(x,y)$ $\Delta\phi = 2\pi f/f_s$

В качестве реализации псевдокода Matlab (реальная реализация находится в 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)]

Затем в каждом такте мы немного вращаем вектор:

%% 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

Это позволяет вычислять генератор только с 4 умножениями за цикл. Тем не менее, я бы беспокоился о фазовой ошибке и стабильности амплитуды. (В простых тестах я был удивлен, что амплитуда не умерла или взорвалась немедленно - возможно, sincosинструкция гарантирует ?). $\sin^2+\cos^2=1$

Как правильно это сделать?

software-implementation oscillator nibot
источник

Ответы:

Вы правы в том, что строго рекурсивный подход уязвим для накопления ошибок при увеличении числа итераций. Еще один надежный способ, как правило, это сделать, это использовать генератор с числовым управлением (NCO) . По сути, у вас есть аккумулятор, который отслеживает мгновенную фазу генератора, обновленный следующим образом:

δ = \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$

Таким образом, в каждый момент времени вы преобразуете накопленную фазу в NCO в желаемые синусоидальные выходы. То, как вы это сделаете, зависит от ваших требований к вычислительной сложности, точности и т. Д. Одним из очевидных способов является просто рассчитать выходные данные как

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

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

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

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

используя любую имеющуюся реализацию синуса / косинуса. В высокопроизводительных и / или встраиваемых системах отображение значений фазы в синус / косинус часто выполняется с помощью таблицы поиска. Размер таблицы поиска (т. Е. Количество квантования, которое вы выполняете для аргумента фазы для синуса и косинуса) можно использовать в качестве компромисса между потреблением памяти и ошибкой аппроксимации. Приятно то, что количество необходимых вычислений, как правило, не зависит от размера таблицы. Кроме того, при необходимости вы можете ограничить размер LUT, используя симметрию, присущую функциям косинуса и синуса; вам действительно нужно хранить только одну четвертую периода выбранной синусоиды.

Если вам нужна более высокая точность, чем может дать LUT разумного размера, тогда вы всегда можете посмотреть на интерполяцию между выборками таблицы (например, линейную или кубическую интерполяцию).

Другое преимущество этого подхода состоит в том, что в эту структуру легко включить частотную или фазовую модуляцию. Частота выходного сигнала может быть модулирована путем изменения соответственно, а фазовая модуляция может быть реализована простым добавлением в напрямую. $\delta$ $\phi[n]$

Джейсон Р
источник

Спасибо за ответ. Как время выполнения sincosсравнивается с несколькими умножениями? Есть ли какие-либо возможные подводные камни, на которые стоит обратить внимание во время modоперации?

нибот

Приятно, что одну и ту же фазу-амплитуду можно использовать для всех генераторов в системе.

нибот

Какова цель мода 2pi? Я также видел реализации, которые делают мод 1.0. Можете ли вы рассказать о том, для чего нужна операция по модулю?

BigBrownBear00

@ BigBrownBear00: операция по модулю - это то, что поддерживает в управляемом диапазоне. На практике, если бы у вас не было модуля по модулю, оно со временем стало бы очень большим положительным или отрицательным числом (общим количеством накопленной фазы). Это может быть плохо по нескольким причинам, включая возможное переполнение или потерю арифметической точности, а также снижение производительности при оценке функций косинуса и синуса. Типичные реализации выполняются быстрее, если им не нужно сначала сокращать аргументы до диапазона .

ϕ [n]

$\phi[n]$

[0, 2 π)

$[0, 2\pi)$

Джейсон Р

Коэффициент против 1,0 - это деталь реализации. Это зависит от области тригонометрических функций вашей платформы. Если они ожидают значение в диапазоне (т. Е. Угол измеряется в циклах), тогда уравнение для будет скорректировано с учетом этой отличающейся единицы. Объяснение в ответе выше предполагает, что используется типичная угловая единица радиан.

2 π

$2\pi$

[0, 1.0)

$[0, 1.0)$

ϕ [n]

$\phi[n]$

Джейсон Р

То, что у вас есть, это очень хороший и эффективный генератор. Проблема потенциального численного дрейфа действительно может быть решена. Ваша переменная состояния v состоит из двух частей, одна из которых является реальной частью, а другая - мнимой частью. Давайте назовем тогда г и я. Мы знаем, что r ^ 2 + i ^ 2 = 1. Со временем это может сместиться вверх и вниз, однако это можно легко исправить путем умножения с использованием коэффициента коррекции усиления, например

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

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

Очевидно, что это очень дорого, однако мы знаем, что коррекция усиления очень близка к единице, и мы можем аппроксимировать это простым разложением Тейлора на

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 )$

Более того, нам не нужно делать это на каждом отдельном образце, но один раз на каждые 100 или 1000 образцов более чем достаточно для поддержания стабильности. Это особенно полезно, если вы делаете обработку на основе кадров. Обновление один раз за кадр просто отлично. Вот быстрый Matlab рассчитывает 10 000 000 образцов.

%% 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
источник

Этот ответ дополнительно объясняется Хилмаром в другом вопросе: dsp.stackexchange.com/a/1087/34576

sircolinton

Вы можете избежать нестабильного смещения величины, если не заставите его рекурсивно обновлять вектор v. Вместо этого поверните вектор-прототип v к текущей фазе вывода. Это все еще требует некоторых функций триггера, но только один раз для каждого буфера.

Нет дрейфа величины и произвольной частоты

псевдокод:

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

Вы можете избавиться от умножения, функций триггера, требуемых для cexp, и от остатка модуля более 2pi, если вы можете допустить квантованное преобразование частоты. например, fs / 1024 для векторного буфера 1024 образца.

Марк Боргердинг
источник