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()
_______________________________________________
time-nuts mailing list -- [email protected]
To unsubscribe send an email to [email protected]
