0
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ ListConvolve[h,x,{1,-1},0] implements the convolution of discrete-time sequences h and x, as defined in standard signal processing textbooks. Regarding visualization of power spectral densities, check the periodogram method. $\endgroup$ – Stelios Mar 2 '16 at 21:08
  • $\begingroup$ The periodogram is interesting, I'll have to look into some examples for that. As for the convolution, I see what you mean. I'm not sure where that leaves me with the filtering however (I'm very new to digital filters, as is obvious..). The way I see it I have some original data, being my white noise in the time domain. These are $N$ amplitudes $y_n$ at times $t_n$. Then I create $K$ components for the filter with amplitudes $u_k$. Using listconvolve leaves me with N-K+1 points. How do I then attribute a time $t$ to these points again? $\endgroup$ – user129412 Mar 2 '16 at 21:22
  • $\begingroup$ It is obvious that you need to work a little on your signals and systems knowledge. You should consider the signal processing stack exchange for this purpose, although, I would advise spending some time on your own studying fundamentals of (digital) signal processing such as filtering, sampling of continuous time signals, etc. $\endgroup$ – Stelios Mar 2 '16 at 21:40
  • $\begingroup$ That seems fair. I'm reading up on some of it now. It seems that I need to do something with zero padding. $\endgroup$ – user129412 Mar 2 '16 at 22:17
  • $\begingroup$ Strictly codewise speaking though, my problem is as I have written above, in that I don't understand what to do with the ListConvolve when I have more time domain samples $u(t)$ than filter amplitudes, in the sense of getting an output for which I understand what point corresponds to what time $t$. Overhang seems to be related to this, but as far as I understand that it just makes the amplitudes cyclic and I'm not really sure if that is a sensible thing to do here. $\endgroup$ – user129412 Mar 2 '16 at 22:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.