Como implementar um oscilador digital?

Eu tenho um sistema de processamento de sinal digital de ponto flutuante que opera a uma taxa de amostragem fixa de amostras por segundo implementada usando um processador x86-64. Supondo que o sistema DSP esteja sincronizado com o que quer que seja importante, qual é a melhor maneira de implementar um oscilador digital em alguma frequência ?$f_s = 32768$$f$

Especificamente, quero gerar o sinal: que para o número da amostra .

Uma idéia é acompanhar um vetor que giramos em um ângulo em cada ciclo do relógio.$(x,y)$$\Delta\phi = 2\pi f/f_s$

Como uma implementação de pseudocódigo do Matlab (a implementação real está em 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)]

Em seguida, em cada ciclo do relógio, rotacionamos o vetor um pouco:

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

Isso permite que o oscilador seja calculado com apenas 4 multiplicações por ciclo. No entanto, eu me preocuparia com erros de fase e estabilidade de amplitude. (Em testes simples, fiquei surpreso que a amplitude não tenha morrido ou explodido imediatamente - talvez a sincosinstrução garanta ?).$\sin^2+\cos^2=1$

Qual é a maneira certa de fazer isso?

Você está certo de que a abordagem estritamente recursiva é vulnerável ao acúmulo de erros à medida que o número de iterações aumenta. Uma maneira mais robusta de fazer isso é usar um oscilador controlado numericamente (NCO) . Basicamente, você tem um acumulador que monitora a fase instantânea do oscilador, atualizada da seguinte forma:

A cada instante, resta a conversão da fase acumulada no NCO para as saídas sinusoidais desejadas. Como você faz isso depende de seus requisitos de complexidade computacional, precisão etc. Uma maneira óbvia é apenas calcular as saídas como

usando qualquer implementação de seno / cosseno que você tiver disponível. Em sistemas de alto rendimento e / ou sistemas embarcados, o mapeamento dos valores de fase para seno / cosseno geralmente é feito por meio de uma tabela de pesquisa. O tamanho da tabela de pesquisa (ou seja, a quantidade de quantização realizada no argumento de fase para seno e cosseno) pode ser usado como uma troca entre consumo de memória e erro de aproximação. O bom é que a quantidade de cálculos necessária geralmente é independente do tamanho da tabela. Além disso, você pode limitar o tamanho da LUT, se necessário, aproveitando a simetria inerente às funções cosseno e seno; você realmente só precisa armazenar um quarto de um período do sinusóide amostrado.

Se você precisar de uma precisão mais alta do que uma LUT de tamanho razoável pode fornecer, sempre poderá analisar a interpolação entre as amostras da tabela (interpolação linear ou cúbica, por exemplo).

Outro benefício dessa abordagem é que é trivial incorporar modulação de frequência ou fase nessa estrutura. A frequência da saída pode ser modulada variando acordo, e a modulação de fase pode ser implementada simplesmente adicionando diretamente a .$\delta$$\phi[n]$

Além disso, não precisamos fazer isso em todas as amostras, mas uma vez a cada 100 ou 1000 amostras é mais que suficiente para manter essa estabilidade. Isso é particularmente útil se você fizer um processamento baseado em quadros. Atualizar uma vez por quadro é bom. Aqui está um rápido Matlab calcula 10.000.000 amostras.

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

Você pode evitar o desvio de magnitude instável se não o fizer atualizar recursivamente o vetor v. Em vez disso, gire o vetor protótipo v para a fase de saída atual. Isso ainda requer algumas funções trigonométricas, mas apenas uma vez por buffer.

Desvio de magnitude e frequência arbitrária

pseudo-código:

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

Você pode eliminar a multiplicação, as funções trigonométricas exigidas pelo cexp e o módulo restante acima de 2pi, se você puder tolerar uma conversão de frequência quantizada. por exemplo, fs / 1024 para um buffer fasor de amostra de 1024.