Улучшение разрешения FFT - масштабирование оси Y

Question

Я пытаюсь построить частотную область эквивалента Бесселя FM-сигнала в Python.

Первая проблема, с которой я столкнулся, заключалась в улучшении разрешения БПФ, чтобы разрешить просмотр меньших полос пропускания. По-моему, я решил это, объединив сумму волновых функций Бесселя с большим, заполненным нулями массивом с двойным размером массива, содержащего волновые данные.

Однако это влияет на масштабирование оси Y - поскольку я увеличиваю размер нулевой прокладки, величина оси Y падает. Какой множитель я могу умножить на ось Y, чтобы отменить это? Эксперимент до сих пор привел меня в тупик ...

Большое спасибо!

import numpy as np 
import matplotlib.pyplot as plt 
import math as mt 
import scipy.special as sc 
def bWave(fc=10000, B=1.25): 
    time = np.linspace(0, 0.5, 20001, True) 
    data = np.zeros(len(time)) 
    retarr = np.column_stack((time,data)) 
    Vc = 6 
    fm = 3 
    for n in range(-5,5): 
     for row in range(len(retarr)): 
      retarr[row][1] += Vc*sc.jv(n,B)*mt.cos(2*mt.pi* 
       (fc+n*fm)*retarr[row][0]) 
    return retarr 
FM_array = bWave() 
# ------------- SIGNAL PLOT OF AM ----------------- 
scaling = 2 #default 2, cuts out symmetry from FFT 
buffer_ratio = 1 
padded_array = 
    np.concatenate((FM_array[:,1],np.zeros(buffer_ratio*len(FM_array[:,1])))) #pad array with zeros 
Y = np.fft.fft(padded_array) #perform FFT 
N = len(Y)/scaling + 1 # get FFT length (avoid reflection) 
T = FM_array[1][0] - FM_array[0][0] #get time interval of FFT 
fa = 1/T #sampling frequency 
Xaxis = np.linspace(0, fa/scaling, N, endpoint=True) # create x axis vector from 0 to nyquist freq. (fa/2) with N values 
plt.plot(Xaxis, (2.0/((N)*scaling)) * np.abs(Y[0:N])) # multiply y axis by 2/N to get actual values 
plt.grid(True) 
plt.show() 

Answer 1

Пара очков:

• Вы уверены, что bWave() функция является правильным. Ваша функция Бесселя не зависит от времени, поэтому вы можете легко найти решение закрытой формы, преобразование Фурье косинуса.
• Не заполняйте нулями, но увеличивайте время сигнала bWave() (см. Код ниже), чтобы увеличить частотное разрешение.
• Использование numpy вместо math функции. Это делает ваш код более читабельным и быстрым.

Следующий код описывает БПФ для разных периодов времени. Чем дольше период времени, тем острее пики станут (преобразование Фурье косинусов):

import numpy as np 
import matplotlib.pyplot as plt 
import scipy.special as sc 


def bWave2(t, fc=10000, B=1.25): 
    """ Useing only numpy """ 
    Vc, fm = 6, 3 
    y = np.zeros(len(t)) 
    for n in range(-5,5): 
     y += Vc*sc.jv(n,B)*np.cos(2*np.pi*(fc+n*fm)*t) 
    return y 


fg, ax = plt.subplots(1, 1) 

fc=10000 
for q in range(0, 5): 
    k = 15001*(1+q) 
    t = np.linspace(0-0.25*q, 0.5+0.25*q, k) 
    y = bWave2(t, fc) 

    Y = np.fft.rfft(y)/N # since y is real, rfft() instead of fft is used 
    f = np.fft.rfftfreq(k, t[1] - t[0]) # frequencies for rfft() 

    ax.plot(f, np.abs(Y), label=u"$\\tau=${}".format(t[-1]-t[0]), alpha=.5) 
ax.set_xlim(fc-50, fc+50) 
ax.grid(True) 
ax.legend(loc="best") 
fg.canvas.draw() 

plt.show() 

Обратите внимание, что /N в rfft(y)/N просто добавляется, чтобы иметь сопоставимые значения БПФ. Так как частота дискретизации постоянна, энергия, т. Е. | Y (f) | ², увеличивается с увеличением периода времени. В вашем коде изменяется частота дискретизации, а также энергия сигнала.