FFT in Python 2.7 merging two codes

Question

I have two codes for an FFT but I need to merge them, because they seems to be working well in some parts. Let me explain. First code:

fft1 = (Bx[51:-14])
fft2 = (By[1:-14])

# NS antena FFT - red
FFTdata = np.sqrt(fft1*fft1)**1
samples = FFTdata.size

# WE antena FFT - blue
FFTdata2 = np.sqrt(fft2*fft2)**1
samples2 = FFTdata2.size

# Adjusting FFT variables
duration = 300 # in seconds
Fs = float(samples)/duration # sampling frequency (sample/sec)
delta_t = 1.0/Fs
t = np.arange(0, samples, 1)*delta_t
FFTdata_freq = np.abs(np.fft.rfft(FFTdata))**1
FFTdata2_freq2 = np.abs(np.fft.rfft(FFTdata2))**1
freq = np.fft.rfftfreq(samples, d=delta_t)
freq2 = np.fft.rfftfreq(samples2, d=delta_t)

# Printing data
plt.semilogy(freq, FFTdata_freq, color='r')
plt.semilogy(freq2, FFTdata2_freq2, color='b')
plt.xticks([0,20,40,60,80,100,120,140,160,180,200,220,240,260,280,300, 
320,340,360,380,400,420,440])

Second code:

# Number of samplepoints
N = 600
# sample spacing
T = 300.0 / 266336.0
x = np.linspace(0.0, N*T, N)
y = np.sin(50.0 * 2.0*np.pi*x) + 0.5*np.sin(80.0 * 2.0*np.pi*x)
yf = fft(y)
xf = np.linspace(0.0, 1.0/(2.0*T), N/2)
plt.plot(xf, 2.0/N * np.abs(yf[0:N/2]))

And now the problem. I have two codes for FFT, but I do not know how to make them work. Well first code shows my data correctly, on the chart you have upper draw, but the scale is wrong, I need that upper chart situated where the bottom draw is placed. I do not know how to do this. Any ideas? fft1 and fft2 are data arrays. Everything happens in 300s = 300000ms.

After changing code thanks to @zck it looks like this from scipy.signal import welch

plt.subplot(212)
plt.title('Fast Fourier Transform')
plt.ylabel('Power [a.u.]')
plt.xlabel('Frequency Hz')
fft1 = (Bx[51:-14])
fft2 = (By[1:-14])

for dataset in [fft1]:
    dataset = np.asarray(dataset)
    psd = np.abs(np.fft.fft(dataset))**2.5
    freq = np.fft.fftfreq(dataset.size, float(300)/dataset.size)
    plt.semilogy(freq[freq>0], psd[freq>0]/dataset.size**2, color='r')

for dataset2 in [fft2]:
    dataset2 = np.asarray(dataset2)
    psd2 = np.abs(np.fft.fft(dataset2))**2.5
    freq2 = np.fft.fftfreq(dataset2.size, float(300)/dataset2.size)
    plt.semilogy(freq2[freq2>0], psd2[freq2>0]/dataset2.size**2, color='b')

I applied some changes. I am only missing Hamming Window, can anybody help, to make from this chart:

That one:

Answer 1

From the quick look it appears that in the upper code snippet you're forgeting to divide by N. It's a math probelm, not a code problem.

In general, if you're after spectra/power spectra use a WOSA(=overlapped segment averaging) method that applies window functions and averages if the sample size allows it. The welch method is included in scipy.signals, you should explore the scipy.signal library as it is extremely convenient in signal analysis.

For your datasets following code should work:

from scipy.signal import welch

plt.figure()
for dataset in [Bx, By]:
    dataset = np.asarray(dataset)
    freq, psd = welch(dataset, fs=dataset.size/300, return_onesided=True)
    plt.semilogy(freq, psd/2)

Note the psd is divided by 2 for return_onesided, if False division by 2 is not required. Hope this helps to produce good looking graphs. Above plots the power spectral density. Pass an argument scaling='spectrum' in case you want power spectra instead of power spectral density.

You can also pass arguments for window functions, the default is hanning but it includes most common windows like blackman, hamming, boxcart etc.

You can find more about this at https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.signal.welch.html