A sound with a rich time-varying spectrum. Listen here.

Natural sounds (specifically speech and vocalizations) vary dynamically along multiple
stimulus dimensions, including time, frequency, and intensity. To understand how our brains process this acoustic information, it is essential to characterize the neural responses to the jointly changing spectral and temporal dimensions. Escabi and Screiner (2002) used the * dynamic moving ripple* to determine linear and non-linear properties of neurons in the inferior colliculus (an auditory midbrain-structure). Basically, they presented sounds that had a randomly time-varying spectrum, and determined the response to these stimuli via spike-triggered averaging. In this tutorial, we will create such a dynamic moving ripple in Matlab (and conveniently ignore how this stimulus can be used to characterize a neuron's auditory response or an organism's behavior).

To create a sound whose amplitude of every frequency component changes as a function of time, we will start creating a sound with no spectrotemporal modulation (the carrier) and later introduce the dynamically-varying envelope. For this we create a broadband complex by simply adding a large number of tones together. In formula-form:

\( s(t) = \sum_{i=1}^N \sin{(2 \pi f_i t + \phi_i)} \)

where \( s(t) \) is the sound, consisting of \( N \) sinusoidal components, each of frequency \( f \), and \( \phi \) is a randomly chosen phase (between \( 0-2\pi \)). More specifically, we will create a 1-sec sound of 128 frequency components, with the lowest frequency starting at 250 Hz, and we will use 20 components per octave (so the largest frequency component will be \( 250 \cdot 2^{128/20} \equal 21112 Hz\). In Matlab:

%% Initialization f0 = 250; % base frequency (Hz) N = 128; % # components fNr = 0:1:N-1; % every frequency step f = f0 * 2.^(fNr/20); % frequency vector (Hz) dur = 1; % sound's duration (s) Fs = 44100; % sample frequency (Hz) ns = round(dur*Fs); % number of time samples t = (1:ns)/Fs; % time vector (s); %% Sum carriers, in a for-loop snd = 0; for i = 1:N phi = 2*pi*rand(1); % random phase between 0-2pi carrier = sin( 2*pi * f(i) * t + phi ); snd = snd+carrier; end snd = snd/N;

Here we a broadband complex of 128 components equally distributed (20/active) from 250 Hz to ~20 kHz.

F0 = 250; % base frequency (Hz) nFreq = 128; % (octaves) FreqNr = 0:1:nFreq-1; % every frequency step Freq = F0 * 2.^(FreqNr/20); % frequency vector (Hz)

The amplitude of every frequency component changes as a function of time, which can be written in mathematical form:

\( S(t,x) = 1+\Delta M \cdot \cos{(2 \pi \omega t + \Omega x)} \)

with \( t \) time (s), \( x \) position of the spectral component in octaves above the lowest frequency, \( \omega \) ripple velocity (Hz), \( \Omega \) rippledensity (cycles\octave), and \( \Delta M \) the modulation depth on a linear scale between 0 and 1. In the next sections, the meaning of the parameters will be explained by generating a ripple in Matlab.

Let's take a sound duration of 1000 ms, a modulation depth \( \Delta M \) of 1 (100%), a ripple velocity \( \omega \) of 4 Hz, and a ripple density \( \Omega \) of 0 octaves/cycle. We will change these parameters later, to see what the effects are.

vel = 4; % omgea (Hz) dens = 0; % Omega (cyc/oct) mod = 100; % Percentage (0-100%) durrip = 1000; %msec Fs = 50000; % sample frequency (Hz) nRip = round( (durrip/1000)*Fs ); % # Samples for Rippled Noise time = ((1:nRip)-1)/Fs; % Time (sec) Oct = FreqNr/20; % octaves above the ground frequency oct = repmat(Oct',1,nTime); % Octave

The next step is to compute the dynamic amplitude modulations of the sound, according to the above formula.

%% Create amplitude modulations completely dynamic in a loop A = NaN(nTime,nFreq); % always initialize a matrix for ii = 1:nTime for jj = 1:nFreq A(ii,jj) = 1 + mod*sin(2*pi*vel*time(ii) + 2*pi*dens*oct(jj)); end end

Finally, we can add this envelope *A* to the carrier sound *snd*, and we will have our first ripple.

% Modulate carrier, in a for-loop snd = 0; for ii = 1:nFreq carr = A(:,ii)'.*sin(2*pi* Freq(ii) .* time + Phi(ii)); snd = snd+carr; end

When we look at the spectogram, we can see that the amplitude is changing only as a function of time. We have just created an amplitude-modulated (AM) sound.

This is exactly what happens if you change only the ripple velocity. Changing ripple density will lead to frequency-modulated (FM) noise. Changing both will lead to dynamic variations as shown in the top-figure. Try this out yourself. A negative density corresponds to an upward direction of the spectral envelopes, a positive density to a downward direction, and \( \Omega = 0\) indicates a pure amplitude modulated (AM) sound.

With the PANDA function* pa_genripple* it is easy to create these ripples and visualize them in a graph. The figure above can also be made with the following Matlab code.

figure t = (1:length(snd))/Fs; subplot(221) plot(t,snd,'k-') ylabel('Amplitude (au)'); ylabel('Time (ms)'); xlim([min(t) max(t)]); axis square; box off; subplot(223) nfft = 2^11; window = 2^7; % resolution noverlap = 2^5; % smoothing spectrogram(snd,window,noverlap,nfft,Fs,'yaxis'); cax = caxis; caxis([0.7*cax(1) 1.1*cax(2)]) ylim([min(Freq) max(Freq)]) set(gca,'YTick',(1:2:20)*1000,'YTickLabel',1:2:20); axis square; subplot(224) pa_getpower(snd,Fs,'orientation','y'); % obtain from PandA ylim([min(Freq) max(Freq)]) ax = axis; xlim(0.6*ax([1 2])); set(gca,'Yscale','linear','YTick',(1:2:20)*1000,'YTickLabel',1:2:20) axis square; box off;