3

I have a sound signal, sampled at 48000 hz. Now I want to know 'how much sound' there is of a specific frequency (or near that frequency) at a specific time. For example, I want to know at 10s how much sound there is near the frequencies $20 \cdot 2^{k/3}$ hz with $k=0,1, \dots, 30$.

I'm trying to find a connection between a sound someone is hearing and an fMRI scan of his brain while he is listening to that sound. This website says that sounds with different frequencies are processed in different areas of the brain, that's why I need to know this.

When doing a bit of research I quickly arrived at short-time Fourier transform (STFT). I don't need an inversible transform, but I can use some of the principles. Let's assume that the signal is continuous. If I understand it right, then I can use the formula $$ y(t, f) = \left| \int_{-\infty}^{\infty} x(\tau)\ w(\tau-t)\ e^{-2 \pi i f \tau} \ \mathrm{d}\tau \right|, $$ where $x(t)$ is the signal and $w(t)$ a window function. I read a Gaussian window function is common: $$ w(t) = e^{-\frac12 (t/\sigma)^2}, \quad \sigma > 0 $$ Should I use this window function or some other window function? If I should use this one, what value should $\sigma$ have, should it be fixed or dependent on $f$? It seems intuitive to me that the width of the window should be proportional to $1/f$, but no one seems to be doing this. Can anyone tell why? Is the first formula the way to go or is there a better way to do what I want?

I'm also thinking whether I want 'near a frequency' at all, maybe I want 'between two frequencies' instead, or maybe something else. I understand you cannot guess what I want exactly, but I hope someone has some ideas anyway.

Paul
  • 2,795
  • 1
    I would suggest the Gabor Transform, which is a specific case of the STFT. The Gabor Transform has the best compromise between time and frequency uncertainty - there is a free parameter which determines whether you want more accuracy in the time domain or the frequency domain. – Ryan Aug 14 '14 at 01:31
  • @Ryan, how should I use this transform exactly? Do you suggest using $\left| \int_{-\infty}^\infty e^{-\pi(\tau-t)^2}e^{-2\pi i f\tau}x(\tau)\ \mathrm{d}\tau \right|$ to calculate the amount of sound at time $t$ near frequency $f$? – Paul Aug 16 '14 at 19:02

1 Answers1

2

The problem has several variables.

You'll need a window (in time) around the time where you want to make the determination. Windows of 125 msec are common in many audio applications.

Once you've selected your sample window, you can choose a window function (you've suggested a Gaussian window). Window functions are mainly used to compensate for spectral artifacts generated by the discrete Fourier transform. A Hamming-type window function is generally a good "middle of the road" window function with reasonable sensitivity for frequencies over the range of magnitudes you're looking at.

The third question is how close to the specific frequencies do you want to measure the energy. Since this experiment is targeting human hearing, you might want to look up the smallest variation in frequency that the human brain can distinguish. Or you might want to make it another independent variable in the experiment.

Given a selection of a sample window length, a window function, and a bandwidth (i.e. the margin of frequencies you consider to be "near" your center frequency), you can then run a discrete Fourier transform (DFT) on your sampled, windowed data. You can then sum the squared amplitudes of all of the spectral components inside your selected bands to obtain the total energy "near" the center frequencies.

Of course, when in doubt, dig into the literature and see what other people in the field are doing. ;)

COTO
  • 585
  • Why the squared amplitudes? If you don't square the amplitudes, it is a linear function of the signal, and I think that is desirable, right? – Paul Aug 15 '14 at 15:25
  • The energy of a signal in both time and frequency domain is obtained by integrating (or adding, in the case of discrete signals) the squared amplitudes of the complex components.

    See http://en.wikipedia.org/wiki/Parseval%27s_theorem

    The specific formula you're interested in falls under "discrete Fourier transform", only in your case you wouldn't be computing the sum over the full range of frequencies, just in the bands that interest you.

     – COTO Aug 15 '14 at 15:50
  • So in formulas, your idea is the following, right? $$ \begin{align} r &= 48000 \ w &= 0.125 \ N &= r\cdot w = 6000 \ x_n &= x\left(t+ \frac{n-(N-1)/2}{r} \right), \qquad n = 0, 1, \dots, N-1 \ k_{min} &= \lceil f_{min} / w \rceil \ k_{max} &= \lfloor f_{max} / w \rfloor \ y &= \sum_{k=k_{min}}^{k_{max}}\ \left| \sum_{n=0}^{N-1} x_n \cdot \left(\frac{25}{46} - \frac{21}{46} \cos \left(\frac{2\pi n}{N-1}\right)\right) \cdot e^{-2 \pi i k n / N} \right|^2 \end{align} $$ – Paul Aug 16 '14 at 18:20
  • Correction: it should be $f_{min/max} \cdot w$ instead of $f_{min/max} / w$ – Paul Aug 18 '14 at 16:46
  • @Paul: It looks right at a glance. Unfortunately I won't have time to give it a thorough looking over until Friday. – COTO Aug 19 '14 at 01:22