Bagaimana cara mengimplementasikan osilator digital?

Saya memiliki sistem pemrosesan sinyal digital floating-point yang beroperasi pada laju sampel tetap sampel per detik diimplementasikan menggunakan prosesor x86-64. Dengan asumsi bahwa sistem DSP terkunci secara sinkron ke masalah apa pun, apa cara terbaik untuk mengimplementasikan osilator digital pada frekuensi tertentu ?$f_s = 32768$$f$

Secara khusus, saya ingin menghasilkan sinyal: mana untuk jumlah sampel .

Satu ide adalah untuk melacak vektor yang kita putar dengan sudut pada setiap siklus clock.$(x,y)$$\Delta\phi = 2\pi f/f_s$

Sebagai implementasi pseudocode Matlab (implementasi sebenarnya adalah dalam 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)]

Kemudian, pada setiap siklus jam, kami memutar vektor sekitar sedikit:

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

Ini memungkinkan osilator untuk dihitung hanya dengan 4 perkalian per siklus. Namun, saya khawatir tentang kesalahan fase dan stabilitas amplitudo. (Dalam tes sederhana saya terkejut bahwa amplitudo tidak mati atau meledak segera - mungkin sincosinstruksi menjamin ?).$\sin^2+\cos^2=1$

Apa cara yang tepat untuk melakukan ini?

Anda benar bahwa pendekatan ketat rekursif rentan terhadap akumulasi kesalahan karena jumlah iterasi meningkat. Satu cara yang lebih kuat yang biasanya dilakukan adalah menggunakan osilator yang dikontrol secara numerik (NCO) . Pada dasarnya, Anda memiliki akumulator yang melacak fase sesaat dari osilator, diperbarui sebagai berikut:

Maka, pada setiap waktu instan, Anda harus mengubah fase terakumulasi dalam NCO ke output sinusoidal yang diinginkan. Cara Anda melakukan ini tergantung pada persyaratan Anda untuk kompleksitas komputasi, akurasi, dll. Salah satu cara yang jelas adalah dengan hanya menghitung output sebagai

menggunakan implementasi sinus / cosine apa pun yang Anda miliki. Dalam throughput tinggi dan / atau sistem embedded, pemetaan dari fase ke nilai sinus / cosinus sering dilakukan melalui tabel pencarian. Ukuran tabel pencarian (yaitu jumlah kuantisasi yang Anda lakukan pada argumen fase menjadi sinus dan kosinus) dapat digunakan sebagai tradeoff antara konsumsi memori dan kesalahan perkiraan. Yang menyenangkan adalah bahwa jumlah perhitungan yang diperlukan biasanya tidak tergantung pada ukuran tabel. Selain itu, Anda dapat membatasi ukuran LUT Anda jika diperlukan dengan mengambil keuntungan dari simetri yang melekat pada fungsi cosinus dan sinus; Anda hanya perlu menyimpan seperempat dari periode sinusoid sampel.

Jika Anda membutuhkan keakuratan yang lebih tinggi daripada yang bisa diberikan LUT berukuran cukup, maka Anda selalu dapat melihat interpolasi antara sampel tabel (misalnya interpolasi linier atau kubik).

Manfaat lain dari pendekatan ini adalah sepele untuk memasukkan frekuensi atau modulasi fase dengan struktur ini. Frekuensi output dapat dimodulasi dengan memvariasikan sesuai, dan modulasi fase dapat diimplementasikan dengan hanya menambahkan ϕ [ n ] secara langsung.$\delta$$\phi[n]$

Selain itu kami tidak perlu melakukan ini pada setiap sampel tunggal, tetapi sekali setiap 100 atau 1000 sampel lebih dari cukup untuk menjaga ini stabil. Ini sangat berguna jika Anda melakukan pemrosesan berbasis bingkai. Memperbarui sekali per frame tidak masalah. Berikut ini adalah Matlab cepat yang menghitung 10.000.000 sampel.

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

Anda dapat menghindari penyimpangan besarnya tidak stabil jika Anda tidak membuatnya memperbarui vektor v secara rekursif. Sebagai gantinya, putar vektor prototipe v Anda ke fase output saat ini. Ini masih memerlukan beberapa fungsi trigonometri, tetapi hanya sekali per buffer.

Tidak ada frekuensi penyimpangan dan frekuensi arbitrer

kodesemu:

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

Anda dapat menghapus multiply, fungsi trigonometri yang diperlukan oleh cexp, dan modulus lebih dari 2pi jika Anda dapat mentolerir terjemahan frekuensi yang dikuantisasi. misal fs / 1024 untuk buffer fasor sampel 1024.