4
$\begingroup$

Let's say I want to measure the peak frequency of a sound. In my understanding, I take the sound snippet (here it has 2585 samples), run in in an fft => 2*2585 data points. I take the module of the first half of the fft (back to 2585 values) and here I have a spectrum.


from scipy import signal
from scipy import fft

spec = fft.fft(audio)
spec = abs(spec)[0:int(len(spec)/2+1)]

enter image description here

What I get is noisy, which I thought was because of the really high resolution of the spectrum (2585 bins) so summed the values 5 by 5 to get only 512 bins.

from scipy.stats import binned_statistic

cc=np.linspace(0, len(spec), len(spec), dtype = "int")
binspec = binned_statistic(cc, spec, 'sum', bins=512)[0]

And there... I end up with an even noisier spectrum, with giant ups and downs. This is not what I expected and this is not how other spectrums I have seen look like. enter image description here What am I doing wrong??

Improve this question
$\endgroup$
7
  • $\begingroup$ Suggest, to show: time-series plot, algorithm you are using, and spectrum you got. We cannot see what your data are, you have done and what you got. So, any answer would be only speculation $\endgroup$
    – WMXZ
    Commented Aug 8 at 9:59
  • 1
    $\begingroup$ I'm not convinced you're doing anything wrong here, and don't see anything obviously problematic with the resulting spectra you've showed. The spectrogram shows little energy below roughly 2 kHz, a high energy band from 2-9 kHz tapering off around 10 kHz and going back to low energy up to the Nyquist frequency. I suspect that if you convert the x-axis of your spectra from bin number to Hz, you'll find exactly that. If the colour scaling on the spectrogram is in dB (log-scale), then that could still correspond to large changes in (linear-scale) amplitude in the spectra. $\endgroup$
    – Brian Miller ♦
    Commented Aug 9 at 0:10
  • $\begingroup$ what about the 6 dB jumps in the bandwidth where there is signal then? This is also happening with white noise, and I can't figure why... $\endgroup$
    – lframond
    Commented Aug 9 at 6:57
  • $\begingroup$ What exactly is the input format you read as sound? Maybe it is like compression artifacts or samples padded with zeros you are detecting there? The distances seem to be almost of the same length, could correspond with a sample length. $\endgroup$
    – Antares
    Commented Aug 10 at 15:27
  • $\begingroup$ There's only a small code snippet here, so the scaling and units for neither the spectrogram intensity nor the y-axis of the spectrum in the two figures are clear. There's no one right way or standard for it, and perhaps I'm mistaken, but in my experience most spectrograms tend to show colours as dB (because linear amplitudes can vary of many orders of magnitude blowing out the color scale). In contrast, the output of abs(fft(x)) would usually be positive linear amplitude (not in dB). Showing the full code for the figures would help clarify whether the plotted spectra here are actually in dB. $\endgroup$
    – Brian Miller ♦
    Commented Aug 11 at 23:53

2 Answers 2

Reset to default
4
$\begingroup$

As far as I can see, the spectrum is typical for a random noise sources. The spectral variation you see is equivalent to the temporal variation of noise sources. If you wanted to obtain an averaged (smoother) spectrum, you should use the welch approach, that uses a smaller window and averages the energy properly (you can interpret the welch also as temporal average of the 'power' spectrogram)

Edit:

To demonstrate, that there is nothing wrong with random spectrum, here is a python code with associated results using the raw data (thank you @lframond):

import numpy as np
import matplotlib.pyplot as plt
from scipy import fft
from scipy.io import wavfile
from scipy import signal
#
fname='c:/users/zimme/downloads/SnippetLframond.wav'
#
def butter_bandpass(lowcut, highcut, fs, order=5):
    return signal.butter(order, [lowcut, highcut], fs=fs, btype='band')
#
def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):
    b, a = butter_bandpass(lowcut, highcut, fs, order=order)
    y = signal.lfilter(b, a, data)
    return y

def mSpectrum(xx,win,nfft,fs):
    # estimate power spectral density (PSD) using rfft directly
    nwin=len(win)
    yy=np.fft.rfft(xx[:nwin]*win,nfft)
    pp=np.abs(yy)**2
    pp[1:] *=2                  # compensate for negative frequencies
    ff=np.arange(int(nfft/2)+1)*fs/nfft
    sc1= fs * (win*win).sum()  # for density
    sc2= win.sum()**2          # for spectrum
    return ff,pp/sc1,pp/sc2

fs,data= wavfile.read(fname)

audio=data[:,0]

audio=butter_bandpass_filter(audio, 2000, 8000, fs, order=5)

ta=np.arange(len(audio))/fs

# fft spectrum
nfft=2<<(len(audio)-1).bit_length()

win=np.hanning(len(audio))
ff,P_dens,P_spec=mSpectrum(audio,win,nfft,fs)

# welch periodogram
nw=256
fw,Pw=signal.welch(audio,fs=fs,nfft=2*nw,window=np.hanning(nw))

#Plots
fig, ax = plt.subplots(3, 1, num=0, clear=True,  figsize=(7,10))
ax[0].plot(ta,audio)
ax[0].set_xlabel('Time [s]')
#
q,f,t,im=ax[1].specgram(audio, Fs=fs)
plt.colorbar(im)
ax[1].set_ylabel('Frequency [kHz]')
ax[1].set_xlabel('Time [s]')
#
#mplot=ax[2].plot
mplot=ax[2].semilogy

mplot(ff/1000,P_dens)
mplot(f/1000,np.abs(q).mean(1))
mplot(fw/1000,Pw,'k--')
ax[2].set_xlim(0,10)
ax[2].set_ylim(1e-5,10000)
ax[2].set_xlabel('Frequency [kHz]')

ax0=ax[0].get_position()
ax1=ax[1].get_position()
ax0.x1=ax1.x1
ax[0].set_position(ax0)
plt.show()

Result: enter image description here Comments:

FFT size was constructed to be longer than data length (data were 5136 samples) and fft size was becoming 16384 (at least 2 times data length to sharpen spectral peaks). Spectrogram was from matplot libray to look similar to OP for welch periodogram use same fft-length as for fft-spectrum

  • spectrum density (here in log scale) is as random as in OP which is simply due to the length of random data used for the fft.
  • welch periodogram and averaged spectrogram are 'smooth' and are identical indicating that periodograms are effectively averaged spectrograms.
  • due to bandpass filter response there is a small VLF component, where welch periodogram and averaged spectrogram differ
Improve this answer
$\endgroup$
0
$\begingroup$

In addition to WMXZ response, you can also compute the spectrum over a longer duration. Assuming the sound has the same properties over time, increasing the duration over which the spectrum is computed will "smooth" the spectrum.

This is not what I expected and this is not how other spectrums I have seen look like.

Maybe the spectrum you are used to see are usually computed over a duration of more than 0.05 s at 96 kHz?

Improve this answer
$\endgroup$
1
  • $\begingroup$ My bird syllables are sadly that short... longer duration will include a lot of noise. However, we got the same kind of gigsaw-spectrums when using long (5 seconds) bursts of noise during tests... and calculating the sum in 300 Hz bins did not solve the problem at all, we still had ~ 6 dB jumps from one bin to the next $\endgroup$
    – lframond
    Commented Aug 22 at 8:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.