Answers.Body: <p>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)</p>
  <p>Edit:</p>
  <p>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):</p>
  <pre><code>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&lt;&lt;(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()
  </code></pre>
  <p>Result:
  <a href="https://i.sstatic.net/AJqP4yC8.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/AJqP4yC8.png" alt="enter image description here" /></a>
  Comments:</p>
  <p>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</p>
  <ul>
  <li>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.</li>
  <li>welch periodogram and averaged spectrogram are 'smooth' and are identical indicating that periodograms are effectively averaged spectrograms.</li>
  <li>due to bandpass filter response there is a small VLF component, where welch periodogram and averaged spectrogram differ</li>
  </ul>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2024-08-09T18:08:29.247Z
Answers.LastActivityDate: 2024-08-27T13:18:59.170Z
Answers.LastEditDate: 2024-08-27T13:18:59.170Z
Answers.UserDisplayName: WMXZ
Answers.Body: <p>In addition to <a href="https://bioacoustics.stackexchange.com/a/1702/73">WMXZ response</a>, 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 &quot;smooth&quot; the spectrum.</p>
  <blockquote>
  <p>This is not what I expected and this is not how other spectrums I have seen look like.</p>
  </blockquote>
  <p>Maybe the spectrum you are used to see are usually computed over a duration of more than 0.05 s at 96 kHz?</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2024-08-22T08:26:45.867Z
Answers.Text: 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
Answers.UserDisplayName: lframond
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2024-08-21T15:28:32.777Z
Answers.LastActivityDate: 2024-08-22T15:17:13.887Z
Answers.LastEditDate: 2024-08-22T15:17:13.887Z
Answers.UserDisplayName: Noil
Body: <p>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 =&gt; 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.</p>
  <pre><code>
  from scipy import signal
  from scipy import fft
  
  spec = fft.fft(audio)
  spec = abs(spec)[0:int(len(spec)/2+1)]
  
  </code></pre>
  <p><a href="https://i.sstatic.net/E1tXYVZP.jpg" rel="nofollow noreferrer"><img src="https://i.sstatic.net/E1tXYVZP.jpg" alt="enter image description here" /></a></p>
  <p>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.</p>
  <pre><code>from scipy.stats import binned_statistic
  
  cc=np.linspace(0, len(spec), len(spec), dtype = &quot;int&quot;)
  binspec = binned_statistic(cc, spec, 'sum', bins=512)[0]
  
  </code></pre>
  <p>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.
  <a href="https://i.sstatic.net/Cls8zFrk.jpg" rel="nofollow noreferrer"><img src="https://i.sstatic.net/Cls8zFrk.jpg" alt="enter image description here" /></a>
  What am I doing wrong??</p>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2024-08-08T09:59:14.900Z
Comments.Text: 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
Comments.UserDisplayName: WMXZ
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2024-08-09T00:10:22.533Z
Comments.Text: 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.
Comments.UserDisplayName: Brian Miller
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2024-08-09T06:57:56.957Z
Comments.Text: 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...
Comments.UserDisplayName: lframond
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2024-08-10T15:27:07.763Z
Comments.Text: 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.
Comments.UserDisplayName: Antares
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2024-08-11T23:53:18.617Z
Comments.Text: 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.
Comments.UserDisplayName: Brian Miller
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2024-08-22T08:32:56.193Z
Comments.Text: I'm using WAV files, so no compression artefacts. Yes, the spectrogram is in dB and the spectrum is not, but the large peaks and troughs also appear in the dB scale (there are two-fold changes in power from one bin to the next in the spectrum)
Comments.UserDisplayName: lframond
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2024-08-24T05:29:45.253Z
Comments.Text: Can you make the timeseries available (2585 samples is not that much)?
Comments.UserDisplayName: WMXZ
Tags: <analysis><fft><spectrum>
Title: Why do my power spectrums look like a jigsaw?
ThreadText: # Why do my power spectrums look like a jigsaw?
  
  (Asked by: username_0 on 2024-08-08)
  
  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)]
  
  
  
  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.
  
  What am I doing wrong??
  
  **Tags:** analysis, fft, spectrum
  
  > **Comments:**
  > **username_1**: 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  
  > **username_2**: 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.  
  > **username_0**: 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...  
  > **username_3**: 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.  
  > **username_2**: 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.  
  > **username_0**: I'm using WAV files, so no compression artefacts. Yes, the spectrogram is in dB and the spectrum is not, but the large peaks and troughs also appear in the dB scale (there are two-fold changes in power from one bin to the next in the spectrum)  
  > **username_1**: Can you make the timeseries available (2585 samples is not that much)?  
  
  ---
  
  ## 2 answers
  
  ### Answer by username_1 on 2024-08-09. Score: 4
  
  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 @username_0):
  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:
  
  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
  
  
  ---
  
  ### Answer by username_4 on 2024-08-21. Score: 0
  
  In addition to [username_1 response](https://bioacoustics.stackexchange.com/a/1702/73), 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?
  
  > **Comments:**
  > **username_0**: 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  
  
  ---
  
  