Generating violet noise with a specific PSD coefficient

Question

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:

1. Creating a white-noise signal array $\mathcal{w}(t)$ with sample frequency $f_s$ and $\sigma = 1$.

2. Performing numerical differentiation on this signal (which is equivalent to multiplying by $f$ in frequency domain).

3. Multiplying by $1/f_s$ to renormalize after differentiating.

4. 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.

Answer 1

Upgraded to full answer.

The diff function implements the difference equation

$$y[n] = x[n]-x[n-1]$$

The transfer function is simply

$$H(z) = 1 - z^{-1}$$ or

$$H(\omega) = 1 - e^{-j\omega}$$

where $\omega$ is the normalized frequency. We can write this as

$$H(\omega) = e^{-j\omega /2} \left( e^{+j\omega/2} - e^{-j\omega/2}\right) = e^{-j\omega /2} \cdot 2 j \cdot \sin(\omega/2) $$

(sorry, I had the factor $2j$ inversed in my original comment).

There are a few things to note here: The linear phase term $ e^{-j\omega/2}$ is equivalent to a half sample delay. That is caused by the fact that the difference is centered around $n = 1/2$ and not around $n = 0$. If you estimate the derivative as $y[n] = x[n+1]-x[n-1]$, that problem would go away (but make other things worse).

For small frequencies we can use $\sin(x) \approx x$ and we'd get

$$H(\omega) \approx j\omega e^{-j\omega/2}$$

which matches the continuous derivative other than the half sample delay.

At higher frequencies you will run into some sort of aliasing. In order to sample a signal without loss, the signal needs to be bandlimited, which is not the case. When you represent a signal in a computer as an array of numbers, it's discrete, and if it's discrete in one domain it's periodic in the other. Hence, it flattens out at the Nyquist frequency: the frequency domain periodicity enforces that (which is exactly what aliasing is).

EDIT: matching the amplitudes

I think your goal is to energy-match the signal before and after the spectral shaping. So if $y[n] = A \cdot (x[n]-x[n-1])$ you have

$$\sum y^2[n] = \sum x^2[n]$$

That's really simple if $x[n]$ is white noise. White noise is uncorrelated with itself other than at a lag of zero and specifically we have $r_{xx}[-1] = 0$. That means you are simply subtracting two uncorrelated sequences and the energy of the sum (or difference) is the sum of the energies. Assuming $x[n]$ has a RMS of 1, than $x[n]-x[n-1]$ has a power of 2 or an RMS of $\sqrt{2}$. So the scale factor simply becomes

$$A = \frac{1}{\sqrt{2}}$$

Your original method doesn't work because you are using a continuous model to solve a discrete problem. It's simply not applicable.

You can do it in the discrete frequency domain. Let's assume a DFT length of $N$ which is sufficiently large with the transform pairs $x[n] \leftrightarrow X[k]$ and $y[n] \leftrightarrow Y[k]$ . We also assume DFT scaling of $1/\sqrt{N}$ in both directions which preserves Perceval's Theorem between discrete time and frequency.

We have $$X[k] = 1 \\ Y[k] = Ae^{-j2 \pi/N \cdot k/2} \cdot 2 j \cdot \sin(2\pi/N \cdot k/2) $$

The magnitude squares (or PSDs) are $$|X[k]|^2 = 1 \\ |Y[k]|^2 = 4A^2\sin^2(2\pi/N \cdot k/2) = 2A^2(1 - \cos (2\pi/N \cdot k)) $$

The integration turns into a sum, so we get $$E_x = \sum_0^{N-1} |X[k]|^2 = N \\ E_y = \sum_0^{N-1} |Y[k]|^2 = 2A^2( \sum_0^{N-1} 1 + \sum_0^{N-1} \cos (2\pi/N \cdot k)) = 2A^2 N $$

Again we see that the one-sample difference of white noise simply doubles the power and that $A = 1/\sqrt(2)$ will match the input power. Your expected PSD becomes

$$|Y[k]|^2 = 2\sin(2\pi/N \cdot k/2) = 1 - \cos(2\pi/N \cdot k)$$

Below is a graph, that shows PSD for both white noise and the diff()'ed version with proper scale. Measurement and expectation match well.

Code:

%% violet noise
fs = 48000; % sample rate
nx = 2^16; % FFT size
% create signals
rng(1); % make it reproducible
x = randn(nx+1,1); % one extra sample for diff()
y = diff(x)/sqrt(2);
x = x(1:nx); % cut down to desired length
% calucalted FFT and PSD 
fx = fft(x)/sqrt(nx);
fy = fft(y)/sqrt(nx);
psdx = fx.*conj(fx); psdx = psdx(1:nx/2+1);
psdy = fy.*conj(fy); psdy = psdy(1:nx/2+1);

% plot it
clf;
k = (0:nx/2)'/nx; % index 0 ... 0.5
fr = k*fs; % frequency vector
semilogx(fr(2:end),10*log10([psdx(2:end) psdy(2:end)]));
hold('on');
grid('on');
xlabel('Frequency in Hz');
ylabel('Level in dB');
set(gca,'ylim',[-100 20]);
set(gca,'xlim',[fr(2) fr(end)]);

% expectation for PSD white is 1 or 0dB
plot(fr,0*fr,'Linewidth',2);
plot(fr,10*log10(1-cos(2*pi*k)),'LineWidth',2);
legend('White actual','Diff actual','White expected','Diff expected', ...
  'Location','SouthEast');
title('white and violet noise, unit power, fs = 48 kHz');