Блог

Я читаю файлы .wav различных записей и пытаюсь построить частотный спектр с помощью БПФ. Результаты не имеют для меня никакого смысла. Я не уверен, что

• У меня ошибка в коде или
• Я неправильно понимаю использование функций библиотеки numpy FFT или
• Я неправильно понимаю математику / обработку сигналов, и я вижу артефакт, вызванный частотой дискретизации.

На этом графике я читаю wav-файл синусоидальной волны 1 кГц, записанной со скоростью 44100 отсчетов в секунду. Я строю только 10 циклов (0,01 секунды) или 442 точки данных.

Я ожидал увидеть один высокоэнергетический пик на частоте 1 кГц. Но на самом деле я вижу эту закономерность, которая для меня имеет мало смысла.

Мой код приведен ниже .. и может содержать ошибки, но я подозреваю, что вижу артефакт, вызванный относительно низким разрешением выборки по сравнению с частотой сигнала.

Любая помощь или идеи о том, что здесь происходит, были бы очень полезны, пожалуйста.

»’

 # -*- coding: utf-8 -*- """ Created on Mon Oct 11 18:58:39 2021  """  import sounddevice as sd import soundfile as sf import numpy as np import matplotlib.pyplot as plt  soundDir = '..\01 - Engine Recordings\'  #soundFileName = soundDir   '500 Hz Sine.wav' soundFileName = soundDir   '1000 Hz Sine.wav'  theSignal, fs = sf.read(soundFileName)  length, numChannels = theSignal.shape  print('The shape of the signal data is = ', theSignal.shape)  # I am only interested in mono data  if (numChannels == 2):  theSignal = np.delete(theSignal,1,1)     # now reduce the length of the file - we only need a few seconds theSignal = theSignal[0:442] # 442 = whole number of cycles  print('Sampling rate = ', fs) print('Number of bins in sample', theSignal.size) print('Therefore, lengh of sample in seconds = ', theSignal.size * 1/fs) xAxis = [(x / fs) for x in range(0, theSignal.size)] # Creat a real-time x axis  lastX = xAxis[len(xAxis) - 1] maxY = np.amax(theSignal) minY = np.amin(theSignal) print('Original Signal : Max amplitude = ', maxY, 'Min ampltude = ', minY)  fig,ax = plt.subplots(2)  # Plot the oscillascope ax[0].set_xlim([0,lastX]) ax[0].set_ylim([minY, maxY]) ax[0].set_xlabel('Time (s)') ax[0].set_ylabel('Amplitude') ax[0].set_title('Oscillascope')  ax[0].scatter(xAxis, theSignal, s=1)   yf = np.fft.fft(theSignal) # Calculate FFT yfAbs = np.abs(yf) # just for plotting xf = np.fft.fftfreq(yf.size, 1/fs) # calculates the frequencies in the center of each bin in the output of fft().   # For PLOTTING PURPOSES ONLY .. we don't need the negative frequency components .. so I can filter these out # Note that we WILL need them if we want to reconstruct the original signal from the Frequency Domain  print('Number of frequency bins = ', xf.size) zeroBinIndex = int(xf.size / 2) print('Therefore, the zero value is at bin', zeroBinIndex, 'With a value of',xf[zeroBinIndex])  # zeroBinIndex = 1200  yfAbs = yfAbs[0:zeroBinIndex] xf = xf[0:zeroBinIndex]  maxX = np.amax(xf) minX = np.amin(xf) print('Highest frequency = ', maxX, 'Lowest frequency = ', minX)  maxyfAbs = np.amax(yfAbs) minyfAbs = np.amin(yfAbs) print('Lowest energy in plot = ', minyfAbs, 'Max energy = ', maxyfAbs)   # Plot the frequency spectrum ax[1].set_xlim([np.amin(xf), np.amax(xf)]) ax[1].set_ylim([np.amin(yfAbs), np.amax(yfAbs)]) ax[1].set_xlabel('Frequency (Hz)') ax[1].set_ylabel('Energy') ax[1].set_title('Spectrum') ax[1].scatter(xf, yfAbs, s=3)  # as a bar chart # xf = np.array(xf).flatten() # yfAbs = np.array(yfAbs).flatten() # ax[1].bar(xf, yfAbs,width = 1)  plt.show()   # sd.play(theSignal, fs) # sd.wait()    

»’