My goal is to create time-spectrum noise signals with spectral densities $S_{xx}(f)$ for pretty much arbitrary spectral densities (think white noise, lorentzian, $1/f$). In finding advice on how to do so I went over to our neighbours on the signal processing stack exchange, and they were very helpful. Their advice was to use that filtering white noise with unitary power spectral density (from independent Gaussian samples with zero mean and unitary standard deviation) with a filter having a frequency response $H(f)$ gives you colored noise with power spectral density $S_{XX}(f) = |H(f)|^2$.
They then went on to suggest that this can be done with Mathematica, using the Frequency sampling-based FIR filter FrequencySamplingFilterKernel, using an evenly spaced sampling of $|H(f)| = \sqrt{S_{XX}(f)}$. Once I have then obtained the FIR filter coefficients, I could use those to filter your white noise with ListConvolve.
However, the trouble is that I'm not pulling this off very well. My initial idea was to simply test if I could transform white noise with $\sigma^2 = S_{XX}(f) = 1$ to white noise with $\sigma^2 = S_{XX}(f) = 4$ using $H(f) = 2$, so lets look into that and see what I'm doing wrong.
I start by generating the Gaussian white noise
times = Range[0, 1/Δf, 1/(2 fmax)];
UnitDataWhite = RandomVariate[NormalDistribution[0, 1], Length[times]];
where I use
K = 100; (* number of frequency components *)
fmax = 100; (*max frequency*)
Δf = fmax/K;(* freq. step *)
in anticipation of having to generate amplitudes for the filter. Using ListPlot this seems fine.
Alright, so how does one then move on..
Time to define some filter amplitudes. Repeating a part of the above
K = 100; (* number of frequency components *)
fmax = 100; (*max frequency*)
Δf = fmax/K;(* freq. step *)
fk = Range[0, fmax, Δf]; (* sampling frequencies *)
FilterAmps=Table[2,{i,1,Length[fk]}];
This should give us some amplitudes for the filter. Introducing the filter itself, we use
FilteredNoise =
ListConvolve[FrequencySamplingFilterKernel[FilterAmps],
UnitDataWhite, {-1, -1}];
Or anyway, I think that is what one should do. I'm a bit confused about the final part however, specifying the overhang. If I leave out {-1, -1} I just get a single value, which is definitely not what we want, but what type of overhang should I choose? It is not clear to me at all, and I suppose it is among the main parts of the question. But perhaps more importantly, how do I know that what I did has worked?
Related to this I suppose is the question of how I can visualize the spectral density to see that it has indeed gone from a constant 1 to a constant 4. I feel like this might be a bit of a challenge. Should I use autocorrelations and fourier transforms?
And finally, I'm still a bit unsure about the parameters used here, but I guess that's more a signal processing matter. Like, I suppose in my code the factor controlling how good the filter is is $K$, is it not? The code seems to slow down drastically with increasing $K$, so accurately including a lot of higher frequencies seems kind of hard to do.
ListConvolve[h,x,{1,-1},0]implements the convolution of discrete-time sequenceshandx, as defined in standard signal processing textbooks. Regarding visualization of power spectral densities, check the periodogram method. – Stelios Mar 02 '16 at 21:08