Dynamic moving ripple

image

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

UNDER CONSTRUCTION

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).

Carrier - fine structure

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.

image

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;