# Generating violet noise with a specific PSD coefficient

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.

• is np.diff really the approximation to differentiation you want to make? Because if I remember correctly, it's just out[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.) Sep 28 '21 at 15:38
• A different way to write it $H(\omega) = \frac{1}{2j} e^{-j\omega T/2}\cdot sin(\omega T/2)$, This way you can see the magnitude directly. Sep 28 '21 at 15:47
• But I don't think that's the problem here. I would take a look at the different scaling option of 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. Sep 28 '21 at 15:49
• Thanks @Hilmar. I believe there should be no issue with 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. Sep 28 '21 at 15:53
• Also, would it be possible for you to point me to a reference where I can see how this differentiation (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. Sep 28 '21 at 16:13

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

• hank you @Hilmar; this makes a lot of sense. Isn't the sampling period missing in the expression though? This has also made me realize a few things: As you rightly pointed out, for 2𝜋𝑓𝑇𝑠 (which is the region I'm interested in) this expression is approx to what I am currently using, which means I should be good. However, this still doesn't answer the factor of $\pi$ discrepancy I have. Sep 28 '21 at 19:29
• Fantastic. Thank you so much for the detailed response. Sep 29 '21 at 14:51
• By the way @Hilmar, is there any chance you can share your code? Thank you! Sep 29 '21 at 17:26
• @EdV: Thanks ! fixed Sep 29 '21 at 18:14