The operating mechanism of spectrogram in audio Science

This paper introduces scipy The version information is as follows

Name: scipy
Version: 1.4.1
Summary: SciPy: Scientific Library for Python
Home-page: https://www.scipy.org
Author: 
Author-email: 
License: BSD
anaconda3/envs/torch1.7.1/lib/python3.7/site-packages
Requires: numpy
Required-by: scikit-learn, resampy, librosa

Pictured , given spectrogram Core calling function in , besides ,spectrogram Function also calls other functions ;

1. Input/output parameter

def spectrogram(x, fs=1.0, window=('tukey', .25), nperseg=None, noverlap=None,
                nfft=None, detrend='constant', return_onesided=True,
                scaling='density', axis=-1, mode='psd'):


   ...
   pass
   
   return freqs, time, Sxx

1.1 Input parameters

 data:  Audio data ,
 fs:    Sampling rate 
 mode:  Plural pattern 

 **kwargs: {
      nperseg: 1000, noverlap: 500, nfft: 5000 }

nperseg: number of per segment;
Stands for every paragraph （ frame ） Point length of , The default value is None, If given nperseg, Then use the given value ;
If window yes Character or tuple , Set to 256,
If window yes It's like an array , Set to window The size of the length ;

noverlap: The overlap between two adjacent frames ;
nfft; For each frame FFT The number of transformed points ;

1.2 Output parameters

    Returns
    -------
    freqs : ndarray
        Array of sample frequencies.
    t : ndarray
        Array of times corresponding to each data segment
    result : ndarray
        Array of output data, contents dependent on *mode* kwarg.

freq:  Frequency domain information  ;   Frequency of samples 
t:     Time domain information ;     Multi frame time information 
spec:    Time spectrum ,  be based on  mode   Choose output ;

1.3 Input , Output calculation

When the input parameter ：

 data= 10001  A numerical ;  fs = 100;  mode = "complex";
 **kwargs: {
     nperseg: 1000, noverlap: 500, nfft: 5000 }

The output parameters are as follows ：

freq =  2501  A numerical ,  Represents frequency ;
t = 19  A numerical ,   representative 19  frame ;
spect = (2501, 19)    The output is a complex matrix ,  because mode = "complex"

2. spectrogram() The operating mechanism in the function

2.1 Choice mode

according to modelist = [‘psd’, ‘complex’, ‘magnitude’, ‘angle’, ‘phase’]

Select the corresponding mode

2.2 call _triage_segments()

window, nperseg = _triage_segments(window, nperseg, input_length=x.shape[axis])

 Input ： window,  nperseg,  input_length()
 Returns the parameter ：  win,  nperseg()

Set up window, nperseg Parameters , be used for spectrogram and _spectral_helper()

Parameter interpretation ：

window: string, tuple, or ndarray There are three types of options ;

• If window The window is through string ,or tuple It specifies , also nperseg Is not specified ,
  be nperseg Set to 256, And returns the length of the specified window ;

• If window It's a arrary_like（） Array like , also nperseg Is not specified ,
  be nperseg Set to the size of the window length ;

• Wrong situation , If the user also provides an array like window , Also provide nperseg The numerical , however nperseg ！= The length of the array window ;
  May be an error ;

When window yes string perhaps tuple Type ：
such as , window = (“tukey”, 0.25), nperseg = 1000;

• Perform the following process
  Judge , nperseg Is it greater than input_length,
  if nperseg > input_length , A warning is issued , It indicates that the length of one frame is set to be greater than the total length of input ,
  Back to the nperseg assignment , The assignment is input_length;

• And then call get_window(window, nperseg) function , obtain win length ;
  return Three parameters ,

  1. A paragraph （ frame ） Sample points in （ The length is nperseg）,
  2. The type of window ;
  3. Whether to sample fftbins,

Final triage_segments() The return is ,

1. window() The number of sample points in a frame , here 1000 Number ;
2. nperseg() The length of a frame sample , Here is a Integer numbers , Values for 1000;

2.3 Set up noverlap

if noverlap is None:
    noverlap = nperseg // 8

2.4 Core call _spectral_helper()

Select according to the mode , mode == stft:

When the mode is STFT perhaps complex when , The two are equivalent ;

freqs, time, Sxx = _spectral_helper(x, x, fs, window, nperseg,
                                            noverlap, nfft, detrend,
                                            return_onesided, scaling, axis,
                                            mode='stft')                                                                                               

Output results freqs,time and Sxx, The following is an analysis of the implementation process ：

1. The judgment mode is psd, stft

2. boundary_func : Whether to carry out boundary extension on both sides of the input signal ,
   The effect of continuation is in data X On the basis of , Extend to both sides （nperseg//2） The length of ,
   The extension methods are (const_ext, even_ext, odd_ext, zero_ext);

3. Judge x, y The two data types are equal

4. Make the output type np.complex64 The plural ;

5. call win, nperseg = _triage_segments(window, nperseg, input_length= x.shape[-1] )
   here win = Previous window , For the same thousand numbers ;
   nperseg Still stands for An integer number 1000 ;

6. Judge nfft The number of points must be >= nperseg, A paragraph （ frame ） Points in ;

7. nstep = nperseg - noverlap

8. according to scaling yes density still spectrum, Find the corresponding scale,
   Here the density Found scale by 0.00344;

9. freqs = sp_fft.rfftfreq(nfft, 1/fs)  
   

   This function takes nfft Medium , unilateral fft Take the real part , therefore rfft The number obtained , yes nfft Half of the ;

   freqs = 2501 ( nfft/2); After execution , freqs= ( 2501,) Is included 2501 Array of Numbers ;

10. _fft_helper()
    The main implementation FFT Calculation of transformation , In this step ;
    For each nperseg The length of the window sequence , Conduct nfft Of length fft Transformation ;

10.1.  Find out  step = nperseg - noverlap = 500
10.2.  Find out  shape = (19, 1000) tuple;
10.3.  Find out  strides = (4000, 8)
 
10.4. result = fun1（x, shape, strides）, ndarray(19, 1000)
    detrend(result),  Each frame of data is de -   Trend effect ;

10.5.  result = win * result;
   result = (19, 1000)
      
10.6. result = result. real  ,  Take out  result  The real part of ,  It is a real number ;

10.7.  Make  func = sp_fft. rfft;

10.8. result = func(result, n= nfft),  
   result = (19, 1000) , nfft = 5000;
   result = (19, 2501)   Is a complex matrix ,  Every number is a plural number ;

11. result *= scale Scale ;

12. time = {ndarray: 19}

13. result = Dimension displacement , resutl: (19, 2501) --> (2501, 19)

2.5 Summary

scipy in spectrogram Function , adopt _spectral_helper() function ;

and _spectral_helper() The core function in is through _fft_helper() Realization ;
Final spectrogram The parameter returned by the function , It's all through _spectral_helper() Function ,spectrogram Three parameters are returned ：

freqs: 2501  Number 
time: 19  Segment frame ,
result: （2501, 19）  A complex matrix ,   Represents the time spectrum ;

3. scipy Medium spectral.py The file is introduced

In this file , Mainly achieved 14 A function ,

among 10 Core functions ：

stft, istft, csd, welch, coherence, periodogram, 
spectrogram, check_COLA, check_NOLA, lombscargle

4 Auxiliary functions ;

_triage_segments()
_spectral_helper()
_fft_helper()
_median_bias()

The rest are function calls in other modules ;

3.1 Call relationship

3.2 Operating mechanism

Specific operating mechanism , Please read here Reference resources spectral.py