In my last article, I mentioned how I was trying to divine means by which I can synthesize the sounds of storms. While I was originally going to write about a photograph, I quickly diverted myself into the synthesis of wind.

When talking about the sounds produced by wind, there are actually several sources to consider. Beyond the direct sound of moving air, we have natural resonances in the environment, the rustle of trees as their branches and leaves strike, and similar interactions with loose features of houses and other structures. Ignoring the interaction of loose objects (which would be a Poisson process), it appears that resonances in the environment are the dominant feature of the sound.

How did I come across this conclusion? Like any good engineer, I recorded myself emulating the wind and analyzed its spectrum. To my astonishment (no, not really), the spectrum exhibited a strong Cauchy distribution, the frequency distribution of physical resonance. Remarkable. For those wondering why the wind doesn't sound like a tuning fork, it's simply that nature is a really crappy resonator.

When synthesizing sounds of this nature, the common techniques is to apply an appropriately-crafted filter to white noise:

Since I'm synthesizing a resonating system, my filter (in frequency space), needs to be of form:

Where f₀ is the frequency of resonance and γ is the full-width at half-maximum (HWHM).

Okay, so where do I find these parameters? In the long term, I'd need to record lots of wind with known properties (namely, wind speed), but for now, I'll steal them from a recording. Like this page. So, after opening up “Wind 2” in MATLAB (whose file name is actually wave1.wav), I perform a simple least-squares fit to its spectrum and claim f₀=760 Hz and γ=290 Hz. In reality, there are likely multiple resonances, but we'll concentrate on the dominant one.

Now what? First, we need to construct the filter. The inverse Fourier transform of H(f) is actually quite simple, being of form:

But, if we are implementing this filter for real (say, in MATLAB), our h(t) won't stretch off into infinity. Instead, we will be building a finite impulse response (FIR) filter. Finite because each sample of noise only effects the output for the finite length of the filter.

Simply cropping h(t) to a limited range won't suffice as this means the real H(f) will be convolved with a sinc function. The sinc function is evil. It makes me angry. The whole point of this exercise is for me to not be angry. To those in which I have confided, this is why I despise moving averages with an unholy passion. Still, any FIR will have an inescapable influence from the sinc function and we can only mitigate its effect.

One common approach is to solve non-linear equations allowing one to minimize the least squares difference between the ideal transfer function and the FIR. We won't be doing that. Instead, we'll simply apply a window function to the cropped FIR. Personally, I'm a big fan of the Hamming function (just like 94% of the engineering community).

Now...how long do we make the FIR? The longer we make the transfer function, the finer the frequency precision (and less severe the sinc), but the more computationally expensive it becomes. If our filter has length N and the signal sampling rate is F_s, the frequency resolution is F_s/N. Clearly, this needs to be finer than our FWHM of 290 Hz. (In practice, we'd probably choose a number of form 2ⁿ+1).

So let's fire up our MATLAB to generate the filter:

Fs = 11025;        % Sample Rate in Hz
g = 290;           % Width of Cauchy distribution (Hz)
f0 = 760;          % Center of Cauchy distribution (Hz)
N = ceil(2*Fs/g);  % Length of the FIR

% Time axis for the FIR
th = (-floor(N/2):floor(N/2))'/Fs;
% Generate the FIR
h = cos(2*pi*f0*th) .* exp(-g*abs(th)) .* hamming(N);

Now, let's generate some wind:

T = 10;            % Number of seconds to generate

n = randn([T*Fs 1]);
y = filter(h, sqrt(sum(h.^2)), n);
sound(y / 4, Fs);

For those lacking a copy of MATLAB, here's a sample: windsynth1.wav

While the basic sound has been accomplished, it needs a bit more. It needs to feel more ‘gusty’ and less uniform. For this, we'll modulate the amplitude of the base sound. Assuming the amplitude is proportional to the wind velocity, and since wind velocity follows a roughly Rayleigh distribution (or so the Internets tell me), the amplitude will also follow a Rayleigh distribution. This is where things get tricky.

Generating a random numbers according to the Rayleigh distribution is easy, but wind changes its velocity very slowly. I clearly need to lowpass filter my sequence; however, there is a problem here: the central limit theorem tells me that if I add up a bunch of random variances (such as would happen in a filter), the resulting distribution will start to term into a normal distribution. Not only that, a really low frequency FIR will be extremely long. Impractically so.

The second part is easy. Instead of a finite filter, we'll use an infinite filter, the infinite impulse response (IIR) filter. What makes it infinite? Simple, the output of the filter is fed back into the filter. Once an input sample makes its way into the output, this feedback means it will never leave so we don't need a long filter for low-frequency filtering nor do we end up with a sinc function to cancel out. Okay, so why didn't we use it before? IIRs are much more constrained in the filter shapes they can produce and have this horrible tendency to oscillate.

So, if we IIR some noise, how do we get something that is Rayleigh distributed? While I was busy kludging together various mathematical functions in heretical fashion, a compatriot pointed out the obvious. Right after I explained said obvious to him. The Rayleigh distribution is a specific case of the chi distribution having two degrees of freedom. What does that mean? I can get my Rayleigh distribution by taking the root-mean-square (RMS) of two independent normal random variables, such as that produced by my IIR. Perfect.

Let's return to the matlab:

% Generate two Gaussians
na = filter(1, [32*Fs 1-32*Fs], randn([T*Fs 1]));
nb = filter(1, [32*Fs 1-32*Fs], randn([T*Fs 1]));
% Take their RMS to generate a Rayleigh
a = sqrt(32*Fs) * sqrt(na .^ 2 + nb .^ 2);              
y2 = a .* y;
sound(y2 / 4, Fs);

Again, for the MATLAB impaired, here's a sample: windsynth2.wav

The IIR used here is simplistic, unrealistic, and complete arbitrary, but until I'm able to collect better wind samples for analysis, it will have to do.

For the curious, I've been using normally-distributed random variables in this example, but uniformly-distributed random variables (such as the rand function in your favorite programming language) can work, too. The filters (by virtue of the central limit theorem) will force anything to become normally distributed. Almost anything. Just make sure it has finite mean and variance. The Cauchy distribution simply won't do.