Comment implémenter un oscillateur numérique?

20

J'ai un système de traitement du signal numérique à virgule flottante qui fonctionne à un taux d'échantillonnage fixe de fs=32768 échantillons par seconde implémenté à l'aide d'un processeur x86-64. En supposant que le système DSP est verrouillé de manière synchrone sur tout ce qui compte, quelle est la meilleure façon de mettre en œuvre un oscillateur numérique à une fréquence f ?

Plus précisément, je veux générer le signal:

y(t)=sin(2πft)
t=n/fs pour le numéro d'échantillon n .

Une idée est de garder une trace d'un vecteur (x,y) que nous faisons tourner d'un angle Δϕ=2πf/fs à chaque cycle d'horloge.

En tant qu'implémentation de pseudocode Matlab (l'implémentation réelle est en 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)]

Ensuite, à chaque cycle d'horloge, nous faisons tourner un peu le vecteur:

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

Cela permet à l'oscillateur d'être calculé avec seulement 4 multiplications par cycle. Cependant, je m'inquiéterais de l'erreur de phase et de la stabilité de l'amplitude. (Dans des tests simples, j'ai été surpris que l'amplitude ne meure pas ou n'explose pas immédiatement - peut-être que l' sincosinstruction garantit sin2+cos2=1 ?).

Quel est le bon moyen de le faire?

nibot
la source

Réponses:

12

Vous avez raison de dire que l'approche strictement récursive est vulnérable à l'accumulation d'erreur à mesure que le nombre d'itérations augmente. Un moyen plus robuste pour cela est généralement d'utiliser un oscillateur à commande numérique (NCO) . Fondamentalement, vous disposez d'un accumulateur qui garde la trace de la phase instantanée de l'oscillateur, mis à jour comme suit:

δ=2πffs

ϕ[n]=(ϕ[n1]+δ)mod2π

À chaque instant, il vous reste donc à convertir la phase accumulée dans le NCO en sorties sinusoïdales souhaitées. La façon dont vous procédez dépend de vos besoins en termes de complexité de calcul, de précision, etc. Une façon évidente consiste à calculer simplement les sorties comme

xc[n]=cos(ϕ[n])

xs[n]=sin(ϕ[n])

en utilisant la mise en œuvre de sinus / cosinus dont vous disposez. Dans les systèmes à haut débit et / ou embarqués, le mappage de la phase aux valeurs sinus / cosinus se fait souvent via une table de recherche. La taille de la table de recherche (c'est-à-dire la quantité de quantification que vous effectuez sur l'argument de phase en sinus et cosinus) peut être utilisée comme compromis entre la consommation de mémoire et l'erreur d'approximation. La bonne chose est que la quantité de calculs requis est généralement indépendante de la taille de la table. De plus, vous pouvez limiter la taille de votre LUT si nécessaire en profitant de la symétrie inhérente aux fonctions cosinus et sinus; il vous suffit de stocker un quart de période de la sinusoïde échantillonnée.

Si vous avez besoin d'une précision supérieure à ce qu'une LUT de taille raisonnable peut vous donner, vous pouvez toujours regarder l'interpolation entre les échantillons de table (interpolation linéaire ou cubique, par exemple).

Un autre avantage de cette approche est qu'il est trivial d'incorporer une modulation de fréquence ou de phase à cette structure. La fréquence de la sortie peut être modulée en variant conséquence, et la modulation de phase peut être implémentée en ajoutant simplement à directement.δϕ[n]

Jason R
la source
2
Merci d'avoir répondu. Comment le temps d'exécution se sincoscompare-t-il à une poignée de multiplications? Y a-t-il des pièges possibles à surveiller lors de l' modopération?
nibot
Il est intéressant de noter que la même LUT phase-amplitude peut être utilisée pour tous les oscillateurs du système.
nibot
Quel est le but du mod 2pi? J'ai également vu des implémentations qui font le mod 1.0. Pouvez-vous développer à quoi sert l'opération modulo?
BigBrownBear00
1
@ BigBrownBear00: L'opération modulo est ce qui maintient dans une plage gérable. En pratique, si vous n'aviez pas le modulo, il deviendrait un très grand nombre positif ou négatif (la quantité totale de phase accumulée) au fil du temps. Cela peut être mauvais pour plusieurs raisons, y compris un éventuel débordement ou une perte de précision arithmétique, et une réduction des performances des évaluations des fonctions cosinus et sinus. Les implémentations typiques sont plus rapides si elles n'ont pas à effectuer d'abord une réduction d'argument dans la plage . ϕ[n][0,2π)
Jason R
1
Le facteur versus 1.0 est un détail d'implémentation. Cela dépend du domaine des fonctions trigonométriques de votre plateforme. S'ils s'attendent à une valeur dans la plage (c'est-à-dire que l'angle est mesuré en cycles), alors l'équation pour serait ajustée pour refléter cette unité différente. L'explication de la réponse ci-dessus suppose que l'unité angulaire typique de radians est utilisée. 2π[0,1.0)ϕ[n]
Jason R
8

Ce que vous avez est un oscillateur très bon et efficace. Le problème potentiel de dérive numérique peut en fait être résolu. Votre variable d'état v comporte deux parties, l'une est essentiellement la partie réelle et l'autre la partie imaginaire. Appelons alors r et i. Nous savons que r ^ 2 + i ^ 2 = 1. Au fil du temps, cela peut dériver de haut en bas, mais cela peut facilement être corrigé par multiplication avec un facteur de correction de gain comme celui-ci

g=1r2+i2

Évidemment, cela coûte très cher, mais nous savons que la correction de gain est très proche de l'unité et nous pouvons l'approcher avec une simple expansion de Taylor à

g=1r2+i212(3(r2+i2))

De plus, nous n'avons pas besoin de le faire sur chaque échantillon, mais une fois tous les 100 ou 1000 échantillons est plus que suffisant pour maintenir cette stabilité. Ceci est particulièrement utile si vous effectuez un traitement basé sur les trames. La mise à jour une fois par image est très bien. Voici un rapide Matlab calcule 10 000 000 d'échantillons.

%% 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
la source
Cette réponse est expliquée plus en détail par Hilmar dans une autre question: dsp.stackexchange.com/a/1087/34576
sircolinton
7

Vous pouvez éviter la dérive d'amplitude instable si vous ne le faites pas mettre à jour récursivement le vecteur v. Au lieu de cela, faites pivoter votre vecteur prototype v vers la phase de sortie actuelle. Cela nécessite toujours certaines fonctions trigonométriques, mais une seule fois par tampon.

Aucune dérive d'amplitude et fréquence arbitraire

pseudocode:

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

Vous pouvez supprimer la multiplication, les fonctions trigonométriques requises par cexp et le module restant sur 2pi si vous pouvez tolérer une translation de fréquence quantifiée. par exemple, fs / 1024 pour un tampon de phaseur de 1024 échantillons.

Mark Borgerding
la source