[time-nuts] Re: Simple simulation model for an OCXO?

[Carsten Andrich]

On 14.05.22 16:58, Matthias Welwarsky wrote:
> I went "window shopping" on Google and found something that would probably fit
> my needs here:
>
> https://ccrma.stanford.edu/~jos/sasp/Example_Synthesis_1_F_Noise.html
>
> Matlab code:
>
> Nx = 2^16;  % number of samples to synthesize
> B = [0.049922035 -0.095993537 0.050612699 -0.004408786];
> A = [1 -2.494956002   2.017265875  -0.522189400];
> nT60 = round(log(1000)/(1-max(abs(roots(A))))); % T60 est.
> v = randn(1,Nx+nT60); % Gaussian white noise: N(0,1)
> x = filter(B,A,v);    % Apply 1/F roll-off to PSD
> x = x(nT60+1:end);    % Skip transient response

Hi Matthias,

the coefficients could come the digital transform of an analog filter 
that approximates a 10 dB/decade slope, see:
https://electronics.stackexchange.com/questions/200583/is-it-possible-to-make-a-weak-analog-low-pass-filter-with-less-than-20-db-decade
https://electronics.stackexchange.com/questions/374247/is-the-roll-off-gain-of-filters-always-20-db-dec
https://m.youtube.com/watch?v=fn2UGyk5DP4

However, even for the 2^16 samples used by the CCRMA snippet, the filter 
slope rolls off too quickly. I've attached its frequency response. It 
exhibits a little wobbly 1/f power slope over 3 orders of magnitude, but 
it's essentially flat over the remaining two orders of mag. The used IIR 
filter is too short to affect the lower frequencies.

If you want a good approximation of 1/f over the full frequency range, I 
personally believe Greenhall's "discrete spectrum algorithm" [1] to be 
the best choice, as per the literature I've consulted so far. It 
realizes an appropriately large FIR filter, i.e., as large as the 
input/output signal, and it's straightforward to implement. Because it's 
generic, it can be used to generate any combination of power-law noises 
or arbitrary characteristics, e.g., phase noise based on measurements, 
in a single invocation.

Best regards,
Carsten

[1] https://apps.dtic.mil/sti/pdfs/ADA485683.pdf#page=5

P.S.: Python Code to plot attached frequency response:

import matplotlib.pyplot as plt
import numpy as np
import scipy.signal as signal

w, h = signal.freqz(
     [0.049922035, -0.095993537, 0.050612699, -0.004408786],
     [1, -2.494956002,   2.017265875,  -0.522189400],
     2**16)

plt.loglog(w/np.pi, abs(h))
plt.plot([1e-3, 1e-0], [0.93, 0.93/10**(0.5*3)], "--")

plt.xlabel("Normalized frequency ($f_\mathrm{Nyquist}$)")
plt.ylabel("Gain")
plt.grid()
-------------- next part --------------
A non-text attachment was scrubbed...
Name: H.png
Type: image/png
Size: 21104 bytes
Desc: not available
URL: <http://febo.com/pipermail/time-nuts_lists.febo.com/attachments/20220514/a8e3aaab/attachment.png>