I am trying to generate a time-domain violet noise signal with the following power spectral density (PSD):
$$ S_n(f) = A^2f^2 $$
Unfortunately, I am having trouble finding the right amplitude coefficient to get the correct value of $A$.
I am generating the signal by:
Creating a white-noise signal array $\mathcal{w}(t)$ with sample frequency $f_s$ and $\sigma = 1$.
Performing numerical differentiation on this signal (which is equivalent to multiplying by $f$ in frequency domain).
Multiplying by $1/f_s$ to renormalize after differentiating.
Multiplying this signal by the root-mean square value of:
$$ \begin{aligned} \bar{\mathcal{v}}_n &= \left(\int_0^{f_s/2} S_n(f) df\right)^{1/2} = \left(\int_0^{f_s/2} A^2f^2 df\right)^{1/2} \\ &= \left(A^2 \frac{1}{3} \left(\frac{f_s}{2}\right)^3 \right)^{1/2} \\ &= \frac{1}{2\sqrt{6}}A{f_s}^{3/2} \end{aligned} $$
so the final expression is:
$$ \mathcal{v}(t) = \bar{\mathcal{v}}_n \frac{1}{f_s}\frac{d\mathcal{w}(t)}{dt} $$
My problem is that the resulting PSD from this signal is off by a factor of $\pi$ (or maybe 3?) with respect to the expected response.
Here is my code in python
:
import numpy as np
from scipy import signal
import allantools as aln
from matplotlib import pyplot as plt
rng = np.random.default_rng()
fs = 10e3 # Sampling freq [Hz]
N = 1e5 # Number of points
A = 1 # Amplitude spectral density coefficient of violet noise [a.u./(Hz^(3/2)]
time = np.arange(N)/fs # time array [s]
vn = np.sqrt(1/3*A**2*(fs/2)**3) # RMS value of signal [a.u.]
# Time-domain violet noise signal
vn_t = vn*np.diff(rng.normal(size=time.shape[0]+1))
# Compute PSD
f, Sn_f = signal.welch(vn_t, fs, nperseg=2048)
plt.loglog(f,Sn_f)
plt.loglog(f,A**2*f**2,'tab:red')
plt.xlabel('frequency [Hz]')
plt.ylabel('PSD [(A.U.)**2/Hz]')
plt.legend(('Simulated PSD','Expected PSD'),loc='lower right')
plt.xlim([1e1,5e3])
plt.grid()
plt.show()
Which results in the following plot:
These are the results if I divide the simulated PSD by $\pi$:
My guess is that I am missing something in the differentiation step, as this is for a Gaussian-distributed random process, but after doing a lot of searching, most of the references I see say that either white-noise signals are non-differentiable, or have just some complicated stochastic differential calculus equations that don't really point to anything practical (like this).
Any help would be greatly appreciated.
Note: I know that I could generate the frequency-domain signal and then perform an ifft to get the time-domain signal of interest, but I am asking this question because I am interested in knowing what would be the correct procedure to generating the time-domain signal directly.
np.diff
really the approximation to differentiation you want to make? Because if I remember correctly, it's justout[i] = in[i]-in[i-1]
, and that might not have the amplitude response you want (i.e., with $f(t)=e^{j\omega t}$, $\left\lvert\frac{\mathrm d}{\mathrm dt} f \right\rvert= \omega$, whereas $\lvert\text{numpy.diff}(f) \rvert= \left\lvert e^{j\omega t}-e^{j\omega t}\cdot e^{-j\omega \Delta t}\right\rvert =\left\lvert e^{j\omega t}\left(1-e^{-j\omega \Delta t}\right)\right\rvert=\left\lvert 1-e^{-j\omega \Delta t}\right\rvert\ne \omega$, so this is far from what it should be.) $\endgroup$signal.welch
and understand the effect of windowing on the scaling. Make sure you get the expected answer for white noise first before trying anything fancy. $\endgroup$signal.welch
as I use it quite often for other purposes (including white noise). Your previous comments are quite interesting though. I think the $\sin(\omega T/2)$ might be where the key to the problem might be, but I need to process your comments carefully before drawing conclusions. Thanks again. $\endgroup$out[i] = in[i]-in[i-1]
) corresponds to this transfer function in frequency domain: $H(\omega) = \frac{1}{2j} e^{-j\omega T/2}\cdot sin(\omega T/2)$? I am fairly inexperienced in signal processing. Thanks again. $\endgroup$