Improving FFT resolution - scaling Y axis

Question

I'm attempting to plot the frequency domain of the Bessel Equivalent of an FM signal in Python.

The first problem I had was improving the resolution of the FFT to enable viewing of the smaller bandwidths - I believe I have solved this by concatenating the sum of Bessel wave functions with a large, zero-filled array with a size double that of the array holding the wave data.

However, this has an effect on the scaling of the Y axis - as I increase the size of the zero padding, the magnitude of the Y axis drops away. What factor can I multiply the Y axis by to reverse this? Experimentation so far has led me to a dead end...

Thanks very much!

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

A couple points:

• Are you sure bWave() function is correct. Your Bessel-function is not time dependent, so you could easily find a closed form solution by Fourier transforming the cosine.
• Do not pad with zeros, but increase the time period of your bWave() signal (see code below) to increase the frequency resolution.
• Use numpy instead of math functions. It makes your code more readable and faster.

The following code plots the FFT for different time periods. The longer the time period becomes, the sharper the peaks will become (the Fourier transform of the cosines):

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()

Note that /N in the rfft(y)/N is just added to have comparable FFT values. Since the sampling rate is constant, the energy, i.e., |Y(f)|², increases with increasing time period. In your code, the sampling rate is changed and with that also the signal's energy.